Polymorphic self-assembly of helical tubules is kinetically controlled

In contrast to most self-assembling synthetic materials, which undergo unbounded growth, many biological self-assembly processes are self-limited. That is, the assembled structures have one or more finite dimensions that are much larger than the size scale of the individual monomers. In many such cases, the finite dimension is selected by a preferred curvature of the monomers, which leads to self-closure of the assembly. In this article, we study an example class of self-closing assemblies: cylindrical tubules that assemble from triangular monomers. By combining kinetic Monte Carlo simulations, free energy calculations, and simple theoretical models, we show that a range of programmable size scales can be targeted by controlling the intricate balance between the preferred curvature of the monomers and their interaction strengths. However, their assembly is kinetically controlled—the tubule morphology is essentially fixed shortly after closure, resulting in a distribution of tubule widths that is significantly broader than the equilibrium distribution. We develop a simple kinetic model based on this observation and the underlying free-energy landscape of assembling tubules that quantitatively describes the distributions. Our results are consistent with recent experimental observations of tubule assembly from triangular DNA origami monomers. The modeling framework elucidates design principles for assembling self-limited structures from synthetic components, such as artificial microtubules that have a desired width and chirality.

Motivated by recent DNA origami experiments in which triangular monomers assembled into icosahedral capsids [1] or helical tubules [2], we consider a model in which triangular monomers assemble into tubules. The tubule structure can be represented as a two-dimensional triangular lattice that is rolled up and closed upon itself in the three-dimensional space. In Fig.  §1, we triangulate a two-dimensional sheet and label the three lattice directions in three different colors. Depending on how the edges where the structure closes match, different tubule structures are obtained. We categorize structures using the naming convention employed for carbon nanotubes, in which a pair of integer numbers (m, n) designates a structure [3]. These integers can be identified by starting from an arbitrary vertex, and determining the number of bonds along two directions for the shortest pathway (labeled in green in Fig.  §1) that travels around the tubule circumference and returns to the starting vertex. m is the number of bonds with primary orientation in the circumferential direction (gray bonds in Fig.  §1) while |n| is the number primarily longitudinal bonds (red or blue bonds). In our convention n is positive for a right-handed helix (if the pathway includes red bonds in Fig.  §1) and negative for a left-handed helix (if the pathway includes blue bonds).

B. Calculating preferred dihedral angles for different tubule geometries
To simplify the model, we assume that all triangular monomers are flat, equilateral, and identical. Unlike a flat triangular lattice, the dihedral angles between a triangular monomer and its neighbors cannot all be 180 • because the tubule has a non-zero principal curvature along the circumferential direction. To find the preferred dihedral angles, we first calculate the position of vertices. In a cylindrical coordinate system with the origin sitting at the longitudinal axis of the tubule, the radial distances ρ are the same for all the vertices, while the azimuth angles φ and heights z vary. As an example, Fig.  §2 shows a schematic of a (6,-2) tubule geometry. We label the diameter of the tubule D ≡ 2ρ, the pitch height h of the spiral line formed by gray bonds, the length of edge l 0 , the azimuth angle difference ∆φ i , and height difference ∆z i between neighbouring vertices connected by different bonds. i = 1, 2, 3 corresponds to blue, red and gray bonds. ∆φ i and ∆z i are the same for any lattice point because of translation invariance. h is positive for the right-handed tubule and negative for left-handed tubule. ∆φ i and ∆z i of an arbitrary tubule geometry (m, n) should satisfy the following equations, m∆φ 3 + max (n, 0)∆φ 1 − min (n, 0)∆φ 2 = 2π m∆z 3 + max (n, 0)∆z 1 − min (n, 0)∆z 2 = 0 ∆φ 1 + ∆φ 2 = ∆φ 3 ∆z 1 + ∆z 2 = ∆z 3 We obtain the first two equations by following the shortest pathway that travels around the circumference of the tubule, the third and the forth equation by circulating around the vertices of a triangular monomer, and the final equation based on the equilateral assumption. With these equations, we solve the positions of the vertices given an arbitrary tubule geometry (m, n) as well as the dihedral angles between adjacent monomers. We define ideal dihedral angles θ respectively as the angles between adjacent monomers that share blue, red and gray bonds in the tubule geometry (m, n). id,2 | decrease. All dihedral angles become closer to 180 degrees, since as the diameter of the tubule structure increases, the cylindrical surface becomes closer to a flat surface. As |n| increases, both D and h increase, θ to the curvature become more asymmetric.

