---
abstract: 'Inspired by recent successes using single-stranded DNA tiles to produce complex structures, we develop a two-step coarse-graining approach that uses detailed thermodynamic calculations with oxDNA, a nucleotide-based model of DNA, to parametrize a coarser kinetic model that can reach the time and length scales needed to study the assembly mechanisms of these structures.  We test the model by performing a detailed study of the assembly pathways for a two-dimensional target structure made up of 334 unique strands each of which are 42 nucleotides long.  Without adjustable parameters, the model reproduces a critical temperature for the formation of the assembly that is close to the temperature at which assembly first occurs in experiments.  Furthermore, the model allows us to investigate in detail the nucleation barriers and the distribution of critical nucleus shapes for the assembly of a single target structure.  The assembly intermediates are compact and highly connected (although not maximally so) and classical nucleation theory provides a good fit to the height and shape of the nucleation barrier at temperatures close to where assembly first occurs.' 
author:
- Pedro Fonseca
- Flavio Romano
- 'John S. Schreck'
- 'Thomas E. Ouldridge'
- 'Jonathan P. K. Doye'
- 'Ard A. Louis'
title: 'Multi-scale coarse-graining for the study of assembly pathways in DNA-brick self-assembly'
---

\[sec:level1\]Introduction
==========================

Multi-component self-assembly has emerged over the past years as a promising route towards the fabrication of supramolecular structures at the nanoscale.  Unlike conventional crystals and many self-assembling systems in soft matter physics where a small number of building block types is used many times to assemble a desired product, in such *addressable assembly*,[@Cadamartiri15; @Frenkel2014; @Jacobs_review_paper] the use of unique building blocks (programmed to occupy a unique “address" in the desired target) offers direct control over shape, size and functionality, thus allowing for a much wider design space. 

One of the most popular ways to achieve addressable multi-component assembly is with short strands of DNA, by virtue of the specificity of its interactions.[@Wang2016; @Rothemund_origami; @PengYin2d] Typically, each DNA strand is divided into various domains, each programmed (through the choice of sequence) to hybridize with a different neighbouring strand in the target structure.  By designing which of the strands hybridize with each of the other strands, the shape and composition of the target structure can be uniquely defined.  Besides offering the necessary specificity to codify interactions, the use of DNA also ensures that bonds between building blocks also have a degree of reversibility, granting self-assembling systems some capacity for self-correction if undesirable bonds eventually form. 

DNA origami,[@Rothemund_origami] where many short strands are designed to bind to a long single-stranded scaffold, has been the most used technique for designing large self-assembling structures in DNA nanotechnology.[@Rothemund_origami; @origami_review] More recently, a new scaffold-free method that uses only short DNA strands has been developed by the Yin group.  The single-stranded tile (or DNA brick) approach involves assembling large three dimensional structures with intricate shapes and complex functionalities from thousands of unique strands, using a one-pot thermal annealing protocol.[@PengYin2d; @PengYin3d; @Peng_Yin_Dephts; @Peng_yin_10nm; @Peng_yin_tubes; @Peng_yin_isothermal; @Peng_yin_design_space; @Peng_yin_reconfigurations; @Peng_yin_10000]

This addressable self-assembly method, may, at first glance, seem straightforward, since each unit is designed to be in just one unique position.  In common with other examples of self assembly of finite objects with highly specified interactions,[@Parent2005; @Wilber_2009; @Hagan_2014aa; @perlmutter2015mechanisms] the relative rates of nucleation and growth need to be carefully controlled.  The nucleation barrier needs to be low enough so that it can be crossed on experimentally accessible time scales, whilst also being large enough that the growth time for structures to reach their equilibrium shape is significantly shorter than the time for the next nucleation event.  If the latter does not hold, the system will rapidly grow many more partially formed nuclei than the maximum number of possible assembled structures, thus using up available monomers (monomer starvation). 

However, there are also additional hurdles to correct assembly that are specific to addressable assembly.  Firstly, the extra entropy associated with the large number of component strands pushes the assembly transition down to lower temperatures compared to an equivalent one-component system.  Secondly, although incorrect binding is a potential problem for any kind of self-assembly, the number of pairs of components between which misbonding could occur is vastly increased.  Thus, the strength of the correct interactions must be strong enough in comparison to undesirable interactions to compensate for the vastly larger number of potential incorrect interactions.  Although DNA provides an ideal material to have the requisite specificity of interactions, it is not inherently obvious that there will be an assembly window where the conditions that there is sufficient thermodynamic driving force and yet misbonding is not prevalent are both satisfied.  Thus, the success of the first demonstrations of the single-stranded tile (SST) method was conceptually somewhat surprising. 

The need to explain the success of the SST assembly method inspired a range of important theoretical work.[@Jacobs_review_paper] For example, Reinhardt and Frenkel, using a lattice model[@lattice_model1] where individual strands of DNA are modelled as patchy particles with rigid patches representing binding domains, were first to demonstrate the successful assembly of a large structure in a model of addressable assembly.  This same model was later used to explore the effect of the coordination number on the nucleation barriers, and hence the assembly success.[@lattice_model2] An off-lattice model,[@Aleks_off_lattice] which provides a somewhat more robust description of the translational and rotational entropy of the single strands, was also found to reproduce the main physical phenomena observed in the simpler lattice model. 

In parallel to these numerical investigations, theoretical calculations by Jacobs *et al.*,[@lattice_model_wang2] helped explain why a time-dependent annealing protocol is necessary to achieve a fully assembled structure.  Due to the large combinatorial entropy for these systems, as temperature is cooled and the critical nucleation barrier is crossed, the free-energy minimum is for a partially formed structure.  Further cooling is needed for the final structure to completely form.  Thus addressable assembly differs from more conventional assembly with a small number of sub-units, where complete assembly can typically occur at a single temperature.  Similar theoretical calculations[@lattice_model_wang] were also used to predict how nucleation barriers depend on the topology of the graph describing how the particles are connected in the assembled state.  More recently, this model was used to investigate strategies to optimise nucleation and growth of unbounded, periodic structures.[@Jacobs_protocol_design] Another strand of research into addressable assembly, also using simplified models,[@Murugan_undesired_usage] has suggested that using carefully chosen non-stoichiometric abundances of assembly components may help to mitigate aggregation and undesirable monomer starvation effects. 

While these calculations have provided important qualitative insights, these minimal models neglect important aspects of these systems, such as the polymeric nature of the constituents, thus preventing quantitative predictions.  More detailed modelling would be beneficial to help guide SST experiments and to work out the relative importance and the interplay of the various mechanisms that affect the addressable assembly of a large structure from a set of individual single strands of DNA. 

![ A typical oxDNA configuration illustrating the fully assembled molecular canvas at $T=25^\circ \textnormal{C}$. This snapshot was taken after allowing an initially fully planar configuration to equilibrate using molecular dynamics simulations with oxDNA2,[@snodin2015introducing] an extension of the original oxDNA that allows for a better description of structural properties of large (kilobase-pair) structures.  Distinct tile species are depicted using different colours.  At this low temperature the entire 334 strand structure is highly prevalent in equilibrium under typical conditions.[]{data-label="fig:membrane"}](membrane.jpg){width="45.00000%"}

The fundamental reaction exploited by the SST assembly method is the hybridization of different domains between separate single strands.  It is natural to model these reactions using oxDNA,[@tom_model_jcp; @oxdna_review] a nucleotide-level model of DNA that has successfully been applied to a diverse set of systems and processes where hybridization plays a key role, ranging from DNA devices such as nanotweezers,[@ouldridge2010dna] and DNA walkers[@ouldridge2013optimizing; @sulc2014simulating] to the dynamics of the displacement reaction[@srinivas2013biophysics; @Machinek:2014aa] to the breaking of duplex DNA under force[@romano2013coarse] to hybridization without[@ouldridge2013dna] and with[@schreck2015dna] secondary structure, and even to the assembly of small origamis.[@snodin2016direct] This success, and the good agreement (where comparison is possible) with experiment, suggests that oxDNA can be used to understand SST assembly. 

On the other hand, the scale of the challenge to model self-assembly for these systems is evident.  For example, the two-dimensional SST rectangle, first presented in Ref.   and illustrated in Fig.  1 using oxDNA, is made up of 334 unique single strands each 42 nucleotides long giving 14,028 nucleotides in total.  Each strand is composed of four concatenated domains of length 10-11-11-10 bases or 11-10-10-11 bases (depending if the strand occupies an odd or even row) that bond with four local neighbours during self assembly.  This rectangle serves as a molecular canvas, from which custom shapes can be created by removing certain strand types.  Using this strategy, over one hundred different shapes were successfully assembled using a one-pot annealing protocol where the temperature is lowered at a fixed rate.[@PengYin2d] While computer simulations with oxDNA can be used to study structural properties of a fully assembled rectangle, as illustrated in Fig.   \[fig:membrane\], it is computationally unfeasible to characterise the spontaneous self-assembly of this structure from a bath of single strands via direct simulations, as the correspondent time scales are far too long even given the speedups that coarse-graining at the level of oxDNA provides.  For example, the fastest experimental protocols that lead to correct assembly are on the order of $1^{\circ}\textnormal{C}$ per hour[@SobczakDNAfolding], while oxDNA time steps are on the order of several fs. 

To make progress, further coarse-graining is necessary.  We are assisted by the observation that the equilibrium ensemble of two complementary DNA strands at ideal conditions is dominated by states in which either most of the base pairs are formed between complementary strands, or no base pairs are formed, separated by a free-energy barrier that emerges mostly due to the sudden reduction of translational freedom when the first base pairs form.  This feature of DNA hybridization is well illustrated by the successes of a two-state model[@SantaLucia1998] in predicting the thermodynamic properties of DNA duplexes and other DNA motifs. 

This time-scale separation suggests that a kinetic model at the level of the chemical master equation, with each state specifying a bonding configuration, may capture the dominant physics of assembly.  We are inspired here by recent success of a kinetic model for DNA origami,[@Origami_plus_model; @dannenberg2015modelling] where excellent agreement with experiment was achieved. 

Our basic coarse-graining strategy for SST assembly of a 2D canvas is to first estimate the rates for key reactions using oxDNA, and then use these rates to parameterise a coarser kinetic model that allows us to work out the full dynamics of assembly.  The ratio of the on to off rates can be estimated using the free energy of hybridization together with a two-state approximation for the kinetics.  While the connectivity network of assembly units seems to be in the natural level of description from both kinetic and structural perspectives, a problem arises because the number of possible non-equivalent connectivity networks grows exponentially with the number of subunits composing the target structure. 

To circumvent this problem of determining transition rates between every pair of configurations we identify a finite set of local changes in the connectivity network associated with the binding/unbinding of a staple, and parameterise all transition rates according to which member of this finite set they correspond.  This approach means that we do not need all theoretically possible transition rates, but instead can treat self-assembly using a limited database of tile association/dissociation rates.  We then use standard kinetic Monte Carlo (KMC) techniques to sample over all possible states. 

The full multi-scale procedure to coarse grain our kinetic model and to calculate the dynamics of assembly is summarised in Fig.  \[fig:cg\_scheme\].  The resulting KMC trajectories can easily reach the time-scales of multiple hours used in the assembly experiments[@PengYin2d; @SobczakDNAfolding] in just a few minutes of computer time.  The only free parameter is the on-rate $k_+^0$ which we assume to be temperature independent, and set by recent experiments of DNA strand displacement reactions.[@winfree_kinetics] This parameter sets time scales, but not free energies (and therefore transition temperatures).  Our calculations, which otherwise have no adjustable parameters, predict a critical temperature that is only a few degrees Celsius different from that measured in experiment. 

We proceed as follows.  In the methods section, we describe our coarse-graining methods in more detail, explaining how we define the key states of our model, how we calculate free-energy differences using oxDNA, how we incorporate these thermodynamic calculations into our kinetic rates model, and how we use enhanced sampling KMC methods.  We also discuss the importance of some effects we neglect in our multi-scale coarse-graining procedure, including a discussion of the role of misbonding.  In our results section, we first consider the assembly (and disassembly) of a single target structure during both temperature ramp and isothermal protocols.  Next we consider the nucleation pathway in more detail, calculating critical nucleus sizes and free-energy profiles as a function of assembly size for different temperatures.  Finally, we present some conclusions, pointing to future applications of our new method. 

![ image](scheme.jpg){width="100.00000%"}

\[sec:level1\]Methods
=====================

\[sec:level2\]Definition of the states
--------------------------------------

In our simulations, we follow the assembly and disassembly of a single structure through monomer association from, and dissociation into, a large bath of monomers.  We define the state of the system as a single vector made up of $N$ components $\vec{x} =(x_{1},x_{2},\cdots ,x_{N})$, where $x_{i}=1$ if strand species $i$ is in the assembly and $x_{i}=0$ otherwise, and $N=334$ is the number of 42-base-long strands in the fully assembled 2D structure that we are studying here.  We can make the simplification that the states only describe strands either fully on or off because of the essential two-state nature of strand association: the subensemble of states in which programmed hydrogen bonds are formed between the incoming strand and its binding domains in the assembly is statistically dominated by states in which all complementary domains are hybridized and most of the base pairs are formed (see details in Section II.B.2). 

The transition rates are denoted by $k(\vec{x}|\vec{y})$ for transitions $\vec{y}\rightarrow \vec{x}$. Possible transitions in the model are those in which a single strand, from now on also called a monomer, either binds to or is removed from the assembly.  For consistency, however, transitions in which the assembly is divided into two disconnected parts due to the removal of a monomer are not allowed as we do not describe the reverse process. 

The kinetic simulations performed throughout this work were initiated (unless stated otherwise) from the state $x_{\mu}=\delta_{\mu \nu}$, where $\nu$ was selected randomly from the set of 334 different tile species $(\nu \in [1,334])$, that corresponds to a single tile.  After the initial step of the simulation is performed (the association of a second tile), the first tile is allowed to dissociate, but we always keep at least one tile in the state vector.  Typically a great number of association and dissociation steps occur before any nucleation events happen, so that memory of the identity of the initial tile is effectively erased. 

By using the state definition described above, we focus our attention on the assembly of a single target structure, while the concentration of strand species is kept constant.  Our model is therefore appropriate to describe the first stages of an experimental assembly setup, when the first structures start to form and there is little variation in the concentration of monomers.  To capture the full complexity of an assembly process, future work needs to describe explicitly the variation with time of the composition of the assembly mixture, which contains multiple assembly fragments and (usually) a fixed number of monomers.  In the present work, we instead focus our attention on the nature of the assembly pathways for single-target assembly, and leave consideration of the collective aspects of the assembly process for future work. 

\[sec:level2\]Transition rate model
-----------------------------------

### \[sec:level3\]Tile association rates

The state space of the model is based on the assumption that the rearrangement of a strand to form secondary domains is fast compared to the initial binding of the first domain when a strand arrives from the dilute, well-mixed solution.  Furthermore, association rates of duplexes are known to be less sensitive to sequence and conditions than dissociation rates.[@Zhang_kinetics; @tom_inc_1; @tom_nm_1; @schreck2015dna] In the spirit of minimal parameterization of the model, therefore, we assume that tile-association rates do not vary with temperature and use the standard value for the second order rate constant of $k_+^0=6\times 10^{6}\textnormal{M}^{-1}\textnormal{s}^{-1}$, that yields the best fit[@winfree_kinetics] to a set of experimental data on second-order hybridization reactions that are part of a displacement reaction system (see details in Section S.X).  Therefore, for a state $\vec{y}$ that is obtained from $\vec{x}$ by adding a single tile, we estimate the tile-association rate as $$k(\vec{y}|\vec{x}) = k_+=k_+^0 [{s}],$$where $[{s}]=100\hspace{1mm}\textnormal{nM}$ is the concentration of strand species $s$ (assumed the same for all strand species) in the systems studied here. 

### \[sec:level3\]Tile dissociation rates

If tile association is treated as simple bimolecular binding reaction, as here, then dissociation rates for the reaction $\vec{y} \rightarrow \vec{x}$ can be related to association rates via $$\frac{k(\vec{x}|\vec{y})}{k(\vec{y}|\vec{x})}= \frac{k_-(\vec{x},\vec{y})}{k_+^0[s]} = e^{\beta \Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})} \times \frac{1M}{[s]}. 
\label{eq:deltaG}$$ Here $  \Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})$ is the standard free-energy change of formation for the reaction $\vec{x} \rightarrow \vec{y}$, and ${k_-(\vec{x},\vec{y})}$ is the rate of dissociation of a strand from $\vec{y}$ to form $\vec{x}$. This quantity takes into account the balance between the entropic cost of bringing a diffusing strand into contact, in the correct orientation, and the enthalpic gain due to the formation of hydrogen bonds and stacking interactions between adjacent bases upon tile association. 

