since watson and crick proposed its 
 three - dimensional structure , much effort 
 has been spent to understand its structure and dynamics .
crick double - helix from separated strands is 
 a biologically very important process observed during replication , 
 translation , or dna repair .
also , recently , the rapidly growing 
 field of computer simulations of biological systems has developed 
 models and tools that can provide some insight into this important 
 process .
all - atom molecular dynamics ( md ) methods offer the 
 deepest insight 
 into the time evolution of dna at the atomic level .
the charmm and amber simulation 
 packages have been used to study the duplex dynamics of dna , mainly 
 in the vicinity of equilibrium .
the highlights of near - equilibrium 
 all - atom dna molecular dynamics simulations include the systematic 
 study of nearest - neighbor effects on various short dna sequences performed 
 by the ascona b - dna consortium , and the 
 microsecond simulations of the drew - dickerson dodecamer or the b dna element .
although all - atom md simulations have proven to be very 
 useful for studying near - equilibrium dynamics , their very high computational 
 cost , related to the large number of degrees of freedom and high - frequency 
 oscillations , forcing a short time - step of integration , currently 
 render them unsuitable for simulating long time - scale nonequilibrium 
 processes , such as formation of duplex dna from separate strands . 
 only recently , those enhanced sampling methods , such as replica exchange , facilitated simulations of formation of very 
 short ( 4 base pair ( bp ) ) duplexes of dna .
the combination of metadynamics and 
 replica exchange , called be - meta method , was applied to slightly larger system , dna hexamer , but the more - important biological process of formation 
 of long helices remains prohibitively expensive for all - atom md simulations 
 and has recently become a target for coarse - graining methodology .
coarse - grained models can be divided into three categories : statistical 
 models , continuum models , and reduced interaction - center models . in 
 the first category ,
the early work of zimm and rice and poland and scheraga neglected 
 all structural and dynamical information and used only the free - energy 
 gain per base - pair formation for computation of the partition function , 
 which , in turn , was transformed into a thermodynamic picture of the 
 process of helix formation .
the continuum models approximate dsdna 
 as a continuum elastic rod and by design can not be applied to the 
 dna formation process .
the last category covers all models in which groups 
 of atoms are replaced by a reduced number of interaction centers connected 
 by elastic and/or rigid virtual bonds .
a properly designed coarse - graining procedure should reduce the 
 number of interaction centers , leading to a significant speedup of 
 computations , but simultaneously , the simplified potential energy 
 function should include the most important interactions that determine 
 the behavior of the real system .
there are two common approaches to 
 the problem of parametrization of a coarse - grained model . in the  top - down 
parametrization 
 method , selected parameters are adjusted to reproduce large - scale 
 properties of the investigated systems . in
 
 parametrization method , interactions are directly determined from 
 either all - atom molecular dynamics simulations or quantum chemistry 
 computations applied to some model systems . in the past decade , 
 many coarse - grained models of dna have been 
 developed , but only several of them managed to form double - stranded 
 dna from separated complementary chains .
a three - bead - per - nucleotide 
 model of de pablo and co - workers was 
 successfully used to study dna melting and the renaturation processes .
the 
 model was designed around the potential energy function , which effectively 
 depends on the single parameter  , which was optimized to reproduce 
 experimental melting curves .
the base  base interactions were 
 approximated by a go-like potential . the model of ouldridge et al .
successfully 
 addresses the phenomena of single - strand stacking , duplex hybridization , 
 and dna hairpin formation processes .
the model was also applied to 
 study dna holliday junctions self - assembly and dna nanotweezers .
although both 
 models successfully fold double - stranded dna from separate complementary 
 chains , they rely on base  base nonbonded potentials , which 
 either distinguish between native and non - native pairs or nearest - neighbors and all remaining pairs .
such an approach should speed up the in silico 
 folding process , because many  unwanted  free energy 
 local minima are eliminated , but in the real physical system , the 
 nonbonded potential energy of two interacting bases does not depend 
 on their sequential position in the polynucleotide chain but rather 
 depends on their relative position and orientation in space .
the third model named nares-2p , which does not incorporate any 
 go-like potentials and successfully folds double - stranded dna 
 from separate strands was recently published by our group .
nares-2p model was developed in parallel to 
 the rigid - body model presented in this paper .
nares-2p is a lower - resolution 
 model with two interaction centers per nucleotide ( located on phosphate 
 groups and bases ) , and it seems to be a minimal physics - based model 
 capable of folding of a double - helix from separate strands .
the rigid - body 
 model presented here has more interaction centers and more internal 
 degrees of freedom than nares-2p , and consequently , it is computationally 
 more demanding . although the conformational space of rigid - body model 
 is larger than in case of nares-2p , it is still capable of locating 
 double - helical structures in simulated annealing process started from 
 separate strands .
the advantage of the dipolar - bead model over recently 
 published nares-2p is its higher resolution .
another model , 
 which is capable of folding dsdna from separate 
 strands , was published recently by cragnolini et al .
the degree of coarse - graining ( 67 interaction sites 
 per nucleotide ) is similar to that of our dipolar - bead model .
however , 
 as opposed to our dipolar - bead model , which places more emphasis on 
 the representation of nucleic - acid bases , the model of cragnolini 
 et al . elaborates on the nucleotide backbone , which is represented 
 by 5 sites per nucleotide than on the nucleic - acid bases , which are 
 represented by 1 or 2 interaction sites , depending on the kind of 
 base ; it also utilizes a classical bonded potential and a many - body 
 nonbonded potential responsible for correct watson 
crick hydrogen - bonding 
 and was parametrized in a  top - down  fashion .
this model 
 does not incorporate any go-like potentials and does not treat 
 the interactions between nearest - neighbor nucleotides differently 
 than those between the other pairs of nucleotides . in this paper , 

we present a physics - based , middle - resolution model 
 of dna , which is capable of folding the double - helical structure from 
 separate complementary strands .
the model was parametrized in the 
  bottom - up  fashion , although some adjustment of parameters 
 was necessary to obtain correct balance of key interactions .
the non - bonded 
 pairwise potentials do not depend on the sequential position of interacting 
 bases in the polynucleotide chain because the native and non - native , 
 nearest - neighbors , and all other pairs rely on the same potential 
 energy function , which includes a lennard - jones and an electrostatic 
 term .
the full , physics - based ( pb model ) version of model keeps stable 
 canonical b dna double helices in a long room - temperature simulations 
 and folds short dna molecules .
the reduced version of the model , in 
 which base  base intrastrand interactions are limited to nearest 
 neighbors ( nn model ) , successfully folds canonical b dna from separate 
 complementary strands for all tested systems .
the coarse - grained model of dna is built 
 of three types of particles:neutral bead , uncharged lennard - jones 
 sphere,charged 
 bead , lennard - jones sphere 
 with an electric charge located in its center,dipolar bead , lennard - jones sphere 
 with an electric dipole located in its center .