II. COMPUTING TUBULE WIDTH FLUCTUATIONS IN THE CONTINUUM LIMIT
We employ a continuum equilibrium model to compute the width distribution of the assembled tubule structures. Consider a tubule with diameter D and length L. The number of monomers N in the tubule is then given by with a 0 the area of the monomer. Using the Helfrich model for the bending elastic energy, the free energy per monomer is given by in which E B is the binding energy, γ is the line tension, T is the temperature, and s is the per-monomer entropy in the tubule structure. The coefficient 3/2 in the binding energy term comes from the fact that a triangular monomer has three edges and each bond is shared by two edges.B is the bending modulus (the continuum limit of the discrete bending modulus B) given by [4] The equilibrium probability density ρ D,L to assemble a tubule structure with diameter D and length L is given by where µ is the chemical potential and N is related to L and D by Eq. (1). Using this distribution, the root-mean-squared diameter fluctuation ∆D is given by, (5) in which the mean of the distribution ⟨D⟩ is given by, The argument of the exponential can be expanded to To simplify the results, we assume that γ L is negligible in the large L limit. Furthermore, at equilibrium the free energy per monomer of the geometry that minimizes the free energy (in this case the target geometry) is approximately equal to the chemical potential µ [5]. Since the bending energy of the target geometry is zero, the equilibrium chemical potential is approximately −( 3E B 2 + T s * ) with s * the entropy per monomer of the target structure. Assuming the entropy is roughly independent of geometry, we have 5 Then, we can perform a saddle point approximation around D = D 0 , setting This results in the mean diameter ⟨D⟩ = D 0 . Using Eq. (3), ∆D is then given by  shows examples of defect-free, defective tubules and structures that fail to close. In the snapshots, boundary edges are labeled in cyan. (C) shows the distribution of the fraction of successful searches for 500 simulations run with the same set of parameters. We perform the search algorithm for all the vertices in the final structure and compute the fraction of successful searches for each simulation. We identify a structure as a defect-free tubule only when over 80 % of the searches are successful. Simulation parameters: EB = 6.4 kBT, B = 25 kBT , the target tubule geometry is (15,0). finsert is the attempt frequency of insert/delete move. finsert = 0.001 for all the figures except for Fig.  §14.
We classify the structures that form in the simulations by searching for the shortest path around the circumference of the structure, moving along the edges of the triangles. Fig.  §3(A) illustrates one such path. We choose a random vertex as the origin. Then, we search for the next vertex connected to the origin by a gray edge, which we call the second vertex. The vector pointing from the origin to the second vertex defines the direction of the search. Then, we search for vertices along either the red or the blue edges and measure the distance between those vertices and the origin. If those vertices are farther from the origin, we move to the next vertex that is connected through a gray edge. If those vertices are closer to the origin, we move along the red/blue edge until we reach the boundary or return to the origin. If the search along red or blue edges reaches the boundary, we move to the next vertex that is connected through a gray edge. This search continues until we return to the origin or terminates if the last vertex connected through a gray edge is on the boundary. If the search can return to the origin, we label it as a success, otherwise we label it as a failure. The lattice numbers m and n are determined by the number of edges traversed on this pathway. m is the number of the gray edges, while |n| is the number of red or blue edges. n is negative if the pathway follows blue edges (a left-hand helix) and positive if the pathway follows red edges (right-hand helix). Fig.  §3(B) shows examples of defect-free tubules, defective tubules, and structures that fail to close. To identify defective tubules and to prevent incorrect identification, we perform the search algorithm described above for every vertex in the final assembly, and then record the successful search rate (i.e., the fraction of starting vertices for which the search algorithm is successful). Fig.  §3(C) shows a distribution of the fraction of successful searches, measured for 500 independent simulations under the same parameter set. By comparing the simulation snapshots and the distribution, we find that the tubes are mostly defect-free around the peak (≈ 0.9 success fraction) and mostly defective when the fraction of successful searches is below 0.8 or multiple tubule geometries are identified in the same structure. That is, we identify a structure as defect-free only if more than 80% of the searches find the same tubule geometry.  The left panel shows the structure of a stress-free tubule. Lmax (solid line) is the distance between two vertices that are farthest apart within the same lap, which is the shortest pathway that traverses the circumference of the structure and returns to the original vertex (dashed line). D is the diameter of the stress-free tubule. RLD is defined as the ratio between Lmax and D for calibration. To get the diameter of an assembled tubule, we pair vertices that are the farthest apart within the same lap and measure the distance between those pairs. We compute ⟨Lmax⟩sim by averaging over all paired vertices within the same lap. The local diameter D local of the assembled tubule in the simulation is computed by ⟨Lmax⟩sim/RLD. We compute the diameter of the tubule as the mean value of D local computed from all laps within the structure. (B) A detailed diameter distribution of 500 defect-free tubule structures assembled under the same parameters set. The colors of the bars represent different tubule geometries. Vertical dashed lines label the diameters of the ideal tubule geometries. Simulation parameters: EB = 6.0 kBT, B = 25 kBT , the target tubule geometry is (10,0), and f fusion = 0.01. 7 Fig.  §4(A) shows how we measure the diameter of the tubule structure we obtained from the simulation. Basically, we perform this measurement by comparing the tubule assembled in the simulation with the stress-free tubule as solved in the previous section I B. The diameter of the stress-free tubule is well defined because all the vertices are on the same cylindrical surface. The diameter of the tubule D is defined as the diameter of the cylindrical surface where vertices sit. On the other hand, we identify the shortest pathway that traverses the circumference of the tubule and returns to the original vertex (the dashed black line on the right panel in Fig.  §4(A)), which we call a lap. Within the lap, we find vertices that are the furthest apart and define the distance between them as L max (labeled by the solid black line). We compute R LD as the ratio between L max and D. R LD is different for different tubule geometries. It increases as the tubule becomes more chiral. We use R LD to calibrate the local diameter of the tubule assembled in the simulation.
To get the diameter of a tubule assembled in the simulation, we first identify all the laps within the tubule structure. Then within each lap, we pair vertices that are the farthest apart and compute ⟨L max ⟩ sim by averaging over all paired vertices. The local diameter D local of the lap is computed by ⟨L max ⟩ sim /R LD . Finally, we compute the diameter D of the structure as the mean value of D local computed from all laps within the structure. Fig.  §4(B) shows a detailed distribution of the measured tubule diameter values D. The counts are labelled by color based on the tubule geometry. We can see that fluctuations of D within the same tubule geometry are small, and the mode of each distribution is close to the diameter of the ideal tubule geometry (labeled in dashed lines). The fluctuations around each peak come from thermal fluctuations of each particular tubule geometry itself. The observed fluctuations of tubule width are thus primarily determined by assembly of different tubule structures, not by our measurement method. In dynamical assembly trajectories, we define the closure size N close as the at structure size N at which an edge fusion move is accepted for the first time in that trajectory. Although it is possible for the fused edge to reopen, the probability of reopening is low and the error this contributes to the closure size measurement is small. Fig.  §5 shows closure size distributions for different target tubule geometries measured from simulations, compared with the kinetic model prediction (Section III D of the main text). The mean closure sizes ⟨N close ⟩ for the target geometries (5,0) and (10,0) are respectively ⟨N close ⟩ = 22 and 120. Note that, in addition to target geometry, the closure size distribution is affected by the binding energy, the bending modulus, and the effective closure frequency f fusion . However, within the parameter ranges discussed in the main text, the variations in ⟨N close ⟩ with different values of these parameters 8 for a given target geometry are typically small (≲ 20%).
The length of a tubule at the point of closure L close can be approximately related to the closure size by L close πD 0 = ⟨N close ⟩a 0 (12) in which D 0 is the diameter of the target geometry and a 0 is the area of the monomer. By comparing L close and D 0 at different target geometries, we find L close typically varies in the range 1.3D 0 ≲ L close ≲ 1.7D 0 . We use 1.5D 0 in the main text for estimating the equilibrium tubule geometry distribution at the closure size. We find that the closure size predicted by the kinetic model is approximately twice as large as the simulation result (Fig.  §5). We hypothesise this is due to the difference between the preclosure geometry we assume in the model and the actual geometry of the structure assembled in the simulation. In the model, we assume that the open structure has a circular shape with a smooth boundary. However, the structure in the simulation has a coarser and more irregular boundary due to the triangular lattice geometry, entropic edge fluctuations, and the nonequilibrium nature of assembly trajectories. The boundary roughness and non-circular shape may favor closure of the structure at a smaller sizes, by bringing distal edges into closer proximity for a given bending energy cost than would be possible for a perfectly circular structure. Nevertheless, despite this discrepancy the kinetic model reproduces the probability distribution measured from the simulation with reasonable accuracy.  Performing the tubule-geometry identification at different time points along simulation trajectories shows that the tubule geometry is essentially fixed after closure. Fig.  §6 shows identified tubule geometries as a function of time for a selection of trajectories at different parameter sets. In all cases, once closure occurs, the tubule geometry is fixed for the remainder of the simulation. This observation is consistent over the full range of parameters that we explore.  Structures that failed to close formed paper roll structures ( Fig.  §3(B)). Results are shown for indicated values of the effective closure rate f fusion , which we controlled by varying the frequency of edge fusion/fission attempts. As f fusion increases by two orders of magnitude, the closure probability increases from 55 % to ∼ 100 %. The kinetic model predictions capture this trend and nearly quantitatively agree with the simulation results.

