Introducing Improved Structural Properties and Salt Dependence into a Coarse-Grained Model of DNA
We introduce an extended version of oxDNA, a coarse-grained model of DNA designed to capture the thermodynamic, structural and mechanical properties of single- and double-stranded DNA. By including explicit major and minor grooves, and by slightly modifying the coaxial stacking and backbone-backbone interactions, we improve the ability of the model to treat large (kilobase-pair) structures such as DNA origami which are sensitive to these geometric features. Further, we extend the model, which was previously parameterised to just one salt concentration (), so that it can be used for a range of salt concentrations including those corresponding to physiological conditions. Finally, we use new experimental data to parameterise the oxDNA potential so that consecutive adenine bases stack with a different strength to consecutive thymine bases, a feature which allows a more accurate treatment of systems where the flexibility of single-stranded regions is important. We illustrate the new possibilities opened up by the updated model, oxDNA2, by presenting results from simulations of the structure of large DNA objects and by using the model to investigate some salt-dependent properties of DNA.
Deoxyribonucleic acid (DNA) performs the crucial function of storing genetic information in living organisms. It is made up of repeating units called nucleotides, each of which consists of a sugar and phosphate backbone plus a base (either Adenine (A), Guanine (G), Thymine (T), or Cytosine (C)) attached to the sugar. Watson-Crick base-pairing, A with T and G with C, along with planar stacking interactions between bases and the constraints of the backbone, leads to the formation of the well-known double-helical structure of DNA.
The Watson-Crick complementarity of DNA also permits the rational design of DNA objects for which the intended structure is the global free-energy minimum, a property which has been exploited to create a wide variety of 2D and 3D nanostructures.Dietz et al. (2009); Han et al. (2013); Ke et al. (2012) Further, these DNA objects can be functionalised,Walsh et al. (2011); Douglas et al. (2012); Kuzyk et al. (2012); Zadegan et al. (2012) with potential applications ranging from nanomedicine to nanoelectronics.Pinheiro et al. (2011)
Theoretical and computational approaches to modelling DNA have been widely exploited to probe the behaviour of the molecule in both a biological and a nanotechnological context. At the finest level of detail, quantum chemistry calculations have been used to study the interactions between nucleotides,Svozil et al. (2010); Mladek et al. (2013); Šponer et al. (2013) although the high computational cost of such an approach limits these methods to interactions between nearest-neighbour base pairs in vacuum. Classical all-atom approaches, where every atom of DNA and the surrounding solvent is modelled as a point particle with effective interactions, have been widely employed to study small DNA motifs,Laughton and Harris (2011); Pérez et al. (2012) and have recently been applied to larger DNA systems.Maffeo et al. (2012); Yoo and Aksimentiev (2013); Maffeo et al. (2014) However, simulating rare-event processes such as the breaking or formation of base pairs remains a challenge with these models, with time scales being the limit of what is currently accessible.Maffeo et al. (2014) At the other end of the scale, theoretical approaches have been developed to understand certain large-scale properties of DNA. These include the wormlike-chain model, which treats DNA as a continuously flexible polymer.Bustamante et al. (1994) While such models can provide useful insights into the physical properties of DNA, they are not detailed enough, by design, to address processes such as duplex formation.
The middle ground between analytical and all-atom approaches is occupied by coarse-grained DNA models. Such models integrate out many of the degrees of freedom of the DNA nucleotide, and often neglect the solvent molecules; these approximations inevitably imply a compromise between accuracy, generality and computational efficiency, so that care must be taken in applying these models to a given problem. However, the simplified picture presented by such an approach can be a strength, as, in addition to greatly increasing the time scale and number of nucleotides that can be studied, it can allow one to understand the generic physics governing the system more easily.
Many coarse-grained DNA models have been developed in recent years;Hinckley et al. (2013); Savelyev and Papoian (2010); Cragnolini et al. (2013); He et al. (2013); Linak et al. (2011); Morriss-Andrews et al. (2010); Araque et al. (2011); Korolev et al. (2014); Maffeo et al. (2014) in this work we improve the model due to Ouldridge et al., oxDNA.Ouldridge et al. (2011); Ouldridge (2011) This nucleotide-level model is designed with a heuristic, “top-down” approach, with a focus on reproducing well-known properties of DNA (such as the helical structure of the B-DNA duplex) and experimental results (such as duplex melting temperatures), rather than, for example, building the model up by integrating out details from an all-atom representation. In the oxDNA model, each nucleotide is a rigid body with three interaction sites that have mutual, highly anisotropic interactions. This treatment is sufficient to obtain good agreement with experimental data on the structural, mechanical and especially the thermodynamic properties of single- and double-stranded DNA. Consequently, the model has provided key insights into many different processes relevant to DNA nanotechnologyOuldridge et al. (2010); Šulc et al. (2013); Ouldridge et al. (2013, 2013); Rovigatti et al. (2014); Schreck et al. (2014, 2014); Srinivas et al. (2013); Mosayebi et al. (2015) and biophysics,Matek et al. (2012); Wang and Pettitt (2014); Matek et al. (2015); Mosayebi et al. (2014); Romano et al. (2013) and importantly, has also been shown to provide direct agreement with experimentally measured properties on a range of systems including DNA overstretching,Romano et al. (2013) a two-footed DNA walker,Ouldridge et al. (2013) and toehold-mediated strand displacement.Srinivas et al. (2013); Machinek et al. (2014)
Despite these achievements, there are some areas where oxDNA can be improved. The model was parameterised to , a high salt concentration similar to that used for many applications in DNA nanotechnology. However, the ability to study DNA behaviour as a function of salt would allow quantitative comparison with a greatly expanded set of experiments, and, in particular, if we wish to apply oxDNA to biological systems we would like to work at physiological salt (). At the same time, the wealth of experimental data for the thermodynamic and mechanical properties of DNA as a function of salt concentrationBosco et al. (2013); Polinkovsky et al. (2014); SantaLucia and Hicks (2004); Taylor and Hagerman (1990) makes fitting the salt-dependent properties of an extended model possible.
While the detailed effect of salt on DNA electrostatics can be highly complex,Mocci and Laaksonen (2012); Potoyan et al. (2013) here we use a simple Debye-Hückel interaction term, as first implemented for oxDNA by Wang and Pettitt in Ref. Wang and Pettitt, 2014, to model how salt screens repulsive interactions. The more coarse-grained description is commensurate with the level of approximation used more generally in the oxDNA framework. We carefully parameterise the model to an extended set of experimental data for melting temperatures and persistence lengths, including the behaviour of single strands.
A second area that merits attention is the performance of the model in simulating the structure of large (kilobase-pair) DNA objects. The model value of the B-DNA pitch did not come under too much scrutiny in the original parameterisation as there is some disagreement in the literature about the precise value of the average pitch.Taylor and Hagerman (1990); Shore and Baldwin (1983); Dawid et al. (2006); Ido et al. (2013); Du et al. (2005); Wang (1979) However, improved experimental techniques in fabrication and imaging of large DNA objectsKato et al. (2009); Dietz et al. (2009); Bai et al. (2012) have presented an opportunity to finely tune this value, because small adjustments to the duplex pitch can result in significant changes to the global twist of a large-scale DNA nanostructure. In addition, simulating these large-scale structures has illustrated the potential importance of nicks and junctions for the effective duplex pitch, so that improving the model’s description of these effects has become a priority. Finally, the original model duplex had grooves with equal widths, whereas B-DNA is known to have a larger major groove and a smaller minor groove. This implies that the positions of the backbone sites in the model, which are directly related to the groove widths, could be more realistic. This detail could be relevant, for example, in origami structure, where the precise backbone positions are known to be important for junction placement.Dietz et al. (2009)
The oxDNA model was previously given sequence-dependent hydrogen-bonding and stacking strengthsŠulc et al. (2012) by fitting to the duplex thermodynamics of the SantaLucia model.SantaLucia and Hicks (2004) As the SantaLucia model gives results at the base-pair step level, one can only extract the average stacking strength for the two stacking interactions present in a given base-pair step. In particular, this means that in oxDNA the AA and TT stacking interactions are the same, whereas it is well known that the AA stacking interaction is significantly stronger than the TT interaction, an important property for DNA nanotechnology, for example, where poly-T single-stranded regions are often used as flexible linkers.Andersen et al. (2008); Tian et al. (2013) To remedy this we use new experimental data to reparameterise the AA and TT stacking interactions in the model.
In the following sections, after briefly introducing the original oxDNA model, we consider each change to the model in turn, and then we highlight some important aspects of the behaviour of the new model.
Ii The original oxDNA model
In the oxDNA model introduced by Ouldridge et al.Ouldridge et al. (2011) strands of DNA are represented by a chain of rigid bodies, with each rigid body representing a nucleotide. The coarse-grained, pairwise potential for the model can be written as a sum of interaction terms
The nearest neighbours (adjacent nucleotides on a DNA strand) interact via , , and , which represent the connectivity between neighbouring backbones, the favourable stacking interactions between neighbouring bases, and excluded volume terms, respectively. All other nucleotide pairs interact with , , and , corresponding to hydrogen bonding between complementary bases, cross-stacking, excluded volume, and coaxial stacking between non-nearest neighbours, respectively. We now highlight aspects of the original model that are relevant to the improvements made in this paper.
Firstly, in its original formulation, the model was parameterised for a sodium ion concentration of , chosen to reflect the high salt concentrations commonly used in DNA nanotechnology. The electrostatic interactions between the negatively charged phosphate groups on the DNA backbone were incorporated into the backbone site’s excluded volume, an approximation which can be justified by the very short Debye screening length at that relatively high ion concentration.
Secondly, the original oxDNA modelOuldridge et al. (2011) represented each nucleotide as a linear rigid body (Fig. 1). The optimal configuration for base-pairing occurs when the two nucleotides point directly at each other. As a consequence, the DNA double helix was symmetric, with the two grooves having equal widths.
Thirdly, in the original oxDNA model introduced in Ref. Ouldridge et al., 2011, all four types of base were treated equally except that only A-T and G-C base pairs could be formed. Later Šulc et al.Šulc et al. (2012) introduced sequence-dependent thermodynamics into the model by making the strengths of the hydrogen-bonding and stacking terms depend on the identities of the interacting nucleotides. The nearest neighbour DNA model of SantaLucia,SantaLucia and Hicks (2004) to which oxDNA was parameterised, does not resolve the difference between AA and TT stacking as it works on the base-pair step level – therefore AA and TT stacking strengths were set to be the same in Ref. Šulc et al., 2012.
In this work we mostly work from the original, sequence-averaged parameterisation of the model rather than the sequence-dependent one, as it is more efficient to fit the thermodynamic parameters to sequence-averaged duplex melting temperatures as given by the SantaLucia model. The exception is the parameterisation of the AA and TT stacking strengths, which did use the sequence-dependent parameters from Ref. Šulc et al., 2012 as a starting point, to allow the best possible comparison between the model and the experimental results that were used for the fitting. After the parameters for oxDNA2 had been obtained, including new values for the sequence-averaged hydrogen bonding and stacking strengths, we then rescaled the sequence-dependent interaction strengths from Ref. Šulc et al., 2012 accordingly for use with the new model.
Iii Introducing different widths for major and minor DNA grooves
B-DNA in the original oxDNA model has equal groove widths, while in reality DNA has a larger major groove and a smaller minor groove. Having realistic widths for the major and minor grooves is equivalent to having appropriately positioned backbone sites in the model, an important feature for the physical properties of many DNA motifs. For example, in DNA origami antiparallel double helices are joined by crossovers, for which the position of the backbone has been shown to be crucial for origami structure.Dietz et al. (2009); Douglas et al. (2009) Another example is anisotropic duplex bending: the duplex can be expected to bend more easily into the major groove than into the minor, if the groove widths are unequal.
The oxDNA nucleotide is composed of three interaction sites: the hydrogen-bonding, stacking and backbone sites. We introduce different groove widths by changing the position of the backbone site while keeping the duplex radius unchanged (Fig. 1), such that, rather than lying on a straight line, the three interaction sites lie in a plane. The new nucleotide shape introduces an additional parameter into the model, the angle between the line from the duplex centre to the backbone site and the line from the duplex centre to the stacking site (Fig. 1(a)). Given the coarse graining of the atoms of the sugar-phosphate DNA backbone into a single interaction site, there is no definitive choice for the precise position of the backbone site and thus the value of the model parameter (Fig. 1(a)). We set , a value which maps onto a full-atom representation of a DNA duplex well by visual inspection, although values of between and would give an equally satisfying visual match.
The backbone site is moved such that the duplex radius is unchanged, and we note that the modification has a negligible added computational cost when simulating the model. However, the thermodynamic and mechanical properties are slightly affected. For the thermodynamics, we found a change of 1-2K in the duplex melting temperatures, and we modified the hydrogen-bonding and stacking strengths using the histogram reweighting method described in Section II C 2 of Ref. Šulc et al., 2014, so that the agreement with experimental melting temperatures was as good as for the original model. The mechanical properties of DNA are less well constrained experimentally and so were not refitted. The mechanical properties for the new model can be found in Section VIII.
One illustration of the importance of the groove widths in oxDNA for the structural properties of DNA assemblies is provided by systems of 3-arm star tiles that are designed to form triangular prismatic polyhedra. We find that modifying the groove widths qualitatively changes the structure of trimers of these tiles (Fig. 2). Specifically, the body of the trimer defines a plane with two distinct faces. Zhang et al. Zhang et al. (2012) found that one of two possible isomers of the polyhedron preferentially formed, implying that the free arms of the trimer systematically pointed in the direction of one of these two faces. We find a consistent result when the groove widths specified by oxDNA2 are used. When equal-sized grooves are used (as in the original oxDNA model) the trimer arms point in the opposite direction.
Iv Effective electrostatic interactions
One major improvement presented in this paper is the introduction of a salt-dependent interaction term in the oxDNA model. Since the original goal of the oxDNA model was to simulate nanotechnology experiments, the thermodynamic and structural parameterisation was carried out at high salt. The very short electrostatic screening length at these conditions allows one to incorporate the electrostatics into a soft excluded volume, somewhat circumventing the necessity of a proper treatment. The original parameterisation of the thermodynamics was carried out at 0.5 M [Na], a high enough value that further increasing it does not significantly change the physics at our level of coarse-graining.
The problem of treating electrostatics properly for DNA in solution is a complicated one. Perhaps the most evident issue is that the typical dimensions of nucleotides are comparable to the Debye length of the solution, rendering a mean-field treatment hard to justify. Also, in some cases the presence of salt ions affects the local structure of nucleic acids, by stabilizing some arrangements or destabilizing others. Ion condensation may also lead to stronger screening of the electrostatic interactions than otherwise expected – this has been incorporated into coarse-grained DNA models through partial effective charges.Hinckley et al. (2014); Maffeo et al. (2010) Thus, in principle, many non-trivial effects must be taken into account when modelling the electrostatic interactions, and the debate on the best way to do so implicitly is still unresolved.Mocci and Laaksonen (2012); Potoyan et al. (2013)
Here, we choose a very simple treatment, based on the Debye-Hückel model for screened electrostatics. This approach has been used previously for other coarse-grained models,Hinckley et al. (2014); Morriss-Andrews et al. (2010) and was first introduced into oxDNA by Wang and Pettitt.Wang and Pettitt (2014) We note that this treatment is consistent with the coarse-grained nature of the model, and we use the same top-down strategy that was used in the original parameterisation to design the effective electrostatic interactions: our goal is to introduce a term in the potential that will reproduce the thermodynamic and mechanical effects of salt concentration on DNA, and thus should be regarded as an effective interaction rather than an attempt to rigorously model the local effects of charges. This should be kept in mind when interpreting the results obtained with oxDNA, in particular at low salt concentration.
Since the modelling of electrostatics is rather crude, we restrict our parameterisation to salt concentrations of M of monovalent salt or greater. This restriction is also due to the fact that we parameterise our thermodynamics to the model of SantaLucia,SantaLucia and Hicks (2004) which was fitted in a similar salt regime. Importantly, physiological conditions fall within this range, which will allow quantitative comparison between simulations of DNA systems with our model and experiments at physiological conditions.
The detailed Debye-Hückel form which is added to the non-bonded interactions in the potential of Eq. 1 takes the following form:
where is the effective charge situated at the backbone site of each of the nucleotides, is the magnitude of the distance between the backbone sites of nucleotides and , is the permittivity of the vacuum, is the relative permittivity of water, is the elementary charge and is a dimensionless effective charge. In principle, depends on ,Mazur and Jernigan (1991) and weakly depends on temperature and salt concentration. However, for oxDNA2 we set to be a constant value, in keeping with the coarse-grained approach taken for the rest of the model. In particular, we choose , the standard value for water. In Eq. 2 we have stressed that the interaction depends on the temperature and on the (monovalent) salt concentration through the Debye length :
where is Avogadro’s number and is Boltzmann’s constant.
To ensure the computational efficiency of simulating the model, we set the interaction to zero at a finite distance; and to allow simulation with molecular dynamics, we introduce a quadratic smoothing potential so that the interaction goes to zero smoothly. The quadratic smoothing, the details of which are reported in Appendix A.1, is introduced after a cutoff , which we choose to be . This cutoff allows us to use all the standard techniques to improve the simulation efficiency via the use of Verlet lists and/or cells.Frenkel and Smit (1996) We have checked that introducing our chosen cutoff, , has a negligible effect on the duplex thermodynamics results used to parameterise the interaction (Fig. 3(a)).
Our representation of DNA uses a single rigid body per nucleotide, and the best choice of where to put the charge is not obvious. All the atoms of the sugar and phosphate groups of the backbone are represented in a single interaction site, and it is thus natural to put the charge, which in real DNA is located on the phosphate, on that interaction site. Importantly, the backbone site of a nucleotide is placed almost in between the phosphate of that nucleotide and the phosphate of the neighbouring one, which could potentially lead to some unphysical effects. Also, we should stress that having a charge at each backbone site means that the DNA has as many charges as nucleotides, which is not always true in real systems: very often, the terminal phosphate at the end is cut off, removing a charge. The absence of this charge can cause measurable effects on the thermodynamics, and indeed the SantaLucia modelSantaLucia and Hicks (2004) requires the presence or absence of the terminal phosphates as an input parameter. In keeping with our coarse-graining approach, we put a half effective charge on each of the terminal nucleotides to incorporate the fact that each charge should be halfway in between our backbone sites, emulating a system with the terminal charge removed, and parameterise to the SantaLucia model in a way consistent with this approach.
The parameter that we have tuned to reproduce the thermodynamics predicted by the model of SantaLucia is the effective charge . To do this, we used thermodynamic integration Vega et al. (2008) to compute the melting temperatures of duplexes of length 5, 6, 7, 8, 10, 12 and 15 as a function of at several salt concentrations, and chose the value that best reproduced the melting temperatures predicted by SantaLucia’s model.
The melting temperature of a duplex is defined as the temperature at which half of the strands in a stoichiometric bulk solution are in the duplex state. We cannot simulate the bulk system with oxDNA. However, at the bulk melting temperature for a given concentration of strands, the bound and unbound free energies obtained for a system of two strands at the same concentration, and , are related by:Ouldridge et al. (2010)
where we have made explicit the dependence of on the effective charge . The constant on the right hand side of Eq. 4 accounts for concentration fluctuations that are present in bulk but not for two strands in a periodic box.Ouldridge et al. (2010) The equation is exact in the dilute limit, which is indeed where the thermodynamics of duplex formation are usually studied and is the relevant limit for oxDNA.
We can take advantage of the relation of the free energies given by Eq. 4 for the thermodynamic integration. For small changes in and , it is possible to write
The last step includes the unusual term . This is because our potential, like many coarse-grained models, depends explicitly on through and .
To compute the change in introduced by a small change in , one can impose the condition in Eq. 4 at the new
Eq. 7 is a differential equation that allows us to follow the change in melting temperature as is changed. We note that Eq. 7 is an extended Clausius-Clapeyron relation,Vega et al. (2008) and the quantities on the right-hand side are readily accessible with separate simulations of the bound and unbound states. Since we are dealing with small systems, all these simulations are very quick and it is thus easy to achieve very accurate results. As a starting point for the thermodynamic integration we use the melting temperatures from SantaLucia, which our original model (equivalent to the current model at ) reproduces within tenths of a Kelvin. Since both and only weakly depend on temperature, we used Eq. 7 assuming to obtain an optimal value for , and checked a posteriori the results of thermodynamic integration with melting simulations. We have found that the melting temperature differences computed with Eq. 7 are accurate to within K.
For a range of salt concentrations and duplex lengths, we performed thermodynamic integration to find the melting temperature as a function of . Some of the results are depicted in Fig. 3. The melting temperature of the duplex decreases with increasing effective charge , but not dramatically so (Fig. 3(a)). This is because the increased interchain repulsion in the duplex state is partially balanced by the lower entropy of the single-stranded state (due to it adopting a slightly more extended state to reduce intrachain repulsion). The data in Fig. 3(b) show that the at which the difference in melting temperature between oxDNA and SantaLucia is minimised depends on the duplex length and lies in the range . The best overall predictions are obtained for an effective charge . It is reassuring that this value does not deviate significantly from (corresponding to the value given by Debye-Hückel theory), and that the best value for varies little with duplex length. In addition, it is not uncommon to use a value of for coarse-grained DNA models.Hinckley et al. (2014); Maffeo et al. (2010) One argument is that ion condensation, which is known to occur for DNA, will screen the phosphate charges more strongly than expected from Debye-Hückel theory. This will lead to a lower effective charge when fitting a model using a Debye-Hückel treatment to experimental results, although such arguments should be applied with caution to a crude mean-field approach such as this one.
We note that introducing this explicit electrostatic term in our model potential will raise the computational expense of each simulation step compared to the original oxDNA model, as the electrostatic term will generally result in oxDNA2 having a larger interaction range than the original model has. This effect increases at lower salt, as the Debye-Hückel term becomes more long-ranged. For example, we find that simulating a 10-bp duplex with MD for a given number of steps takes as long with oxDNA2 at [Na]=0.5 M and as long with oxDNA2 at [Na]=0.1 M as it did for the original oxDNA.
V Improving structure prediction for large-scale DNA objects
The model was originally parameterised for small single- and double-stranded DNA structures. If we wish to study large-scale structures, we need to ensure that we can reproduce the existing experimental data for these larger constructs. A good test case is provided by the work of Dietz et al.,Dietz et al. (2009) who measured the global twist of three different DNA origami structures, described as 10-by-6 helix bundles (Fig. 4). We denote the three origami structures as L-, N- and R-type. In the experiment, they were designed to impose a pitch of 10, 10.5, and 11 base pairs per turn (bp/turn) on the constituent DNA double helices and these different designs exhibited left-handed, no, and right-handed global twist, respectively, when multimerised to form ribbons and visualised with transmission electron microscopy. One might think that this result implies that DNA has a natural pitch of 10.5 bp/turn, since the design with that inherent periodicity did not result in a globally twisted system. However, it is not that straightforward.
|Helix Bundle Type||oxDNA2||experiment|
|L||-22.2||-31 +/- 5|
|R||26.9||26 +/- 5|
Unsurprisingly, when simulated with the original version of oxDNA, the N-type origami still showed a significant right-handed global twist, chiefly because the duplex pitch for the model was 10.36 bp/turn. However, even when simulated with a version of oxDNA modified so that the model duplex pitch was 10.5 bp/turn, the N-type origami displayed a right-handed twist. One reason for this is the following: in oxDNA, the helical twist across nicks and junctions is larger by about and degrees respectively than for a normal duplex step (we call this an overtwist), so that fractionally fewer base pairs are required for a helix turn than would otherwise be expected. Although these differences are small, they add up constructively to create a global twist on a structure that would otherwise have no such twist. For a nicked duplex in oxDNA, the overtwist occurs because, opposite the nick, the non-nicked strand prefers a larger twist than for duplex DNA (just as stacked single strands do in the model), and the nicked backbones lack the FENE spring that would usually oppose this tendency. A similar argument applies for strands with junctions.
Although there is no direct evidence available to show whether this overtwist is physically realistic, there are multiple lines of experimental evidence that suggest that the pitch of duplex DNA is close to 10.5 bp/turn.Taylor and Hagerman (1990); Shore and Baldwin (1983) Given the evidence from Dietz et al.Dietz et al. (2009) that the effective pitch in DNA origami structures is also close to 10.5 bp/turn, we decided to reduce the overtwisting in oxDNA as much as is possible. To achieve this, we modify the coaxial stacking term of the potential, , so that the overtwist is for a junction and for a nick (see Appendix A.3 for details of the changes to the potential). We have verified that this change has very little effect on the other features of the model.
Even when the overtwist at nicks and junctions has largely been removed, we require an intrinsic duplex pitch of bp/turn in order for the oxDNA N-type helix bundle to have zero global twist, and we set the pitch to this value (at [Na]=0.5 M) for oxDNA2 by modifying . We found that this modification to the backbone potential changed the duplex melting temperatures by around 2 K, and we refitted the model’s hydrogen bonding and stacking strengths to correct for this using the same method as was used for the thermodynamics refitting when implementing unequal helical groove widths. The requirement for a pitch of bp/turn, rather than the bp/turn that one might expect, is due to a few subtle effects that are currently being investigated further.
The global twist measured for the helix bundles with the new oxDNA2 model (Appendix B.1 describes how the global twist was measured from simulations) is compared to the experimental results in Table 1, while typical simulation snapshots are shown in Fig. 4. The slight modifications to used to set the model pitch are given in Appendix A.2. We note that the modifications to changed the duplex melting temperatures in the model by 1-2 K, which were then refitted using histogram reweighting (as described in Section II C 2 of Ref. Šulc et al., 2014) to give an agreement with experimental melting temperatures that was as good as for the original model.
Vi AA/TT sequence dependence
As a further improvement, the oxDNA2 model incorporates a more realistic sequence-dependent stacking interaction, which is achieved by differentiating between the AA and TT stacking interaction strengths. In the previous parameterisation of oxDNA,Šulc et al. (2012) the sequence-dependent base-pairing and stacking interaction strengths were obtained by fitting the oxDNA duplex melting temperatures to the SantaLucia model, a nearest-neighbor model that is able to predict experimental duplex melting temperatures very well.SantaLucia and Hicks (2004) The SantaLucia model is designed at the level of base-pair steps, where each base-pair step consists of four bases, with a free-energy difference between its single-stranded state and its duplex state. In total there are 10 unique base-pair steps in the SantaLucia model: AA/TT, AT/AT, TA/TA, GC/GC, CG/CG, GG/CC, GA/TC, AG/CT, TG/CA, GT/AC (for example, AG/CT refers to a base-pair step with complementary bases AG on one strand and CT on the other, both specified in the to direction). Therefore, fitting oxDNA to the SantaLucia model only allows one to find the sum of the strengths of the two stacking interactions between the nucleotides that are within a base-pair step. For example, it is only possible to find the average of the AA and TT stacking strengths from the SantaLucia model.
However, experimental evidence from sequence-dependent measurements of the mechanical properties of single-stranded DNA,Uzawa et al. (2013); Bosco et al. (2013) as well as hairpin stabilities and closing rates, Goddard et al. (2000); Bonnet et al. (1998) have revealed that sequences of A bases are much stiffer than equal length sequences of T bases. It has been argued that this is evidence that consecutive AA bases stack much more strongly than TT bases do. Here we use original experimental data on the stabilities of hairpins with either a poly-A or a poly-T loop, to differentiate between the AA and TT stacking strengths. All hairpins have 6-bp stems and loops of either or bases. Details of the experimental systems and the sequences are given in Appendix C.
To parameterise our model to reproduce the experimental data, our procedure is to vary the AA stacking interaction strength, , while fixing the sum of the AA and TT stacking interaction strengths to the value that is obtained by fitting oxDNA’s duplex melting temperatures to the SantaLucia model, i.e.
where is the TT stacking strength, is the strength for both AA and TT obtained from fitting oxDNA to the SantaLucia model. We vary the AA stacking interaction strength to match the experimental differences of the thermal stabilities of the poly-A-loop and poly-T-loop hairpins, specifically , where is the difference in free energy between the bound (hairpin) state and the unbound (open) state for a hairpin with a poly-X loop.
We obtain the bound and unbound free energies , where , by using thermodynamic integration. We start by calculating , the free energy of a reference state for which the effective charge on the backbone site, , is set to 0, and . , the free energy of the state with salt concentration , stacking strengths and , and effective backbone charge , is then obtained by solving the following integrals
where measures the deviation of the AA and TT strengths from , such that Eq. 8 is satisfied, i.e.
Taking the derivatives of the free energy in Eq. 9 we obtain,
where is the average stacking energy of all AA(TT) nucleotides in a strand in state with salt concentration , defined as in Eq. 10, and backbone charge . The terms appearing in the integrands of Eq. VI can be obtained by running short simulations of bound and unbound states separately and therefore the free energies can be calculated in an efficient manner. The results for the difference in between oxDNA and experiment as a function of for the hairpins with 21- and 31-base loops are shown in Fig. 5; the value for the AA stacking strength that minimises this difference is found to be (corresponding to ), which is not too dissimilar to the preliminary value suggested in Ref. Šulc et al., 2012. We note that this value gives satisfactory predictions for for a wide range of salt concentrations down to 0.05 M.
Vii The oxDNA2 model
In summary, the oxDNA2 potential can be written as
where a indicates that the term is either modified for oxDNA2, or, in the case of , new in oxDNA2. The modified parameters for oxDNA2 are compared with those for oxDNA in Table A1, and a full account of the changes is given in Appendix A. All other parameters remain the same as in the original model.
We emphasise that, after all the relevant changes to the potential were made, the hydrogen bonding and stacking parameters were modified to ensure that the close agreement to experimental duplex melting temperatures achieved for the original model was retained with oxDNA2 (this was done immediately after the changes to , as described at the end of Section V).
Viii Physical properties of oxDNA2
The structural, mechanical and thermodynamic properties of DNA for the new version of oxDNA presented in this paper, which we call oxDNA2, are slightly different from the properties for the original oxDNA model. We briefly highlight the most important of these changes here; in addition these properties are given as a function of salt concentration as this is now possible with the new model. We note that all of the results described in this section were computed using the final version of the oxDNA2 model as summarized in Section VII. The details of the simulations used to compute these results are given briefly in the figure captions and in detail in Appendix B.3 and Table A3.
The structural properties of the new model, specifically the pitch and rise of double-stranded DNA, are presented in Fig. 6. As might be expected, the rise increases with decreasing salt concentration, due to the greater repulsion between backbone sites. The pitch also increases with decreasing salt, consistent with the measured increase in rise and slight decrease in neighbouring backbone-backbone distance measured as salt concentration is decreased (the duplex radius remains approximately constant). Although there is not much experimental evidence to compare this with, there is some indication that the pitch is roughly constant for low salt concentrations (0.162 M and below).Taylor and Hagerman (1990) As mentioned earlier, the pitch is chosen (by modifying the bonded neighbour backbone-backbone interaction) so that the global twist of origami structures agrees with experimental measurements. Specifically, we set the backbone-backbone interaction so that the helix bundle designed to have no global twist has no global twist in the model at [Na]=0.5 M. This results in a pitch of roughly bp/turn at this salt concentration, compared to bp/turn in the original model, and experimental values of around 10.45 suggested by cyclisation experiments,Taylor and Hagerman (1990); Shore and Baldwin (1983) albeit in the presence of some divalent salt.
The thermodynamics of duplex formation are shown in Fig. 7. The transition width for the yield curve of a 10-bp duplex at 0.5 M [Na] in oxDNA2 is largely unchanged from the original oxDNA, and the transition widths depend weakly if at all on salt in oxDNA2. The free-energy profiles for duplex formation in oxDNA2 show that the free-energy cost of forming the first base pair decreases with increasing salt, presumably due to the reduced energetic cost of bringing the two single strands close together as the electrostatic repulsion between backbone sites becomes more short-ranged. The slope of the bound region of the free-energy profile also decreases with increasing salt, indicating a reduced free-energy cost of forming subsequent base pairs.
In Fig. 8, the duplex melting temperatures as a function of salt and for different duplex lengths in oxDNA2 are compared to the results from the SantaLucia model, to which oxDNA2 was parameterised. As expected, we find good agreement.
The mechanical properties of the model double-stranded DNA as a function of salt, specifically the persistence length and torsional stiffness, are shown in Fig. 9. The persistence length was calculated by computing the correlation of the helix axis as described in Ref. Ouldridge et al., 2011. The effective torsional stiffness, , was computed from the linear regime of the torque response curve of a 60-bp duplex under tension (see Appendix B.2 for details). The simulations were carried out under a linear force of 30pN, a high force regime where we expect that approximates the true torsional stiffness .Moroz and Nelson (1997); Matek et al. (2015) The persistence length at [Na]=0.5M, about 123bp, is very slightly lower than for the original model, 125bp, and the persistence length at different salt concentrations is consistent with the rather broad range of values suggested by experiment.Savelyev (2012); Herrero-GalÃ¡n et al. (2013) The decrease in persistence length with increasing salt is expected due to the decrease in repulsion between the duplex’s backbone sites, which makes the duplex less stiff, in agreement with experiment.Savelyev (2012) The slight decrease in torsional stiffness with increasing salt can be rationalised in the same way. Regardless of the salt concentration, the torsional stiffness measured for oxDNA2 ( fJfm) is lower than for the original modelMatek et al. (2015) ( fJfm). Single-molecule twisting experiments give a value of fJfm for a pulling force of 3.5 pNJanssen et al. (2012) and a pulling force of 45 pN,Bryant et al. (2003) at around [Na]=0.1 M. There is limited experimental data on the salt dependence of torsional stiffness in DNA; however the torsional stiffness has been reported to be roughly constant in the range [Na]=0-0.162 M.Taylor and Hagerman (1990)
The effects of various motifs on duplex melting temperatures in oxDNA2 are compared to results from the original oxDNA and the SantaLucia model in Table 2. These effects are either barely changed in oxDNA2, or are now closer to the SantaLucia values than the original oxDNA values were.
The melting temperatures of hairpins at 0.5 M [Na] are shown in Fig. 10. The hairpin melting temperatures are lower than for the original oxDNA by about 2 K, which in turn had hairpin melting temperatures about 3 K lower than those predicted by SantaLucia. The gap between the oxDNA2 hairpin melting temperatures and those given by the SantaLucia model widens somewhat as the salt concentration is lowered (Fig. 11), indicating a stronger salt dependence in oxDNA2, with a typical underestimate of around 7K at a salt concentration of 100mM for a relatively short loop length of 6 bases. For longer loops, the difference in predicted melting temperatures between oxDNA2 and SantaLucia widens further. This difference is unsurprising, as in oxDNA2 (and physical DNA) ssDNA becomes stiffer at lower salt concentrations, making the formation of a hairpin less favorable,Mosayebi et al. (2014) whereas in the SantaLucia model the loops’ contribution to hairpin stability is salt-independent.SantaLucia and Hicks (2004)
It seems plausible that oxDNA2’s performance is better than implied by the salt-independent loop contribution in the SantaLucia model. In Fig. 12 we consider hairpins with a short 5-bp stem and a long 31-base loop, making the hairpin thermodynamics particularly sensitive to the change in loop stiffness with salt. The hairpin stability (i.e. the free energy difference between the bound and unbound state) as predicted by oxDNA2 and the SantaLucia model is compared to experimental data determined using FRET (experimental details in Appendix C). As expected we find the hairpin stabilities for oxDNA2 have a much stronger salt dependence (steeper gradient) than those of SantaLucia (shallow gradient). At higher salt the oxDNA2 results show a similar slope to the experimental curve, but at lower salt (0.2 M and below) oxDNA2 shows a stronger destabilization of the hairpins with decreasing salt (a steeper gradient) than experiment. Thus, although oxDNA2 does a reasonable job of capturing these physical effects, the comparison suggests that the single strands may be experiencing too much repulsion at low salt. We also note that the SantaLucia model is unable to predict the large difference in stability of the two hairpins due to their different loop sequences, because the SantaLucia model is insensitive to the loop sequence except for the two bases adjacent to the stem.
|Bulge (1-base bulge in 8-bp duplex)|
|Bulge (2-base bulge in 8-bp duplex)|
|Terminal mismatch (1-bp mismatch in 5-bp duplex)|
|Internal mismatch (2-bp mismatch in 8-bp duplex)|
We examine the flexibility of single strands as a function of salt in oxDNA2 by measuring the radius of gyration, , of a single strand of DNA at different salt concentrations, sequences and strand lengths, and comparing to experimental results (Fig. 13). The oxDNA2 model reproduces the overall trend of increasing with salt, and is in agreement with the somewhat noisy experimental data. We note that oxDNA2 is able to capture two important effects: the greater stiffness (and hence larger ) of the more strongly-stacking poly-A strands; and the greater salt dependence of the poly-T strands’ compared to that of the poly-A strands, caused by the weaker stacking of the poly-T strands, which means that electrostatic repulsion makes a greater relative contribution to its stiffness.
OxDNA has been applied to a wider variety of systems than any other coarse-grained model of DNA. The modifications and extensions to oxDNA presented here open up a variety of new potential applications for the model. With the introduction of an explicit salt-dependent term in the potential, the model can be used to simulate systems under physiological conditions, and to investigate the salt-dependent behaviour of DNA. Also, the parametrization presented here allows a quantitative comparison to experiments run in a wide range of salt concentrations rather than just in the high-salt limit. The introduction of major-minor grooving adds detail to the model, and, combined with other small modifications, allows the use of oxDNA2 to accurately characterise the structural properties of large DNA nanostructures, such as DNA origami. In particular the helical pitch and the twist angles at nicks and junctions were fine-tuned to obtain a correct global twist for a test-case 3D origami structure. In the absence of definitive experimental values for these individual parameters we chose a combination that we deemed physically reasonable and that produces equilibrium structures that compare well with experimental ones. Finally, the sequence dependence in the model has been extended by introducing different interaction strengths for the AA and TT stacking. This change will be particularly useful for studying the effects of stacking in single strands or single-stranded sections (e.g. hairpin loops) as poly-A and poly-T sections are often used as paradigmatic examples of strongly and weakly stacking sequences, respectively, and for modelling DNA nanostructures where poly-T loops are often used as flexible linkers.
There are some limitations associated with the new features of the model. The Debye-Hückel treatment of the electrostatics is perhaps an oversimplification, and we should not expect it to capture all of the complexities associated with the electrostatics for DNA, but it is not straightforward to think of a different approach that would be consistent with our level of coarse graining. In particular, the fact that hairpin stability is reduced faster than in experiment suggests that single strands may experience too much repulsion in oxDNA2 at low salt, and care should be taken in making predictions based on our models at low salt concentration if the system under investigation depends crucially on the thermodynamics of long single-stranded sections.
We are currently exploiting the improvements in the model’s structural prediction to carry out investigations on DNA origami structures as well as on structures composed of multi-arm tilesZhang et al. (2009). At the same time, we are studying the salt dependent thermodynamics of a diverse set of systems, which would not have been possible without the introduction of salt dependence into the potential. The introduction of different strengths for the AA and TT stacking in the model allows us to better capture these effects, and in addition gives us more accurate sequence-dependent hairpin thermodynamics and kinetics, which are currently being exploited to further study hairpins. We note that a simulation code implementing both the original oxDNA model and the new oxDNA2 model, including an implementation using GPUs,Rovigatti et al. (2015) is available for download from dna.physics.ox.ac.uk.
The authors are grateful to the Engineering and Physical Sciences Research Council for financial support. M. M. was supported by the Swiss National Science Foundation (Grant No.PBEZP2-145981). The authors thank Oxford’s Advanced Research Computing, E-infrastructure South and the PolyHub virtual organization for computer time.
- Dietz et al. (2009) H. Dietz, S. M. Douglas, and W. M. Shih, Science, 325, 725 (2009).
- Han et al. (2013) D. Han, S. Pal, Y. Yang, S. Jiang, J. Nangreave, Y. Liu, and H. Yan, Science, 339, 1412 (2013).
- Ke et al. (2012) Y. Ke, L. L. Ong, W. M. Shih, and P. Yin, Science, 338, 1177 (2012).
- Walsh et al. (2011) A. S. Walsh, H. Yin, C. M. Erben, M. J. A. Wood, and A. J. Turberfield, ACS Nano, 5, 5427 (2011).
- Douglas et al. (2012) S. M. Douglas, I. Bachelet, and G. M. Church, Science, 335, 831 (2012).
- Kuzyk et al. (2012) A. Kuzyk, R. Schreiber, Z. Fan, G. Pardatscher, E.-M. Roller, A. Hogele, F. C. Simmel, A. O. Govorov, and T. Liedl, Nature, 483, 311 (2012).
- Zadegan et al. (2012) R. M. Zadegan, M. D. E. Jepsen, K. E. Thomsen, A. H. Okholm, D. H. Schaffert, E. S. Andersen, V. Birkedal, and J. Kjems, ACS Nano, 6, 10050 (2012).
- Pinheiro et al. (2011) A. V. Pinheiro, D. Han, W. M. Shih, and H. Yan, Nat. Nano, 6, 763 (2011).
- Svozil et al. (2010) D. Svozil, P. Hobza, and J. Šponer, J. Phys. Chem. B, 114, 1191 (2010).
- Mladek et al. (2013) A. Mladek, M. Krepl, D. Svozil, P. Cech, M. Otyepka, P. Banas, M. Zgarbova, P. Jurecka, and J. Šponer, Phys. Chem. Chem. Phys., 15, 7295 (2013).
- Šponer et al. (2013) J. Šponer, J. E. Šponer, A. MlÃ¡dek, P. BanÃ¡Å¡, P. JureÄka, and M. Otyepka, Methods, 64, 3 (2013).
- Laughton and Harris (2011) C. A. Laughton and S. A. Harris, Wiley Interdisc. Rev. Comput. Mol. Sc., 1, 590 (2011).
- Pérez et al. (2012) A. Pérez, F. J. Luque, and M. Orozco, Acc. Chem. Res., 45, 196 (2012).
- Maffeo et al. (2012) C. Maffeo, B. Luan, and A. Aksimentiev, Nucleic Acids Res., 40, 3812 (2012).
- Yoo and Aksimentiev (2013) J. Yoo and A. Aksimentiev, Proc. Natl. Acad. Sci. USA, 110, 20099 (2013).
- Maffeo et al. (2014) C. Maffeo, J. Yoo, J. Comer, D. B. Wells, B. Luan, and A. Aksimentiev, Journal of Physics: Condensed Matter, 26, 413101 (2014a).
- Bustamante et al. (1994) C. Bustamante, J. F. Marko, E. D. Siggia, and S. Smith, Science, 265, 1599 (1994).
- Hinckley et al. (2013) D. M. Hinckley, G. S. Freeman, J. K. Whitmer, and J. J. de Pablo, J. Chem. Phys., 139, 144903 (2013).
- Savelyev and Papoian (2010) A. Savelyev and G. A. Papoian, Proc. Natl. Acad. Sci. USA, 107, 20340 (2010).
- Cragnolini et al. (2013) T. Cragnolini, P. Derreumaux, and S. Pasquali, J. Phys. Chem. B, 117, 8047 (2013).
- He et al. (2013) Y. He, M. Maciejczyk, S. Ołdziej, H. A. Scheraga, and A. Liwo, Phys. Rev. Lett., 110, 098101 (2013).
- Linak et al. (2011) M. C. Linak, R. Tourdot, and K. D. Dorfman, J. Chem. Phys., 135, 205102 (2011).
- Morriss-Andrews et al. (2010) A. Morriss-Andrews, J. Rottler, and S. S. Plotkin, J. Chem. Phys., 132, 035105 (2010).
- Araque et al. (2011) J. C. Araque, A. Z. Panagiotopoulos, and M. A. Robert, J. Chem. Phys., 134, 165103 (2011).
- Korolev et al. (2014) N. Korolev, D. Luo, A. P. Lyubartsev, and L. Nordenskiöld, Polymers, 6, 1655 (2014).
- Maffeo et al. (2014) C. Maffeo, T. T. M. Ngo, T. Ha, and A. Aksimentiev, J. Chem. Theory Comput., 10, 2891 (2014b).
- Ouldridge et al. (2011) T. E. Ouldridge, A. A. Louis, and J. Doye, J. Chem. Phys., 134, 085101 (2011).
- Ouldridge (2011) T. E. Ouldridge, PhD Thesis, University of Oxford (2011).
- Ouldridge et al. (2010) T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, Phys. Rev. Lett., 104, 178101 (2010a).
- Šulc et al. (2013) P. Šulc, T. E. Ouldridge, F. Romano, J. P. K. Doye, and A. A. Louis, Natural Computing (2013).
- Ouldridge et al. (2013) T. E. Ouldridge, P. Šulc, F. Romano, J. P. K. Doye, and A. A. Louis, Nucleic Acids Res., 41, 8886 (2013a).
- Ouldridge et al. (2013) T. E. Ouldridge, R. L. Hoare, A. A. Louis, J. P. K. Doye, J. Bath, and A. J. Turberfield, ACS Nano, 7, 2479 (2013b).
- Rovigatti et al. (2014) L. Rovigatti, F. Smallenburg, F. Romano, and F. Sciortino, ACS Nano, 8, 3567 (2014).
- Schreck et al. (2014) J. S. Schreck, T. E. Ouldridge, F. Romano, A. A. Louis, and J. P. K. Doye, ArXiv e-prints (2014a), arXiv:1412.6309 .
- Schreck et al. (2014) J. S. Schreck, T. E. Ouldridge, F. Romano, P. Šulc, L. Shaw, A. A. Louis, and J. P. K. Doye, ArXiv e-prints (2014b), arXiv:1408.4401 .
- Srinivas et al. (2013) N. Srinivas, T. E. Ouldridge, P. Šulc, J. M. Schaeffer, B. Yurke, A. A. Louis, J. P. K. Doye, and E. Winfree, Nucleic Acids Res., 41, 10641 (2013).
- Mosayebi et al. (2015) M. Mosayebi, A. A. Louis, J. P. K. Doye, and T. E. Ouldridge, arXiv preprint arXiv:1502.03623 (2015).
- Matek et al. (2012) C. Matek, T. E. Ouldridge, A. Levy, J. P. K. Doye, and A. A. Louis, J. Phys. Chem. B, 116, 11616 (2012).
- Wang and Pettitt (2014) Q. Wang and B. M. Pettitt, Biophys. J., 106, 1182 (2014).
- Matek et al. (2015) C. Matek, T. E. Ouldridge, J. P. K. Doye, and A. A. Louis, Sci. Rep., 5, 7655 (2015).
- Mosayebi et al. (2014) M. Mosayebi, F. Romano, T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, J. Phys. Chem. B, 118, 14326 (2014).
- Romano et al. (2013) F. Romano, D. Chakraborty, J. P. K. Doye, T. E. Ouldridge, and A. A. Louis, J. Chem. Phys., 138, 085101 (2013).
- Machinek et al. (2014) R. R. F. Machinek, T. E. Ouldridge, N. E. C. Haley, J. Bath, and A. J. Turberfield, Nat. Comm., 5, 5324 (2014).
- Bosco et al. (2013) A. Bosco, J. Camunas-Soler, and F. Ritort, Nucleic Acids Res., 42, 2064 (2013).
- Polinkovsky et al. (2014) M. E. Polinkovsky, Y. Gambin, P. R. Banerjee, M. J. Erickstad, A. Groisman, and A. A. Deniz, Nat. Comm., 5, 5737 (2014).
- SantaLucia and Hicks (2004) J. SantaLucia, Jr. and D. Hicks, Annu. Rev. Biophys. Biomol. Struct., 33, 415 (2004).
- Taylor and Hagerman (1990) W. H. Taylor and P. J. Hagerman, J. Mol. Biol., 212, 363 (1990).
- Mocci and Laaksonen (2012) F. Mocci and A. Laaksonen, Soft Matter, 8, 9268 (2012).
- Potoyan et al. (2013) D. A. Potoyan, A. Savelyev, and G. A. Papoian, WIREs Comput. Mol. Sci, 3, 69 (2013).
- Shore and Baldwin (1983) D. Shore and R. L. Baldwin, J. Mol. Biol., 170, 983 (1983).
- Dawid et al. (2006) A. Dawid, F. Guillemot, C. Brème, V. Croquette, and F. Heslot, Phys. Rev. Lett., 96, 188102 (2006).
- Ido et al. (2013) S. Ido, K. Kimura, N. Oyabu, K. Kobayashi, M. Tsukada, K. Matsushige, and H. Yamada, ACS Nano, 7, 1817 (2013).
- Du et al. (2005) Q. Du, C. Smith, N. Shiffeldrim, M. Vologodskaia, and A. Vologodskii, Proc. Nat. Acad. Sci. USA, 102, 5397 (2005).
- Wang (1979) J. C. Wang, Proc. Nat. Acad. Sci. USA, 76, 200 (1979).
- Kato et al. (2009) T. Kato, R. P. Goodman, C. M. Erben, A. J. Turberfield, and K. Namba, Nano Lett., 9, 2747 (2009).
- Bai et al. (2012) X.-C. Bai, T. G. Martin, S. H. W. Scheres, and H. Dietz, Proc. Nat. Acad. Sci. USA, 109, 20012 (2012).
- Šulc et al. (2012) P. Šulc, F. Romano, T. E. Ouldridge, L. Rovigatti, J. P. K. Doye, and A. A. Louis, J. Chem. Phys., 137, 135101 (2012).
- Andersen et al. (2008) F. F. Andersen, B. Knudsen, C. L. P. Oliveira, R. F. Frohlich, D. Kruger, J. Bungert, M. Agbandje-McKenna, R. McKenna, S. Juul, C. Veigaard, J. Koch, J. L. Rubinstein, B. Guldbrandtsen, M. S. Hede, G. Karlsson, A. H. Andersen, J. S. Pedersen, and B. R. Knudsen, Nucleic Acids Res., 36, 1113 (2008).
- Tian et al. (2013) C. Tian, C. Zhang, X. Li, C. Hao, S. Ye, and C. Mao, Langmuir, 30, 5859 (2013).
- Douglas et al. (2009) S. M. Douglas, H. Dietz, T. Liedl, B. Högberg, F. Graf, and W. M. Shih, Nature, 459, 414 (2009).
- Šulc et al. (2014) P. Šulc, F. Romano, T. E. Ouldridge, J. P. K. Doye, and A. A. Louis, J. Chem. Phys., 140, 235102 (2014).
- Zhang et al. (2012) C. Zhang, W. Wu, X. Li, C. Tian, H. Qian, G. Wang, W. Jiang, and C. Mao, Angew. Chem. Int. Edit., 51, 7999 (2012).
- Hinckley et al. (2014) D. M. Hinckley, J. P. Lequieu, and J. J. de Pablo, J. Chem. Phys., 141, 035102 (2014).
- Maffeo et al. (2010) C. Maffeo, R. Schöpflin, H. Brutzer, R. Stehr, A. Aksimentiev, G. Wedemann, and R. Seidel, Phys. Rev. Lett., 105, 158101 (2010).
- Mazur and Jernigan (1991) J. Mazur and R. L. Jernigan, Biopolymers, 31, 1615 (1991).
- Frenkel and Smit (1996) D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications (Academic Press, 1996).
- Vega et al. (2008) C. Vega, E. Sanz, J. Abascal, and E. Noya, J. Phys.: Condens. Matter, 20, 153101 (2008).
- Ouldridge et al. (2010) T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, J. Phys.: Condens. Matter, 22, 104102 (2010b).
- Uzawa et al. (2013) T. Uzawa, T. Isoshima, Y. Ito, K. Ishimori, D. E. Makarov, and K. W. Plaxco, Biophys. J., 104, 2485 (2013).
- Goddard et al. (2000) N. L. Goddard, G. Bonnet, O. Krichevsky, and A. Libchaber, Phys. Rev. Lett., 85, 2400 (2000).
- Bonnet et al. (1998) G. Bonnet, O. Krichevsky, and A. Libchaber, Proc. Nat. Acad. Sci. USA, 95, 8602 (1998).
- Whitelam et al. (2009) S. Whitelam, E. H. Feng, M. F. Hagan, and P. L. Geissler, Soft Matter, 5, 1251 (2009).
- Kumar et al. (1992) S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, J. Comput. Chem., 13, 1011 (1992).
- Moroz and Nelson (1997) J. D. Moroz and P. Nelson, 94, 14418 (1997).
- Savelyev (2012) A. Savelyev, Phys. Chem. Chem. Phys., 14, 2250 (2012).
- Herrero-GalÃ¡n et al. (2013) E. Herrero-GalÃ¡n, M. E. Fuentes-Perez, C. Carrasco, J. M. Valpuesta, J. L. Carrascosa, F. Moreno-Herrero, and J. R. Arias-Gonzalez, J. Am. Chem. Soc., 135, 122 (2013).
- Janssen et al. (2012) X. J. A. Janssen, J. Lipfert, T. Jager, R. Daudey, J. Beekman, and N. H. Dekker, Nano Lett., 12, 3634 (2012).
- Bryant et al. (2003) Z. Bryant, M. D. Stone, J. Gore, S. B. Smith, N. R. Cozzarelli, and C. Bustamante, Nature, 424, 338 (2003).
- Sim et al. (2012) A. Y. L. Sim, J. Lipfert, D. Herschlag, and S. Doniach, Phys. Rev. E, 86, 021901 (2012).
- Zhang et al. (2009) C. Zhang, Y. He, M. Su, S. H. Ko, T. Ye, Y. Leng, X. Sun, A. E. Ribbe, W. Jiang, and C. Mao, Faraday Discuss., 143, 221 (2009).
- Rovigatti et al. (2015) L. Rovigatti, P. Šulc, I. Z. Reguly, and F. Romano, J. Comput. Chem., 36, 1 (2015).
- De Michele et al. (2012) C. De Michele, L. Rovigatti, T. Bellini, and F. Sciortino, Soft Matter, 8, 8388 (2012).
- Yakovchuk et al. (2006) P. Yakovchuk, E. Protozanova, and M. D. Frank-Kamenetskii, Nucleic Acids Res., 34, 564 (2006).
- Tsukanov et al. (2013) R. Tsukanov, T. E. Tomov, Y. Berger, M. Liber, and E. Nir, J. Phys. Chem. B, 117, 16105 (2013a).
- Tsukanov et al. (2013) R. Tsukanov, T. E. Tomov, R. Masoud, H. Drory, N. Plavner, M. Liber, and E. Nir, J. Phys. Chem. B, 117, 11932 (2013b).
- Nir et al. (2006) E. Nir, X. Michalet, K. M. Hamadani, T. A. Laurence, D. Neuhauser, Y. Kovchegov, and S. Weiss, J. Phys. Chem. B, 110, 22103 (2006).
Appendix A Modifications to the oxDNA potential in oxDNA2
a.1 Introducing an electrostatic term into the model potential
We introduce an explicit electrostatic term into the model using Debye-Hückel theory. To aid computational efficiency when simulating the model, we ensure that the interaction goes to zero within a finite distance , and we introduce a quadratic smoothing function so that the interaction goes to zero smoothly. The form of the electrostatic term is then
where is the temperature, is the salt concentration, and , the distance above which smoothing is introduced, is chosen to be equal to ( is defined below). The Debye-Hückel-like term, , is given by
where is the effective charge (which is equal to 1 in Debye-Hückel theory but which we choose to be 0.815 after fitting to experimental data), is the charge of an electron, is the permittivity of a vacuum, is the relative permittivity of water (which we choose to be 80), and is the distance between the pair of interacting backbone sites. is the Debye screening length, and is given by
where is Boltzmann’s constant and is Avogadro’s number. The quadratic smoothing term, , is given by
and and are constants chosen such that is smooth and differentiable. Note that the above definition implies that the cutoff radius for the Debye-Hückel term, , depends on the salt concentration and the temperature .
a.2 Modifying the average pitch of the model duplex
The model duplex pitch can be controlled by modifying the position of the minimum in the FENE bonded backbone potential, , which has a small effect on the thermodynamics (which was refitted using the histogram reweighting method described in Section II C 2 of Ref. Šulc et al., 2014) and little effect on the other properties of the model DNA. For oxDNA2, was set to give a global twist of zero for the N-type helix bundle due to Dietz et al.;Dietz et al. (2009) the new and old values for are given in Table A1.
a.3 Modifying the average coaxial stacking term in the oxDNA potential
We modified the coaxial stacking term in the oxDNA potential in order to change the twist angle across a nick and across an origami junction, while changing the thermodynamic properties of a nicked duplex and other motifs as little as possible. This was achieved by moving the position of the minimum with respect to the angle (shown in Fig. A1) of the coaxial stacking term, , and by increasing the well depth, , to make the minimum more energetically favourable. Increasing the well depth was necessary because it is less favourable for a pair of coaxially stacking nucleotides in a nicked duplex to stack with this tighter angle; a stronger overall coaxial stacking interaction compensates for this effect. The changes to the model are shown in Table A1.
The altered optimal angle requires an alteration of the modulation function involving for coaxial stacking. In the original oxDNA, the modulation is given by , where is defined in Ref. Ouldridge, 2011. The second term was introduced to ensure that the potential is differentiable at , where the gradients of the two terms cancel. In oxDNA, this additional term had a very small effect on the numerical value of the potential for . Using the same trick in oxDNA2, however, is problematic, because is much closer to . The second term, , would then be substantial for and would distort the modulating function in an undesired fashion. To avoid this, we replace the second term, , with a new function :
The factor in the coaxial stacking term is then replaced by . The values for and are chosen to ensure that the coaxial stacking potential is smooth and differentiable; they are specified in Table A1.
For oxDNA2, we made an additional change to the coaxial stacking potential: we removed the and terms. These terms had allowed only right-handed blunt-ended stacking, disallowing left-handed blunt-ended stacking; however, since the development of the original oxDNA model, this constraint has been deemed unnecessary, as the blunt-ended coaxial stacking interaction was found to be likely to be too weak compared to experiment,De Michele et al. (2012) and there is no experimental evidence indicating that left-handed blunt-ended coaxial stacking is not possible.
Incorporating the changes described above, the final form for the oxDNA2 coaxial stacking term is
Table A2 gives a comparison of the thermodynamic properties of the motifs relevant to the coaxial stacking term in oxDNA and oxDNA2. The free-energy change upon stacking across a nick, , is slightly more negative in oxDNA2, indicating that a nicked duplex stacks slightly more strongly across the nick in oxDNA2 than in the original oxDNA. In both cases the stacking is stronger than is seen in experiment ( in Ref. Yakovchuk et al., 2006 – a discussion comparing the value of measured for the original oxDNA with experimental results is given in Ref. Ouldridge, 2011). The stabilising effect of a hairpin stem on duplex hybridisation (see Fig. A2, , is found to be slightly stronger than for the original oxDNA, and is on average in better agreement with the stability predicted by the SantaLucia model. Finally, dimerisation of two duplexes through blunt-ended coaxial stacking is found to be less favourable in oxDNA2 than in oxDNA. This is despite the mitigating effect of removing the terms in the coaxial stacking interaction, which strengthens blunt-ended stacking. A previous studyDe Michele et al. (2012) has shown that oxDNA already underestimated the stability of blunt-ended coaxial stacking.
a.4 Modifying the hydrogen bonding and stacking interaction strengths
Model duplexes were found to fray slightly less often when different widths for the major and minor duplex grooves were introduced compared to the equal-groove-width case, which led to a slight raising of the duplex melting temperatures. In addition, the change to the position of the minimum in the FENE bonded backbone potential, , was found to destabilise the duplex state, slightly lowering the duplex melting temperatures. In the process of the reparameterisation of oxDNA to incorporate these two changes, the duplex melting temperatures were reset on the introduction of each change by making small modifications to the stacking and hydrogen bonding strengths using the histogram reweighting method described in Section II C 2 of Ref. Šulc et al., 2014, giving an agreement with experimental melting temperatures that was just as good as before. The new parameters are specified, and compared to the old parameters, in Table A1, while Fig. 7 gives a comparison of the duplex thermodynamics with the old and new versions of the model.
Appendix B Simulation methods
b.1 Measuring the global twist of an origami
We measure the global twist of the origami due to Dietz et al.Dietz et al. (2009) in the following way: The origami (for example, the L-type, although this applies equally well to all three types) can be thought of as a set of antiparallel double helices on a hexagonal lattice joined together by crossovers. The overall shape of this structure is roughly that of a (possibly twisted) rectangular prism. To measure the twist, we consider the two faces of the rectangular prism which are opposite each other and approximately normal to the axes of the double helices. For each face, we define a vector along its long edge, and we superimpose these two vectors into a plane perpendicular to the average of the helix axes. We then define the global twist as the angle between these two vectors.
To avoid end effects, the top and bottom rows of helices (3 helices in the top row and 3 helices in the bottom row) were excluded for this analysis, as were approximately 22 base pairs from each end of every double helix.
b.2 Measuring the torsional stiffness of a duplex in the oxDNA model
The torsional stiffness of a duplex in oxDNA2 was measured according to the scheme outlined by Matek et al.Matek et al. (2015) In summary, this method involves pulling a duplex while also twisting it, and measuring the torque required to impose a given twist. In this case, a 60-bp duplex was pulled at a constant force of 30 pN and held with virtual traps which imposed a helical pitch different from the equilibrium pitch. For each salt concentration, simulations were run with imposed pitches of 9.8, 10.1, 10.4, 10.7, 11.0, and 11.3 bp/turn, with each simulation run for at least MD steps. The torque observed was plotted as a function of imposed twist (this is called a torque response curve), and the gradient of the linear region used to approximate the effective torsional stiffness.
The torsional stiffness is computed as
where is the change in torque, and the change in superhelical twist density, measured in the linear regime of the torque response curve, is the double helical rise for a relaxed duplex and is the twist angle across a base-pair step for a relaxed duplex, measured in radians. According to Moroz and Nelson,Moroz and Nelson (1997) the effective torsional stiffness measured by twisting a duplex under tension can be related to the true torsional stiffness by
where is the linear force and is the bending stiffness. Our simulations are in the high-force regime, where should be a good approximation for .
b.3 Details of the simulation protocols
The simulation results presented in this paper were obtained using either a virtual move Monte Carlo (VMMC) algorithm or a molecular dynamics (MD) algorithm with a Brownian thermostat. Note that, where VMMC steps are reported in the main text, they refer to the number of accepted moves during the simulation. Details of the simulation parameters used are given in Table A3.
For the duplex free-energy profiles, yields and melting temperatures (Fig. 7 and Fig. 8) we used umbrella sampling with an order parameter defined in terms of the number of native base pairs (0 to where is the number of base pairs in a fully formed duplex). A base pair was defined as being formed if the potential’s hydrogen-bonding energy term for that pair was lower than . For the hairpin melting temperature calculations (Fig. 10, Fig. 11 and Fig. 12) we used an analogous order parameter for the native base pairs in the hairpin stem.
The oxDNA2 results reported in Fig. 12 were computed using thermodynamic integration as described in Section VI. For each oxDNA2 data point shown in the figure, first an umbrella sampling simulation of each hairpin (sequences shown in Table A4) was carried out for a version of the model with , where is given by Eq. 10, and in order to find the reference free energy ( in Eq. VI). The simulations were run for VMMC steps and used the parameters specified in Table A3. The integrals in Eq. VI were discretised into 15 gridpoints for the integral in and 10 gridpoints for the integral in . Then simulations to compute these integrals were run for roughly VMMC steps for each grid point, so that each data point in the figure required 25 grid point simulations, with one umbrella sampling simulation per hairpin.
|Parameter||Standard||Bulge||Nick stack.||Hairpin stab.||Blunt coax||Hairpin||Helix bundle twist|
|max. cluster size||-||-||-||-||-||-||-||-|
Appendix C Experimental Methods
c.1 DNA Hairpin Design
The hairpin constructs consist of a loop with different lengths (21 or 31 nucleotides) and sequences (adenosine or thymine), and six base pairs in the stem, as shown in Table A4.
A 35-bp long, double-stranded DNA was used to keep the hairpin section away from the coverslip and the origami surfaces. The hairpins were labeled with donor and acceptor fluorophores (ATTO-550 and ATTO-647N, respectively), positioned such that the hairpin open state yielded low FRET values () and the closed state yielded high values (Figure A3).
c.2 Annealing Procedures and Experimental Setup
The top and bottom strands were annealed at C (1.5 M in 10 L TE-NaCl buffer [10 mM Tris, pH 8.0, 1 mM EDTA], and 100 mM NaCl) and then gradually cooled (30 min) to room temperature (using PCR).
The measurements were performed on a diffusion-based single-molecule FRET-ALEX setup, as described elsewhere.Tsukanov et al. (2013) In brief, a green CW laser beam (532 nm, CL532-025-L, Crystal Laser, Reno, NV, USA) was aligned/misaligned into a single-mode fiber using an acousto-optic modulator (AOM; R23080-2-LTD, Neos Technologies, Melbourne, FL, USA), alternating with a red diode laser (640 nm, 1069417, Coherent, Auburn, CA, USA) that was electronically switched on/off. The AOM and the red laser were computer controlled with a 12.4-s on-time, a 12.6-s off-time, and a phase shift of 12.5 s. The intensity rise and fall times were less than 50 ns, and there was no time overlap between the lasers. The laser beams were combined by a dichroic mirror (Z532RDC, Chroma, Bellows Falls, VT, USA) and coupled into a single-mode fiber (P1-460A-FC-2, Thorlabs, Newton, NJ, USA). The laser intensities were tuned such that the doubly labeled species would yield a stoichiometry, , of roughly 0.5 (90 W for the green laser and 35 W for the red laser, measured after the fiber while alternating). After collimation (objective PLCN10Ã/0.25, Olympus America, Melville, NY, USA), the combined green and red beams were introduced into a commercial inverted microscope (IX71, Olympus America) and focused about 70 m inside the sample solution with a water-immersion objective (NA 1.2, 60, Olympus America). The emitted fluorescence was separated from the excitation light by a dichroic mirror (ZT532/638RPC, Chroma), focused into a 100-m pinhole (P100S, Thorlabs), re-collimated, split by a second dichroic mirror (FF650-Di01, Semrock, Lake Forest, IL, USA), passed through filters (band-pass filter, Semrock FF01-580/60, for the donor channel and a long-pass filter, Semrock BLP01-635R, for the acceptor channel), and focused into two single-photon avalanche photodiodes (SPAD; SPCM-AQRH-13, Perkin-Elmer Optoelectronics, Fremont, CA, USA). The TTL signals of the two SPADs were recorded as a function of time by a 12.5-ns resolution counting board (PCI-6602, National Instruments, Austin, TX, USA) and in-house prepared LabView acquisition software.
c.3 Measurement Solution and Conditions
Measurements were carried out on diluted samples (40 L of 1 pM hairpins) that were placed on a coverslip and sealed with silicone and an upper coverslip. The bottom coverslip was KOH-treated to prevent sticking to the surface, sonicated for 15 min in 1 M KOH solution, thoroughly washed with distilled water, and dried in air. The measurement buffer comprised 10 mM Tris pH 8.0, 1 mM EDTA, 10 g/mL bovine serum albumin (BSA, Sigma-Aldrich) to reduce sample sticking, 2 mM Trolox (Sigma-Aldrich) to reduce photo-bleaching and photo-blinking of the dyes, and indicated NaCl concentrations. Data were collected for over 30 minutes, and temperature was maintained at C.
c.4 Data analysis
Data analysis was performed with the in-house written LabView software as described elsewhere.Tsukanov et al. (2013) The beginnings and endings of bursts were determined by the all-photons-burst-search (APBS, parameters: = 200, = 10, and = 500 s. For each burst, and were calculated according, respectively, binned (0.01 bin size), and plotted on one-dimensional and histograms and on a two-dimensional histogram.Tsukanov et al. (2013) Donor-only and acceptor-only populations were rejected by ALEX. The Obtained -histogram showed two distinct peaks and was assumed to have quasi two-state dynamicsTsukanov et al. (2013) (Figure A4).
c.5 Experimental Data
The results of the hairpin kinetics experiments are shown in Fig. A5.