neutral bead , uncharged lennard - jones 
 sphere , charged 
 bead , lennard - jones sphere 
 with an electric charge located in its center , dipolar bead , lennard - jones sphere 
 with an electric dipole located in its center .
the dna chain is built of six types of chemical units : 
 phosphate group ( p ) , deoxyribose ( s ) , adenine ( a ) , thymine ( t ) , guanine 
 ( g ) , and cytosine ( c ) . in the coarse - grained approximation ,
atoms 
 that constitute each unit are replaced by the three types of beads 
 defined above .
the deoxyribose ring is replaced by one neutral bead , 
 the phosphate group is replaced by one charged bead and each base 
 is replaced by a set of dipolar beads .
a schematic 
 representation of the coarse - grained units overlapped on their atomic 
 representation is shown in figure 1 .
coarse - grained representation of fragments of the dna chain superimposed 
 on the atomic representation .
four basic building blocks a , t , g , 
 c attached to small fragments of backbones are shown .
charged beads , 
 which replace phosphate groups , are marked red ; neutral beads , which 
 replace deoxyriboses , are marked dark blue .
dipolar beads have different 
 colors for different bases : a , pink ; g , green ; t , cyan ; and c , yellow . 

the internal 
 degrees of freedom used for the definition of the bonded part of the 
 potential energy function ( eq 2 ) are defined 
 in tables 14 . for clarity ,
symbols of beads used for definition of bonded potentials 
 of the backbone are shown for thymine only .
the longer fragment of 
 coarse - grained strand of dna with atgc sequence overlapped on the 
 atomic model is shown in the right panel .
sequential numbers of consecutive 
 building blocks are assigned . each chemical unit in the coarse - grained approximation constitutes 
 a rigid body .
the origin of the local coordinate frame of each rigid 
 body is located at its center of mass , and its axes are aligned with 
 the principal axes of the moment of inertia tensor . the position and 
 orientation of each rigid body with respect to the global coordinate 
 frame is described by a vector 
tensor pair ( r , q ) , where r is the position of the center of 
 mass and q is the rotation matrix . in the model ,
spherical ( p and 
 s ) and planar ( a , t , g , c ) .
the spherical rigid bodies interact with 
 other particles only by central forces ( see  potential energy function  ) , and their orientation 
 is irrelevant ; therefore , a constant rotation matrix qsphere = const is assigned to them .
the orientation of 
 each planar rigid body is described by a 3  2 rotation matrix qplanar and varies as the rigid body changes its 
 orientation during a simulation run .
the geometry of the chain is 
 defined in vector  tensor space ( ri , qi )
where i = 1, ... ,n and n is the 
 number of rigid bodies .
the relations between global coordinates ( ri , qi ) and distances between various beads , used for definitions 
 of nonbonded interactions , are shown in figure 2 .
relations between global coordinates ( r , q ) and distances between beads for interactions of ( a ) two spherical 
 rigid - bodies , where there is only one distance rij = |rij| ; ( b ) spherical and planar rigid body , where distances between interacting 
 beads are functions of global coordinates given by eq 9 ; and ( c ) two planar rigid - bodies , where distances between 
 interacting beads are given by eq 13 .
rigid bodies are connected by 
 elastic virtual bonds and form a 
 nucleic acid chain as shown in figure 1 .
three 
 rigid bodies connected by two virtual bonds form a virtual - bond angle 
 and four rigid bodies connected by three virtual bonds form a virtual - bond 
 dihedral angle .
all virtual bonds , virtual - bond angles , and virtual - bond 
 dihedral angles are listed in tables 1 , 2 , and 3 , respectively .
x 
 and s  b  y 
 are virtual - bond angles between the s  b virtual bond and the 
 x and y axes of the moment of inertia tensor of a rigid body .
figure 1 ) p5s  b  x 
 virtual - bond dihedral angle describes rotation of the base around 
 the s 
the algorithm of leimkuhler and 
 reich , which is based on the rattle constraint method applied to rotation matrices , 
 was used for the description of rotation of rigid bodies .
a slightly 
 different version of this algorithm , based on the shake ( rather than the rattle ) constraint method , 
 showed superior stability compared to 
 the other commonly used quaternion - based integration scheme .
however , 
 this algorithm ( based on shake ) was not implemented in our work .
the time - evolution of the system is defined in the vector - tensor 
 phase space ( ri , pi , qi , pi ) , where pi and pi are the linear momentum and generalized angular 
 momentum of the i - th particle , respectively . for 
 a planar rigid body ,
the rotation matrix is a 3  2 tensor , subject 
 to the orthogonality condition qq = i2 , where i2 is a 2  2 unit diagonal matrix .
the leapfrog ( verlet ) 
 algorithm was used for propagation of translations , 
 and the r - rattle ( rotational rattle ) discretization 
 scheme was applied to the rigid - body rotations . as in 
 our unres model of 
 polypeptide chains developed earlier
, the effective energy function of our nucleic - acid model stems from 
 the potential of mean force ( pmf ) of polynucleotide systems immersed 
 in water .
the pmf is then to be expanded into a cluster - cumulant series , 
 which , in general , gives rise to the neo - classical local and two - body 
 and nonclassical multibody terms .
however , 
 because we keep more explicit degrees of freedom than in the unres 
 model , the multibody terms do not seem to be necessary .
in particular , 
 in contrast to the unres and nares-2p models , we do not average over 
 any degrees of freedom of the nucleic - acid bases , the interactions 
 between which are highly directional .
in the unres model , we did average 
 over the angles of rotation (  ) of the peptide groups about 
 the c ... c axes ; hence , the multibody terms in the unres force 
 field are important . in our
highly reduced 
 nares-2p model of nucleic acids , averaging 
 is carried out over the rotation of nucleic bases about their long 
 axes .
the potential energy function of the system is a sum of 
 the bonded and nonbonded parts:1although both the bonded and nonbonded parts 
 depend on the cartesian - rigid - body degrees of freedom ( r , q ) , it is more convenient to define the former as a 
 function of the internal degrees of freedom defined in tables 14 ( the internal degrees of freedom are functions of the ( r , q ) vector ) .
the bonded part has a commonly 
 used form2where:3a3b3c3dboth bond - stretching 
 and angle - bending terms 
 are approximated by harmonic potentials .
kd and di0 are harmonic force constants and equilibrium distances for bond 
 stretching interactions ; k and 
 i0 are harmonic force constants 
 and equilibrium angles for the bond - angle bending potential .
the general 
 form of the dihedral - angle torsion potentials is a fourth - order fourier 
 series , although such a long expansion was used only for the backbone 
 dihedral angles ( see details in the parametrization section ) .
the torsional 
 terms are the only higher ( second ) order 
 terms that are present in the energy function ; however , they still 
 belong to the neo - classical term set .