F. Tubule growth rates
To describe the growth dynamics, we define the growth rate k grow as the steady-state time-derivative of the structure size N . Fig.  §8(A) shows how we measure k grow . We see that once the tubule structure has closed and filled gaps in the structure to achieve a roughly steady-state morphology, the net growth rate fluctuates around a constant value. The fluctuations arise due to stochastic association events and because monomers can detach from the structure (in particular from the two open ends where monomers have fewer contacts) due to the relatively weak binding energy E B required for assembly of well-formed tubules. To estimate the net growth rate, we choose a time point τ 0 and perform a local time-average by performing linear regression on N as a function of time τ within a range of 2 × 10 5 time steps centered at τ 0 , and compute the slope of the linear fit. Designating the structure size at τ 0 as N 0 , the slope of this linear fit provides one measurement of the growth rate at size N 0 . We compute the growth rate by averaging over 500 measurements under the same parameter for each value of N or N boundary . We select 2 × 10 5 time steps as the fitting range because this is typically the smallest time interval for which there is a significant change in N (> 10) beyond the critical nucleus size N C within the parameter space that we simulated. By choosing the smallest such time interval, we avoid over-smoothing fine features on the N, τ curve. Fig.  §8(B) shows k grow as a function of N for different assembled tubule geometries (colored symbols), as well as points along the trajectories before tubules have closed (open symbols). We observe that k grow increases with N before closure and plateaus after closure once the tubule geometry reaches a steady state. The steady-state growth rate increases with the width of the assembled tubule geometry. The inset of Fig.  §8(C) shows that the growth rate of closed tubules is proportional to m, which determines the tubule diameter and thus the mean value of N boundary for a closed structure. To understand these trends, we measured k grow as a function of the number of boundary edges N boundary ( Fig.  §8(C)). We find that the growth rate increases roughly linearly with N boundary , but with a non-zero intercept of N boundary ≈ 7. The non-zero intercept is reasonable since the growth rate should approach zero as N approaches the critical nucleus size N C . For this parameter set, N C = 5 monomers, and the number of boundary edges at the critical nucleus size ranges from 7 ≲ N boundary ≲ 12, consistent with the observed intercept. We observe that the steady-state growth rate for closed structures is somewhat larger than that for open structures, but is roughly linearly related to the number of boundary edges N boundary . This trend is consistent with stable monomer attachment events at each open edge, as would be expected for conditions under which monomer addition is significantly biased over monomer detachment. However, the net growth rate is somewhat smaller for open structures compared to closed structures at the same value of N boundary . This trend can be understood from the fact that addition of a monomer with only a single contact is relatively unstable at these parameter values due to binding entropy penalties. Thus, stable monomer addition requires formation of at least two contacts. Formation of two (or more) contacts occurs more readily for a closed structure compared to an open structure due to the larger curvature of the rim in the latter case.
Despite the geometrical differences between open and closed tubule structures, k grow depends linearly on N boundary during growth in both phases. Below closure, considering the limit that the critical nucleus size is small and N boundary scales with √ N before closure, k grow can be approximated by in which k 0 grow is a factor that can be measured from the simulation, which depends on B, E B and the attempt frequency of the monomer insertion/deletion moves. Fig.  §8(D) shows the growth rate k grow of assembled tubules with geometry (10,0) measured at structure size N = 150 (k grow becomes independent of size shortly after the tubule closes). k grow increases as E B increases and B decreases. Except when the growing disk has a size commensurate with closure of a stress-free tubule (that is, the diameter of the disk is related to the preferred diameter of the tubule by D disk = πD 0 so that when the edges of the disk touch each other the structure has the stress-free-curvature), there will be a free energy barrier to closure arising from the bending energy required to adopt a curvature (i.e. monomer-monomer binding angles) that differ from the preferred geometry. In particular, the disk must curve in such a way that its edges approach each other closely enough to bind. Once these contacts are made and the monomers bind, then the additional binding energy can compensate the unfavorable bending energy. Thus, an estimate of the free energy barrier to closure can be obtained by computing the bending energy difference between the disk at the point just before closure occurs, and in its stress-free state. Here, we are making the good approximation that we can neglect configurational entropy differences between these states. Fig.  §9(A) shows a schematic of the closure process. Depending on the size of the disk, different tubule geometries can be formed by closure, and the free energy barrier is a function of both the disk size and the geometry being formed.
We have used several approaches to compute the closure free energy barrier. First, we measured the closure rates for an ensemble of unclosed proto-tubule configurations, constructed to have well-formed triangular lattices such as the ones that occur during assembly dynamics under productive assembly conditions ( Fig.  §9(B)). The pre-closure structure that we choose resembles a regular hexagonal geometry if flattend onto the two dimensional plane, yet with a 'notch' for accepting the edge fusion move. The reason to choose this spatial configuration is that we have identical binding energies for all three edges and the growth is nearly isotropic before closure. We can use a pair of numbers [a, b] to describe the geometry of the pre-closure structure, in which a is the number of triangles along the long (horizontal) axis and b is the number of triangles along the short (vertical) axis of the flattend pre-closure structure. The pre-closure structures we use for the closure rate measurement are [8,8], [9,8], [10,10] and [11,10], which assembles (8,0), (9,0), (10,0), (11,0) tubules when they close. The aspect ratio of these pre-closure structures is roughly held the same. The target tubule geometry is set to (10,0) for all the closure rate measurements, however, the pre-closure structure limits the tubule geometries that can assemble via closure.
First, we generate each pre-closure structure as a two dimensional triangular lattice using a home-made MATLAB code and use that as the initial state in our dynamical simulations. Then, we perform only vertex moves to equilibrate the configuration at fixed size and topology with the preferred dihedral angles favouring a (10,0) tubule. We performed 10 4 sweeps to make sure that the pre-closure structure equilibrated. Finally, we allowed both vertex and edge fission/fusion moves (so that the structure can close), and measured the fraction of structures that remain open as a function of elapsed time. We performed this measurement for various diameters (the long axis of the open structure, which equals the circumference of the closed structure) and B values for a target geometry of (10,0).
The bottom plot in Fig.  §9(B) shows that the fraction of open structures decays exponentially as a function of time. By performing an exponential fit with decay time τ c for each data set, we extract the closure rate k close = 1/τ c . We observe that k close decreases as the bending modulus B increases or the assembled tubule geometry deviates further from the target tubule geometry.
Because the measured closure time distributions are exponential, we assume that the closure rate follows the Arrhenius mechanism, and we estimate the closure free energy barrier ∆G with k 0 as the closure rate when the diameter of the unclosed structure equals the circumference of the target tubule geometry. bend . Considering that all monomers in the open structure have dihedral angles that favor the tubule geometry (m, n), the bending elastic energy difference between such any configuration and the stress-free configuration ∆E (m,n) bend is given by We observe that ∆G bend with 0.3 ≲ α ≲ 0.5 as discussed in the main text. Therefore, the closure rate can be approximated by, In the inset of Fig.  §9(C), we also compare the elastic energy difference measured in the simulation ∆E (m,n) el to the bending energy estimate ∆E (m,n) bend . We measured ∆E (m,n) el as the difference between the elastic energy just before the structure closed and the average elastic energy of the structure at equilibrium (with fixed size and topology). Notice that we observe the same scaling, ∆E bend , which supports the assumption above that the closure free energy barrier arises primarily due to the elastic energy difference between the closed and open configurations. Nevertheless, the measured elastic energy difference is only ∼ 30 % of the elastic energy if the entire unclosed structure adopts a curvature commensurate with the closed structure (inset of Fig.  §9(C)). We hypothesize that this discrepancy arises because the bending is nonuniform and only a few contacts are required to initially stabilize the closed structure. For the kinetic model, we use α = 0.3 to account for these effects.

