A Single Nucleotide Resolution Model for Large-Scale Simulations of Double Stranded DNA

The computational modelling of DNA is becoming crucial in light of new advances in DNA nanotechnology, single-molecule experiments and in vivo DNA tampering. Here we present a mesoscopic model for double stranded DNA (dsDNA) at the single nucleotide level which retains the characteristic helical structure, while being able to simulate large molecules -- up to a million base pairs -- for time-scales which are relevant to physiological processes. This is made possible by an efficient and highly-parallelised implementation of the model which we discuss here. We compare the behaviour of our model with single molecule experiments where dsDNA is manipulated by external forces or torques. We also present some results on the kinetics of denaturation of linear DNA.


Introduction
Since the discovery of the structure of the deoxyribonucleic acid (DNA) [1][2][3] , the geometry of the double-helix and its topological implications have engaged and fascinated the scientific community 4,5 .It is becoming more and more evident that not only is the genetic information encoded in the DNA sequence of primary importance, but also that changes in its three-dimensional structure can alter crucial biological functions, such as gene expression and replication [6][7][8][9][10] .At the same time, the rapid improvement of techniques using DNA functionalised colloids 11,12 , DNA-origami 13 and, more generally, supra-molecular DNA assembly 14 is setting new standards for DNA-based nano-technology.This has farreaching applications, ranging from materials science (to create new DNA-based and possibly biomimetic materials), to medicine (to be used in, e.g., gene-therapy and drug delivery).
In light of this, the formulation of accurate theoretical and computational models that can efficiently capture the behaviour of DNA, either in vivo or in vitro, is of great importance in order to understand a number of outstanding biological problems, and also to assist the advance of DNA-based nanotechnology.
Several fully atomistic models for double-stranded (ds) DNA are available in the literature [15][16][17] .While these give an accurate description of the dynamics of DNA molecules and their interaction with single proteins, the complexity of the all-atom approach places severe limits on the size (up to about a hundred base-pairs) and time scales (of the order of µs) which can be probed 18 .Coarsegraining, where large collections of atoms or molecules are represented by single units, allows larger systems to be simulated for longer at the expense of molecular detail.One of the most challenging aspects in designing a computational model is to retain the key microscopic details necessary to answer a given question while "trimming" the rest.At the large scale limit, entire eukaryotic chromosomes can be modelled using simple beadand-spring polymer models 7,19 , where each monomer can represent up to 3000 base-pairs (bp) and the simulated time can reach time-scales spanning minutes 19 or even days 20 ; similar chains of beads can also be used to model naked DNA, though clearly such an approach neglects microscopic details such as the base-pair specificity or the double-stranded structure.While in some cases these models can still capture the essential physics 21 , in others they are only a crude approximation of the real systems.Several successful mesoscopic models have recently been proposed which aim at bridging the gap between the "all-atom" and "bead-spring" limits [22][23][24][25] .Nevertheless, a coarse grained model able to retain the necessary physical microscopic details while allowing simulations of the several tens or hundreds of kilo-base pairs that would be needed to address many biologically relevant questions, is still currently needed.
Examples of biological processes for which such a mesoscopic approach would be highly valuable can be classified in two broad categories: processes where DNA is mechanically manipulated by enzymatic machines (for example during replication or transcription which require opening of the double-helical structure), or processes where interactions between DNA and proteins depend more subtly on the topological and geometrical properties of the double-helix.An example of the latter class of problems is the so-called "linking number paradox", where it has been observed that the unbinding of DNA from a nucleosome releases only one unit of writhe, rather than the 1.7 which were stored 26,27 ; the resolution of the paradox is that the nucleosome also stores some twist (the terms twist and writhe are explained below).To complicate the picture even more, there are several proteins which operate to alter the DNA topology, whose collective actions may sometimes trigger complex feedback mechanisms that are crucial for biological functions 28,29 .For a model to be applicable to such problems, it must possess both a good accuracy in mimicking the geometry of the double-helix, and the ability to consider long molecules on which many proteins may act simultaneously, so that cooperative effects can be investigated.
Motivated by this goal, in this paper we introduce a single nucleotide resolution coarse-grained model for ds-DNA which retains several biologically-relevant DNA features, while being capable of delivering large-scale simulations.The model is implemented in the LAMMPS molecular dynamics engine 30 which allows us to comfortably study molecules on the order of thousands of bp (kbp).Because the code is fully parallel and highly scalable, it is portable to supercomputers to reach the length and time scales needed for some of the biological applications just mentioned.The scope of this work is to present the construction of the model, starting from the known geometry of DNA 4 (Sec.2), and to discuss the validation of its main physical features, i.e. helical pitch, persistence length and torsional rigidity (Sec.3).These properties are traditionally addressed via singlemolecule experiments 31,32 in vitro, and we here provide an indirect validation via simulated single-molecule experiments, obtaining a remarkably good agreement with the experimentally observed values (Sec.4).Finally, we present an application of this model to the dynamics of DNA denaturation, and discuss further future applications.These range from the study of DNA denaturation to that of supercoil dynamics in the presence of topological proteins (Sec.6).The flexibility of the model and the scalability provided by the LAMMPS engine means it provides a solid framework on which to base further studies of the topological properties of DNA and DNAprotein interactions.