the harmonic uimproper potential was introduced to maintain the correct 
 position of the center of mass of the base with respect to the sugar 
 ring and surrounding phosphate groups .
k and i0 are harmonic force constant 
 and equilibrium values of the improper dihedral angles . for a detailed 
 geometric definition of the internal degrees of freedom related to 
 the bonded potentials see figure 1 and tables 14 .
the force constants 
 and equilibrium positions were determined by fitting the analytical 
 expressions to the potentials of mean force for model systems .
the 
 nonbonded energy is a sum of van der waals and electrostatic 
 energies and is given by4where the p , b , 
 and s symbols denote the phosphate 
 group , base , and sugar ring , respectively .
phosphate energy5a5binclude 
 only the lennard - jones term and depend 
 only on the distance rij between two beads .
the lennard - jones potential has a commonly used 
 form6where , alj and blj are lennard - jones 
 coefficients . the phosphate 
charge repulsion term7the three remaining energy terms of eq 4 involve interactions with the planar rigid - bodies 
 ( bases ) .
the sugar  base energy is a sum of excluded volume 
 interactions between a neutral bead and the dipolar beads .
this kind 
 of all - repulsive term is also used in the unres force field to account 
 for excluded - volume interactions between the peptide groups and side 
 chains8where nj is the number of dipolar 
 beads of the j - th 
 base and rijk is the distance 
 between the neutral bead of the i - th sugar and the k - th dipolar bead of the j - th base .
it 
 should be noted that rijk distances 
 are functions of the distance rij between the sugar and the base and the orientation qj of the base ( see figure 2).9where rjk is 
 the position of the k - th dipolar bead in its 
 rigid - body local coordinate frame .
the excluded volume potential has 
 a form10where the parameters aex and bex depend 
 on the type of 
 interacting beads .
the phosphate - base interaction additionally includes 
 the charge - dipole terms11 the final term of eq 4 has the most complicated 
 form and describes the interaction between the dipolar beads of two 
 bases12where ni and nj are the 
 number of dipolar beads located on the i - th and j - th bases , respectively .
the distances rijkl between dipolar beads are functions of 
 the relative position rij and orientations qi , qj of the rigid bodies13where ri(j)k(l ) is the position of the k - th ( l - th ) dipolar bead in the local coordinate 
 frame of i - th ( j - th ) base .
hckel model14a14b14cwhere q = eq is the negative unit charge located 
 on phosphate group , n = r / r is a unit vector pointing from the first to the second interaction 
 center , and pi(j ) are electric dipole vectors described in the global coordinate 
 frame:15where pi(j ) is an electric dipole vector in the local 
 coordinate frame of the i - th ( j - th ) 
 base .
four functions16a16b16c16dapproximate the screening of the 
 electrostatic 
 interactions by ions and solvent , where  is the debye screening length . for charge 
charge interactions 
 distance - dependent dielectric constant of the following form was applied:17where inf = 78.0 is a dielectric 
 constant of bulk water and int is a dielectric constant 
 of the helix 
the switching distances r0 and r1 , which determine boundary between unscreened 
 and screened electrostatic interactions , were set to the values of 
 4  and 13  .
the  parameter 
 was adjusted to keep continuity of dielectric constant . because 
 the interactions between the nearest neighbors , that is , 
 particles separated by less than three virtual bonds , are included 
 in the bonded part of the potential , they were excluded from the nonbonded 
 interactions .
the nonbonded interactions and 
 solvation effects are very important for dna folding and its conformational 
 stability and were approximated by the excluded - volume potential ( eq 10 ) for p  b and s  b , the lennard - jones 
 potential for the remaining nonbonded interactions ( eq 6 ) , and a multipole - multipole potential energy function with 
 debye  hckel screening .
the optimization procedure developed 
 previously was applied to the extended 
 system , which now includes the phosphate group , deoxyribose , and the 
 four bases .
the parameters of the b  b interactions , including 
 the positions r of the dipolar beads in the local 
 coordinate frame , the electric dipole vectors p 
 in the local coordinate frame , and the lennard - jones parameters alj , blj , were retained 
 from the 3445 model of ref ( 47 ) . for convenience
, they can also be found in tables a - ii 
 to a - v of supporting information .
the missing 
 parameters for the extended system are the charges located on neutral 
 and charged beads , their positions in the local coordinate frame , 
 and the lennard - jones parameters of p 
p , p  s , s  s , 
 s  b and p  b interactions .
the charge of a neutral bead 
 was , by definition , set to zero , the charge of a charged bead was 
 set to |eq| , 
 which is equal to the charge of the phosphate group in the dna backbone . 

the positions of the charges in both the neutral and charged beads 
 were fixed at the center of mass of the corresponding unit .
consequently , 
 the only parameters , which must be determined by the optimization 
 procedure were the lennard - jones alj and blj for p 
the reference energy grids for missing 
 lennard - jones interactions were computed in the same manner as the 
 reference grids in ref ( 47 ) and the optimization procedure was 
 applied to the extended system . although it was initially planned 
 that the lennard - jones potential would be applied to all pairwise 
 interactions in the system , extensive testing of dsdna dynamics showed 
 that its stability was improved when the lennard - jones potential of 
 p  b and s  b interactions was replaced by the excluded - volume 
 one .
therefore , the potential was switched to an excluded volume of 
 the form of eq 10 , and the alj and blj parameters were 
 retained from the fitting procedure , that is , aex = alj , bex = blj .
the values of the computed 
 lenard - jones and the excluded volume parameters , presented in the 
 form of  = ( a / b ) and  = b/(4a ) , are collected in table a - i of the supporting 
 information .
the debye screening length  was set 
 to 8  , which corresponds to the physiological ionic concentration . in this work ,
the 
 stability of the helix and folding efficiency was tested with several 
 values of int and various distances r0 and r1 of the dielectric 
 constant function ( eq 17 ) ( see conformational stability in the results and 
 discussion section ) .
the virtual bond , virtual - bond 
 angle , and virtual - bond dihedral angle parameters were derived to 
 reproduce the behavior of model systems in the all - atom representation . 

the potentials of mean force ( pmf ) of each internal 
 degree of freedom of the coarse - grained model were calculated by using 
 the umbrella sampling method .
analytical functions of the form 
 of eqs 3a3d were 
 then fitted to the numerical pmf s obtained from the wham procedure .
the parameters of all virtual bond , virtual - bond angle , and virtual - bond 
 dihedral angle degrees of freedom were computed by using four model 
 systems , which were obtained from three - nucleotide backbones with 
 methyl end groups .
the base of the middle nucleotide was either adenine , 
 cytosine , guanine , or thymine , for a total of four different model 
 systems .
one such system with a cytosine base in the middle nucleotide 
 is shown in figure 3 .
this is a minimal size 
 molecule for which every internal degree of freedom of our coarse - grained 
 model can be defined .