![ Examples of local environments for tile association.  The red tiles denote the assembly fragment that contains the first neighbours of the binding site (coloured grey).  The sections in the assembly fragment that are hybridized are shown in blue.  The numbers 10 and 11 indicate the length (in bases) of the binding domains positioned above.  The examples involve (a) one, (b)-(d) two, (e) three and (f) four binding domains.  The full list of 32 local environments used can be found in the supplementary Section S.VI.[]{data-label="fig:le_examples"}](local_env_d.jpg){width="45.00000%"}

The major biophysical input into our model is the use of oxDNA simulations to estimate $ \Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})$, and hence ${k_-(\vec{x},\vec{y})}$, for each tile dissociation reaction.  While oxDNA can be used to directly calculate relative rates using either direct molecular dynamics calculations or using techniques such as forward flux sampling (FFS),[@hairpin_stacking; @Machinek:2014aa; @schreck2015dna; @srinivas2013biophysics; @design_knotted_DNA] for the relatively simple reaction rates that we need to calculate, an equilibrium calculation of free-energy difference in Eq.  \[eq:deltaG\], combined with an estimate of hybridization rates from the literature, is more practical. 

In principle oxDNA could be directly used for all possible tile-association events, but since the assembly can proceed in a very large number of different ways there are far too many possible transitions to make this direct approach feasible.  Additionally, for larger tile-fragments, a lot of the computational load would involve degrees of freedom that are far from the site where hybridization is occurring.  To simplify our calculations, we apply a *local environment approximation* in which we assume that the enthalpic and entropic factors that govern $\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})$ can be approximated by considering the minimal assembly fragment that contains the first neighbours of the attaching site in the assembly (see Figs.  \[fig:cg\_scheme\].I and \[fig:cg\_scheme\].II). 