13
H. The dependence of assembled tubule geometries on binding energy Fig.  §10 shows the geometry distributions of tubules assembled in simulations with target structure (10,0) and different values of the binding energy E B . As E B increases, the fraction of defect-free tubules decreases, while the fraction of defective tubules and open structures increases. This is consistent with observations in other self-assembly systems that the probability of defective structures increases with binding energy because the rate of annealing defects decreases exponentially with increasing binding energy [6,7]. However, the geometry distribution within the defectfree population does not change significantly. This can be understood as follows, based on the results of the kinetic model (section III of the main text). The tubule closure size is substantially larger than the critical nucleus size for the parameter regime that we focus on, so the growth dynamics is strongly forward-biased and thus the growth rate depends only weakly on E B (or µ). Therefore, although the kinetic model predicts that the width distribution will shift to larger values with increasing growth rate, this effect will be very weak with increasing E B when the closure size is large compared to the critical nucleus size.    §12 shows the schematic of our kinetic model. Prior to closure, the assembling structure is modeled as a curved circular disk. The diameter of the flattened disk is D disk . The pre-closure structure grows with a growth rate k grow , and in the mean time, attempts closure with a closure rate is k close . The transitions that are marked with red 'x' symbols in the diagram are considered low probability, which means that the closed structure does not reopen. D is the diameter of the closed assembled tubule structure, while D 0 is the diameter of the target tubule structure. close (N ) is a function that indicates whether a structure of size N is geometrically compatible with closing into a structure (m, n). We compute I (m,n) close (N ) in the following simplified manner. Assuming a circular disk forms a cylinder geometry by bending around a single axis and closing where the two outermost points meet, the curvature κ of the closed structure is However, we must account for the fact that only discrete tubule geometries are allowed. In general, closure to many geometries will require fluctuations from the perfect circular disk that we have assumed thus far. In SI section V C 2 and V D, we perform a calculation in which we explicitly compute the Boltzmann-weighted probability of disk shape fluctuations in the continuum limit. That calculation and analysis of aspect ratios of unclosed disks within simulation trajectories show that the magnitude of shape fluctuations is generally small. Therefore, we simplify the calculation here by neglecting shape fluctuations, but assuming that closure only occurs for geometries that have curvatures close to the circular disk curvature: κ ∈ [(1 − δ)κ(N ), (1 + δ)κ(N )]. δ is the cutoff threshold δ, which we determined from the extent of aspect ratio fluctuations of open structures within simulation trajectories (See Fig. §13 and SI section V D for details about the aspect ratio fluctuation measurement. Within the parameter space that we have studied, the surface energy is primarily determined by the binding energy, and the shape distribution of disks near closure size does not vary dramatically with B. Thus, for simplicity we use the same value of δ for all parameters with the same target geometry. We use δ = 0.05 for the (10, 0) target.
The accessibility indicator function is then given by with Θ the Heaviside function.

C. Continuum limit
In the main text, we introduce a discrete kinetic model that reproduces the distribution measured in the simulation. Here, we perform an analogous calculation in the continuum limit. We find that the resulting model qualitatively matches the simulation results, but is not as quantitatively accurate as the discrete model.

Model without shape fluctuations
First, we assume that the preclosure geometry of the structure is always circular, without shape fluctuations. If such structure closes, the diameter D of the closed tubule is given by where N is the number of monomers in the assembled structure, a 0 is the area of each monomer, and D disk is the diameter of the flattened disk. If the disk closes, D disk will equal the circumference of the closed tubule. Thus, the diameter D of the closed tubule is related to the structure size N by In terms of the kinetic rates, we again approximate k grow by where k 0 grow is a kinetic factor that can be measured from simulations. The closure rate k close is approximated by We use the Helfrich energy model to estimate the bending energy differences between tubules with different diameters [8].
where D is the diameter of the tubule that the open structure is about to close into, D 0 is the diameter of the target tubule geometry, andB is the effective bending modulus at the continuum limit, which is related to the bending modulus B of the discrete triangular lattice byB = √ 3B/2 [4]. The closure rate k close can be further simplified to D 0 and N 0 are related by, where k 0 close is a kinetic factor that can be measured from simulations. Finally, we evaluate the closure probability as a function of time. In the absence of closure, the structure remains at each size N for a time ∆t(N ) = 1/k grow (N ), and the time at which a structure first grows to size N is t N = (N ). The probability that a structure with size N remains open at a time t N ≤ t < t N +1 is given by By summing over the smaller sizes, we can compute the probability that a structure has stayed open until size N as The probability for the structure to close at size N is then given by Since D and N are related by Eq.20, the width distribution is equivalent to the closure size distribution.

Model with shape fluctuations
To incorporate shape fluctuations of the pre-closure structure into the model, we consider the structure before closure as elliptical rather than circular. The aspect ratio ψ is defined as the ratio between the lengths of the major and minor axes. The orientation angle ϕ is defined as the angle between the major axis and the direction of the non-zero principle curvature (the circumferential direction of the tubule once it closes). We assume that ϕ follows a uniform distribution from 0 to π. Assuming the system is quasi-equilibrated before closure, the distribution f of ψ, ϕ should follow Boltzmann distribution. According to classical nucleation theory (CNT), the free energy barrier G to forming an open structure from unassembled monomers can be written as where ∆g is the per-monomer bulk free energy difference between the structure and the fluid phase, σ is the permonomer surface free energy, N is the number of monomers in the structure, and N boundary is the number of monomers at the structure boundary, which depends on ψ.
Assuming ∆g does not change with ψ, the bulk term depends only on N and the shape distribution is completely determined by the surface energy of the open structure. The distribution f of ψ and ϕ at size N is given by where Z ψ is the partition function obtained by summing over ψ and ϕ, assuming that the latter follows a uniform distribution, which contributes a factor of 1/π. l 0 is the stress-free length of the edge and L(N, ψ) is the perimeter of the open elliptical geometry, which is given by where a 0 is the area of each monomer, a is half the length along the major axis, and b is half the length along the minor axis. The first equation enforces that the area of the structure is conserved for fixed N . The second equation gives the perimeter of the elliptical geometry. The surface energy density σ can be measured from thermodynamic integration as introduced later in SI section VII. In terms of the kinetic rates, we again approximate k grow by where k 0 grow is a kinetic factor that can be measured from simulations. The closure rate k close is approximated by We use the Helfrich energy model to estimate the bending energy differences between tubules with different diameters [8].
where D 0 is the diameter of the target tubule geometry andB is the effective bending modulus at the continuum limit defined above. D is the diameter of the tubule that the open structure is about to close into, and is given by This equation is obtained by assuming that the circumference of the closed tubule geometry is equal to the length of the elliptical geometry along the direction of the non-zero principal curvature. Finally, we evaluate the closure probability as a function of time. In the absence of closure, the structure remains at each size N for a time ∆t(N ) = 1/k grow (N ), and the time at which a structure first grows to size N is t N = (N ). The probability that a structure with a certain geometry remains open at a time t N ≤ t < t N +1 is given by By summing over ψ and ϕ, we can compute the probability that a structure with size N remains open at a time t N ≤ t < t N +1 as By summing over the smaller sizes, we can compute the probability that a structure has stayed open until size N as The probability to assemble a tubule geometry with diameter D is then computed by summing over all sizes n that can close into the diameter D, where S D is the subset in parameter space of ψ and ϕ that enables the structure to close with diameter D.  32)). The structure size is N = 120, which is also the average closure size of the (10,0) tubule. To measure the aspect ratio distribution in simulations, we first performed simulations in which the size was allowed to grow to N = 120. Then we changed the preferred dihedral angles to zero and increased the bending rigidity to B = 400 k B T so that the structures flattened. Then, we rotated the flattened structure to minimize deviations of the vertex positions from the X-Y plane, and fit the positions of the vertices on the structure boundary to an ellipse. The fitting was performed by minimizing a distance function that computes the sum of the distances from vertices on the structure boundary to surface of the trial ellipse. The statistics were collected over 1000 independent simulation trajectories. In the model prediction, the surface energy density σ is obtained from fitting to the G − N free energy curve introduced later in SI section VII, which gives σ ∼ 4.5 k B T per monomer. For the simulation measurements, the mean of the aspect ratio distribution is 1.20, which is close to the aspect ratio of a perfect hexagon ( ∼ 1.15). The standard deviation of the distribution is 0.49. For the model prediction, the mean and standard deviation of the distribution are 1.12 and 0.13 respectively. Comparing these two distributions, we find the spread of the model prediction is narrower and the peak is closer 1 (a perfect circle). We hypothesize this is because the actual pre-closure geometry is a triangulated lattice and it does not have a smooth boundary, contrary to the model assumption.

