Self-assembly of two-dimensional binary quasicrystals: a possible route to a DNA quasicrystal

Patchy-particle models have been used extensively to study a wide range of behaviours in computer simulations [30, 31, 65, 69–83], including self-assembly and crystallisation, and represent one of the simplest types of 'toy model' which can account for the complexity of behaviour seen in experiments on a number of colloidal systems.

In our simulations, we use the Metropolis Monte Carlo scheme [84] with volume moves [85, 86] and periodic boundary conditions. In simulations with multiple particle types, we furthermore allow moves in which two particles of distinct types are exchanged with one another in order to help facilitate equilibration.

We model particles with attractive patches using a simple angular modulation of the attractive part of the Lennard-Jones potential,

$$\begin{eqnarray}V\left(\boldsymbol{r}_{{ij}},{{\varphi}_{i}},{{\varphi}_{j}}\right)=\left\{\begin{array}{*{35}{l}} {{V}^{\text{LJ}}}\left({{r}_{ij}}\right) & {{r}_{ij}}<{{\sigma}_{\text{LJ}}}, \\ {{V}^{\text{LJ}}}\left({{r}_{ij}}\right){{V}^{\text{A}}}\left({{\widehat{\boldsymbol{r}}}_{ij}},{{\varphi}_{i}},{{\varphi}_{j}}\right) & {{\sigma}_{\text{LJ}}}\leqslant {{r}_{ij}}, \end{array}\right. \end{eqnarray} \tag{ 1 }$$

where $\boldsymbol{r}_{{ij}}$ is the interparticle vector connecting the particles i and j, r_ij is the magnitude of this vector, and ${{\varphi}_{i}}$ and ${{\varphi}_{j}}$ are the orientations of the particles i and j, respectively. The Lennard-Jones potential is given by

$$\begin{eqnarray}&&{{V}^{\text{LJ}}}\left({{r}_{ij}}\right)=4\varepsilon \left[{{\left(\frac{{{\sigma}_{\text{LJ}}}}{{{r}_{ij}}}\right)}^{12}}-{{\left(\frac{{{\sigma}_{\text{LJ}}}}{{{r}_{ij}}}\right)}^{6}}\right],\end{eqnarray} \tag{ 2 }$$

and we use a potential cutoff of ${{r}_{\text{cut}}}=3{{\sigma}_{\text{LJ}}}$ and shift the potential so that it equals zero at r_cut. The angular modulation term in the potential is given by a product of gaussian functions,

$$\begin{eqnarray}&&{{V}^{\text{A}}}\left({{\widehat{\boldsymbol{r}}}_{ij}},\,{{\varphi}_{i}},\,{{\varphi}_{j}}\right)=\underset{k,\,l}{{\max}}\,\,\left\{\exp \left[\frac{-\theta _{kij}^{2}}{2\,\sigma _{\text{pw}}^{2}}\right]\exp \left[\frac{-\theta _{lji}^{2}}{2\,\sigma _{\text{pw}}^{2}}\right]\right\},\end{eqnarray} \tag{ 3 }$$

where ${{\sigma}_{\text{pw}}}$ is a parameter reflecting the patch 'width' and ${{\theta}_{kij}}$ is the angle between the patch vector of patch k on particle i and the interparticle vector $\boldsymbol{r}_{{ij}}$. The product of gaussian functions is evaluated over all possible patch pairs $\left\{k,\,l\right\}$, and the optimum combination is chosen, i.e. a pair of particles can only interact via that pair of patches which is most energetically favourable.

In simulations with multiple patch types, the angular modulation of equation (3) is modified to include a prefactor that depends on the interaction matrix of the two patches k and l considered. In all simulations considered here, the matrix elements of this matrix are either zero or unity, i.e. all patches that interact have the same strength, and patches that do not interact do not contribute at all to the energy in the attractive part of the Lennard-Jones potential. The most energetically favourable pair of patches is still chosen in the computation of the angular modulation in equation (3).

In order to characterise the structures we observe in our simulations, we classify each particle according to its nearest-neighbour environment [30, 31]. To do this, we determine the neighbours of each particle, using a simple spherical cutoff of $1.38\,{{\sigma}_{\text{LJ}}}$ [30], and then determine how many neighbours each neighbouring particle shares with the particle we are classifying. We classify particles into three distinct types of environment, σ, H and Z, with common neighbour signatures of {21111}, {22110} and {222222} respectively, as illustrated in figure 1. This labelling corresponds to the equivalent Frank–Kasper phases [87, 88]. In all simulation snapshots shown in the next section, particles are coloured according to their classification following the colour-coding shown in figure 1, namely cyan (σ), violet (H) and red (Z), with particles whose environments give any common neighbour signature not listed above depicted in green. In simulations with multiple particle types, the base particle colour corresponding to figure 1 is mixed with varying amounts of black for particles of types A and C (as defined below) in order to help to distinguish them.

Zoom In Zoom Out Reset image size

In all simulations reported here, the patch width was chosen to be ${{\sigma}_{\text{pw}}}=0.3$ rad. For the patchy potential introduced in equation (1), this is quite a narrow patch width [65], making interactions very angularly dependent. Such a narrow patch width allows us to account for the relatively highly directional nature of DNA multi-arm motif structures and prevent competition between environments of different valencies. However, for precisely the same reason, this choice of patch width leaves us in a region of parameter space where, for pentavalent particles we considered previously [30, 31], the quasicrystal was not thermodynamically stable, as the patches are so narrow that it is not possible for six neighbours to bond competitively: instead, the σ phase was stable for pentavalent particles under these conditions [31].

We note that in the quasicrystals we studied previously [30, 31], the most common structural feature was a series of edge-sharing dodecagonal motifs; one of the most common such motifs is shown in figure 2(a). In the 'unit cell' of the approximant crystal corresponding to this motif, also illustrated in figure 2(a), there are two particles in a hexagonal (Z) environment and 24 particles in a σ environment. Since the former correspond to particles having six neighbours and the latter to only five, we can surmise that in order to achieve our goal of assembling quasicrystals with particles with a narrow patch width, including both hexa- and pentavalent patchy particles in a simulation box in a ratio of 1:12 would be a sensible choice.

Zoom In Zoom Out Reset image size

We have run simulations with a mixture of hexa- and pentavalent patches in this ratio (figure 3). A quasicrystal phase forms spontaneously when the liquid phase is cooled (figure 3(c)); its quasicrystallinity is confirmed by the diffraction pattern shown in figure 3(c), which clearly exhibits dodecagonal symmetry: since this dodecagonal order is coherent through the whole of the simulation box, this is a single quasicrystal. There are several features of this quasicrystalline configuration worth noting. The large majority of hexavalent patches are at the centres of dodecagons, and many of these dodecagons are arranged locally in the motif of figure 2(a). However, a common alternative motif for a quasicrystalline approximant is depicted in figure 2(b), based on overlapping, rather than edge-sharing, dodecagonal motifs. We can see that there are sections of the structure in figure 3(c) which are locally like both the edge-sharing and the overlapping approximants, both of which form a triangular lattice with longer and shorter distances between the 'vertices' of the lattice. However, it is worth noting that there are also a number of other motifs of dodecagon centres, such as rectangular and isosceles triangular ones, which feature in figure 3(c).

Zoom In Zoom Out Reset image size

As we decrease the temperature further still, the quasicrystalline approximant is expected to become the stable phase, as the configurational entropy of the quasicrystal becomes relatively less important compared to the enthalpic stability of the approximant [31]. In order to check that the quasicrystal is thermodynamically stable at intermediate temperatures, rather than simply a kinetic product, we have computed the free energies of the competing phases. To find the free energy of the edge-sharing approximant (figure 3(b)), we used Frenkel–Ladd integration [42] at ${{k}_{\text{B}}}T/\varepsilon =0.1$ and $\sigma _{\text{LJ}}^{2}\beta P=0.5$ (where $\beta =1/{{k}_{\text{B}}}T$), and, for consistency checking, also at $\sigma _{\text{LJ}}^{2}\beta P=1.5$, and then integrated this free energy along iso-$\left(\beta P\right)$ curves [43]. We found the free energy of the fluid phase by integrating from the ideal gas, bearing in mind that it is a multicomponent gas. The free energy of the quasicrystal was then set by equating it to the free energy of the fluid at the point at which the quasicrystal and the fluid phase are at equilibrium; we determined this condition by direct-coexistence simulations [31]. By finding where the free-energy curves of the approximant and the quasicrystal cross, we can obtain the relevant coexistence curve between the two phases. Finally, at very high pressures, the hexagonal (Z) plastic crystal phase dominates because its density is larger [31], whether or not all neighbour–neighbour interactions can be satisfied. However, bearing in mind that the corresponding DNA systems are very dilute solutions, if we work at reasonable pressures, the Z phase does not need to be considered further.

At such reasonable pressures (i.e. for $\sigma _{\text{LJ}}^{2}\beta P\lesssim 10$), the quasicrystal is thermodynamically stable relative to the fluid at temperatures below about ${{k}_{\text{B}}}T/\varepsilon =0.2$. The approximant crystal takes over below about ${{k}_{\text{B}}}T/\varepsilon =0.1$. However, while the quasicrystal–fluid coexistence point is relatively straightforward to determine, the coexistence point with the approximant crystal is less so, and we are likely to be overestimating the approximant's stability to some extent. The reason for this is that it is considerably more difficult to equilibrate the quasicrystal at low temperatures when there is a mixture of hexa- and pentavalent components in the system than in the previously considered work, and our calculation of the coexistence point between the approximant and the quasicrystal depends on how well equilibrated the quasicrystal is. The lower temperature limit we have determined thus gives the minimum region of stability of the quasicrystal, and its true region of stability is expected to be somewhat larger still, since, unlike for pure pentavalent particles, the binary quasicrystal could largely be fully bonded, and so we expect little difference in energy between such a configuration and the approximant. An approximate phase diagram is shown in figure 4⁵. Interestingly, the quasicrystalline phase is stable over a wider range of temperatures than for the pure pentavalent patchy particle system we previously studied; [31] this is because the energy difference between the quasicrystal and the approximant in the current system is smaller than that between the quasicrystal and the σ phase for the pentavalent particles.

Zoom In Zoom Out Reset image size

Unlike in the single-component phase diagram considered in [31], the low-temperature phase, at which the configurational entropy afforded by the quasicrystal is no longer as important as it is to maximise the bonding, is the approximant crystal rather than the σ phase, since the presence of hexavalent particles means that bonding cannot be maximised in the σ geometry for all particles. However, it would alternatively be possible that the mixture may phase separate into a σ and a Z phase. We have confirmed that the zero-temperature enthalpy of the approximant phase is lower than the enthalpy of a 1:12 mixture of Z and σ phases comprising solely hexa- and pentavalent particles respectively at all pressures considered, even with no interfacial energy penalty imposed. The approximant crystal is therefore stable with respect to phase separation at zero temperature.

The stability of the quasicrystal over a range of conditions demonstrates that it is possible to set up a competition between penta- and hexavalent co-ordination by means other than having a large patch width. In particular, in previous work, the quasicrystal was never found to be stable for patch widths below ${{\sigma}_{\text{pw}}}\approx 0.45$ [31]. By contrast, we have shown here that a quasicrystal phase is thermodynamically stable even if the patch width is considerably narrower (i.e. ${{\sigma}_{\text{pw}}}=0.3$). This is very encouraging if our aim is to construct quasicrystals using DNA multi-arm motifs. However, at this stage, the set-up we have considered completely ignores one of the most important reasons why DNA is so popular in experiment: using DNA makes it very easy to design mutually orthogonal interactions. Indeed, the multi-arm DNA tiles equivalent to the patchy particles that we have considered thus far, where all the sticky ends at the end of the arms have a favourable interaction with all other sticky ends, might be rather more difficult to self-assemble than envisaged because growth might be arrested as multiple undesigned bonds are formed, and so it is prudent to investigate whether it remains possible to form quasicrystals if the patches are made to be more specific. However, while the kinetics of self-assembly may be less frustrated as the interactions are made more specific, it is important to bear in mind that it is in the nature of quasicrystals that they are not completely ordered (indeed, this is what affords them their additional entropic stability): it is not possible, by construction, to have a set of interactions for which the fully bonded configuration is uniquely determined to be the quasicrystal.

As a first step in exploring the factors that govern quasicrystal self-assembly in binary mixtures corresponding to DNA star tiles, we consider the interaction set required to form two of the quasicrystalline approximants considered above (figure 2), before considering how these interactions could be relaxed to allow the variety of environments typical of our target dodecagonal quasicrystal to form. Let us first consider how we can make every distinct type of interaction that can be identified in the unit cell of figure 2(a) different. This can be achieved by making the pentavalent particles of two different types, as depicted in figure 5(a). In this set-up, patches only interact with complementary patches denoted by the same letter and an asterisk; other patch pairs do not interact at all. The unit cell of the edge-sharing approximant of figure 2(a) can readily be identified in the approximant shown in figure 5(b). However, whilst the approximant is stable in roughly the same conditions as it was before, the quasicrystal is not expected to form with such specific interactions. Even when cooled very slowly, kinetic products such as that shown in figure 5(c) are obtained. Whilst the underlying approximant crystal ordering can certainly be identified in this figure, the fact that these dodecagonal motifs are not orientated in the same way throughout the simulation box means that large gaps must be left in order to reduce the strain in the system, and since the bonding is so specific, no particles can be used to 'glue' the different regions together. The resulting structure is therefore, unsurprisingly, full of defects: it is not a single crystal, but has several crystalline domains with grain boundaries between them. The fact that there are multiple crystallites is confirmed by the smeared out diffraction pattern that involves a superposition of patterns from the different crystallites (figure 5(c)).

Zoom In Zoom Out Reset image size

The alternative approximant motif of figure 2(b) is based on overlapping rather than edge-sharing dodecagonal motifs. In the 'unit cell' of this alternative approximant crystal, there are 2 particles in a hexagonal (Z) environment and 12 particles in a σ environment, necessitating a ratio of 1:6 of hexa- and pentavalent patchy particles in the simulation box. If the bonding interactions are made fully specific with this alternative set-up, as shown in figure 6(a), we can again stabilise the approximant crystal (figure 6(b)). However, similarly to the edge-sharing fully specific system, on cooling, multiple nucleation events occur, leading to multiple crystallites of the overlapping approximant separated by highly defective grain boundaries.

Zoom In Zoom Out Reset image size

In order to improve the kinetics whilst retaining enough plasticity in the interactions to allow quasicrystals, rather than just crystals, to form, we can aim to strike a balance between the full specificity considered in simulations illustrated by figures 5 and 6 on the one hand and the completely non-specific bonding of figure 3. To this end, we have considered an alternative set-up in which we allow a competition between the overlapping and edge-sharing dodecagonal motifs to be set up. As illustrated in figure 7(a), by permitting the outlying 'type C' particles to bond freely with two of the patches of 'type B' particles, we can assemble either structures analogous to those of figure 5, involving all three types of particle, or of figure 6, involving only particles of types A and B, with the excess of particles of type C forming regions of σ-environments that can fill the gaps between 'approximant' motifs that are orientated in different ways. We have chosen the composition of particle types in this set-up to be the same as the one we considered for the edge-sharing approximant system of figure 5.

Zoom In Zoom Out Reset image size

With this set-up, cooling a liquid again results in a quasicrystalline phase, such as that shown in figure 7(c), with the corresponding diffraction pattern confirming its dodecagonal symmetry. As in the non-specific case of figure 3, there is again a variety of ways in which the dodecagons pack in addition to the two triangular lattice patterns of the approximants of figure 2. To verify that the quasicrystal is thermodynamically stable, we have repeated the calculations considered above for the non-specific case at $\sigma _{\text{LJ}}^{2}\beta P=1.5$. In particular, we have computed the free energies of the edge-sharing and overlapping approximants as well as the σ phase using Frenkel–Ladd integration. The free energy of the edge-sharing approximant matches the free energy of a system combining 7 parts of the overlapping approximant and 6 parts of the σ phase; this ratio accounts for the excess of 'type C' particles when an overlapping approximant is formed at the considered composition of particles of types A, B and C (see figure 7(a)). At temperatures above ${{k}_{\text{B}}}T/\varepsilon \approx 0.1$, the quasicrystal's free energy is lower than that of the approximants, confirming that it is the thermodynamically stable phase across a range of temperatures.

The fact that the quasicrystalline phase is thermodynamically stable for this system of 'intermediate' specificity suggests that a DNA star tile system with interactions chosen in this way would be the most likely to result in a two-dimensional soft DNA-based quasicrystal. However, it is certainly the case that we have ignored a number of considerations when abstracting the system to the toy-model level considered here. For example, it is not at all clear a priori that dodecagonal motifs comprising DNA star tile would be sufficiently planar to permit the growth of suitably large two-dimensional quasicrystalline structures, although this consideration may be mitigated to a large extent by performing the self-assembly on a surface. In order to address this point, we turn briefly to a more realistic potential of DNA molecules themselves.

OxDNA [68, 89] is a coarse-grained DNA model at the nucleotide level that allows the simulation of systems of large numbers of nucleotides. The nucleotides are modelled as rigid bodies interacting with a series of effective interactions (the solvent is not explicitly modelled) that account for hydrogen bonding between Watson–Crick base pairs, stacking between bases, electrostatic repulsion between the phosphates, excluded volume and chain connectivity. These interactions have been fitted to reproduce the thermodynamics of hybridisation and the structural and mechanical properties of both single-stranded and double-stranded DNA. Here, we use the second version of the model ('oxDNA2') that includes fine-tuned properties to reproduce better the properties of large nanostructures—in particular DNA origamis [68].

The oxDNA model is the most widely used coarse-grained model of DNA at the nucleotide level, and has been used to study the biophysical properties of DNA, a wide variety of DNA nanotechnology systems and applications in soft matter materials [90]. These applications have confirmed the model's robustness and quantitative accuracy (e.g. [91]'s reproduction of six orders of magnitude variation in the kinetics of strand displacement). Particularly relevant to the current application is the model's ability to account for the structural properties of both bulged duplexes [92], a motif that is crucial to the properties of DNA star tiles, and polyhedral assemblies of star tiles [93], and to rationalise how the kinetics of star tile self-assembly can be controlled through the size of the bulges they contain [94].