terminal bases were removed to avoid implicitly 
 including electrostatic and van der waals interactions of bases since 
 our model has the nonbonded parameters derived separately .
model systems for parametrization were created from three - nucleotide 
 chains by removing the bases from the nucleotides at the 5- 
 and 3-end and replacing them with methyl groups .
the middle 
 nucleotide was left unchanged and can be cytosine ( in the figure ) , 
 thymine , adenine , or guanine .
umbrella sampling simulations , described in the supporting information , of different model systems 
 produced 
 results that are similar for each internal degree of freedom .
figure 
 a - i in the supporting information shows 
 the individual and average pmfs for several virtual bonds , virtual - bond 
 angles , and virtual - bond dihedral angles .
therefore , it was decided to average the degrees 
 of freedom for purines ( model systems with adenine or guanine ) and 
 pyrimidines ( model systems with cytosine or thymine ) , except for the 
 degrees of freedom such as the s 
p3 
 virtual - bond dihedral angle , which contains two deoxyriboses , and 
 therefore , the base type can not be uniquely assigned .
for these angles , 
 the averaging was done across all four model systems ( bases ) .
some 
 degrees of freedom , which influence orientation of bases with respect 
 to glycosidic bond , were assigned separately for each type of base 
 ( see tables 2 and 3 ) . for each degree of freedom , the distributions of data from each 
 window and for appropriate model systems were combined , the biasing 
 potential was removed , and the pmf along that coordinate was calculated 
 by using the wham procedure implemented in a program by dr . 

the analytical expressions 
 of the form of eqs 3a3d were then fitted to pmf s . 

the fitting was done using a trust - region nonlinear least - squares fit algorithm as implemented in matlab software . 
 figure a - ii shows the comparison of several pmfs derived from the 
 simulations with the potential calculated using analytical expressions .
although both the force constants and the equilibrium positions 
 for the virtual bonds and virtual - bond angles were determined by fitting 
 to the pmfs of model systems , it was decided to shift them slightly 
 to match the equilibrium values of the  ideal  b - dna 
 conformation , in the coarse - grained model . 

b 
 with arbitrary force constant and an equilibrium values set to those 
 of the  ideal  b - dna were 
 incorporated .
the rotation of the base around the virtual glycosidic 
 bond ( s  b ) was restricted by a cosine - shaped potential with 
 the minimum corresponding to the anti conformation 
 of the base and the barrier height was preliminary set to 7.0 kcal / mol . 

the internal degrees of freedom for which the bonded terms were determined 
 and used in the model are listed in tables 14 .
we note that a model system 
 and bonded degrees of freedom almost 
 identical to those presented here were used in work by morriss - andrews 
 et al . however , the predicted values 
 for the equilibrium values and force constants were somewhat different . 

there are several reasons for this , first , morriss - andrews et al . 

use statistics from free simulations and charm27 force field to calculate 
 the degrees of freedom while this work used umbrella sampling along 
 the degrees of freedom and amber s ff99bsc0 force field .
most 
 importantly though , the position of the sugar in our model is at the 
 center of mass of sugar ring while they set it at the c1 atom . 

this affects not only the equilibrium values but the force constants 
 as well , which makes the direct comparison between bonded degrees 
 of freedom of two models difficult .
the model was tested on the following 
 tree systems , for which the structure was determined experimentally:drew 
dickerson 
 dodecamer , 
 12 base - pair dsdna with the sequence 5-cgcgaattcgc-3 and its 
 compliment , which has exactly the same sequence ( pdb code 1bna),narayana 
weiss hexadecamer , 
 16 base - pair dsdna with the sequence 5-actacaatgttgcaat-3 
 and its compliment 5-attgcaacattgtag-3 ( pdb 
 code 3bse).21 base - pair dsdna 
 with the sequence 
 5-acagcttatcatcgatcacgt-3 and its compliment 5-acgtgatcgatgataagctgt-3 
 ( pdb code 2jyk ) . drew 
dickerson 
 dodecamer , 
 12 base - pair dsdna with the sequence 5-cgcgaattcgc-3 and its 
 compliment , which has exactly the same sequence ( pdb code 1bna ) , narayana 
weiss hexadecamer , 
 16 base - pair dsdna with the sequence 5-actacaatgttgcaat-3 
 and its compliment 5-attgcaacattgtag-3 ( pdb 
 code 3bse ) .
21 base - pair dsdna 
 with the sequence 
 5-acagcttatcatcgatcacgt-3 and its compliment 5-acgtgatcgatgataagctgt-3 
 ( pdb code 2jyk ) .
the following 60-nuclotide sequence 
 of ref ( 34 ) was used 
 for computation 
 of persistence lengths of single and double stranded dna:hetss = d(catcctcgacaatcggaaccaggaagcgccccgcaactctgccgcgatcggtgttcgcct)hetds = 
 hetss + complementary strand hetss = d(catcctcgacaatcggaaccaggaagcgccccgcaactctgccgcgatcggtgttcgcct ) hetds = 
 hetss + complementary strand for canonical ensemble simulations , berendsen s 
 thermal - bath weak - coupling algorithm was 
 implemented .
the stability of the integration algorithm was checked by simulations 
 of the model system in the microcanonical ensemble .
the energy of 
 the system was minimized and initial velocities were assigned from 
 the maxwell distribution ; then , the integration of the equations of 
 motion was performed without coupling of the system to the thermal 
 bath .
the ability of the model to fold double - stranded dna from 
 separated 
 chains was tested with a two - stage simulated - annealing procedure .
two complementary strands of each molecule were separated 
 by a random distance and randomly rotated with respect to each other . 

the separation distance ensured that two strands did not overlap for 
 any random orientation and was different for chains of various lengths . 

for 1bna , the 
 separation distance between centers of mass was picked randomly from 
 40 to 45  interval and for 3bse and 2jyk separation intervals were 5055 
  and 6065  , respectively . in the first stage , 
 the initial structure was heated up to 600 k ( 620 for 2jyk ) with velocities 
 randomly assigned from maxwell distribution for 0.1 ns and then cooled 
 to 300 k ( 320 for 2jyk ) with a temperature step t = 50 k. the final 
 temperature was set approximately around 20 k below melting temperature 
 of each molecule .
the simulation at each temperature was 0.1 ns long . 
 at the end of each simulation ,
a contact was assumed 
 to be present when the centers of mass of any two bases from separate 
 chains came within a distances of at most 6.5  .
if at least 
 one contact was present , the simulation was extended by 1.5 ns for 
 drew  dickerson ( 2.0 ns and 2.5 ns for narayana 
weiss 
 and 2jyk , respectively ) 
 at 300 k in the second stage of the canonical simulation ; otherwise , 
 the cycle was terminated and the final structure was excluded from 
 the analysis .
the total simulation time of a two - stage simulated annealing 
 cycle was 2.23.2 ns ( depending on the chain length ) , and each 
 one was followed by a short energy minimization with powell s 
 algorithm .