E. Accounting for shape fluctuations in the discrete model
For the discrete model in the main text, we introduce a cutoff threshold δ to approximately account for shape fluctuations. We use δ = 0.05 for the system with target geometry (10,0). If we assume thet the geometry of the open structure is elliptical, the aspect ratio is 1.1 at the threshold we choose, which is close to the mean of the aspect ratio distribution given by the model prediction. Furthermore, the fraction of structures that fail to close is sensitive to δ. Setting δ = 0.05 leads to a close match between the fraction of structures that fail to close between the discrete model and simulations (see Fig.  §7).
F. Comparing the models with and without shape fluctuations to simulation measurements Figure. §14 compares the tubule diameter distributions estimated from the continuum models with (blue dashed bars) and without (red solid bars) shape fluctuations against simulation measurements ('•' symbols). We observe that predictions from both of the models accurately estimate the fraction of the distribution near the mean. However, the model which accounts for shape fluctuations more accurately describes the skew measured in simulations, in particular that the skew is biased toward structures that are larger than the mean. We hypothesize that this is because the shape fluctuations allow the pre-closure structures to close to tubule geometries with larger diameters than are accessible to a perfect circle. Therefore, when the size of a pre-closure structure is small, the shape fluctuation would drive the structure to close at a diameter closer to the target diameter D 0 because the closure energy barrier is lower, resulting in a larger mean tubule width when shape fluctuations are accounted for. However, this effect is small and quantitative; the two models give qualitatively similar predictions throughout the parameter space that we have studied. The mean diameterD as a function of B. D0 is the diameter of the ideal (10,0) tubule, the target tubule geometry is (10,0), and EB = 6 kBT . Figure.  §15(A) shows that the effective closure rate prefactor (not accounting for the free energy barrier in Eq.(13) of the main text, Eq. (15) of the SI) increases with the bending modulus increases at fixed f fusion . In principle, the faster closure rates favor the tubule width distribution to shift to smaller widths (since the tubule rarely reopens once it closes). However, Figure.  §15(B) shows that the mean diameterD increases as B increases, although the effect is very small. This effect arises because increasing B increases the free energy barrier, thus reducing the closure rate for tubules with diameters that deviate from D 0 , thus reducing the probability for premature closure. This elastic energy effect dominates over the prefactor, and shifts the diameter distribution closer to D 0 with increasing B.

VII. THERMODYNAMIC INTEGRATION
We use an adapted thermodynamic integration method to compute the free energies of different tubule geometries. Our primary goal is to estimate the bulk and surface free energy densities that are required for the CNT calculation of the free energy profile as a function of assembly size.

Initial Configuration
Evolution Trajectory The motions of particles that interact with the magenta, blue, and orange vertices (in the Einstein sold) are constrained to eliminate rigid body motions of the entire structure. The particle interacting with the magenta site is fixed at its ideal lattice position.
The particle corresponding to the blue vertex is only allowed to move along the direction that connects between the magenta and the blue vertices. The particle corresponding to the orange vertex is only allowed to move in the plane that contains the magenta, blue, and orange vertices.

A. Thermodynamic integration implementation
We implemented thermodynamic integration as follows. First, we construct an ideal solid lattice corresponding to the helical tubule geometry that our target states are based on. The spatial coordinates of the vertices of the lattice are computed by the method we introduce in the SI section I B. Then, we apply the following Hamiltonian to the system, Here λ is a parameter that increases from 0 to 1 during each thermodynamic integration run so that the Hamiltonian changes from U 0 to U 1 . U 0 is an Einstein solid potential, in which each particle is connected to the ideal position of the lattice vertex with a spring. ⃗ r i is the instantaneous position of the ith particle, while ⃗ r 0,i is the ideal position of the lattice vertex that the ith particle is connected to. K is the spring constant of the Einstein solid potential. We performed thermodynamic integration computations with K = 100 k B T /l 2 0 . U 1 is the Hamiltonian of our system. In each thermodynamic integration run, we perform Monte Carlo sampling with only vertex moves allowed to equilibrate the system at fixed topology. We perform independent runs with λ ranging from 0 to 0.9 with an increment of 0.05 and from 0.9 to 1.0 with an increment of 0.01. For each λ, we equilibrate the structure by performing 10 6 sweeps, with the number of vertex moves in a sweep equal to the number of vertices N .

21
The free energy difference between our system and the Einstein solid model is given by The integrand is the sample average of the difference between U 0 and U 1 when the system is evolved under the Hamiltonian H(λ). U 0 and U 1 are computed by applying the corresponding potential function to the instantaneous spatial configuration in the simulation. One challenge is that, as λ → 1, the Einstein solid potential weakens and the structure undergoes rigid body translations and rotations, which affect U 0 since the coordinates of the ideal lattice are fixed in the lab frame. Therefore, we use the method introduced in the Ref. [9] to prevent rigid body motions of the structure. We constrain the motion of three vertices in the structure as shown in Fig.  §16. The particle that interacts with the magenta vertex (in the Einstein solid model) is fixed. The particle that interacts with the blue vertex is only allowed to move along the direction that connects between the magenta and the blue vertices. The particle that interacts with the orange vertex is only allowed to move in the plane that contains the magenta, blue, and orange vertices.
The free energy of the constrained Einstein solid model is given by where N v is the number of vertices in the structure.