Our aim here is to use oxDNA to explore the structural stability of quasicrystalline arrays made out of DNA star tiles. To achieve this, we first need to generate starting initial structures for these arrays. We do this by first designing the arrays using vhelix [95], a recent DNA nanostructure design programme that allows free placement of the component DNA helices in space rather than on a lattice [96]. We then convert the vhelix design into a starting geometry for the oxDNA model. This geometry is not yet a suitable starting point for molecular dynamics simulations, as it may have particle overlaps or extended bonds that give rise to unreasonably large energies and very large forces. We therefore first relax the structure using a steepest-descent-like minimisation technique. The details of these procedures will be described elsewhere [97].

The structures are then simulated with a molecular dynamics algorithm employing an Andersen-like thermostat [98] both to keep the temperature constant and to generate diffusive motion of the nucleotides, as is appropriate for molecules in solution. As the systems we study contain thousands of nucleotides (the largest has 50 tiles and 22 704 nucleotides), to make the simulations feasible on a reasonable time scale, they are run on GPUs using a specially developed code [99]. We consider systems of tiles both when free in solution, as is typical during the assembly process, and when adsorbed on a surface, as is the case when visualised by some types of microscopy (e.g. AFM). The interaction of the nucleotides with the surface is modelled with a simple one-dimensional Lennard-Jones interaction that depends only on the distance of a nucleotide from the surface⁶.