(a) (d) (b) (c)
Fig. 1 (a) The level of coarse-graining of the model is here summarised by encapsulating the atoms forming one nucleotide into one bead-patch complex.The small yellow sphere represents the position of the phosphate with respect to the complex, while the pink sphere denotes the position of the hydrogen bond between bases.The blue sphere approximates the excluded volume of the nucleotide.(b) This panel shows the main interaction sites between consecutive beads in the same strand.The equilibrium distance between patches (E-F) is set to 0.34 nm while the one between beads centres (A-B) to 0.46 nm.This leads to an equilibrium distance of 0.7 nm between the external edge of the backbone (C-D).These distances are set so that the correct pitch of 10 bp is recovered.(c) Two nucleotides are bonded via a breakable harmonic spring.Their distance is set so that the full chain thickness is around 2 nm, as that of B-DNA.(d) Representation of the double-stranded DNA model.The red chain also shows the beads which interact sterically (solid red) as well as the phantom beads (solid grey).The faded red spheres represent the steric interaction volume of the red beads.Neither the interacting beads nor the ghost beads along the blue chain are shown to ease visualisation.

The model
We start by considering a complex made of two spherical monomers (see Fig. 1(a)), one of which represents the sugar-phosphate backbone ("bead" hereafter, shown in blue), while the other represents the nitrogenous base ("patch" hereafter, shown in violet), and is placed at a distance of 0.5 nm from the bead centre.Beads have an excluded volume so that they cannot overlap, whereas patches have no associated excluded volume.In order to see the resolution of the system, a fictitious nucleotide structure lying inside the bead is shown in black in Fig. 1(a).Although a phosphate group is not directly included in the current version of the model, is marked in yellow in Fig. 1(a) for clarity; this is sitting 0.5 nm from the bead centre but slightly away from the antipodal point to the patch.Each bead-patch complex represents a single nucleotide, and acts as a rigid body; we connect a chain of these bodies via FENE bonds of length d bp = 0.46 nm between the beads to represent one strand of DNA.We set the distance between two consecutive patches along the strand (E-F in Fig. 1 Nucleotides belonging to different strands are bonded together with breakable harmonic springs between two patches, representing hydrogen bonds (see Fig. 1(c)).The equilibrium bond distance is set to zero; if the extent of the bond increases beyond a critical value r c = 0.3 nm, the bond breaks, modelling denaturation.
While the pitch of the chain is set by the ratio of the base pairing distance and the distance between successive phosphate groups on a DNA strand, the right-handedness is imposed using a dihedral potential between the quadruplets of monomers forming two consecutive nucleotides (A, E, F and B in Fig. 1(b)).This potential regulates the angle between the planes A-E-F and E-F-B.The minimum of this potential is arbitrarily set at 36°, as this matches the geometry of a regular dsDNA helix.
In order to limit the splay of consecutive nucleotides (also called "roll" 4 ) we used a stiff harmonic potential so as to keep the angle between particles E, F and B (two patches and one bead) at 90°(Fig.1(b)).This interaction imposes the planarity between consecutive bases in the same strand.Finally, the last ingredient of this model is a Kratky-Porod potential regulating the angle between three consecutive patches along one strand.This allows us to finely regulate the chain stiffness.
The excluded volume around each bead depicted in Fig. 1(d) (faded red spheres) has diameter 1 nm.Since we use spherical beads rather than asymmetrically shaped ones (this is important for the speed of the algorithm), the geometry of the double-strand depicted in Fig. 1(b) and (d) would involve a large degree of overlapping which would lead to a large steric repulsion.To avoid this we consider two types of beads in each strand: sterically interacting beads (shown as small solid red spheres for one strand in Fig. 1(d)) are intercalated by two ghost beads (depicted as small grey spheres) which don't interact sterically along the same strand but they do interact with all the beads on the complementary strand with an excluded volume of 0.5 nm.This choice ensures that only non-overlapping beads sterically interact with one another.In addition, this allows us to preserve the correct thickness of the chain (2 nm for B-DNA), to maintain the desired distance between contiguous nucleotides and avoid the strands crossing through one another.
This model is based on few crucial geometric constraints of double-stranded DNA while the aim of delivering large-scale simulations is achieved by using spherical monomers that interact via standard potentials.These are efficiently implemented in LAMMPS and deliver a highly scalable performance in large scale parallel simulations (see Appendix A for more details).In Fig. 2 we show a typical equilibrated configuration using the presented model for a 1000 bp molecule.

Parameterisation
Our model has several parameters which can be varied to control the pitch, bending and torsional properties of the simulated DNA molecule.Nonetheless, we are interested in modelling the B form of dsDNA, of which two main physical properties are: the persistence length l p = 50 nm 150 bp, and the torsional rigidity C/k B T 60 − 80 nm 177 − 235 bp [33][34][35] .Due to the interplay between the potentials presented in the previous Section, there is no simple mapping between individual simulation parameters and the resulting physical properties; instead we obtain a simulated molecule with the correct values of l p and C via a systematic tuning of the parameters.In this Section we measure these properties from the microscopic positions of the beads in equilibrated DNA molecule configurations.Then in the following Section, we use the parametrised force field to simulate single-molecule experiments, showing that the DNA molecules show the correct macroscopic response to mechanical manipulations.

Persistence Length
The persistence length of dsDNA is a well-studied physical property that plays an important role in the wrapping of dsDNA around histone octamers to form the chromatin fibre, as well as in many other biological processes.In physical terms it gives a measure of the length-scale over which the direction of the chain is no longer correlated with itself.Following the description of an elastic rod by Moroz and Nelson 36 one can define the bending rigidity via the elastic energy functional where l p is the bending persistence length, s the arclength parameter and t(s) = dr/ds the tangent to the chain (at s) whose location in space is described by r(s).This quantity can also be readily measured by computing the tangent-tangent correlator: In our model, we use the position of the patches to extract the centreline of the dsDNA molecule, where the tangent at the nth patch at position r(n One can compute the tangent-tangent correlator along this curve and obtain the persistence length by extracting the exponent of the exponential decay.In order to avoid finite-size effects due to the presence of ends, we neglect the two terminal segments (∼5 bp at each end).The resulting curve is Fig. 3 The tangent-tangent correlator t(n) • t(n ) computed for a chain 300 bp long; it shows an exponential decay as in Eq. ( 2) with a decorrelation length lp = 143 ± 7 bp.Points show correlations measured from the simulations (average over time), and the line shows a fit to Eq. 2. Error bars give the standard error in the mean.
shown in Fig. 3.The exponential fit returns a persistence length l p 143 ± 7 bp, in agreement with experimentally observed values.

Torsional Rigidity
The behaviour of DNA when twisted is regulated by its torsional rigidity.There are several well known examples in which this property is crucial for important biological processes, such as transcription and gene expression 28,29 .Furthermore, the high torsional stiffness of DNA molecules implies that, when placed under torsion, they preferentially bend, thereby creating writhe and plectonemes 4 .In order to take this feature correctly into account, it is therefore crucial to accurately model the competition between bending and torsional rigidities 35 .
Following Moroz and Nelson 36 once again, we first define the torsional stiffness of an elastic rod C via the elastic energy functional where Ω 3 (s) is the rate of rotation of a local reference frame along the curve around the tangent t(s), defined as in the previous section.Analogous to the measurement for the bending persistence length via the tangent-tangent correlator, we here measure the torsional persistence length by computing the decorrelation of the twist angle.This correlator can be quantified by defining a local reference frame for each base pair, and tracking the transformation of the frames from one base pair to the next via the Euler angles.Each local frame is specified by the tangent vector t(n) as defined above, a normal vector f (n), defined as the projection of the vector connecting two beads in a base-pair, onto the plane perpendicular to t(n), and a third vector vector v(n) = t(n) × f (n), perpendicular to both t(n) and f (n).
The Euler angles between the frames at n and n + 1 can be used to obtain the twist increment between those base-pairs, and the correlation function for the total twist between m consecutive base-pairs Ω(m) calculated.Since the DNA has an equilibrium twist angle θ 0 = 36°per bp, we subtract this out, and calculate the correlation for the residual twist ∆Ω(m) = Ω(m) − mθ 0 .It can be in fact shown 37 (see also Appendix B) that the average cosine of the residual total twist between any two reference frames separated by m bases exhibits an exponential decay as: where we define l τ = 2C the characteristic torsional correlation length.We obtained cos ∆Ω(m) from a 300 bp long DNA molecule and averaged it over time.The curve obtained is shown in Fig. 4 on top of which we show the fitted exponential which has a characteristic decay length l τ = 512 ± 18 bp 174 ± 6 nm, which is consistent with experimental estimates valid for the B-form of dsDNA.

F F
Fig. 5 In order to simulate single-molecule experiments the model for dsDNA is anchored to a surface at the bottom end while being stretched with a constant force F from the top end.We then monitor the end-to-end elongation along the z-direction, Rz, and report its equilibrium value for a given force in Fig. 6.

Validation through single molecule experiments
Many cellular processes, such as replication and transcription, are carried out by proteins acting on single DNA segments.In light of this, recent years have seen an increasing interest in experimental techniques such as optical tweezers and atomic force microscopy, that can probe the response of DNA to external stresses (modelling the effect of DNA-binding enzymes) at the singlemolecule level.In particular, the stretching and twisting behaviour of DNA under external forces and torques has been thoroughly investigated 31,[33][34][35][38][39][40] .
In this section we aim at reproducing the conditions of two different experiments, in order to test the response of our model DNA to stretching and twisting.This also provides us with an independent method to evaluate its persistence length and torsional rigidity.In the following, we therefore keep the parameters fixed at the values used in the previous Section, and do not further tune them to achieve the experimentally known behaviours.

Response to Stretching
The classic elastic response of DNA to an external stretching force F is that of an entropic spring with relaxed length R 0 ∼ N ν with ν = 0.588 for a self-avoiding polymer.The force required to induce an end-to-end distance R z = [r(L) − r(0)] • e z for a chain of length L and persistence length l p can be approximated using the worm-like chain (WLC) result 41,42 : where excluded volume effects are neglected (a good approximation when L is not much larger than l p , as in our case).In order to test this result we performed simulations in which a constant pulling force directed along e z and acting on the last base pair of the dsDNA was applied, while the other end of the molecule was anchored at a surface (see Fig. 5).
The force-extension curve 31 observed for a chain 300 bp long is reported in Fig. 6 as data points, while the solid curve is the fit to Eq. ( 5).The fitting results in values for both L and l p , that we can compare with the values set in our model.In particular for a 300 bp chain we obtain L = 100.3± 1.7 nm (which gives a bp step size of 0.33 ± 0.01 nm) and l p = 47 ± 2 nm 140 ± 7 bp.When l p is measured from the tangent-tangent correlation for the same chain without applied force a value of l p = 49 nm was obtained.The results are therefore in good agreement with the calculation and the tuning of the persistence length performed in the previous Section.

Response to Twisting
The torsional stiffness of DNA can be calculated by computing the twist response of dsDNA to an imposed external torque, for instance applied by a magnetically con-Fig.7 The model DNA is anchored to a surface at the bottom end while being stretched with a constant force F, and a torque Γ is applied at the top end.We then monitor the linking number and report its equilibrium value for a given torque.With this information is possible to compute the superhelical density.
trolled macroscopic bead 35,43 (see Fig. 7).For different magnitudes of the applied torque, |Γ|, we compute the superhelical density, σ.The level of supercoiling is determined by the linking number Lk, the number of times one DNA strand wraps round the other.
Since a dsDNA chain has a preferred equilibrium linking number Lk 0 , the superhelical density is defined as σ = (Lk − Lk 0 )/Lk 0 .The well-known White-Fuller theorem 44 Lk = T w + W r.
relates the linking number of the edges of a ribbon (Lk) to the twist (T w), i.e. the extent of rotation of the two ribbon edges about the axis, and the writhe (W r), i.e. the self-crossing of the ribbon centreline.Although the chain we use is not closed into a loop, and therefore it is not possible to formally define a linking between the strands, it is possible to compute the linking number between two "artificially" closed strands 45,46 which follow the paths of the DNA strands along the chain backbone and then join the respective ends far away from the molecule (see Appendix C).By applying a force to the molecule, we keep it straight, consequently imposing null writhing and, in turn, ensuring that the twist is equal to the computed linking number.By measuring the deviation of twist ∆T w from the equilibrium value T w 0 defined as the num- Fig. 8 Response to torque experiment.Here we show the linear regime for small |σ|.Fitting the data points gives the torsional rigidity C using Eq. ( 7).
ber of base-pairs divided by the pitch p = 10 bp, we can readily obtain σ.
With this information it is possible to recover the response curve of the molecule to an external torque.A feature of this is a linear regime for small |σ| which we recover (see Fig. 8).The torsional rigidity, C, can finally be calculated as 36 where a 0 = 0.34 nm is the double helical rise for a relaxed dsDNA and θ 0 is the equilibrium twist angle across a base-pair step in the relaxed case.The data points shown in Fig. 8 are obtained from simulations of a 600 bp long chain anchored at a surface to one end, while the other end was pulled by a constant force of 16 pN and different torques, Γ = Γ • e z , were applied.From the fit we get the value of torsional persistence length C = 88 nm 260 bp in good agreement with experimental results [47][48][49] .One can finally use the relation between the torsional persistence length l τ and the torsional stiffness C obtained from the twistable worm-like chain theory 37 , which gives l τ = 2C = 176 nm, very close to the measurement performed in the previous section (l τ = 174 ± 6 nm).

DNA Denaturation and Supercoiling
DNA denaturation is the separation and unwinding of the two strands, transforming a DNA duplex into two isolated and unbound single strands 50 .This process can be driven by heating a solution of dsDNA molecules, and a critical "melting" temperature T m can be defined as the temperature at which 50% of a long dsDNA molecule is denatured.This critical temperature commonly depends on the genetic sequence, pH and salt concentration 47,51,52 .Localised, temporary, and dynamic denatured segments are often referred to as "bubbles".
It is well known that local denaturation has several biological implications such as favouring transcription initiation, DNA repair or recombination 28,53,54 and that the dynamics of these bubbles can be affected by torsional stress, which is itself often regulated by enzymes, such as RNA polymerases [55][56][57] .This fascinating interplay between the elasticity and biology of DNA has received much theoretical and experimental attention 47,54,[57][58][59][60][61] , but there have been remarkably few attempts to tackle it from a computational point of view 62,63 .Whereas theoretical models can capture the thermodynamics of a "stress-induced DNA-duplex destabilisation" (SIDD) 64 , elucidating the kinetics of such a process, under both equilibrium and out-of-equilibrium conditions, is an important question that can be addressed using numerical investigations.
In this Section we show that our model can readily recapitulate DNA denaturation upon decreasing the stiffness, K 2 , of the spring connecting patches in the two strands (U hb ).While the most common strategy is increasing the solution temperature, here we focus on a pathway that more closely mimics a change in salt concentration 52 or solution pH.
In Fig. 9 we show the fraction of denatured basepairs as a function of time for three different choices of K 2 .As the energy of the bond is decreased, we observe the unbinding of two strands nucleating from the ends of the chain, as observed experimentally 65 .We then observed that the denaturation spreads to the middle of the molecule, finally melting the whole chain when K 2 1.2k B T and producing two single strands.
Single stranded DNA (ssDNA) is much more flexible than its bound counterpart.In order to mimic this behaviour in our model, we eliminate both the dihedral and the Kratky-Porod interactions between nucleotides which are part of a "bubble" larger than two base-pairs.This results in single strands with a persistence length of around 2 bp which are extremely flexible, as one can appreciate from the snapshots in Fig. 9.
It is worth stressing that while our model can show the reverse of partial denaturation, by for instance increasing K 2 back to higher values, it cannot create hybridised molecules in which bases pair with partners other than those which they started with (i.e.secondary structures cannot form).This is a limit of the current model which we aim to improve in the future.It is also worth highlighting that this single-nucleotide resolution model can, in principle, readily incorporate sequence specificity.This can be done, for instance, by defining two types of harmonic bonds connecting patches in the complementary strands and by using springs with different stiffness such that K 2 (AT ) < K 2 (CG).In light of this, we expect that this model, thanks to its high scalability when run in parallel, will be of use to investigate the dynamics of denaturation in long dsDNA molecules, whether torsionally relaxed or supercoiled.
As a preliminary step to show that our model can readily take into account supercoiling, in Fig. 10 we give an example of a simulation for supercoiled DNA.A model dsDNA ring of contour length equal to 500 bp is initialised with a linking number deficit of ∆Lk = Lk 0 − Lk = −3 (47 turns instead of the usual 50 for a pitch of 10 bp).In a linear molecule this deficit would be quickly washed out by the free motion of the ends, whereas in a closed molecule, the difference creates a negative supercoiling σ = ∆Lk/Lk 0 −0.06 which is conserved throughout the dynamics.The amount of supercoiling can then be distributed into the torsional or bending degrees of freedom as long as the White-Fuller theorem 44 is satisfied (see Eq. ( 6)).Since the torsional stiffness of DNA is bigger than the bending rigidity, much of the twist is quickly converted to writhe, as can be readily seen in Fig. 10.

Discussion
The interplay between the physics and biology of DNA is one of the most intriguing topics in biophysics.While computational models can strongly aid the understanding of this fascinating open problem, the computational resources for such an expensive task have traditionally been limited.Researchers often use either very detailed and accurate all-atoms models, which can only cover short time and length scales, or coarse-grained models, which can follow the evolution of the system for much longer times, but at the expense of neglecting key physical properties of dsDNA.Mesoscopic models have been recently proposed to fill in the gap between these two approaches 22 , but they have not yet been exported to a highly efficient and parallel environment.Here, we have proposed a mesoscopic model that can be readily implemented at minimal cost into LAMMPS, one of the most popular molecular dynamics codes for atomistic and mesoscopic simulation.
Our model aims at bridging the gap between all-atoms and coarse-grained models for dsDNA; while it is currently less sophisticated than other mesoscopic models, most notably in the treatment of sequence specificity or hybridisation, this model exploits the scalability of LAMMPS, and is ideally suited to study problems such as DNA-protein interactions, or the denaturation of supercoiled DNA, where it is essential to consider long molecules, as well as to simultaneously model doublestranded and denatured regions.
This model can also be extended to include basepair specificity, and variable salt or pH concentration, while allowing the user to reach biologically relevant time and length scales.In this paper we have shown that this model is capable of reproducing DNA melting and, more importantly, of tracking the dynamics of supercoiled molecules ∼ 1000 bp long for up to ∼ 2 ms.In the near future, we aim to use this model to investigate further the interplay between denaturation and supercoiling, especially in light of its connection to gene expression 28,29 .
We should also highlight that the presented model has several limitations which arise from the compromise between accuracy and scalability.For instance, our model lacks the ability of reproducing realistic hybridisation events where distant parts of the chains can become bonded forming an intermediate hairpin.The choice of neglecting the modelling of such events allows us to employ short-ranged neighbour calculations -i.e. the algorithm does not include O(N 2 ) loops -which markedly improves the computational speed-up.
Furthermore, we also extensively tested the scalability of the model (Fig. 11).It features very good speed-up up to hundreds of processes when deployed in parallel.These results are for so-called "strong scaling" where the number of processes is increased while the total problem size, in our case the number of nucleotides, is kept constant.The scaling tests were performed on ARCHER, a Cray XC30 supercomputer with 4920 compute nodes, each consisting of two 2.7 GHz 12-core Intel Ivy Bridge processors and Aries Interconnect (Dragonfly topology).Two different benchmarks were investigated.They consisted both of linear, double-stranded DNA strands of a length of 600 bp each.The strands were initialised as a regular array of 10 × 10 or 40 × 40 strands, respectively to form a total system of 60 kbp and 960 kbp.The daily simulation times were derived from the loop timings of runs with 30, 000 timesteps (60 kbp) and 10, 000 timesteps (960 kbp) and were compared with those of a run with 24 processes (MPI-tasks), corresponding to one fully occupied node on ARCHER.We made use of the "shift" load-balancing algorithm in LAMMPS, which repositions the cutting planes between the single processes in order to mitigate a potential load imbalance between the individual processes (further details and full input files are available upon request).
For the smaller problem size of 60 kbp we observe a parallel efficiency of about 50% at 512 MPI-tasks, allowing to run for about 2 ms per day.More processes do not lead to a further speedup and the parallel efficiency decreases rapidly due to the relatively small number of "atoms" per process (LAMMPS requires several hundred atoms per process to show good scaling behaviour).The Two different benchmarks were used, a small one with 60 kbp and a 16 times larger one with 960 kbp.The results are compared with the timings of a run with 24 processes for each benchmark, corresponding to one fully occupied node on the ARCHER XC30 architecture.This leads to parallel efficiencies (see inset) in excess of 100% for 1 and 8 processes.
Total simulation times of up to 2 ms per day are feasible.
larger benchmark of 960 kbp shows a parallel efficiency of about 50% at 2048 MPI-tasks, which permits simulation times of about 0.4 ms per day.Compared to the smaller benchmark the performance degrades more slowly in this case, making simulation times of up to 1 ms per day at 8192 MPI-tasks feasible.These results strongly encourage its use on a larger scale.
Other existing models 22,23 might therefore be more suitable for studies of DNA-DNA hybridisation leading to DNA origami and synthetic DNA assemblies.The model we presented here might instead be more apt to study denaturation, supercoiling and DNA-protein interactions as previously discussed.
Finally, exploiting the ability of LAMMPS to function as a library coupled to external programs, we aim to design systems in which ATP-driven proteins interact with the model dsDNA.This paves the way to the attractive avenue of molecular dynamics simulation of large-scale out-of-equilibrium and biologically inspired systems, which are appealing to a broad range of researchers.

Conclusions
In summary, we have introduced a coarse-grained singlenucleotide model for dsDNA, which can be readily implemented in computationally efficient and parallelised engines.We tuned the model in order to reproduce the crucial physical features of dsDNA such as bending and torsional rigidities.We then tested our model by simulating single-molecule experiments so as to independently check the parameterisation and the response of our model to external manipulation.Finally, we studied denaturation and the dynamics of supercoiled DNA.We have shown that this implementation can comfortably reach length and time scales that are relevant to both single molecule and biological experiments, therefore making our model interesting for applications.In the future we intend to refine this model and to extend it in order to study biologically-inspired out-of-equilibrium scenarios.
As described in Sec. 1 the minimum of this potential is set to α 0 = 90°.
Finally, the bending rigidity is given by a potential on the angle θ formed by three consecutive patches that reads The parameters for each potential are reported in simulation units in Table 1.

Interaction Parameters
Backbone: Table 1 Parameter values in the model and expressed in simulation units.

A.3 Simulation units
the simulation units to physical ones can be done by setting the fundamental units: distance, energy and time.These are shown in Table 2.The chosen system of reference is a bath at room temperature T = 300 K and with the viscosity of water η = 1 cP.

Parameter
Experimental units Table 2 Mapping between simulation and physical units.

Distance (σ
Finally, the numerical integration is performed in an NVT ensemble by a standard velocity-Verlet algorithm with integration time-step ∆t = 0.005τ Br . (16)

Appendix B Computing the torsional persistence length
To obtain the torsional properties of the DNA molecule described in Sec.3.2 we consider a discrete elastic rod.As described in Ref. 37 , Eq. 3 is an integral over the rate of rotation of the Darboux frame (or material frame) of reference with respect to the distance along the rod.We first find the discrete approximation to this in terms of the Euler angles α n , β n , γ n which describe the rotation which generates the frame at segment n + 1 from that at segment n.To do this we make the approximation that the step size between segments is constant and denote it a; this gives where twist angle between the frames is given by α n +γ n , so the total angle between m consecutive beads is given by Ω(m) = m n=1 (α n + γ n ).An appropriate measure of the thermal fluctuations about the equilibrium twist is given by the mean of the cosine of this angle; since this quantity will decrease with m, and we identify the decay constant as the torsional persistence length cos Ω(m) = e −ma/lτ .The ensemble average is found in the usual way by taking the integral over the phase space of the system; in the small a limit this gives l τ = 2C.This is the case for an elastic rod; for the DNA molecule, the non-zero equilibrium twist between each base-pair will appear in the energy functional, so this must be subtracted from Ω(m) so that the ensemble average is a simple exponential decay.The Darboux frame at each DNA base-pair is given by the tangent vector, and the normal vector defined as the projection of the vector connecting the two beads onto the plane perpendicular to the tangent.The "closure" procedure can be performed on a pair of linear open curves to construct a closed pair whose linking number can be formally defined through the Gauss' integral (see eq. ( 17)).In this case the curves are linked once.See text for further details.
the dsDNA segment and extending the curves away from the pair of strands.Reached a certain distance by following, for instance, t 1R and t 2R , one can close the contour by defining a vector f R that joins the two new terminal beads (see Fig. 12).By following this procedure one can finally construct a pair of closed oriented curves γ R and γ B , for instance "stitching" C R , t 1R ,f R ,−t 2R , and similarly for the blue curve.Their linking number can be computed through the numerical evaluation of the double integral where r R and r B are the vectors defining the position of the segments along the curves γ R and γ B , respectively.If the centreline running through the pair of curves has no self-intersections (null writhe) then the linking number is equal to the twist.It is also worth mentioning that tightly wound curves, such as those obtained from dsDNA configurations, can lead to imprecise numerical evaluation of the integrals in eq. ( 17).In fact, the computation of Lk can become unreliable when |r R − r B | dr R dr B .The numerical evaluation can be arbitrarily improved by replacing the DNA backbones by contours more finely interspersed with points, i.e. enhancing the resolution of the integral by decreasing the infinitesimal element dr.Clearly, this can slow down the computation of Lk.We found a good compromise between and speed by adding three intermediate points every pair of beads for which we consistently measured the correct linking number during topology-preserving simulations (for instance by considering circular dsDNA).
(b)) at 0.34 nm by means of a Morse potential; the difference between the lengths A-B and E-F implies that the distance between the implicit phosphates at the external edge of the beads (C-D in Fig. 1(b)) is d ph = 0.7 nm.The ratio between d bp and d ph is well known to crucially regulate the correct pitch of the chain 4 (for details about the potentials used see Appendix A).

Fig. 2
Fig. 2 An example of an equilibrated configuration of a 1000 bp double-stranded DNA molecule, as simulated with the model presented in Sec. 1.

Fig. 4
Fig.4The average of the cosine of the total twist angle ∆Ω(m) is computed for a chain 300 bp long; in this figure we show the correlator to decay exponentially as in Eq. (4) with a characteristic length lτ = 512 ± 18 bp.Data points are obtained from simulations while the line is an exponential fit with f (d) = e −d/lτ .

Fig. 6
Fig.6Force-extension curve from the simulation (data points) and fitted by the function in Eq. (5) (solid line).The free parameters for the fitting are the total polymer length L and the persistence length lp, both of which are in agreement with the fixed parameters of the model (see text).

Fig. 9
Fig. 9 This figure shows the fraction of denatured base pairs f as a function of time and for different bond energies connecting the patches of paired bases.Snapshots from simulations are also shown.The energies used range between K2 = 1.0 kBT and K2 = 1.4 kBT .We always observe that in linear dsDNA the denaturation process nucleates from the ends, as argued from experiments 65 .

3 Fig. 10
Fig. 10 This figure shows the relaxation of a negatively supercoiled circular dsDNA.(left) The molecule (500 bp long) is initialised as a perfect ring from which three full turns are removed.(right) As the system evolves, the lack in twist is converted into writhe, and the molecule assumes stable buckled configurations.This behaviour is expected for a real dsDNA molecule because the torsional stiffness is larger than the bending rigidity.

Fig. 11
Fig.11This plot analyses the strong scaling behaviour of the model.The figure shows the simulation time in microseconds per day as a number of processes (MPI-tasks).Two different benchmarks were used, a small one with 60 kbp and a 16 times larger one with 960 kbp.The results are compared with the timings of a run with 24 processes for each benchmark, corresponding to one fully occupied node on the ARCHER XC30 architecture.This leads to parallel efficiencies (see inset) in excess of 100% for 1 and 8 processes.Total simulation times of up to 2 ms per day are feasible.

Fig. 12
Fig.12The "closure" procedure can be performed on a pair of linear open curves to construct a closed pair whose linking number can be formally defined through the Gauss' integral (see eq. (17)).In this case the curves are linked once.See text for further details.