B. Thermodynamic integration results: Free energies of different tubule geometries
We used our thermodynamic integration implementation to compute the free energy ∆G of different tubule geometries assembled with the target geometry (10,0) as a function of the structure size N . Since the number of surface bonds is fixed and finite, the G − N curve should follow a straight line according to CNT. The bulk free energy density at finite temperature g (m,n) T finite can then be extracted from the slope of the linear fitting to the G − N curve. As a comparison against the thermodynamic integration results, we estimate the bulk energy density at zero temperature g (m,n) T0 as the contribution from the elastic energy, given by in which E B is the binding energy, i represents different edges of each triangular monomer, θ (m,n) id,i is the ideal dihedral angle that favours a (m, n) tubule geometry, and θ 0 id,i is the ideal dihedral angle of the target tubule geometry. The thermodynamic integration results indicate that the free energy difference between different tubule geometries primarily arises due to differences in bending elastic energy. Fig.  §17 shows that the free energies g (m,n) T finite scale linearly with the elastic energy estimates g (m,n) T0 , with slope 1, for all structures that we have considered. As B increases, g (m,n) T finite for all tubule geometries generally increases because the entropy loss increases by adding a monomer to the structure. The difference in g (m,n) T0 between different tubule geometries also increases with B because the difference in bending energy increases according to Eq. (46). The inset plots the difference between g (m,n) T finite and g 0 T finite with respect to the diameter D of the tubule, with B = 20 k B T . The solid curve is the estimation of the bulk energy density difference ∆g D between a tubule with arbitrary diameter D and the tubule with the same diameter as the (10,0) tubule based on the Helfrich energy model [8]: in which a 0 is the area of the monomer and D 0 is the diameter of the target tubule. This model only considers the bending energy, yet the free energy density difference g (m,n) T finite − g 0 T finite for most of the discrete tubule geometries falls onto the curved predicted by the model, except for the chiral tubule geometries such as (9,±1) and (10,±1). We hypothesis that this effect reflects the inconsistency between the triangular lattice direction and the curvature direction of the tubule structure itself that arises in chiral structures.
Despite the small geometric effect observed for chiral lattices, the thermodynamic integration results show that the primary factor determining the tubule free energy density is the bending elastic energy difference between a given structure and the target geometry as computed from the continuum Helfrich energy. This result also suggests that the bulk entropy density is similar for different tubule geometries under the same bending rigidity.   Table. I presents the parameter values that we used for the simulations shown in the main text. f insert , f wedge fusion and f crack fusion are the attempt frequency of the corresponding kinetic moves. (See the following SI Section IX and X for details about the algorithm and the moves).

IX. COMPUTATIONAL MODEL AND KINETIC MONTE CARLO DETAILS
In this section we provide additional details about the model and Monte Carlo simulations that we used to generate the results in the main text. The model and algorithm are closely based on ref. [10] and identical to the one used in ref. [11]. In particular, we consider flexible triangular subunits which can bind to each other along edges with a set of preferred dihedral angles that set the preferred curvatures of the assembling sheet. Monte Carlo simulations are performed in the grand canonical ensemble at fixed µV T , with µ the chemical potential of subunits in the bath. Each Monte Carlo simulation involves a single cluster undergoing assembly and disassembly, with subunits taken from or returned to the bath respectively, as well as structural relaxation moves.

Energies
In the tubule model, each three edges of the triangular subunits are of different types, t(p) = 1, 2, 3, for edge index p and each edge can only bind to an edge of the same type on a neighboring subunit.
The total energy of the system is given by where the first sum goes over all edges, with n s the number of subunits in the cluster. The second sum only goes over bound edges (i.e. non-boundary, adjacent edges, so there are 2n b terms in the sum, with n b as the number of bonds). The 1/2 factor corrects for double counting. The stretching energy is defined as: where k s is the stretching modulus, l p is the instantaneous length, and l 0 is the stress-free (rest) length of an edge. For the tubule model we set the stretching modulus and rest length equal for all edges. The bending energy is quadratic in deviations from the preferred dihedral angle: with p and q adjacent edges and t(p), t(q) the edge types. B is the bending modulus and is set equal for all edge types. θ t(p)t(q) 0 is the preferred dihedral angle between edges with types t(p) and t(q). Since only edges of the same types are allowed to bind to each other, t(p) = t(q) ≡ t for all adjacent edge pairs pq, and θ t(p)t(q) 0 ≡ θ t 0 . By adjusting θ t 0 independently for all three edge types t = 1, 2, 3 one can tune between different (m, n) target geometries. The binding energy between two edges p and q (with the same type t(p) = t(q) = t) is given by Binding energies corresponding to all edges are set equal to In addition to the above terms, each subunit has at its center of mass a spherical excluder of radius 0.2l 0 to prevent subunit overlaps. Finally, to prevent extreme distortions of subunits, maximum edge length fluctuations are limited to l 0 /2 < l < 3l 0 /2.

Coarse-graining
Our model is motivated by the triangular DNA origami subunits developed in Sigl et al. [1], in which subunits bind through lock-and-key 'patches' along subunit edges in which attractive interactions are generated through bluntend stacking of unsatisfied nucleotides. Therefore, in our model we define attractive bonds along subunit edges. In particular, attractive bonds occur at each shared pair of subunit edges with the same type. Because the interactions in the experimental system are driven by nucleotide stacking, they are extremely short-ranged in comparison to the subunit size (the subunit edge lengths are approximately 60 nm). Therefore, in our simulations we avoid resolving the short length scale fluctuations in separation distance between bound edges and their associated vertices by coarsegraining as follows.
A microstate i is defined as the position of all the 3n s vertices of n s subunits: i → (⃗ x 1 , ⃗ x 2 , ...⃗ x 3ns ) The grand canonical probability density of finding the system around state i is where µ is the chemical potential and λ 3 is the standard state volume. This probability density has the dimensions of 1/volume 3ns corresponding to all the 3n s vertices of the subunits. Due to bonds, however, some pairs of vertices are confined within a binding volume v a . We consider a square-well potential so that the binding energy is constant within this volume. Analogous to Ref. [10], we can then coarse-grain to avoid resolving intra-bond fluctuations. We assume that fluctuations of bound edges are sufficiently small that each pair of vertices at either end of a bound edge pair are constrained within a binding volume v a . Note that we constrain vertices rather than edges so that the coarsegrained microstate can be represented in terms of positions of vertices rather than edges, which is easier to implement computationally. In the coarse-grained system, a coarse microstate is specified by the coordinates corresponding to the independent vertex degrees of freedom (with 1 degree of freedom for each bound vertex group and unbound vertex): where n v is the number of independent bound vertex groups and free vertices. The probability of such a coarse-grained state is given by the net weight of all the corresponding fine-grained microstates: where n VB is the number of vertex-bonds and is given by n VB = 3n s − n v . For simplicity, we take the limit in which 3 √ v a is small in comparison to the length scale over which the elastic energy varies, so that the energy is constant within the bound volume v a . Then f (i) is a constant, and the probability density is given by where E Γ is the total energy of state Γ (including stretching, bending and binding energies). The coarse graining process is illustrated in Fig. §18.