Our aim here is use oxDNA to check whether there might be any structural reasons why the dodecagonal quasicrystals that we have seen for the above patchy particles might not be realisable using DNA star tiles. Previous experiments do not suggest any obvious hindrances. For example, both five-arm and six-arm tiles have been produced and found to assemble into two-dimensional crystalline arrays, forming σ [59] and hexagonal [58] crystals, respectively. Furthermore, mixtures of three- and four-arm tiles with specifically designed interactions have been shown to be able to produce more complex crystal structures [60, 61]. One possible complication is that, if the tiles are not flat, but rather the arms possess an intrinsic preference to bend in a given direction, then if all tiles face in the same direction, this has the potential to lead to curvature in the resulting structure that could hinder the assembly of an extended two-dimensional structure. One solution is to design the tiles so that they alternate in orientation, and any curvature cancels out [57], but this approach is not available for structures that possess polygons with an odd number of edges. However, in the examples above [58–61], conditions and designs were still found for which the assembly of extended two-dimensional arrays dominated over the formation of finite closed objects (e.g. icosahedra for five-arm star tiles [59]). Furthermore, self-assembly of DNA star tiles on a surface has also been shown to be possible [100].

Example five- and six-arm star tiles are illustrated in figures 8(a) and (b). The tiles consist of a long central strand, five or six medium-length strands that bridge two arms and five or six short strands which bind at the ends of each arm. The bulges on the long strand between the arms provide flexibility, allowing the arms to bend back on themselves. In the current examples, there are four nucleotides in the bulges, the same as was used experimentally to produce extended structures with five- and six-arm tiles [58, 59]. The lengths of the arms of the tiles are also the same as in these experiments.

