Tiling a tubule: how increasing complexity improves the yield of self-limited assembly

We consider the assembly of cylindrical tubules from triangular subunits, as illustrated in figure 1(a). The subunits are flat equilateral triangles, which bind edge to edge. Since a tubule is an object with non-zero curvature, that curvature must be encoded in some way into the subunit. We accomplish this goal by considering subunits that form a specified dihedral angle between their edges when they bind. Therefore, each triangle has associated with it three dihedral angles, one for each side. We consider a case where the interactions are specific and all three sides of the triangle are distinguishable. The simplest set of interactions that prescribe the assembly of a tubule is side 1 binds to side 1, side 2 to side 2, and side 3 to side 3. The specificity of these interactions is important in order to preserve the same local curvature everywhere across the assembly and to ensure that the pattern is deterministic. A schematic of such specific interactions and dihedral angles are shown in figure 1(b).

Zoom In Zoom Out Reset image size

Any tubule formed from triangular subunits can be classified uniquely by a pair of indices. Consider a tubule as a rolled-up triangular lattice or a sheet with periodic boundary conditions along two parallel lines (figure 1(c)). Within this conceptual framework, a tubule can be constructed by choosing any two points on this plane and enforcing them to be periodic with one another. In other words, these two points overlap when the sheet is rolled up into a tubule, as in figure 1(c). To classify the tubule that forms, we count how many lattice edges need to be traversed between the periodic vertices. Taking the unit vectors of this tiling to be m and n, as shown in figure 1(c), we can go between the vertices in m and n steps in the respective lattice directions, giving a total displacement of w = m m + n n. The corresponding tubule would be identified as (m, n) and have a width of $w=\sqrt{{m}^{2}+{n}^{2}+mn}$. Here, the width w refers to the width of the strip on the triangular lattice (figure 1(c)), which corresponds to the circumference of the closed tubule.

For any given tubule, there is a unique set of dihedral angles between adjacent triangular edges that will yield that desired structure. These values can be obtained by constraining the vertices of the lattice to lie on the surface of a cylinder and finding the angle between adjacent faces. Note that in this construction, the specified dihedral angles are constant along a given lattice direction, denoted by indices 1, 2, 3 in figure 1(c). This constraint imposes an orientational order for the triangular subunits, which we will return to later.

We explore the assembly outcomes using grand canonical Monte Carlo simulations. Specifically, we use the event-driven Monte Carlo algorithm developed by Rotskoff and Geissler [39] and Li et al [40], and adapted to tubules in reference [41], in which an assembled structure exchanges subunits with a bath at fixed chemical potential. Each triangular subunit is modeled by three vertices and three straight edges connecting the vertices. The Hamiltonian of the system is given by,

$$\begin{align}\hfill H& =\frac{1}{2}\sum\limits _{\text{Bound}\;\text{Edges}}{E}_{\mathrm{B}}+\sum\limits _{\text{Edges}}\frac{1}{2}S{(l-{l}_{0})}^{2}\hfill \\ \hfill & \quad +\sum\limits _{\text{Adjacent}\;\text{Faces}}\frac{1}{4}B{(\theta -{\theta }_{0})}^{2}.\hfill \end{align} \tag{ 1 }$$