Implementation and data structure
The simulation is implemented on top of the OpenMesh library [12]. Subunits are implemented as triangular mesh elements. OpenMesh uses the halfedge data structure which is suitable to implement triangles with directed normals (Fig. §19). The directed halfedges allow for a clockwise iteration through the boundary of a triangle, which makes the two faces of the triangles distinguishable. Only halfedges with opposite orientations can bind together, making it impossible to form a Mobius strip, for example. The data structure and the resulting iterators in OpenMesh allow for an easy and efficient iteration over the neighborhood of mesh elements (vertices, edges and faces). The implementation of mesh element rearrangements is less straightforward, but we implemented it via the insertion and removal of virtual triangles. In addition, OpenMesh allows for the storage of various properties on mesh elements, allowing storage of edge types and face types stored on the elements. To improve readability in the upcoming sections, we will not represent halfedges separately.

X. THE MONTE CARLO MOVES
In this section we detail the Monte Carlo moves of the simulation. Our algorithm has 11 moves: vertex displacement, simple subunit insertion/deletion, wedge insertion/deletion, wedge fusion/fission, crack fusion/fission, and edge Figure 18. Coarse-graining of an example cluster configuration. In this configuration, the number of subunits is ns = 5, the number of initial (before coarse-graining) vertices is 3ns = 15, and the number of vertices after coarse-graining is nv = 7. The red circles indicate bound vertex groups, and the number of vertex degrees of freedom that have been eliminated by coarsegraining in this configuration is nVB = 1 + 3 + 2 + 2 = 8 = 3ns − nv. Motivated by DNA origami subunits in Sigl et al. [1], the attractive interactions (i.e. 'bonds') in this model occur along edge-pairs of the same type shared by two subunits. In this configuration there are n b = 4 bonds. Figure 19. The halfedge data structure used by OpenMesh. Each edge is represented by two directed edges. Boundary edges are no exception and thus are represented by a non-boundary halfedge and a boundary halfedge (in green). This latter is irrelevant for our model. Directed edges allow for the unambiguous definition of face normals, for efficient iterations of the element's neighborhood as well as boundary iterations.

fusion/fission.
Detailed balance. For the transition between state Γ and Γ ′ detailed balance corresponds to [10,13]: where α(Γ → Γ ′ ) is the probability of generating a Γ → Γ ′ move attempt (trial), p acc (Γ → Γ ′ ) is the probability of accepting the move, and P (Γ) = ρ(Γ)d nv(Γ) ⃗ x is the equilibrium probability of finding a system in a voxel of volume d nv(Γ) ⃗ x.
Next, we use Eq. (55) to define the acceptance criteria for each MC move. The acceptance criteria are derived in detail for the wedge fusion/fission move; the steps to follow are the same for all other moves.

A. Vertex displacement
In this move, a vertex is randomly selected, a random uniform displacement is drawn, and the vertex is displaced to its new position according to: with d max the maximum displacement. The move is accepted with a probability p acc = exp(−∆E/k B T ) where ∆E is the (bending plus stretching) energy change due to the displacement. The parameter d max can be adjusted during a burn-in period to optimize convergence to equilibrium. Generally optimal values are on the order of the typical length scale of thermal fluctuations dictated by the elastic energy, leading to acceptance probabilities on the order of 50%. In our simulations typical values are between d max = [0.01l 0 , 0.1l 0 ]. The vertex displacement move is illustrated in Fig. §20: the number of subunits n s , number of vertices n v , number of vertex bonds n VB and number of bonds n b remains unchanged during this move. Throughout the article, we define a Monte Carlo sweep as the number of steps required for the system to have undergone (on average) n v vertex moves, with n v the number of vertices in the structure. One sweep is one time step in our simulation. The other kinetic moves are attempted less frequently than the vertex moves. Denoting f v as the vertex move attempt frequency, monomer insertion/deletion moves are attempted with frequency f i = 10 −3 f v , while wedge fission/fusion and crack fission/fusion moves are attempted with frequency f wf = f cf = 1˛0 −4 f v . Edge fission/fusion moves are attempted as specified in the main text, controlled by the parameter f fusion . The attempt frequency of other moves is fixed throughout all the simulations reported in this work. While the acceptance probabilities are different for different types of moves, they all guarantee detailed balance as shown below.

Simple insertion
In this move, an edge is randomly selected from the set of all boundary edges, where a new subunit will be attached. The number of such boundary edges is n e . Subunits can be inserted in n r different rotations, where n r is the number of distinct rotational states for a subunit which has one edge aligned with the edge of a neighboring subunit. For our triangular subunits with three distinct edge types, n r = 3. In our algorithm, during insertion of a subunit its rotational state is chosen randomly from the set of three possibilities. If the aligned edge is not complementary to the type of the boundary edge, then the move is rejected. In this work, the two edges must be of the same type to be complementary.

27
The positions of two of the new subunit's vertices (those at either end of the edge being bound) are set equal to the positions of the corresponding vertices of the boundary edge to which it is binding. The third vertex position is randomly chosen from within a volume v add centered at the equilibrium position of the new vertex.
Thus, the attempt probability for a simple insertion is given by: = n e f i τ n r × 1 n e n r (v add /d⃗ x) .
Then, applying Eq. (55) and the attempt probability for the reverse move (simple deletion, presented next, Eq. (61)), the acceptance probabilities for a simple insertion is ∆E i→j is the energy change due to the move and includes the stretching energy of the newly inserted subunit, its bending energy along the shared edge, and the binding energy due to the creation of an extra bond. During this move, one new (edge) bond and two new vertex bonds are created; i.e. n b → n b + 1 and n VB → n VB + 2. Moreover, the number of vertices in the structure increases by one, n v → n v + 1.

Simple removal
The reverse move to simple insertion is simple removal. Subunits that can be deleted with this move are those with two boundary edges. The number of simply removable subunits is n sr . One of these is selected randomly, so the attempt probability is and, using Eq. (55) and Eq. (59), the acceptance probability is During this move, the structure loses one (edge) bond and two vertex bonds; n b → n b − 1 and n VB → n VB − 2. The number of vertices in the structure decreases by one, n v → n v − 1.
If there are multiple species with chemical potentials µ k , detailed balance must be satisfied for each species, individually. Moreover, each species can have different insertion rates f k i . To keep α < 1, we ensure that the insertion rate f i constrained by n e f i τ n r < 1 (63) n sr f i τ < 1 In equilibrium, one can use adaptive rates, i.e. reduce k i on the run if the above condition is not satisfied. In that case, sampling is not taken for the ensuing several time steps. Alternatively, the rates may be set to a low enough value from the beginning and only tested on the run to ensure that the α < 1 condition is satisfied. This latter technique is appropriate for dynamical runs as it keeps the rates constant throughout the simulation.
Moreover, we must ensure that v add is large enough so that the vertex does not leave the v add volume during structural relaxation moves; otherwise the insertion/deletion moves would not be reversible and the detailed balance would be violated. For a better convergence, one could choose a gaussian distribution N (⃗ r) for the position of the new vertex instead of a uniform distribution 1/v add . In this case, this distribution has to be accounted for in the acceptance probabilities p acc (i → j) and p acc (j → i).