Zoom In Zoom Out Reset image size

We have constructed DNA analogues of three example motifs that are important for the quasicrystalline structures observed in our patchy-particle simulations, namely a dodecagon, three overlapping dodecagons, and three edge-sharing dodecagons. The simulations of these structures showed that they are all stable at room temperature with the correct topology of the network maintained throughout the simulation. Example configurations that have been adsorbed on a surface (figures 8(c)–(e)) clearly show the expected structures. Due to the flexibility of the star tiles, the quadrilaterals in the network need not be perfectly square. Furthermore, bending is not always localised to the bulge regions of the tiles, and can sometimes occur at the four-way junctions in the arms leading to further distortions from the idealised geometries, even for the triangles.

By contrast, when free in solution, although the topology of the network is retained, the motifs are highly fluxional and no longer look anything like the idealised two-dimensional target structure (figure 8(f)). This is both because of the inherent flexibility of the star tiles and because all tiles face in the same direction, so any tendencies for the arms to bend away from the plane in a preferred direction are additive. The non-planarity is probably also exacerbated by the relative small size of the motifs, as consequently a large number of the arms (over 18% even for the largest motif) are on the edge of the motif and lack the constraint of being connected to another tile.

We should emphasise that the flexibility of the structures and the large fluctuations away from planarity do not affect the stability of the networks, nor do they mean that further self-assembly of the networks is necessarily hindered. When a new tile binds to a free arm on the edge of the motif, further binding is probably still most likely to occur in the intended way. However, the non-planarity may also make allowed, but not intended, arm binding more likely than when in a planar geometry because the relevant arms have been brought closer together by the non-planar fluctuations. This again emphasises the importance of annealing to facilitate the melting away of incorrect bindings.