the number of sa cycles performed for each system 
 is shown in table 5 .
dissociation of chains 
 was prevented by application of a flat - bottom restraint on the distance 
 between their centers of mass with a force constant kcm 
cm = 1 kcal/(mol  ) and a minimum 
 distance of d0 = 40  for 1bna ( 50  for 3bse and 60  for 2jyk ) , at which the attractive 
 harmonic force was switched on .
the artificial restraint force was 
 switched off when the distance between the centers of mass dropped 
 below 40  for 1bna ( 50  for 3bse and 60  for 2jyk ) .
the first row shows 
 the total numbers 
 of simulated annealing cycles computed for each molecule .
the second 
 row shows the number of trajectories for which the cooling cycle ended 
 up with at least one contact between bases .
statistics of contacts 
 for these trajectories , which were extended by 1.5 ns to 2.5 ns canonical 
 simulation , is shown in the next 5 rows .
numbers of long trajectories 
 with various percentage of native contacts are presented in rows 37 .
the final structures of the 
 simulated annealing cycle were compared 
 to the experimental reference by computing the all - bead rmsd .
the 
 persistence length lp is a measure of 
 flexibility of a polymer and can be computed from decay of correlation 
 function of unit vectors tangent to helical axis18where s is the position along 
 the length of the chain . for ssdna
the 
 tangent vector is defined as a normalized vector between consecutive 
 neutral beads ( centers of mass of deoxyriboses ) . for dsdna tangent 
 vector
sugar 
 distances of complementary bases , separated by 5 base pairs . a new software package for this coarse - grained 
 model , composed of rigid - bodies , was developed .
the software is written 
 in c / c++ programming language and currently can perform the following 
 tasks : read the pdb - format structure and transform it into the coarse - grained 
 representation , optimize the geometry of the molecule ( by energy minimization ) , 
 perform microcanonical and canonical ensemble simulations with berendsen s 
 weak thermal - bath coupling , and perform simulations of heating and 
 simulated annealing cycles .
the 
 efficiency of the parallelization up to 12 cores is shown in figure 4 .
the stability 
 of the verlet - r - rattle algorithm was tested on the 
 drew  dickerson double - helix .
microcanonical simulations of 
 10 ns length were run with various time - steps : dt = 10 , 12 , 14 , 16 , 18 , 20 , 22 , 24 , 26 , 28 fs .
the algorithm 
 maintains good stability up to dt = 22 fs resulting 
 in an order of magnitude speedup compared to the commonly used all - atom 
 time - step of 2 fs ( with shake applied on hydrogens ) , and shows superior 
 stability compared to the quaternion - based rigid - body integration 
 scheme , applied previously to the model 
 of a dna chain .
the drift of the energy over a 20 ns trajectory with 
 dt = 22 fs was only 0.017 kcal / mol , compared to 10 
 kcal / mol over 1 ns simulation with the same time - step using a quaternion 
 integration scheme ( see figure 5 in ref ( 19 ) ) .
the drifts of the energy 
 and the rms fluctuations for various time - steps are collected in table 6 .
all canonical simulations and simulated annealing 
 cycles utilized time - step dt = 10 fs .
stability of the verlet - r - rattle 
 integration algorithm over 0.1 
 s microcanonical simulation of the drew 
drift was calculated as 
 a difference between the averages of the last and the first 1 ns of 
 each trajectory .
the efficiency 
 of the model was compared to that of all - atom simulations 
 performed with the amber 11 package .
dickerson 
 dodecamer immersed in a truncated octahedron box filled with 5099 
 tip3p water molecules .
a 10  cutoff and the particle - mesh ewald 
 summation method were applied for electrostatic interactions .
the 
 simulation of 1 ns took 17476 s on 8-cores using intel xeon e5645 
 2.4 ghz cpus of the computing cluster of the chemistry research computing 
 facility of baker laboratory of chemistry and chemical biology , cornell 
 university .
the same simulation with our coarse - grained model and 
 dt = 20 fs took 235 s with the same cpu , leading 
 to almost 2 orders of magnitude speedup .
it should be noted that almost 
 2 orders of magnitude speedup was achieved without any cutoffs applied 
 to the nonbonded interactions .
implementation of cutoffs , especially 
 to the short - range lennard - jones ( u  r ) and the dipole 
dipole ( u  r ) interactions 
 will further increase the speed of computation , especially for long 
 chains .
the initial 
 crystal structure ( 3bna , 3bsa , and 2jyk ) was heated from 
 100 to 300 k with t = 50 k. the simulation 
 at each temperature was 0.1 ns long , and the production run at 300 
 k was 0.1 s long .
snapshots of trajectories were recorded every 
 1 ps . the initial simulations of the drew 
dickerson dodecamer 
 led to a collapse of the helical structure , caused by overestimation 
 of the strength of the lennard - jones interactions with phosphates 
 and sugars , which were fitted to reproduce the amber potential energy 
 in vacuum .