Wedge insertion
Wedges are positions in the structure where a triangle can be inserted via attaching to two edges (Fig.  §22). In a wedge move, we pick randomly from the set of available wedge positions in the structure, and pick a random orientation for the new subunit. Denoting the number of wedge positions in a given structure as n w , the attempt probability for a wedge move is In contrast to the simple insertion move, there is no need for random vertex displacement in a wedge move because all three vertices of the new subunit are fixed by the three vertices of the wedge position. Combining Eq. (65) and the attempt probability for wedge removal (Eq. (67)), The acceptance probability for a wedge insertion is During a wedge insertion, two edge bonds and three vertex bonds are created; i.e., n b → n b + 2 and n VB → n VB + 3, but the number of vertices is unchanged, n v → n v . ∆E i→j includes the binding energy of the two newly formed bonds, the two bending energies along the two newly bound edges and the stretching energy of the newly inserted subunit.

Wedge removal
The reverse move of wedge insertion is wedge removal. In a wedge removal, we randomly choose one of the removable wedges from a given structure. With the number of removable wedges as n wr , the attempt probability is Using Eq. (65), the acceptance probability for a wedge removal is then We have the following constraints on rates f i for wedge insertion/removal: As for simple insertion and removal, in the case of multiple species, detailed balance is satisfied for each species separately for wedge insertion/removal.

Wedge fusion
In this move, a fusable wedge is closed, without inserting a new subunit (Fig. §23). That is, the two vertices on either side of the wedge opening are merged into a single vertex. Fusable wedges are vertex pairs that i) form a wedge (as in the case of wedge insertion) and ii) are within a separation distance of l fuse . Denoting the number of fusible wedge positions as n w , in each MC step, a wedge fusion is attempted with probability n w f wf τ , where f wf is an adjustable parameter controlling the relative probability of attempting wedge fusion. Then, a wedge position is selected randomly from the set of all n w fusible wedges. The attempt probability is thus Using Eqs. (71) and (73), the acceptance probability for fusion moves is where v fuse = (4π/3)(l fuse /2) 3 is the volume of a sphere with diameter l fuse , and ∆E i→j is the energy change due to the fusion, including changes in bending, stretching, and binding energies. A fusion move increases the number of edge bonds and vertex bonds by one, n b → n b + 1 and n VB → n VB + 1; the factor of v a appears in Eq. (72) to account for the latter.

Wedge fission
Wedge fission, in which a wedge is opened, is the reverse of the wedge fusion move. Fissionable edges are those edges that can be split along their boundary vertex to obtain a wedge. Denoting the number of such edges as n f , the probability of attempting a wedge fission move during an MC step is n f f wf τ . If a fission move is attempted, then an edge is selected randomly from the n f fissionable edges. The position of one of the new vertices is selected randomly within the sphere of volume v fuse centered at the original position of the merged vertices, and the other new vertex is placed in the opposite direction from the original position, at the same distance. Thus, the attempt generation probability is and the acceptance probability is We verify that detailed balance holds between wedge fusion and fission as follows. There are two cases to consider: In this case, p acc (i → j) = 1 and p acc (j → i) = (v fuse /v a ) exp(−∆E j→i /k B T ). Then and Using: ∆E j→i = E i −E j , n s,i = n s,j (because the move leaves the subunit number unchanged), n VB,i = n VB,j −1 (one vertex bond is broken upon fission) and n v,i = n v,j + 1 (an extra vertex is being born upon fission), we see that the two are equal and detailed balance holds.
In this case, p acc (i → j) = (v a /v fuse ) exp(−∆E i→j /k B T ) and p acc (j → i) = 1. Then and P j × α(j → i) × p acc (j → i) = 1 Z Ω v n VB,j a exp[−(E j − µn s,j )/k B T ] 1 λ 9ns,j × d nv,j ⃗ x (80) Using again ∆E j→i = E i − E j , n s,i = n s,j , n V B,i = n VB,j − 1 and n v,i = n v,j + 1, detailed balance holds. Note that detailed balance is satisfied regardless of the values of f wf τ or v fuse , but as with all of the move frequencies these parameters can be optimized during burn-in to accelerate convergence to the equilibrium distribution P (i). In our simulations, we find that the optimal value of v fuse is on the order of the optimal value of d max for analogous reasons: if v fuse is too small there will be very few vertex pairs identified as fusable, so n w will be low. If v fuse is too large, there will be many fusion candidates but most fusion attempts will be rejected due to the large elastic energy change necessary for the merging deformation.
Most importantly, we note the constraint on the parameters f wf τ to ensure that generation probabilities do not become larger than unity. Because each attempt is generated as a three step process, using three probabilities, one has to ensure that all those probabilities are less than 1. Specifically, n w f wf τ < 1 (82) n f f wf τ < 1.
E. Crack fusion / fission

Crack fusion
Crack fusion closes a crack within the structure; i.e., two adjacent pairs of edges are merged (Fig. §24). Cracks are identified as 4-edge-length holes inside the structure. If the vertices of the hole are labeled A, B, C, D then the polygon ABCD forms a closed loop (see Fig. §24). The crack can be closed by either merging vertices A and C (and correspondingly edges CD to DA and AB to BC) or by merging vertices B and D (and correspondingly edges AD to AB and CD to CB). Each 4-edge-length loop thus defines two potential fusable cracks. However, an additional condition for a crack to be fusable is that its merging vertices must be within a distance l fuse (A and C or D and B in this example). In this work, we have set the crack fusion volume to be the same as that for wedge fusion to reduce the number of parameters, but it is not necessary that they be the same.and the acceptance probability is There are two edge bonds and one vertex bond formed during a crack fusion.

Crack fission
The reverse move for crack fusion is crack fission. With the number of potential cracks as n cf : As for the case of wedge fusion/fission, the crack fusion attempt frequency parameter f cf is constrained by the conditions maintaining probabilities smaller than unity: F. Edge fusion / fission

Edge fusion
During this move two non-neighbor edges are fused (Fig. §25). Fusable edges are non-neighboring edge pairs whose corresponding vertices are within a separation distance l fuse . Since edges are directed, they can only fuse such that, after fusion, they point in the opposite direction. Assuming the edges to be fused are A → B and C → D (see Fig. §25), vertex A will merge into vertex D and vertex B will merge into vertex C. Edges are counted as fusable if A is within a distance l fuse to D and B is also within a distance l fuse to C. The attempt probability is analogous to that for wedge and crack fusion/fission, with n e the number of fusable edges and f fusion the edge fusion frequency parameter. The acceptance probability is During edge fusion, one edge bond and two vertex bonds are created.

Edge fission
Edge fission is the reverse move of edge fusion. n ef is the number of breakable edges, that is, those edges that have both vertices on the boundary and which would not result in breaking the structure apart. α(j → i) = n ef f fusion τ × 1 n ef (v fuse /d⃗ x) 2 (92) The factor 1/(v fuse ) 2 arises because we must select a random position for each pair of vertices, independently. The acceptance probability is then To maintain probabilities within unity, the edge fusion frequency parameter f fusion is constrained by n e f fusion τ < 1 (94) n ef f fusion τ < 1.