Using this approximation, each possible transition in the model falls within a relatively small number of distinct cases (local environments), each specified by the number, size, and position of binding domains that hybridize during tile association, and the number and position of bases adjacent to the binding domains on the assembly that coaxially stack.  Several examples of local environments are illustrated in Fig.  \[fig:le\_examples\], and the full set of 32 local geometries that we use is given in the supplementary information (Section S.VI). 

For each local environment, we determine the free energy of tile association by simulating the minimal assembly fragment that contains all first neighbours of the binding site, and the incoming tile with the correct sequence, in a periodic box using oxDNA.  In the simulations we use the average-base version of the model,[@tom_model_jcp] which allows complementary base-pairing between AT and CG bases, but averages over the hydrogen-bonding and stacking parameters.  Although tile-binding energies are expected to depend on the tile’s sequence, we use this approximation in order to focus our attention on generic mechanisms of tile association rather than the effect of using a particular sequence.  In addition, the simulations are set up so that the formation of secondary structure for single strands is prevented.  The advantage of these simplifications is that only one calculation is needed for each local geometry.  Sequence dependence can easily be included into oxDNA,[@vsulc2012sequence] but for simplicity we have not done this here.  Future work may include sequence dependence.  This may shine light on conflicting predictions in the literature, with some arguing that adding sequence heterogeneity stabilises assembly intermediates,[@lattice_model_wang] while others predict that it hinders assembly.[@Brenner_SA]

The free-energy profiles were first sampled as a function of the number of native base pairs between fragment and incoming strand, using the virtual-move Monte Carlo (VMMC) algorithm of Whitelam and Geissler[@Whitelam_VMMC_2] that attempts moves of dynamically selected clusters of nucleotides, greatly accelerating sampling in oxDNA over Metropolis Monte Carlo (see details in Section S.II.A).  In addition, umbrella sampling [@Torrie_Umbrella] was used in order to improve the sampling of rare transitions between free-energy minima (see details in Section S.II.B).  In the simulations, a base pair is considered bonded if the respective hydrogen bonding energy falls bellow $-0.596\hspace{1mm} \textnormal{kcal}\hspace{0.7mm}\textnormal{mol}^{-1}$ (which corresponds to about 15% of the typical hydrogen bond energy).  The number of base pairs between designed complementary domains is defined as the order parameter and used to divide the configurational space into regions biased by the umbrella potential.  Base pairs that were not designed to be part of the assembled structure are not allowed to form in the simulations.  Due to the complexity of reactions in which a tile associates to the assembly through more than one binding domain, the order parameter was augmented with other variables, that specify the distance between the centre of mass of each binding domain of the tile and the respective complementary sequence in the assembly fragment.  To improve the efficiency of the calculations, the order parameter space was split into sampling windows, and the statistics were later combined using the weighted histogram analysis method (WHAM).[@kumar_wham] The simulations were performed at $T=50^\circ \textnormal{C}$, which is within the range where the transition temperature was expected to be, and the statistics were extrapolated to other temperatures using single histogram reweighting [@vsulc2012sequence] (see details in Section S.II.C). 

For each local environment, the tile-association free energy $\Delta G_{\textnormal{sim}}$ is defined as the free-energy difference between the subensemble of states in which at least one base pair is formed between complementary domains $(\Omega _c)$, and the subensemble of diffusive states, in which no base pairs are formed between the assembly and tile $(\Omega _d)$: $\beta \Delta G_{\textnormal{sim}}=-\textnormal{ln}[p(\Omega_c)/p(\Omega_d)]$. As an example, for the local environment depicted in Fig.  \[fig:cg\_scheme\], $\Omega _c$ corresponds to the subensemble of states satisfying $N_{10}+N_{11}>0$ where $N_{10}$ and $N_{11}$ denote the number of base pairs formed between complementary domains that are 10 (between green and blue tiles) and 11 (between green and red tiles) nucleotides long, respectively, and $\Omega _d$ corresponds to the subensemble of states satisfying $N_{10}+N_{11}=0$. To explore the accuracy of the two-state approximation, we also consider free-energy landscapes in which the bound state is further subdivided into states with different degrees of bonding. 

The simulations were performed in a small volume $v$, implying a higher effective concentration $1/v$ than in experiments, so as to reduce diffusion times.  We made sure that the systems were large enough so that periodic images could not interact.  Afterwards, the tile-association free energies and free-energy profiles were extrapolated to concentrations relevant to experiment by adding the factor $\Delta \Delta G_{\mathrm {vol}}=k_{B}T \ln(u/u_0)$ to the free energies of states with base pairs.  Here, $u_0$ is the desired tile concentration and $u$ the tile concentration used in the simulations, which varied between $4.38\times 10^{-6}\textnormal{M}$ (for the smallest systems), and $1.72\times 10^{-7}\textnormal{M}$ (for the largest ones).  This extrapolation provides the best description for the thermodynamics of the binding process under the relevant conditions.  Finally, we determined the effective tile-association free energies at standard concentration ($\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}$), as employed in the kinetic model.  To do so, we employed the method in Ref.   that relates bonding probabilities observed in a high concentration, single-replica simulation to $\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}$.

In Fig.  \[fig:cg\_scheme\].II we show a typical equilibrium free-energy profile of tile association obtained from the simulations, determined here for local environment 7 (see Section S.VI).  This same environment is used throughout Fig.  \[fig:cg\_scheme\].  We also present the respective probability profile, projected along a single reaction coordinate that specifies the total number of base pairs formed at $50^\circ \text{C}$. The probability profile in Fig.  \[fig:cg\_scheme\].II helps illustrate that the ensemble is statistically dominated by states in which either most base pairs are formed (and all complementary domains hybridized) or no base pairs are formed.  This feature is most pronounced in cases in which the tile binds with the assembly fragment through 3 and 4 binding domains (see Section S.VI). 

Given the simple two-state nature of the transitions, one might ask whether a simpler approach could be used to estimate the free energy of tile association.  For example, the nearest-neighbour (NN) model of Santa Lucia,[@SantaLucia1998] the predictions of which were used to fit the oxDNA model, provides an accurate description of the thermodynamics of hybridization.  This model could be used for binding events where a tile is only able to form a single domain (e.g. Fig.  \[fig:le\_examples\]a).  However, the NN model is unable to provide any estimate of the free-energy barrier to initiating a second (Fig.  \[fig:cg\_scheme\].II) or subsequent domains (see Section S.VI).  The accurate capturing of these barriers, which arise from the loss of conformational entropy of the strands involved on initiating the formation of these additional binding domains, is the key extra ingredient that using the oxDNA model provides.  We find that this barrier is typically around $8 kT$.