the stability of the system was improved by scaling of 
 lennard - jones interaction energy between neutral and charged beads 
 ( p  p , s 
the stability of a double - helix was also tested with varying r0 and r1 parameters 
 of eq 17 .
table 7 shows 
 the average number of native contacts as a function of distances r0 and r1 for drew - dickerson 
 dodecamer with int = 8.0 .
the largest average number 
 of native contacts is located in the lower right corner of table 7 , which means that sufficient electrostatic repulsion 
 of phosphate groups of interacting chains is necessary for stability 
 of double - helix .
smaller switching distances r0 and r1 causes instability of 
 the helix termini , which can even destabilize helix interior . on the 
 other hand
the reduced ( and less physical ) model , in which intrastrand 
 base  base interactions were limited only to nearest - neighbors , 
 was tested .
the nn model generates very stable double helices even 
 for short switching distances r0 = 4  
 and r1 = 13  , as applied in the 
 ref ( 40 ) .
it seems 
 that stronger electrostatic repulsion of phosphate groups is necessary 
 to prevent deformations of helix termini caused by base
, stronger electrostatic 
 repulsion affects folding rate in the fast simulated annealing procedure 
 ( see section folding of dsdna from separate chains ) .
table 8 collects mean helical parameters 
 obtained from trajectories for three model systems for full physics - based 
 model and for the simplified nearest - neighbor model .
these parameters 
 were calculated from coarse - grained representations of experimental 
 reference structures and the structures obtained in simulations , respectively . 

a precise determination of all helical parameters would require conversion 
 of the coarse - grained structurs to an all - atom representation ; however , 
 the values calculated by using a coarse - grained representation are 
 sufficient for comparison of the simulated structures with the respective 
 experimental structures . for all model systems ,
the average number 
 of native contacts is larger for the nn model . for shorter molecules , 1bna and 3bse ,
the radius of gyration 
 within nn approximation is close to the value obtained for experimental 
 structure . for the pb model
applied to 1bna and 3bse the radius of gyration is slightly too 
 short , which is caused by deformation of helix termini . for the longest 
 molecule both models overestimates the radius of gyration .
both models 
 underestimate helix radius and overestimate helical rise for all dsdnas , 
 although these discrepancies are relatively small .
larger deviations 
 are observed for the size of minor groove , which is significantly 
 overestimated by both models .
pb and nn denotes the full physics - based and the nearest - neighbor 
 model , respectively .
the 
 helical rise is approximated 
 by the distance between consecutive centaral points of lines connecting 
 centers of mass of watson  crick - paired bases .
persistence lengths 
 of both hetss and hetds were computed from 0.1 
 s canonical 
 simulations at 300 k with 145 mm salt concentration for both full 
 pb and simplified nn model .
pb model produces smaller persistence 
 lengths than nn model , although the relative difference is much bigger 
 for ssdna .
this behavior is expected , because switching off some of 
 the intrastrand base  base interactions should decrease bending 
 of the chain , especially for ssdna , in which bases are exposed to 
 solvent and can easily form intrastrand hairpins with complementary 
 bases .
figure 6 presents typical configurations 
 of het60ss and het60ds molecules taken from t = 300 k canonical simulation .
the conformation of ssdna 
 taken from full pb model simulation differs significantly from the 
 one taken from nn approximation .
the ssdna obtained from full pb model 
 bends and forms hairpins , which can act as kinetic traps in the folding 
 process of dsdna .
typical configurations of het60ds ( left panel ) , 
 het60ss within nearest - neighbor approximation ( middle panel ) 
 and 
 het60ss for full physics - based model ( right panel ) .
persistence lengths were obtained 
 by fitting of eq 18 the 
 the correlation functions obtained from canonical ensemble trajectory 
 computed at 300 k with 145 mm salt concentration .
the value of persistence length 
 of ssdna obtained with nn model 
 agrees reasonably well with experimental values .
the persistence length of dsdna appears to be around three 
 times too small compared to experimental values , and it is close to the value 
 computed by early version of 3spn model .
although the agreement with experimental data is not perfect it 
 appears to be quite satisfactory for  bottom - up-parametrized 
 physics - based coarse - grained model and this discrepancy will be corrected 
 in next versions of the model .
the 
 influence of mismatches was tested on 
 the 30-bp system of 10 g / g pairs flanked by 10 a / t pairs at both 5 
 and 3-end . during the 0.1 s canonical simulation at 
 300 k , opening of mismatched g / g pairs
although 
 interaction of complementary bases is different than that between 
 mismatched bases , the model does not seem to be very specific with 
 respect to base  base interactions .
this problem may be result 
 of exaggeration of base  base lennard - jones interactions with 
 respect to electrostatics and will be addressed in the next version 
 of the model .
the simulated 
 annealing ( sa ) procedure described in the methods section was applied to some model systems .
only 3 out of 2100 trajectories led 
 to the final helical structures with more than 10 native base 
much better folding rate was achieved with nn model 
 in which intrastrand kinetic traps ( hairpins ) , resulting from intrastrand 
 long - range base  base interactions , were eliminated .
the statistics 
 of the folding process of three model systems within nn approximation 
 are presented in table 5 .
after the majority 
 of simulated annealing cooling cycles ( first stage ) , the chains did 
 not come into contact because of electrostatic repulsion of the highly 
 charged chains and the relatively short time of the simulation .
the 
 frequency of making at least one contact after a cooling cycle is 
 different for each simulated system and changes from 50% for the dodecamer 
 to 6% for the hexadecamer , and to 8% for the 21-bp system .
this phenomenon 
 can be explained by increasing volume of an accessible space for longer 
 chains , which decreases probability of making contact during fast 
 cooling process .
the effect is enhanced by increased electrostatic 
 repulsion of the longer negatively charged backbones .
after 
 the second stage of simulated annealing , native - like structures were 
 obtained .
the percentage of the native - like structures , based on the 
 number of native contacts , is shown in table 5 .
it should be noted that the numbers summarized in table 5 were obtained from sa calculations which do not 
 reflect the actual amount of the native - like structures that a force - field 
 can produce .
nevertheless 24% of extended 
 trajectories led to double - helical 
 structure with more than 10 native contacts for 1bna molecule .
roughly 
 the same number of structures broke up into separate chains or formed 
 misfolded structures without native contacts .
41% trajectories ended 
 up in misfolded structures with 13 native contacts , which 
 were mostly parallel dsdna s , in which 5-end of one 
 chain makes nonative contacts with 5-end of the other chain . 

figure 7 shows the potential energy of final 
 structures as a function of all - bead rmsd computed with respect to 
 experimental structure .
the one located in the lower left corner of figure 7 is a cluster of correctly folded antiparallel dsdna s 
 and the other one represents misfolded parallel dsdna s , in 
 which 3-ends ( and 5-ends ) of chains are paired .
the 
 important observation that comes from this graph is that lowest energy 
 structure is clearly located in the correctly folded group and its 
 energy is around 10 kcal / mol lower than misfolded parallel dsdna , 
 although the lowest energy structure is not the structure with lowest 
 rmsd .
the mean potential energy of antiparallel dsdna cluster is around 
 8 kcal / mol lower than the mean energy of the parallel dna cluster , 
 which clearly favors correctly folded structures , although this difference 
 for real systems might be larger ( experimental structures of the parallel 
 drew  dickerson dodecamer were not observed ) and further refinement 
 of the electrostatics and/or lennard - jones balance of the model is 
 required to increase this gap and improve specificity of base  base 
 interactions .
two major 
 clusters , corresponding to antiparallel and parallel dsdna , and the 
 energy - minimized experimental structure are pointed by arrows .
the reference structure marked 
 by bigger dot in figure 7 is the energy - minimized 
 experimental structure . 

the optimization of experimental structure was necessary to remove 
 small steric clashes in original structure caused by imperfectness 
 of the coarse - grained model .
the potential energy of original 1bna 
 molecule is more than 100 kcal / mol higher than that of the optimized 
 one , but the all - bead rmsd between the optimized and original structure 
 is only 0.4  .
this effect can be expected for the coarse - grained 
 model with force - field parameters developed in the  bottom - up  
 fashion .
the potential energy of the reference structure is almost 
 exactly the same ( < 1 kcal / mol ) as the lowest - energy structure obtained 
 from simulated annealing procedure .
the folding efficiency of 
 longer polymers decreases as the dimensionality 
 of phases spaces increases ( see table 5 ) .
16% 
 of long trajectories of 3bse molecule led to final structures with more than 75% 
 of native contacts ( 1316 contacts ) .
for the 2jyk molecule , the folding 
 rate is almost the same as that for 3bse and equals 15% .
nevertheless , the energy 
 vs rmsd graphs ( figures 8 and 9 ) are qualitatively similar to these obtained for 1bna molecule 
 ( figure 7 ) . in both cases , two clusters of 
 antiparallel and parallel dsdnas are clearly visible and correct 
antiparallel 
 cluster has mean potential energy around 20 kcal / mol lower than this 
 of incorrect 
these values should effectively 
 discriminate misfolded structures in the real- , long - time folding 
 process .
this mean - energy difference increases to 22 kcal / mol for 
 the longest 2jyk molecule .
the energy - minimized reference structures , marked by bigger 
 dot in figures 8 and 9 , have rmsd only 0.6  away from original structure . in both 
 cases ,
the reference structures have potential energy around 10 kcal / mol 
 lower than these obtained from simulated annealing procedure .
this 
 observation suggests that longer cooling process should lead to improvements 
 in geometry of final structures , as the minimum of potential energy 
 was not reached by any simulated annealing cycle for longer model 
 systems .
figure 10 presents the three lowest - energy 
 structures of tested molecules superimposed on their experimental 
 structure .
the all - bead rmsds of 1bna , 3bse , and 2jyk molecules with respect to experimental 
 references are 2.1  , 3.1  , and 4.2  , respectively , 
 and as expected , they increase with the length of polymers .
all structures 
 clearly form double - helices , although some discrepancies of the geometry 
 are visible , especially for the longest molecule , for which overestimation 
 of the size of minor groove is clearly visible .
lowest - energy structures 
 of 1bna , 3bse , and 2jyk molecules obtained 
 from extended simulated 
 annealing procedure , superimposed on experimental structures .
their 
 all - bead rmsds are 2.1  , 3.1  , and 4.2  , respectively . 

the backbone of the computed 
 structures is shown in blue and the bases are a ( cyan ) , t ( pink ) , 
 g ( green ) , c ( yellow ) .
analysis of 158 trajectories of 1bna , which led to structures with 1012 
 native contacts revealed two qualitatively different mechanisms of 
 folding .
the representative of the first most common mechanism 
 ( 93 of 158 
 trajectories ) is shown in figure 11 .
two separated 
 chains ( 0 ns ) were heated to 600 k and evolved into random chains 
 ( 0.11 ns , t = 550 k ) . at 0.56 ns ( t = 350 k )
the 5-gcgc - end established native contacts with 
 the complementary 3-cgcg - end , and since then , more native 
 contacts were formed , and the dsdna molecule was zipped . at 0.68 ns , 
 all native watson 
snapshots of trajectory of a zippering mechanism of the 
 folding 
 process of the drew  dickerson dodecamer . a different and less common folding process ( 65 of 158 trajectories ) 

it was started from 
 two separated ssdna chains , which were heated to 600 k and , after 
 0.3 ns ( t = 500 k ) , evolved into a completely random 
 conformation .
first , contacts between the bases of the 5-ends 
 of the chains were established after 0.49 ns . it should be noted that 
 these are non - native watson  crick g / c contacts .
after that , 
 the chains started to slide with respect to each other along a watson 
after 0.89 ns , two g / c contacts were estabilished but also 
 a / c and a / g mismatched pairs were formed .
a 1.4 ns snapshot shows 
 a molecule with 10 contacts established and with distorted termini . 

snapshots of trajectory of a slithering mechanism of the 
 folding 
 process of the drew 
as can be expected , the zippering mechanism was also observed 
 for 
 longer molecules , 3bse and 2jyk . 
 nevertheless , some trajectories with a slithering folding mechanism 
 were also recorded . at first , this observation is surprising because , 
 as opposed to the 1bna molecule , the initial mis - alignment of the two strands creates a 
 number of mismatched base 
therefore , the fact that 
 the simulated folding of 3bse and 2jyk proceeded almost as frequently by the zippering and by the slithering 
 mechanism probably resulted from low base  base pairing specificity 
 of the model . the snapshot from helix initiation process for 3bse molecule
the 3-end of a complementary 
 strand made contacts with 3-end of the leading strand . besides 
 three nonative watson 
crick contacts ( two a / t and one g / c ) 
 and two mismatched contacts
three 
 watson  crick pairs and two mismatched pairs form initial structure 
 from which monomers start slithering with respect to each other until 
 final proper double - helix is formed .
the helix nucleation occurred 
 between nucleotides of 3-ends but also the mechanisms where 
 3-end ( or 5-end ) of one chain made a contact with 
 midsection of the other chain were observed .
examples of slithering 
 pathways of dna hybridization of ( a ) 3bse and ( b ) 2jyk molecules . for the longest chain ,
the slithering mechanism was observed 
 only 
 for two out of 22 folding trajectories . in one case ,
nucleation started 
 from formation of contacts of 3-end of one chain with middle 
 section of the other chain ( see figure 13b ) . 

in the second slithering trajectory , the 5-end made initial 
 contact with middle section of the other chain .
no trajectory where 
 two 3-ends ( or 5-ends ) made initial contact was recorded . at this point ,
the analysis of dna duplex hybridization is reduced 
 only to qualitative observations presented above .
further development 
 of the model is necessary , especially further reduction of computational 
 cost via implementation of cutoffs , to obtain large set of trajectories , 
 which is required for proper analysis of thermodynamics and kinetics 
 of dna hybridization process .
nevertheless , the mechanism of 
 hybridization can be qualitatively 
 compared to hybridization mechanisms proposed by de pablo and doye groups . 

the former group showed two mechanisms also called slithering and 
 zippering , which follows double - helix nucleation process .
the former mechanism is dominant for repetitive 
 sequences and zippering is characteristic of random sequences .
both 
 mechanisms seem to be similar to folding mechanisms obtained with 
 our dipolar bead model , although it should be noticed that sequences 
 of simulated systems are not exactly the same . the drew 
dickerson 
 dodecamer is symmetric and might have pathways of folding that are 
 similar to slithering mechanism observed for 3spn model .
our model 
 shows that slithering and zippering mechanisms for 1bna are almost equally 
 probable .
the 3bse and 2jyk molecules 
 have a sequence that can be described as random , and their dominant 
 folding mechanism is molecular zippering . the zippering mechanism 
 was also observed for the model of ouldridge 
 et al .
they showed that there is a bias 
 toward initial structures with contacts between ends of the strands , 
 an observation that can be confirmed by our model .
one of 
 the mechanisms of duplex formation observed for repetitive sequences , 
 so call inchworm mechanism , was not observed in our simulations , although 
 we did not test polymers with exactly the same sequence .
it should also be noticed that tendency for slithering , 
 instead of inchworm - like behavior , in our model may result from insufficient 
 specificity of base  base interactions .
finally , it 
 should be stressed that these results were obtained 
 in unrestrained simulations started from separated chains .
this mode 
 of simulations is different from that of many other works in which 
 simulations are started from the experimental structure . unlike many 
 other force - fields developed for dna , our force - field , therefore , 

, it was shown that dihedral 
 angle bonded potentials are not prerequisites for double - helix formation , 
 which is rather driven by base  base dipole interaction .
also , the model developed by ouldridge et al . 
 does not require backbone dihedral angle potentials to obtain dna 
 duplex hybridization .
a question can , 
 therefore , be asked : are the backbone potentials a necessary factor 
 for double - helix formation , or is this process solely dependent on 
 the base stacking and pairing ?
p3s3 , 
 s5p5s  p3 , s5p5s  b , and 
 s3p3s  b dihedral angle potentials switched 
 off .
dickerson dodecamer showed 
 the same efficiency of double - helix hybridization process ( 24% ) .
the 
 lowest energy structure shown in figure 14a 
 is only 3.2  away from the experimental structure .
next , improper 
 rotamer - like potentials p5p3b  s , p5s  b , 
 and p3s  b were additionally switched off .
the lowest 
 energy structures have a ladder - like shape , as shown in figure 14b .
it seems that backbone dihedral angle potentials 
 are not really necessary for dna duplex hybridization .
switching off 
 rotamer - like potentials ( p5s  b , p3s  b , 
 and p5p3s  b ) causes formation of ladder - like 
 structures instead of double - helices .
this effect might be related 
 to the lack of specificity of base  base interactions , which 
 might be the result of overestimation of base  base van der 
 waals interactions with respect to dipole  dipole electrostatic 
 interactions .
nevertheless , the dipolar - bead model suggests that only 
 backbone bond - stretching , angle - bending , and rotamer - like potentials 
 are necessary for double stranded dna hybridization .
lowest energy structures 
 obtained from simulated annealing procedure 
 with the nn model with some of bonded potentials switched off .
( b ) additionally , rotamer - like potentials p5s  b , 
 p3s  b , and p5p3s  b switched 
 off .
a 
 physics - based rigid - body coarse - grained model of dna is proposed . 

the bonded interactions are fitted to the potentials of mean force 
 of the model system .
the equilibrium virtual - bond lengths and virtual - bond 
 angles were set to the values of the ideal b - dna double - helix .
the 
 nonbonded interactions are approximated by lennard - jones , excluded 
 volume and electrostatic interactions of charges and dipoles .
the model does not incorporate any go-like potentials , 
 which force the structure toward the experimental one . in the full 
 physics - based model , nonbonded potentials do not distinguish between 
 near - neighbors and any other base pairs .
a full physics - based version 
 of the model the nearest - neighbor approximation , in which intrastrand 
 base  base interactions are limited to nearest - neighbors only , 
 was tested .
it should be stressed that in the reduced model all interstrand 
 base  base interactions are present . the r - rattle rigid - body integration algorithm , 
 implemented in our software package , demonstrated superior stability 
 compared to a quaternion - based algorithm , with low drift of the total 
 energy even for a time - step as long as 20 fs .
the total speedup of 
 computations compared to all - atom simulations is close to 2 orders 
 of magnitude .
stable double - helix trajectories 
 were obtained for both pb and 
 nn models after adjustment of several force - field parameters as discussed 
 in the conformational stability section . 

the geometry of two experimental dna structures is preserved during 
 0.1 s simulations at 300 k. the persistence lengths of ssdna 
 and dsdna are underestimated by the factor of 23 compared 
 to experimental data , discrepancy that seems to be reasonable for 
 physics - based model designed in the  bottom - up  fashion . 

the other drawback of the model , which seems to have the same origin 
 as discrepancy in persistence lengths , is insufficient specificity 
 of base  base interactions , which affects simulations of dsdna 
 with mismatches .
opening of mismatched pairs was not observed during 
 0.1 s canonical simulations at 300 k. nevertheless , the 
 model successfully folds short ( 1221 
 bp ) dnas from separated strands .
the full pb model is only able to 
 fold shortest double - stranded dna dodecamer from separated chains 
 with relatively low efficiency caused by formation of intrastrand 
 hairpins  kinetic traps in the folding process .
the reduced 
 nn model removes kinetic traps and greatly improves the folding efficiency . 

all three tested molecules with lengths varying from 12 to 21 bp were 
 folded from separate strands by nn model , although the efficiency 
 of folding drops with increasing size of the simulated system .
the 
 possible explanation is that the cooling process is too fast to sample 
 conformational space effectively and find the global minimum on the 
 complicated multi - minima free energy hyper - surface .
application of 
 a replica - exchange method , which will be implemented in our software 
 package , should help to explain the cause of the low helical populations 
 produced by the simulated annealing procedure .
the more 
 frequent process called zippering starts from formation of native 
 base  base contacts between 3 and 5-ends of 
 two chains and is followed by helix propagation . in the second , less 
 frequent process called slithering , initial non - native watson
 crick 
 contacts ( also mismatched contacts for 3bse and 2jyk ) are established between 5-ends 
 ( or 3-ends ) of two strands ; then , chains slide with respect 
 to each other along the watson 
this process was observed for all model systems 
 although it is much more frequent for the dodecamer , which has specific 
 sequence , than for other molecules .
the slithering process might also 
 be an artifact caused by insufficient specificity of base 
base 
 interactions . finally , the influence of backbone dihedral angle 
 and rotamer - like 
 potentials on a double - helix formation process was tested .
the backbone 
 dihedral potentials do not seem to be a necessary prerequisite for 
 dsdna hybridization , but without rotamer - like potentials , ladder - like 
 structures instead of double helices were observed .
base 
 interactions , underestimation of persistence lengths of both single 
 and double - stranded dna , and small influence of mismatched pairs on 
 double - helix dynamics .
the origin of mentioned problems most probably 
 arises from overestimation of lennard - jones base  base interactions , 
 and we will address these problems in the next evaluations of the 
 dipolar bead model .
 a middle - resolution 
 coarse - grained model of dna is proposed . 
 the 
 dna chain is built of spherical and planar rigid bodies connected 
 by elastic virtual bonds . 
 the bonded part of the potential energy 
 function is fit to potentials of mean force of model systems . 
 the 
 rigid bodies are sets of neutral , charged , and dipolar beads . 
 electrostatic 
 and van der waals interactions are parametrized by our recently developed 
 procedure [ maciejczyk , m. ; spasic , a. ; liwo , a. ; scheraga , h.a . j. comp . 
 chem.2010 , 31 , 1644 ] . 
 
 interactions with the solvent and an ionic cloud are approximated 
 by a multipole  
 multipole debye  hckel model . 
 
 a very efficient r - rattle algorithm , for integrating 
 the movement of rigid bodies , is implemented . 
 it is the first coarse - grained 
 model , in which both bonded and nonbonded interactions were parametrized 
 ab initio and which folds stable double helices from separated complementary 
 strands , with the final conformation close to the geometry of experimentally 
 determined structures . 
Introduction
Methods
Results and Discussion
Conclusions