The first term is the binding energy. E_B is the energy difference between a pair of bound and unbound edges, and is set to be the same for all favorable interactions. The second term is the stretching energy. S is the stretching modulus of the edge, l is the instantaneous length of the edge, and l₀ is the stress-free length of the edge. The third term is the bending energy. B is the bending modulus of the edge pairs and is again set to be the same for all edges. θ is the instantaneous dihedral angle between two subunits and θ₀ is the preferred dihedral angle for a given tubule structure and type of edge-pair. See supporting information (https://stacks.iop.org/JPCM/34/134003/mmedia) (SI) section 1 for a detailed description of our computational methods.

For a given set of input parameters—the bending modulus, the lattice numbers m and n of the target, and the number of unique species—we perform one thousand independent simulations and analyze the distribution of tubule types that form. We prescribe the equilibrium dihedral angles, θ₀, to favor (m, 0) and enforce the binding rules specified previously. We tune the various energies in the Hamiltonian to keep the supersaturation low enough that the structure can nucleate, grow, and close near to equilibrium in a reasonable time scale. The simulation starts with a single subunit and grows by adding subunits onto the pre-existing structure, ending once the tubule has a length roughly three times its circumference. Finally, we determine the tubule type of the end state and compute the distribution of tubule structures for each condition.

We only examine defect-free tubules in the following analysis. We consider a tubule to be defect-free if it has the same tubule structure, characterized by the same pair of indices (m, n) along its entire length. While it is possible for defective tubules to form in our simulations, we find that the likelihood of forming a defective tubule does not depend on the number of subunit species used for assembly. However, the fraction of defect-free tubules does vary between 95%–63% depending on the targeted tubule width (see SI figure S3). We hypothesize that defect creation is related to the degree of supersaturation and the kinetics of monomer addition near closure. This hypothesis is consistent with the type of defective tubules that we observe, which tend to have a central region that does not close properly, resulting in a tubule with a heterogeneous tubule type along its length. See SI section 2 for a description of the algorithm that we use to determine the tubule type and for details on the defect rates.

The results of our simulations show that assembly yields a broad distribution of tubule types. Figure 2(a) shows an example of such a distribution. The bevel angles for this example were chosen to target a (10, 0) structure. While we do see that the majority of tubules formed the target state, there is a broad distribution of adjacent states that have also formed, extending one to two lattice steps in both m and n directions. The probability for assembling off-target states falls off further from the target state as expected, since the dihedral angle differences, and hence the elastic energy cost, become more significant.

Zoom In Zoom Out Reset image size

To understand the origin of this distribution, we consider the process by which assembly occurs and hypothesize that the distribution is dictated by the mechanism of tubule closure. Assembly begins when subunits come together to nucleate a curved sheet. Once the extent of the sheet is large enough, it can close to form a short tubule, which then grows by the addition of subunits to its free ends. In the pre-closure state, thermal fluctuations can cause the sheet to close at larger or smaller widths around the target state, sampling the energies of different closed states. We hypothesize that once the tubule closes, it is highly unlikely to open back into a sheet because doing so would require rupturing multiple edge–edge interactions simultaneously. Therefore, once the assembly has gone through closure and begun to extend, the tubule type is essentially fixed, even if it is not the global free-energy minimum. While the energy difference between off-target states and the target structure continues to increase as the tubule grows longer, without the ability to open again, the system cannot equilibrate to the designed global energy minimum. Indeed, in simulations, we find that the tubule type does not change after a structure has closed (see figure S2(c)).

In addition to running simulations, we also estimate the distribution of structures that may arise during assembly from calculations of the elastic energy at closure. The penalty for forming off-target states comes from the elastic energy cost of not abiding by the prescribed dihedral angles. Recall that the bending energy, E_θ , for any binding site is given by ${E}_{\theta }=\frac{1}{2}B{(\theta -{\theta }_{0})}^{2}$. For the case of our triangular subunits, the elastic energy per subunit is half the sum over all three binding sites. When closure happens, it occurs for a finite number of subunits in the assembly, N_C, all of which contribute to the elastic energy. The full elastic energy cost at closure is therefore

$$\begin{equation}{E}_{\theta }=\frac{1}{4}{N}_{\mathrm{C}}\sum\limits _{i=1,2,3}B{({\theta }_{i}-{\theta }_{i,0})}^{2},\end{equation} \tag{ 2 }$$

where the extra factor of $\frac{1}{2}$ comes from the fact that a single binding site is shared between two subunits. To estimate the closure size, we assume that subunit addition occurs isotropically, forming a circular disk with a diameter of the tubule width, w. Taking the ratio of the disk area to the area of a subunit gives ${N}_{\mathrm{C}}=\pi {w}^{2}/\sqrt{3}$. To get a more accurate accounting of the free-energy cost of forming misassembled structures, one would need to include the surface energies and the various entropic contributions, but we find that just considering the bulk elastic terms reproduces the scaling that we find in simulation.

We see the elastic energy at closure is lowest at the target state and increases by a few k_B T at the nearest off-target states, suggesting that these states are likely accessible in a system at finite temperature (figure 2(b)). To get an estimate of the tubule type distribution, we compute the probability of a structure s according to $P(s)={\text{e}}^{-\alpha E(s)/{k}_{\mathrm{B}}T}/Z$, where Z is the partition function and α is a factor that adjusts the height of the free-energy barrier between the pre-closure and the closed states to account for missing contributions to the closure rate (e.g. deviations of pre-closed tubules from the assumed circular shape and other entropic effects). Supplementary simulations found a value of α = 0.3 and are described in SI section 3. Without this additional factor the energetics predict narrower distributions. Figure 2(c) shows the probabilities according to the energy landscape in figure 2(b), which bears a similar shape and extent to the simulated distribution in figure 2(a).

Estimates of the number of accessible states follow power-law scalings with the tubule width and bending rigidity. From the energetics calculation, we make an estimate of the number of different structures, N_types, by counting how many states have a probability greater than 0.25%. From simulations, we look at the final tubule structures that form in four hundred different runs. We find that the number of different structures that assemble is sensitive to both the target width and the bending rigidity: the number of accessible states increases with increasing width and decreases with increasing bending rigidity (figure 2(d)). We find that these same data collapse to a single curve when rescaled by $w/\sqrt{B}$ (figure 2(e)). In the inset of figure 2(e) we show the same data but using α = 1.0 for the energetics, finding that the number of types predicted is lower than what we see in simulation. This implies that α is capturing some aspects of the kinetic processes during closure.

The scaling for the breadth of the tubule type distribution comes from a balance between the closure size and the fluctuations of the curvature of the sheet before closure. We consider the Helfrich energy of a curved sheet [42], $E=\frac{1}{2}BA{({\Delta}\kappa )}^{2}$, where A is the area of the sheet and Δκ is the deviation of its curvature away from the ideal curvature. We can approximate this curvature as (Δw/w²), where Δw is the fluctuations of the width and w² is the area of the sheet. Looking at the size of these fluctuations on an energy scale of k_B T shows that ${\Delta}w\sim w/\sqrt{B}$. Therefore, the scaling that we find in figure 2(e), N_types ∼ Δw², arises from the fact that thermal fluctuations populate a region of vertices around the target vertex with an area of Δw².

These results illustrate a fundamental hurdle for self-limited assembly: in a thermal system, it is difficult to achieve specificity of a target state when the self-limited length scale is large compared to the subunit size. Small fluctuations of the dihedral angles between subunits become amplified as the number of subunits in the self-limited length scale increases. Even though the rigidity of individual dihedral angles may remain the same, the fluctuations of the self-limited length scale grow proportionally to the self-limited lengthscale itself, as we saw from the Helfrich energy. Compounding this effect, the process of irreversible closure prevents the assembly from visiting different states at later times to further relax. Therefore, the breadth of the distribution is driven by a kinetic process, which yields a larger variety of states than would be expected in equilibrium. It is this bottleneck to high-yield of the specific target that must be engineered around, either by making ever stiffer subunits, circumventing closure-control by seeding nucleation with specific geometries, or by altering the energy landscape near closure. Here, we will explore this last direction by considering how multiple types of interacting subunits limit the accessible states at closure, and thereby prune off-target states from the distribution of final structures.

To proceed, we extend our framework to allow for multiple species of triangles, where a species of triangle is a subunit with a distinct set of specific interactions encoded in its edges. The first task is to identify allowed patterns of multiple species that still satisfy the requirements imposed by the tubule geometry. This challenge involves finding periodic patterns in a triangular lattice with multiple colors of triangles. There are many ways that one can imagine periodically coloring a triangular lattice, but not all of these tilings will necessarily preserve the physical rules for tubule assembly. Here, we introduce three rules for multispecies tubule assembly. First, to have a deterministic assembly, there need to be unique interactions between subunits: each side can only bind to one other side within the mixture of many subunit species. In the context of a tiling, this constraint means that once the colors of tiles adjacent to any subunit have been specified, there cannot be different neighbors at any other location in the tiling. Second, to impose constant dihedral angles along specific lattice directions, only rotations of 180 degrees for a given subunit type are allowed throughout the tiling (figure 3(a)). Third, if we view all interactions between subunits as creating allowed dimers with specific orientations, then we also must be able to construct closed vertices on the plane using either three or six of our specified interactions (figure 3(b)); this last constraint will force the system to assemble into a deterministic pattern that is the same everywhere in the tiling.

Zoom In Zoom Out Reset image size

The concept of an interaction matrix helps to make these rules for subunit interactions more concrete. Elements of an interaction matrix will either be non-zero, prescribing an allowed attraction between specific sides of triangles, or zero, meaning no binding is allowed (figure 3(c)). In SI section 4, we describe how the restrictions mentioned above can be translated to allowed constructions of interaction matrices. Following this construction, we enumerated all allowed tilings of a triangular sheet for up to ten different species of subunits. Figure 3 shows one allowed tiling for four species, which are illustrated as four different colors.

Looking at the different tilings that can be wrapped into tubules, we notice that some of them will be more useful for restricting neighboring states than others. The most important feature of a tiling in this regard is its set of primitive vectors (figure 5). For any tiling we can define a pair of vectors, a₁ and a₂, that go between similar vertices in the tiling, have a minimal length, and are maximally orthogonal. These primitive vectors can then be used to identify which tubule types can be formed from a certain pattern. Recall that for tubule type (m, n) there is an associated displacement vector w = m m + n n. If w can be made from a linear combination of a₁ and a₂, then (m, n) is an allowed state for that tiling (figure 4(a)). See the SI for patterns and interaction matrices.

Zoom In Zoom Out Reset image size

When the two primitive vectors are large, the similar vertices are farther apart, leading to greater distances between allowed states. Furthermore, the closer in magnitude the primitive vectors are to one another, the more uniform the restriction of states will be around the target state. If there is a difference in the length of the primitive vectors then the tiling will be anisotropic and there will be an additional orientation dependence of the restricted states with respect to the orientation of the triangles. Two patterns made from three species of triangles are shown in figures 4(b) and (c). From these two patterns we can see that figure 4(b) shows isotropic distances to nearby vertices while figure 4(c) shows anisotropic distances. Depending upon how w is aligned with respect to the anisotropic primitive vectors, the assembly outcome changes. For instance, if w is aligned along a shorter primitive vector, the density of states in that direction will be larger and will result in larger fluctuations of the widths of tubules that form. In the other case, when w is aligned along the longer primitive vector, there will be larger fluctuations in the chirality of tubules that form. To exemplify some of these points, we will look in detail at the patterns formed from two colors.

To illustrate how multispecies tilings can limit the number of accessible off-target states, we first consider the allowed tilings composed of two species. We find that there are only th-ree unique patterns that satisfy the restrictions for assembling tubules (figure 5), with each pattern having a unique interaction matrix required for assembly. For each pattern, we also compute the allowed tubule geometries. Interestingly, we find that one of the three patterns (figure 5(a)) does not give additional restrictions compared to a single-species tiling. In contrast, the other two tilings reduce the available states by half, albeit with the same restrictions as one another (figure 5(b)). Note that even though some patterns have the same restrictions, all three patterns have different interaction matrices, and thus there is only one designed ground-state tiling for each pattern.

Zoom In Zoom Out Reset image size

As mentioned above, we rationalize the observation that some tilings restrict the allowed states while others do not by looking at the primitive vectors of the tilings. For the pattern in figure 5(a), we see that every vertex of the tiling is the same, meaning that all points in m, n-space are allowed tubule types. This occurs because the length of its primitive vectors both match the subunit length. In contrast, the two patterns in figure 5(b) have two distinct vertices that cannot overlap with one another on a tubule. Because the two vertex types appear with equal frequency, each vertex can only bind with half of the total vertices in the pattern, thereby reducing the number of accessible states by a factor of two. Additionally, since the primitive vectors for the patterns in figure 5(b) are not equal length, the restrictions imposed on the allowed states are anisotropic. For this pattern, any value of m is allowed, while only even values of n are permitted. Depending on whether or not one wants to more tightly restrict the chirality of the structure or the available widths, one can adjust the orientation of the target displacement vector, w, with respect to the two primitive vectors.

A subtle point about the two patterns in figure 5(b) is that even though they produce the same restrictions on allowed states, they require different numbers of specific interactions to encode their patterns. For example, we see that there are five and four unique matrix elements for the top and bottom patterns in figure 5(b), respectively. We hypothesize that reducing the number of interaction types to restrict a greater number of states will be an important design criteria for multiple-species experiments, since there will inevitably be a finite capacity for specific interactions beyond which undesired crosstalk between edges becomes non-negligible.

Now that we understand how a tiling with multiple species can change the number of accessible states, we explore how to design a system that targets a single state with a high yield. As we saw in figure 2(d), the number of off-target states, N_types, increases as the area of fluctuations around the target vertex, Δw², increases. To compensate for this effect, we expect that the required primitive vector length of a tiling should be comparable to the fluctuations of the self-limited length scale to achieve high yield of the target state. Figure 6(a) shows how a distribution of tubule types relates to the spatial distribution of allowed states around the target vertex. We see that the accessible vertices are nearly uniformly distributed around the target vertex, with a slightly larger variation in the w direction. Therefore, we expect that tilings that have primitive vectors of similar lengths will eliminate off-target states most effectively. Figure 6(b) shows examples of two tilings and their respective excluded regions. The first pattern shows an isotropic tiling with the nearest accessible vertices shown connected by the dotted line. The second pattern shows an extreme case in which the excluded region is highly anisotropic. In the short direction, smaller fluctuations would be needed to access another allowed vertex. Going forward, we restrict ourselves to the use of isotropic or near-isotropic tilings; the patterns used for the different numbers of species are shown in figure 6(c).

Zoom In Zoom Out Reset image size

To see how increasing the number of species, N_colors affects the assembly specificity, we perform both simulations and energy calculations as before. Specifically, we calculate the number of accessible states for tubules of different target widths, w_target, and bending rigidity, which are assembled from different numbers of subunit species.

We find that the assembly yield of the target state increases as the number of species increases. Figure 6(d) shows the number of tubule types that form, N_types (left), as well as the probability of forming the target state, P_target (right), as we change the number of species, rescaled by w²/B. Recalling that the number of tubule types for a single species, N_types, grows as w²/B (figure 2(e)), we consider the quantity w²/B as the area of the fluctuations around the target vertex. For isotropic tilings the primitive vectors have a length that grows as $\sqrt{{N}_{\text{colors}}}$, meaning the excluded area around a target vertex will grow linearly with the number of species. When N_colors B/w² ≈ 1 the area of the disallowed region and the area of fluctuations become comparable. At this point we observe full specificity in both N_types and P_target. Both the energetics and the simulation show the same scaling.

To clearly illustrate how adding additional species impacts the tubule distributions, we show examples of three different target widths for N_colors = 1, 3, 9 in figure 6(e). As the number of species is increased there is a reduction in the number of types of structures that form, with a corresponding increase in the fraction that form the target state. Note that even though the specificity of the target is increased, the extent of the off-target states is not impacted, as seen most clearly in the three-species case. These results further illustrate that the minimum number of subunit species to achieve full specificity is proportional to the size of the fluctuations of the system.