Similarly, this is the key limitation of previous developed patchy-particle models for SST assembly[@lattice_model1; @lattice_model2] that use the NN model to provide the strength of the patch-patch interactions.  These models do exhibit a free-energy barrier for the formation of a second domain, but this arises from the loss in entropy of the patchy particle assembly (e.g. a patchy particle with only one bond is able to freely rotate about this bond, but this degree of freedom becomes constrained on the formation of a second bond) and not from a realistic description of DNA’s polymeric degrees of freedom.  In Section III.C, we will highlight some of the effects that our more accurate description of the thermodynamics of multiple domain formation has on the observed assembly pathways. 

Simulation of the coarse-grained kinetic model
----------------------------------------------

### Kinetic Monte Carlo algorithm

Having developed a database of rates, the task shifts to determining the time evolution of the assembly.  For that purpose we implemented a standard KMC algorithm.[@gillespie_kmc] Starting at assembly state $\vec{x}$, the algorithm determines the nearest neighbourhood of tiles composing the fragment and binding sites at the interface.  For tiles composing the fragment, the algorithm attributes the relevant dissociation rate from the kinetic database determined in advance.  If the fragment is composed by only two tiles, a single dissociation rate is attributed that corresponds to the dissociation of the dimer.  For sites at the interface the algorithm attributes an association rate (which we assumed to be independent of local environment and temperature).  Next, the algorithm determines the cumulative rate function $R(i)=\sum_{j=1}^{i}k(\vec{x_j} |\vec{x})$ of transitions $\vec{x}\rightarrow \vec{x_i}$ allowed, and chooses which transition to carry out by generating a random number $u\in [0,1)$ and finding the transition $i$ that satisfies $R(i-1)< uR(n_{trans}) \leq R(i)$, where $R(n_{trans})$ is the cumulative rate of leaving state $\vec{x}$: $R(n_{trans})=\sum_{i=1}^{n_{trans}}k(\vec{x_i} |\vec{x})$. In the case the assembly state is a dimer, the rate for dimer break-up is only counted once and the two possible destination states are selected at random if dissociation is picked by the algorithm.  Time is then updated, such that $\delta t$ follows the probability distribution for the time of first transition leaving $\vec{x}$: $P(\vec{x},t+\delta t)=\exp ( -R(n_{trans})\delta t )$, which is satisfied by choosing $\delta t=-R(n_{trans})^{-1}\log v$ where $v\in [0,1)$ is a random number. 

### The simulated ensemble

If the simulation is run for an arbitrarily long time under constant conditions, it reaches a steady state.  Here we discuss the nature of that steady state, and its relation to the equilibrium yield of partial assemblies.  In an ideal (dilute) solution, in which partial assemblies evolve independently in a large bath of tiles at fixed concentration $[s]$, the steady-state concentration of a cluster $\vec{x}$ is given by the equilibrium condition $${[\vec{x}]} = [s] \left(\frac{[s]}{1{\textrm M}} \right) ^{N_{\vec{x}}-1}   \exp  \left(-\beta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x})\right)
\label{eq:conc}$$ in which $G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x})$ is the free energy of assembly of $\vec{x}$ from its constituent strands at standard concentrations, and $N_{\vec{x}} $ is the number of strands in assembly state $\vec{x}$. Perfectly accurate values of $\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})$ would satisfy $$\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y}) = G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{y}) - G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x}). 
\label{eq:thermoSC}$$

In principle, one could study the assembly process given by a well-defined $\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})$, $k_+^0$ and $[s]$ in a fixed volume, and infer the equilibrium concentration of partial assemblies from the frequency with which they are observed in the simulation volume.[@bulk_extrapolation] Instead of simulating a fixed volume, however, we pursue an “assembly-eye view" of the same process in which we simply track the growth and shrinkage of a single assembly.  In this case, the steady-state probability of the system occupying a given state $\vec{x}$ in equilibrium is given by $$p^{\textnormal{eq}}(\vec{x}) = 
\begin{cases}
\frac{1}{2Z}, \hspace{5mm} N_{\vec{x}} =1\\
\frac{1}{Z} \left(\frac{[s]}{1 \mathrm M}\right)^{N_{\vec{x}}-1} \exp \left( -\beta \Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x})\right), \hspace{5mm} N_{\vec{x}} >1,
\end{cases}
\label{eq:ss}$$ where $Z$ is a normalising factor.  Eq.  \[eq:ss\] can be verified by observing that, given these probabilities, the net flux of trajectories between all $\vec{x}$ is zero under the simulation procedure (yielding a detailed-balanced steady state).  We note that $p^{\textnormal{eq}}(\vec{x})/p^{\textnormal{eq}}(\vec{y})$ is consistent with the equilibrium concentrations in Eq.  \[eq:conc\] unless either $N_{\vec{x}} =1$ or $N_{\vec{y}}=1$ (in this case, the ratio deviates by a factor of two).  In order for the equilibrium probabilities in the “cluster’s-eye view" to map directly into concentrations, we correct, after sampling, the equilibrium probability of finding the system in a state $\vec{x}$ with $N_{\vec{x}}=1$ by a multiplicative factor $2$.

In fact, due to the local approximation to $\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x},\vec{y})$, and the finite sampling accuracy of oxDNA simulations, Eq.  \[eq:thermoSC\] is only approximately satisfied.  Detailed balance is therefore marginally violated in simulations; we perform the following calculation to gauge whether deviations from detailed balance are significant. 

Let $\vec{x}_i$ and $\vec{x}_j$ be two assembly states, and $\wp_{ij}\equiv \vec{x}_i \rightarrow \vec{x}_{i+1} \rightarrow \vec{x}_{i+2} \rightarrow \cdots \rightarrow \vec{x}_{j-1}  \rightarrow \vec{x}_j$ be a reaction pathway that connects $\vec{x}_i$ and $\vec{x}_j$. If detailed balance is satisfied exactly, the total free-energy change $$G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x}_j) - G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{x_i}) = \sum_{ \vec{y} \rightarrow \vec{z} \in \wp_{ij}}   \Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}(\vec{y},\vec{z}),
\label{eq:G_balance}$$ should be pathway-independent.  In Section S.VIII we consider a few examples of assembly states $\vec{x}_i$ and $\vec{x}_j$, and determine the factor on the RHS of Eq.  \[eq:G\_balance\] for different paths $\wp_{ij}$. We find the variation with path to be small for the examples considered (see Section S.VIII).  While differences between paths imply that detailed balance is not strictly satisfied, we find these differences to be relatively small and assume that they can be neglected as a first approximation for the remainder of this manuscript.  A similar approach was taken in a recent model of DNA origami folding.[@Origami_plus_model; @dannenberg2015modelling]

We further gauge the effect of detailed balance deviations by determining free-energy profiles as a function of assembly size using the biasing method introduced in Section II.C.3, which is based upon the assumption of detailed balance, and comparing the resulting profiles with those obtained through unbiased sampling (see Section S.IX).  We find good agreement between both methods which further supports our assumption that deviations from detailed balance can be neglected as a good approximation. 

### \[sec:level2\]Enhanced KMC sampling

Using a KMC algorithm to determine equilibrium properties can be very inefficient if states of very low probability need to be crossed frequently to obtain good statistics.  Here, we address this issue by devising a simple biasing scheme, similar in spirit to umbrella sampling, that allows us to relate equilibrium properties computed using biased transition rates with equilibrium properties of the reference system.  Supposing that the state space of the system is divided in windows according to a certain order parameter $\boldsymbol{Q}$, if transition rates are modified according to $$k(\vec{y}|\vec{x})\rightarrow \frac{W[\boldsymbol{Q}(\vec{y})]}{W[\boldsymbol{Q}(\vec{x})]}k(\vec{y}|\vec{x}),
\label{eq:4}$$ then, assuming detailed balance in the underlying rates $k(\vec{y}|\vec{x})$, the probability distribution $P'(\boldsymbol{Q})$ of the modified system relates to those of the reference system $P(\boldsymbol{Q})$ through $$\frac{P(\boldsymbol{Q}_a)}{P(\boldsymbol{Q}_b)}=\frac{P'(\boldsymbol{Q}_a)}{P'(\boldsymbol{Q}_b)}\left[ \frac{W(\boldsymbol{Q}_b)}{W(\boldsymbol{Q}_a)}\right]^{2},
\label{eq:enhanced_kmc}$$ where $\boldsymbol{Q}_a$ and $\boldsymbol{Q}_b$ denote any two points of the order parameter space.  In practice, the weights $W(\boldsymbol{Q})$ were first adjusted so that every value of the order parameter could be observed in the simulations, regardless of the initial state.  From the first estimate of the probability distribution $P(\boldsymbol{Q})$ a new set of weights were derived from Eq.  \[eq:enhanced\_kmc\] by fixing $P'(\boldsymbol{Q})=\textnormal{const}$. The simulations were performed again using the new set of weights and the method repeated until a final set of weights $W(\boldsymbol{Q})$ is found, such that any two values of the distribution $P'(\boldsymbol{Q})$ differ by no more than 10%. 

\[sec:level2\] Neglected effects
--------------------------------

In our coarse-graining scheme we ignore several effects that we discuss below. 

### Long-range bridging effects

Firstly, [*long-rage bridging effects*]{} can be relevant if a tile either incorporates or dissociates from a position in the assembly that acts as a bridge connecting two distant sections (an example is presented in Section S.VII).  To gauge how important these cases are, we performed independent simulations using two different versions of the model.  In the first version, the rates of addition or removal of “bridging-tiles" were assumed to be the same as for an otherwise equivalent environment in which both sections of the assembly are constrained, and the local approximation can be applied.  As an example, the rates of tile association or tile removal for the local environment in Section S.VII were assumed to be the same as the rates for the environment in Fig.  \[fig:le\_examples\]b.  In the second version of the model, the addition (and removal) of “bridging-tiles" was prohibited altogether.  All the simulations performed were repeated using both versions of the model.  As the entropic penalty for bringing together both sections of the structure is not accounted for in the first version of the model, bridge formation is artificially favored compared to a (hypothetical) more complete version of the model where non-local effects are treated correctly.  Nevertheless, in the simulations performed with the first version of the model bridge formation occurred so sporadically that the overall time scales and thermodynamic properties barely varied with relation to those observed using the second version of the model.  Thus, it was concluded that bridge formation is rare enough that it can be ignored.  The results presented throughout this work were obtained using the second version of the model, where bridge formation was prohibited. 

We note that throughout this work we study assembly under moderate supercooling conditions, where we find growth to occur through a series of compact intermediates with a high number of bonds between components.  However, at higher supercooling the shape of intermediates is more ramified as tile association tends to be more irreversible.  In this regime, bridge formation might become more relevant to the assembly pathways, and thus require an explicit treatment. 

### Multiple copies of the same tile species in a structure

Secondly, we *ignore states in which multiple copies of the same tile species are present in a structure*, each copy bonded through a different set of (programmed) binding domains.  During tile association, as the first contacts begin to form, we expect further bonding to occur relatively quickly (it is a localised first-order process) compared to the association of a second tile, which is a second-order process whose rate we expect to be much slower at strand concentrations typical of experiment.  Of course, at sufficiently high strand concentrations this approximation is likely to begin to break down, but we expect this to only occur well above the range relevant to experiments. 

### Local environment approximation

Thirdly, we tested the [*local environment approximation*]{} for a few simple cases with 1 or 2 binding domains where we included a much larger number of strands in the assembly.  The differences in the respective free-energy profiles were very small, and so we expect this to be a good approximation. 

### Misbonding

Fourthly, we ignore [*misbonds*]{} in our kinetic model, assuming interactions respect the domain-level description, and only occur between domains that are intended to be complementary.  This approximation greatly simplifies our calculations.  In Section S.III we sketch an upper bound on the average effect of misbonds.  We briefly summarise the calculation below. 

Firstly, we attempt to identify an upper bound for the typical binding probability of an unintended strand to an available “site" within an assembly.  We assume that this site has 42 available bases (in general, an overestimate) and further treat these bases as if they were presented by a single contiguous strand, rather than on separate tiles (again, overestimating the stability of bonding by neglecting entropic penalties due to the complexity of the environment).  We generate a representative set of 42-base sequences for the site by using the sequences of the 334 tiles themselves. 

For each of these representative sites, we consider the total binding free energy of strands that are not intended to bind to any of the domains in the site.  We perform this calculation using the NN model,[@SantaLucia1998] considering the contribution to the ensemble of all possible base-pairing configurations with either one or two base-paired sections of 2 base pairs or more. 

At the relevant temperatures, we find that the binding free energies for a given strand and a given site are significantly higher than if the strand is designed to bind with that site, with a few extreme sequences with free-energies of association comparable to an $N_b=1$ local binding domain (see Fig.  S2).  Averaging over sequences, we can compare the relative probabilities that an available site will bind to a correct strand or an incorrect one in equilibrium.  We find that at the range of temperatures where self-assembly takes place, our upper bound on the misbond probability falls below that of an $N_b=1$ local binding domain (see Fig.  S3), which is the weakest of the correct binding motifs.  Because the binding probability at equilibrium of tiles with a higher connectivity $N_b>1$ is much higher than for tiles with a single binding domain (Fig.  S3), we conclude that misbonding has a negligible effect on assembly at equilibrium. 

We also estimate the effect of misbonding on dynamics by employing a simplified variant of our assembly model where correctly bonded tiles are not allowed to dissociate from the structure once bonded, but incorrectly bonded tiles can associate and dissociate with rates that were derived using our upper bound on the average probability for misbonding.  At higher temperatures ($45^\circ\text{C}-50^\circ\text{C}$), we find assembly times during an isothermal assembly protocol to be largely insensitive to incorrect bonding (see Fig.  S4).  At lower temperatures (typically $T<40^\circ\text{C}$) we find the assembly dynamics for an isothermal protocol to be increasingly slowed down due to misbonding. 

However, it should be kept in mind that at these low temperatures, our approach of studying growth of a single cluster is likely to be inappropriate, and the rapid nucleation of many clusters that are unable to grow to completion because of monomer starvation is likely to have a more dominant negative effect on the assembly dynamics than misbonding. 

Interestingly, the upper bound on misbonding does not appreciably change if we simply take a set of random strands (see Fig.  S2).  It should be kept in mind that typical DNA design spaces are enormous[@louis2016contingency] – for this system there are $1.9 \times 10^{25}$ different possible sequences of length $42$ – and so it is not hard to see that the probability of finding two random sequences that have a significant amount of complementarity is small.  For example, the probability of two random sequences of length $42$ having a complementary section 10-bases long, the typical length of the bonding domains in the assembly, is $6.2\times 10^{-4}$ (see Section S.V).  Most complementary motifs are short, and it is rare to find one more than about 5 bases long between any two strands in a set of 334 randomly generated tiles.  In the experiments [@PengYin2d] it was also found that randomly picking sequences (that is not trying to design out misbonding) did not make the SST system more likely to aggregate, which accords with our findings above.  For the SST assembly system, it appears to be the case that the exponentially large design space means that misbonding, which in the past was often thought to be a very significant barrier to the success of an addressable assembly method, can be avoided quite straightforwardly. 

![ Annealing and melting curves using different cooling/heating rates.  Also shown in the figure (dashed black curve) is the assembly size at equilibrium as a function of temperature.  Inset: derivatives with respect to temperature of the curves in the main plot.  The assembly peaks (inset) have maxima at $51.0^{\circ}\text{C}$ for $3\text{h}/^\circ\text{C}$, $50.7^{\circ}\text{C}$ for $1\text{h}/^\circ\text{C}$, $49.4^{\circ}\text{C}$ for $5\text{min}/^\circ\text{C}$ and $51.64^{\circ}\text{C}$ for equilibrium.  The melting peaks (inset) have maxima at $51.9^{\circ}\text{C}$ for $3\text{h}/^\circ\text{C}$, $52.1^{\circ}\text{C}$ for $1\text{h}/^\circ\text{C}$ and $53.4^{\circ}\text{C}$ for $5\text{min}/^\circ\text{C}$. []{data-label="fig:assembly_rate"}](assembly_rate.jpg){width="45.00000%"}

![ image](frames.jpg){width="100.00000%"}

\[sec:level1\]Results
=====================

Assembly with a thermal ramp protocol
-------------------------------------

The oxDNA model predicts a critical temperature of $T_\text{c}=51.64^\circ\text{C}$, defined as the temperature where it first becomes thermodynamically favourable for a large assembly to form (the details of how $T_c$ is calculated are given in Section III.C). 

In the original experiments to which we are comparing  [@PengYin2d; @SobczakDNAfolding] a thermal annealing protocol with steps of $1^\circ\text{C}$ was used.  For a cooling protocol of $3 \text{h}/^\circ\text{C}$, the maximum assembly rate was measured at $53^\circ\text{C}$ by Sobzak [*et al.*]{}.[@SobczakDNAfolding] Given the resolution of the $1^\circ\text{C}$ steps, we estimate that for experiments, the critical temperature is somewhat above $53^\circ\text{C}$, possibly close to $54^\circ\text{C}$, where structure growth was first observed in the assembly mixture.  A full comparison with the dynamics of assembly in the experiments would need to reflect the correct ensemble, which means a fixed number of tiles, so that the change in concentration of tiles with assembly growth can be taken into account.  Nevertheless, the critical temperature observed in single target simulations corresponds to the onset temperature where it first becomes thermodynamically favourable for structures to form in an assembly scenario, and so should be accurately represented by our single-target system.  The agreement with our calculated critical temperature is remarkable given that no adjustable parameters are used in the model that affect its thermodynamic predictions (note that the forward reaction constant $k_+^0$ only affects the time scales).  In Kelvin, the appropriate temperature scale for statistical mechanical models, this agreement is better than 1$\%$. Note that this agreement also implies that the oxDNA model is accurately capturing the magnitude of the barriers to multiple domain formation. 

Keeping in mind that our current single-target assembly approach does not take into account the effects of a fixed number of monomers as in experiments, it is nevertheless instructive to study the effect of cooling and heating protocols.  Starting at $T=60^\circ\text{C}$, we simulate an annealing process by decreasing the temperature in discrete steps (assumed to be instantaneous) of $\delta T=0.1^\circ\text{C}$ until the temperature reaches $T=40^\circ\text{C}$. This is then followed by a melting protocol where the temperature is increased back to $T=60^\circ\text{C}$ using the same temperature step.  For each cooling/heating rate used, we produced $10^{4}$ independent trajectories, each initiated from the assembly state $x_{\mu}=\delta_{\mu \nu}$, where $\nu$ was selected randomly from the set of 334 different tile species $(\nu \in [1,334])$. In Fig.  \[fig:assembly\_rate\] we present the average assembly sizes observed at the end of each temperature step for both annealing and melting.  Also shown in Fig.  \[fig:assembly\_rate\] are the (finite) derivatives of both annealing and melting curves (inset), and the equilibrium assembly size.  For each of the cooling/heating rates used ($15\text{min}/^\circ\text{C}-3\text{h}/^\circ\text{C}$), we find differences between the peaks for assembly and melting, i.e. the system is not at equilibrium and exhibits hysteresis.  Interestingly, melting is found to occur closer to $T_c$ than assembly.  Also, we find the melting peaks to be sharper than the assembly peaks (Fig.  \[fig:assembly\_rate\] - inset), which suggests that melting has a higher degree of cooperativity than assembly. 

![ Assembly size versus time curves from isothermal simulations at different temperatures. []{ data-label="fig:trajectories"}](free_trajectories.jpg){width="45.00000%"}

![ Equilibrium assembly size (blue diamonds) as a function of temperature, and the respective size fluctuations.  The vertical bars indicate the range of sizes that are observed at equilibrium with $95 \%$ probability.  Note that the fluctuations in size are not symmetric, fluctuations are larger to smaller than average sizes.[]{data-label="fig:assembly_sizes"}](eq_sizes.jpg){width="45.00000%"}

Isothermal assembly
-------------------

We also simulated isothermal assembly by producing $10^{3}$ uncorrelated trajectories at several (constant) temperature values in the range $45^{\circ}\textnormal{C}-51^{\circ}\textnormal{C}$. In Fig.  \[fig:snapshots\] we show typical oxDNA representations of states of the assembly taken at different time slices from one of the trajectories at $51^{\circ}\textnormal{C}$. The time evolution of assembly yield, presented in Fig.  \[fig:trajectories\], shows a number of interesting features.  Firstly, at lower temperatures the size of assembly increases quickly after the trajectories are started.  However, at higher temperatures, the system exhibits increasingly longer waiting times in which the assembly does not grow beyond a few tiles, followed by a relatively quick growth towards large stable intermediates.  Note also the large range of waiting times in Fig.  \[fig:trajectories\].  For example, at the critical temperature, we measure a waiting time on the order of several years.  A small drop in temperature to $51^{\circ}\textnormal{C}$ yields much faster assembly on the order of hours.  By $49^{\circ}\textnormal{C}$ the waiting time is reduced to just a few minutes.  All the above is consistent with a nucleation barrier to assembly that rapidly decreases in size with decreasing temperature. 

At the critical temperature of $51.64^{\circ}\textnormal{C}$, the final structure has an average size of around 260 tiles, with large fluctuations (See Figs.  \[fig:trajectories\] and \[fig:assembly\_sizes\]).  As the temperature is dropped, the assembly grows to a larger final state, with smaller fluctuations.  Fully assembled structures involving all 334 strands only appear with a probability greater than 50$\%$ at temperatures of $36^\circ$C and below. 

It should be kept in mind that these calculations were performed for single target assembly in a constant background concentration of unbound tiles, and thus, structure growth is not hindered by the consumption of monomers by competing nuclei.  However, we can infer from our isothermal curves (Fig.  \[fig:trajectories\]) that there is a narrow temperature window, between $\sim 49^{\circ}\textnormal{C}$ and $51.64^{\circ}\textnormal{C}$, where monomer starvation will be avoided in an isothermal assembly scenario, as there is a clear timescale separation between waiting times and growth. 

Free energy of assembly and nucleation barriers
-----------------------------------------------

In order to further understand the behaviour observed in Figs.  \[fig:trajectories\] and \[fig:assembly\_sizes\], free-energy profiles were obtained using the biased KMC simulation method described in Section II.C. In the calculations, we used the size of the assembly $(N)$ as the biasing parameter (see Section II.C).  The results (Figs.  \[fig:free\_energy\_large\] and \[fig:CNT\_fit\]) confirm the existence of a nucleation barrier that increases quickly as the temperature approaches the critical temperature from below.  The critical nucleus size for each temperature was confirmed by using unbiased simulations to verify that for trajectories initiated at the critical nucleus $50\%$ go on towards the formation of a large assembly, whilst in the other $50\%$ the nucleus shrinks back down. 

The nucleation behaviour changes rapidly over just a few degrees.  The largest critical nucleus is $N^*=50$ at the critical temperature of $T=51.64^\circ\text{C}$, but by $T=51^\circ\text{C}$ it has already dropped to $N^*=16$. Just a few degrees lower, at $49^\circ\text{C}$, the critical nucleus is only $N^*=3$, so that the nucleation time will be very fast and no longer significantly limit the rate of assembly. 

The free energy exhibits a minimum at an assembly size that is less than that of the full 334-strand structure, and the size at this minimum grows with decreasing temperature.  This behaviour arises due to the finite size of the target structure.  Close to completion the number of possible combinations of tiles for a given assembly size decreases rapidly, reaching $1$ for the complete assembly.  This entropic effect is characteristic of finite-size assemblies and is absent, for example, in conventional crystal nucleation, where growth is solely limited by supersaturation. 

This feature of the free-energy curves explains the final assembly sizes observed in Figs.  \[fig:trajectories\] and \[fig:assembly\_sizes\].  The changing location of this minimum also explains why a thermal cooling ramp is needed for full assembly of the SST system.[@lattice_model_wang2] At temperatures where the nucleation barrier is large enough that a separation of timescales emerges between nucleation (slow) and further growth (fast), the final equilibrium structure is smaller than the fully formed structure.  At lower temperatures, where a fully formed structure is thermodynamically stable, the nucleation barrier is so small that nucleation will be fast compared with growth, leading to many partially formed structures and to monomer starvation if assembly is conducted under isothermal conditions.  Thus, an annealing process is needed that enables controlled nucleation and growth towards near completion at high temperature, followed by the growth towards completion as the temperature is decreased. 

![ Free-energy profiles as function of assembly size for different temperatures. []{ data-label="fig:free_energy_large"}](free_energy1.jpg){width="45.00000%"}

![  Higher resolution plot that shows the nucleation barriers for different temperatures, and the respective CNT fits (obtained using the free energy values for $N<200$).  The dashed vertical lines indicate the critical nucleus sizes for the nucleation barriers simulated: $N^*=3,7,16,50$ for $T=49^\circ\text{C},50^\circ\text{C},51^\circ\text{C},51.64^\circ\text{C}$, respectively. []{ data-label="fig:CNT_fit"}](free_energy1a.jpg){width="45.00000%"}

To characterize the assembly reaction pathways in more detail, biased simulations were performed at several temperatures to determine the free-energy landscape as a function of the size $(N)$ and connectivity $(E)$ of an assembly, the latter defined as the number of bonded domains in the structure.  In the calculations, we used the size of the assembly $(N)$ and the connectivity $(E)$ as the biasing parameters (see Section II.C).  As can be seen in Fig.  \[fig:free\_energy\_N\_E\], for a fixed size, higher connectivity is generally favoured, as this means a more favourable enthalpy.  But entropy also plays a role and disfavours structures with full connectivity, so the most likely connectivity typically lies just below the maximum.  In Fig.  \[fig:P\_E\_51\], we illustrate this entropic effect for the $N=16$ critical nucleus size at $T=51^\circ\text{C}$, and for the $N^*=50$ critical nucleus size at the critical temperature ($T_c=51.64^\circ\text{C}$) in Fig.  S.IV.  This entropic effect is increasingly relevant for larger nuclei, because the number of possible alternative structures with connectivity just below that of the maximum increases with size. 

For example, if we look closely at the free-energy profiles in Fig.  \[fig:CNT\_fit\], we see some small features at small nucleus size with nuclei with $4$, $6$ and $9$ strands showing deviations to lower free energy, but with the curves quickly becoming smooth beyond this point.  These sizes correspond to complete $2\times 2$, $2\times 3$ and $3\times 3$ rhomboids.  For example, as illustrated in Fig.  \[fig:cooperative\], the extra stability of the $N=6$ cluster is because the last strand that binds to form this cluster gains two binding domains in contrast to the previous or next strand to add to the growing assembly, which can only bind with one domain. 

In summary, these results show that, for all but the smallest nuclei, the most likely assembly pathways are through a set of highly connected structures with a near rhomboidal shape, but which are not the the most connected structures. 

This behaviour contrasts with previous theoretical work,[@lattice_model_wang] using a more coarse-grained approach where it was found that maximally connected structures have the lowest free energy and are the most likely to be formed, and where the free-energy profiles show much stronger features associated with sizes where cycles of connections can be completed, that persist to significantly larger sizes.  The reason for these differences with the current work is the comparative size of the free-energy barriers associated with initiating the second and subsequent domains.  In the previous patchy particle models, these barriers are significantly lower, leading to a greater stabilization for those strands that bind with two domains and hence to a stronger preference both for nuclei that are fully-connected, and for sizes for which compact, “closed” motifs are possible.  By contrast, when the full entropic cost of initiating the second and subsequent domains is properly accounted for, as we do here, the free-energetic advantage of strands that can bind with multiple domains is significantly diminished. 

![ Free-energy landscape for assembly, as a function of the size ($N$) and number of bonded domains ($E$) between tiles, at $T = 51^\circ\text{C}$. The bubbles show typical examples of assembly states representative of the respective macrostate ($N,E$).  The critical nucleus size at $T = 51^\circ\text{C}$ is $N=16$ (dashed vertical line).[]{data-label="fig:free_energy_N_E"}](free_energy2.jpg){width="45.00000%"}

![ Probability of a nucleus as a function of the number of bonded domains ($E$) between tiles for structures with $N=16$ components (critical size at $T = 51^\circ\text{C}$).  Above the curve are depicted examples of assembly states, representative of states $(N=16,E)$ with different connectivities.[]{data-label="fig:P_E_51"}](connectivity_16.jpg){width="45.00000%"}

![ Cooperative growth of a closed motif.  Closed motifs are structures where each tile is connected with the remaining through at least two binding domains, and the binding sites at the interface have a single domain available for binding.  In order for a closed motif to grow, one tile first needs to bind through a single binding domain $(i)$, that creates one binding site with two domains available for binding ($ii$-dark grey), that have a chance of being occupied if a second tile further associates, with extra stabilization effect $(iii)$.[]{data-label="fig:cooperative"}](cooperative.jpg){width="45.00000%"}

We next consider whether the main features of the nucleation barrier can be captured by classical nucleation theory (CNT), which models the competition between the addition of monomers, which lowers the free energy and the formation of a surface (or line in 2D), which costs free energy.  For a circular structure free of anisotropies, this leads to a simple expression

$$\Delta G(N)=- |\Delta \mu |N+ \Lambda N^{1/2},
 \label{eq:CNT_fit}$$

where $N$ is the number of monomers forming the structure, $\Delta \mu <0$ is the chemical potential difference between monomers forming the structure and monomers in the reservoir, $\Lambda=2\pi^{1/2}\rho_{s}^{-1/2}\lambda $ where $\lambda $ is the line tension, and $\rho_{s}$ is the number density of monomers forming the structure. 

To estimate the number density of monomers we simulate a fully assembled structure (Fig.  \[fig:membrane\]) using molecular dynamics simulations with oxDNA2,[@snodin2015introducing] an extension of the original oxDNA model that allows for a better description of the structural properties of large (kilobase-pair) DNA nanostructures.  The simulations were performed at $T = 25^\circ\text{C}$ and after a period of equilibration, we measured the average distance between the centres of mass of contiguous tiles.  We then estimate the area occupied by each tile by assuming that contiguous tiles can be considered to be contained in the same plane.  Using this method, we obtain the value $\rho_{s}=5.82\times 10^{16}\, \textnormal{m}^{-2}$ for the surface density of monomers. 

We use the CNT expression to fit our free-energy curves for assemblies with up to 200 monomers; we exclude larger sizes as CNT is not able to capture the entropic effects that lead to a minimum in the free energy near to the fully-assembled state.  Note also that in order for the CNT fit to coincide with our free-energy curves at $N=1$ (where we defined $G(1)=0$ for all temperatures) we use the fitting variable $N'=N+1$ in Eq.   \[eq:CNT\_fit\]. 

In Fig.  \[fig:CNT\_fit\], we show that the nucleation barrier for several temperatures is fit well by the classical CNT expression with the values of the parameters $|\Delta \mu |$ and $\lambda$ given in Table \[table:1\].  Of course this simple expression will not capture the size-specific features observed for very small nuclei, and the minimum at larger sizes, but it is nevertheless remarkable how well it fits the overall curve.  As is found for many other systems, the surface (line) tension does not change much with temperature.  Instead, the rapid drop in the barrier height with temperature is mainly caused by the rapid change in $\beta \Delta \mu$.

  ---------------------- ----------------------- ------------------------------
   $T (^\circ\text{C})$   $\beta |\Delta \mu |$   $\lambda (\textnormal{pN}) $
      \[0.5ex\] $49$              1.26                        1.42
           $50$                   0.99                        1.61
           $51$                   0.67                        1.51
     $51.64\, (T_c)$              0.36                        1.58
         \[1ex\]                                 
  ---------------------- ----------------------- ------------------------------

  : CNT fit parameters

\[table:1\]

Although experimental measurements on the line tension of DNA assemblies have not to our knowledge been made, it is interesting that predictions are fairly similar to typical experimental measurements on various lipid bilayers which range from ($\lambda =1\, \textnormal{pN}$[@line_tension_exp1] to $\lambda =4\, \textnormal{pN}$[@line_tension_exp2]), as well as predictions using Monte-Carlo simulations for a coarse grained lipid model ($\lambda =3\pm 2\, \textnormal{pN}$[@line_tension_mod1]). 

\[sec:level1\]Conclusions
=========================

In this work, we develop a two-step coarse-graining approach that uses extensive thermodynamic calculations with a nucleotide-level model of DNA, to parametrise a coarser kinetic model that is able to reach the time and length scales needed to describe explicitly the self-assembly of a two-dimensional $334$-strand SST structure first produced experimentally by the Yin group.[@PengYin2d] The resulting kinetic model provides a detailed description of the assembly trajectories and, without using adjustable parameters, it predicts a transition temperature which is consistent with previous experimental observations.  The model used to calculate the thermodynamics of tile association, oxDNA, describes explicitly the polymeric degrees of freedom of DNA, which makes it more suitable to capture the entropic penalties associated with initiation of a second and subsequent domains during tile association than previous patchy-particles models.[@lattice_model1; @lattice_model2]

For the SST system, we show that there exists a narrow temperature regime where the nucleation barrier is low enough to be crossable on experimental time-scales, but large enough to allow a time-scale separation between growth and nucleation.  At the same time, the exponentially large design space for DNA means that misbonds do not overwhelm correct assembly at experimentally relevant temperatures. 

We find that the critical nuclei are made up of an ensemble of structures that are highly, but typically not maximally, connected.  The nucleation barrier near the temperature where assembly first occurs is surprisingly well described by CNT, with the fits providing reasonable estimates for the surface (line) tension and chemical potential difference.  At lower temperatures and for smaller nuclei, the full model exhibits structure due to cooperative effects not included in CNT.  In contrast to systems with only a few components, the only way to achieve full assembly is to use a cooling protocol, where the system first crosses the nucleation barrier to a partially formed assembly, and further cooling stabilises the full assembly.  While our calculations are quantitative for a specific 334 strand SST system, a number of the basic qualitative results above are also found with simpler models[@lattice_model1; @lattice_model2; @Aleks_off_lattice; @lattice_model_wang2; @lattice_model_wang; @Jacobs_protocol_design] and are likely to generically hold for a wider range of target structures[@PengYin3d; @Peng_Yin_Dephts; @Peng_yin_10nm; @Peng_yin_tubes; @Peng_yin_isothermal; @Peng_yin_design_space; @Peng_yin_reconfigurations] that use the addressable SST assembly method. 

There are a number of possible extensions of our calculations.  A particular important direction will be to perform simulations in an NVT ensemble which would allow for a more careful comparison to collective properties (such as yield) that are measured in experiments, in particular the role of monomer starvation in reducing the yield at higher cooling rates could be quantitatively elucidated.  Extending our kinetic model to consider this situation is in principle straightforward.  Future tests to our approach may include a careful compassion between the distribution of structure sizes observed in the experimental samples and that predicted by the kinetic model for a similar background concentration of monomers. 

Making different shapes, as done in Ref.   by leaving out certain tiles, would also be straightforward to investigate, and it would be interesting to see how the nucleation pathway might change for a more complex shape.  Extending the method to consider SST structures in 3 dimensions is also feasible, but would of course require the detailed thermodynamics of the different local environments relevant to this new geometry to be characterized. 

In this paper we used an average-base model that enables us to look at the generic behaviour of these SST systems.  However, it should be relatively uncomplicated to extend the model to investigate a specific set of sequences by using the NN model to estimate how the relevant $\Delta G^{{\ensuremath{-\kern-4pt{\ominus}\kern-4pt-}}}$ would differ from that calculated for the average-base model.  The model could then be used to study how sequence might affect the nucleation pathway - one might imagine that nucleation is more likely to occur in regions with higher GC content - and how it could be used to optimise the assembly of particular structures. 

Our methodology, where oxDNA is used to parameterise the rates for a limited set of geometries, which are then fed into a KMC scheme, can also be extended to other DNA tile systems such as 2D double-crossover tiles[@winfree1998design] which are widely used for studying molecular computation. 

Lastly, we point out that our modelling approach is based upon the assumption that interactions relevant for tile-incorporation are typically local, and to a good approximation only involve components in the local neighbourhood of the binding site; as a consequence our modelling strategy may not be adequate to address self-assembling systems where relevant interactions are typically long ranged.  In particular, the description of DNA origamis, in which strand incorporation typically involves bringing into close proximity distant sections of the circular scaffold strand (and is further complicated by other strands already incorporated in the structure[@dannenberg2015modelling]) is at the present beyond the scope of our multi-scale modelling approach. 

supplementary material {#supplementary-material .unnumbered}
======================

See supplementary material for further details on the simulations performed with oxDNA (S.I and S.II), misbond formation calculations (S.III and S.V), the structure of the critical nucleus at $T_c=51.64^\circ\text{C}$ (S.IV), the full set of 32 local geometries and respective free energy profiles (S.VI), an example of a bridging tile (S.VII), the tests performed to gauge the effect of detailed balance deviations (S.VIII and S.IX), and additional details on the experimental system considered to parametrise $k_+^0$ (S.X). 

The authors wish to thank the EPSRC for financial support (EP/I001352/1 and EP/J019445/1).  PF thanks FCT-Portugal for funding through grant SFRH/BD/94405/2013.  TEO is supported by a Royal Society University Research Fellowship.  We acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out some of this work (doi:10.5281/zenodo.22558). 

[60]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty [****,  ()](http://dx.doi.org/10.1038/nmat4184) [****,  ()](http://dx.doi.org/10.1038/nmat4178) [****,  ()](http://dx.doi.org/10.1021/jacs.5b11918) [****,  ()](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC5027514/) [****,  ()](http://dx.doi.org/10.1038/nature04586) [****,  ()](http://dx.doi.org/10.1038/nature11075) [****,  ()](http://dx.doi.org/10.1016/j.chempr.2017.02.009) [****,  ()](http://science.sciencemag.org/content/338/6111/1177.abstract) [****,  ()](http://dx.doi.org/10.1038/nchem.2083) [****,  ()](http://dx.doi.org/10.1002/smll.201501370) [****,  ()](http://science.sciencemag.org/content/321/5890/824.abstract) [****,  ()](http://dx.doi.org/10.1021/nl4019512) [****,  ()](http://pubs.acs.org/doi/abs/10.1021/ja4062294) [****, ()](http://dx.doi.org/10.1002/ange.201402437) [****,  ()](http://dx.doi.org/10.1038/nature24648) [****,  ()](http://www.sciencedirect.com/science/article/pii/S0042682205003429) [****,  ()](https://doi.org/10.1063/1.3243581) [****,  ()](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4318123/) [****,  ()](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4382372/) [****,  ()](http://link.aps.org/doi/10.1103/PhysRevLett.112.238103) [****, ()](http://dx.doi.org/10.1039/C5FD00135H) [****,  ()](http://dx.doi.org/10.1039/C6SM01031H) [****,  ()](http://www.pnas.org/content/112/20/6313.abstract) [****,  ()](http://scitation.aip.org/content/aip/journal/jcp/142/2/10.1063/1.4905670) [****,  ()](http://dx.doi.org/10.1039/C5SM01841B) [****,  ()](http://dx.doi.org/10.1038/ncomms7203) [****,  ()](http://aip.scitation.org/doi/full/10.1063/1.4921957?crawler=true) [****,  ()](http://aip.scitation.org/doi/full/10.1063/1.3552946) [****,  ()](http://dx.doi.org/10.1039/C3CP53545B) [****,  ()](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.104.178101) [****,  ()](http://pubs.acs.org/doi/abs/10.1021/nn3058483) [****,  ()](https://link.springer.com/article/10.1007/s11047-013-9391-8) [****,  ()](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3905871/) [****,  ()](https://www.nature.com/articles/ncomms6324) [****, ()](http://aip.scitation.org/doi/full/10.1063/1.4792252) [****,  ()](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3799446/) [****,  ()](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4513862/) [****,  ()](http://pubs.acs.org/doi/abs/10.1021/acsnano.5b05865) [****,  ()](http://science.sciencemag.org/content/338/6113/1458.abstract) [****,  ()](http://www.pnas.org/content/95/4/1460.abstract) [****,  ()](http://dx.doi.org/10.1038/nature14860) [****, ()](http://aip.scitation.org/doi/full/10.1063/1.4933426) [****,  ()](http://pubs.acs.org/doi/abs/10.1021/ja906987s) [ ()](http://dx.doi.org/10.1038/nchem.2877) [****,  ()](http://pubs.acs.org/doi/abs/10.1021/bi00063a022) [****,  ()](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1888818/) [****, ()](http://pubs.acs.org/doi/abs/10.1021/jp510061f) [****,  ()](https://www.nature.com/articles/ncomms10803) [****, ()](http://aip.scitation.org/doi/full/10.1063/1.4754132) [****,  ()](http://www.pnas.org/content/108/13/5193.abstract) [****, ()](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3516813/) [****,  ()](http://www.sciencedirect.com/science/article/pii/0021999177901218) [****,  ()](http://dx.doi.org/10.1002/jcc.540130812) [****,  ()](http://stacks.iop.org/0953-8984/22/i=10/a=104102) [****,  ()](http://dx.doi.org/10.1021/j100540a008) [****,  ()](https://doi.org/10.1016/j.shpsc.2015.12.014) [****,  ()](http://dx.doi.org/10.1021/jp001299a) [****,  ()](http://www.sciencedirect.com/science/article/pii/S0006349508705216) [****,  ()](http://www.sciencedirect.com/science/article/pii/S000930840500054X) [****,  ()](http://www.nature.com/nature/journal/v394/n6693/abs/394539a0.html)
