In Silico Single-Molecule Manipulation of DNA with Rigid Body Dynamics
[['The mechanical and topological properties of DNA and protein-DNA assemblies are of primary importance in many biological processes , including transcription , replication , chromatin organization and remodeling .', 'Since techniques have become available enabling the manipulation of single-molecules [1] , [2] , a large amount of experimental data have been accumulated on the mechanical response of DNA and protein-DNA assemblies under stretching forces and twisting torsions , in particular from optical and magnetic tweezers experiments [2]–[5] .', 'In magnetic tweezers experiments , a DNA molecule is grafted at one end to a coverslip and at the other end to a magnetic bead .', 'The bead is trapped in the magnetic field of a pair of magnets that may be translated , thus exerting a varying force on the bead .', 'Moreover the pair of magnets may be rotated at a certain number of turns , thus constraining the linking number of the DNA molecule .', 'After the stretching force and the number of turns are applied to the bead , the only physical variable that can be directly measured is the DNA extension , i . e . the distance between its two ends .', 'Therefore the interpretation of the experimental results requires an important modeling effort , particularly in the more complex cases where DNA is associated with proteins , as for instance in chromatin assemblies [6] , [7] .', 'Although theoretical approaches may be successful in some cases [8]–[10] , simulations are often crucial tests of the proposed model validity , when they are not the unique possible way of dealing with the system complexity .', 'In this spirit , we aim to develop an efficient tool to manipulate single-molecules in silico reproducing optical and magnetic tweezers experiments .', 'This task is challenging since the DNA model should have precise specifications to reproduce the behavior of DNA accurately .', 'We need to:', '( i ) model a polymer , i . e . an articulated chain;', '( ii ) reproduce the effective diameter of DNA ( depending on electrostatic conditions ) and , when proteins are present , have the possibility to model their shape and steric hindrance;', '( iii ) deal with collisions , especially in order to reproduce DNA supercoiled structures ( plectonemes ) and steric effects in DNA-protein assemblies;', '( iv ) reproduce DNA twisting and bending elasticities;', '( v ) include statistical mechanics features to account for temperature and thermal motion .', 'Beside these essential points , we also wish to simulate the system dynamics , which may be important in some cases , e . g . when hysteresis is observed under magnetic tweezers [7] or for in vivo chromosome dynamics experiments in the cell nucleus [11] , [12] .', 'This ambitious list of specifications is beyond the reach of traditional simulation approaches where particles interact through 2-body potentials ( as in Molecular Dynamics or Monte Carlo simulations [13] , [14] with a given force field ) .', 'The need to deal with frozen degrees of freedom in coarse grained modeling may be addressed through holonomic constraints , as in the SHAKE algorithm [15] , [16] , where an iterative approach is adopted .', 'However , collision detection and steric hindrance may only be accounted for in this scheme by introducing additional steps .', 'More recently , non iterative algorithms have been developed [17] , that subsequently led to the development of new powerful tools , called “physics engines” .', 'These have been designed by the engineering and robotics communities to reproduce the dynamic behaviour of articulated systems of rigid bodies .', 'Physics engines are acquiring an increasing importance , notably in the fields of computer graphics and video games , where they are now widely used to simulate rigid body motion under realistic conditions and in real-time .', 'Open Dynamics Engine ( ODE ) is one of the most popular rigid-body dynamics open source library for robotics simulation applications [18] .', 'As other physics engines , ODE simulates the kinematics of articulated systems by using permanent joints that impose holonomic constraints , instead of bond potentials .', 'The same method is used to manage collisions: when overlapping between bodies is detected , a temporary joint is locally created that reproduces the action of the contact forces , without the need for explicit permanent 2-body interaction potentials ( see section “Materials and Methods” for details on how ODE manages joints and collisions ) .', "These extremely efficient simulators haven't , up to now , been used in statistical mechanics .", 'Although well adapted to mechanical simulations , physics engines lack coupling to a thermal bath .', 'The main novelty of our approach is the implementation of Langevin-Euler equation in the ODE software .', 'Moreover we improve the simulation efficiency of this Langevin dynamics by extending the “global thermostat” algorithm designed by Bussi and Parinello in 2008 [19] to physics engines .', 'This algorithm allows much faster yet unbiased sampling of the phase-space .', 'As a first step toward simulating DNA-protein assemblies , we focus here on bare DNA and show how to perform in silico single molecule manipulation of DNA .'], ['In rigid body dynamics simulations run with ODE , the state of a system consisting of rigid bodies is described by the positions of their centres of mass , a quaternion representation of their orientations , and their linear and angular velocities and respectively .', 'These velocities are collected in the column vector .', 'We use the superscript T to denote the transpose of a vector or a matrix everywhere in this article .', 'The vector then collects all linear and angular momenta , where is a block diagonal matrix whose elements are the mass matrices and inertia matrices of the N bodies ( with the identity matrix ) .', 'The Newtonian dynamics equation then reads where the generalized force is a vector collecting forces and torques applied to the system .', 'These forces and torques may be external , due for example to gravity or magnetic fields , or internal , as a consequence of the mechanical constraints between the rigid bodies that make up the system .', 'Most notably , in articulated systems , as is the case of polymers , rigid bodies are connected by mechanical joints .', 'A joint is a relationship that is enforced between two bodies so that they can have only certain positions and orientations relative to each other , and ODE provides different types of joints according to the kind of articulation that has to be implemented , e . g . ball-and-socket , hinge , slider or universal .', 'Mathematically a joint imposes some holonomic constraint between both connected bodies .', 'Such a constraint is an equation that reads where is the distance between both joint bearings , e . g . the center of the ball of one body and the the center of the socket of the other one .', 'The constrained distance is purely geometrical , depending only on the relative position and orientation of both jointed bodies .', 'The position and orientation of each of the bodies the articulated system is composed of depend on time .', 'Therefore the constraints of the articulated system can be derived with respect to time to get the kinematic constraints in the form where we introduce the jacobian matrix of constraints ( see subsection “Exact solution for when there are no collisions” for a detailed example ) .', 'This velocity-based description is used in ODE as in most game/physics engines .', 'So , mechanical joints exert reaction forces and torques on the joint bearings .', 'These internal mechanical constraints can be collected into a generalized constraint force which , by virtue of the principle of virtual work , reads where is a vector of Lagrange multipliers that precisely accounts for the reaction forces and torques coming from the joint bearings [16] .', 'The Newtonian dynamics equation therefore reads where and stand for the external and internal contributions to the generalized force respectively .', 'As the constraint force reads , the Newtonian dynamics equation becomes an equation for in the form: .', 'Solving this equation for should moreover satisfy the holonomic constraints at every timestep .', 'However the discretization used in the numerical calculation results in errors on so that is generally not equal to 0 . Then , in order to have at the next timestep , the kinematic constraint should be adapted accordingly .', 'Indeed according to the Euler semi-implicit integration scheme which is used in velocity-based algorithms .', 'Hence .', 'But then this implies that the kinematic constraint is not equal to zero at time , i . e . , so that the joint bearings will continue to move apart afterwards .', 'In order to keep both and close to zero at every timestep , ODE introduces an error reduction parameter in the kinematic equation [18] .', 'This parameter has to be adjusted to some optimal value between 0 ( no correction at all ) and 1 ( complete correction of in one timestep ) .', 'However setting is not recommended since , as said above , this would imply that the joint bearings will continue to move apart afterwards with maximal velocities .', 'ODE recommends values between and .', 'In addition to , ODE introduces a second ingredient to soften the rigid constraints by allowing the violation of the constraint equation by an amount proportional to the restoring force .', 'More explicitly , a “constraint force mixing” diagonal matrix is defined , such that ( implicit integration ) [18] .', 'This is equivalent to introducing a spring-damper system ( spring constant of and damping constant of ) with implicit integrator between the joint bearings; this can be understood as analogous to a bead-spring model .', 'Nevertheless there is a major difference between this effective spring and a regular spring: the term constrains the velocity whereas a regular spring constrains the acceleration .', 'As a result , no energy is stored in this effective spring , at odds with a regular spring which stores an averaged energy ( see below subsection “Preliminary tests of validity and performance of the global thermostat” along with the histograms of energy in Figure 1 ) .', 'In particular , ODE uses a powerful software called libccd [20] to detect collisions between two convex shapes .', 'Whenever overlapping is detected between two rigid bodies , ODE attaches a temporary joint between them called a “contact joint” .', 'Defining vector ( resp . ) that connects the center of mass of body 1 ( resp . 2 ) to the contact point and denoting the common normal to both bodies at the contact point ( directed from 2 to 1 ) , the kinematic constraint imposed by the contact joint would read in the perfect case when the holonomic constraint imposed by the contact joint reads exactly .', 'However , in practice is not equal to 0 because of discretization errors , hence the kinematic constraint imposed by the contact joint actually reads: ( 1 ) with ( 2 ) ( 3 ) The right hand side of Eq .', '( 1 ) deals with the already existing overlapping of the two bodies in contact at time when collision is first detected , or with their residual overlapping while the contact joint exists .', 'By inserting the constraint force into the equation of motion and taking the first-order discretisation of this equation , one can easily get the following expression to be solved for : ( 4 ) This equation is of the form .', 'Importantly , the addition of the term to each diagonal term of the matrix provides a symmetric positive definite matrix , thus greatly increasing the solution accuracy of Eq .', '( 4 ) .', 'From this equation , the vector of Lagrange multipliers , hence the constraint force , can be determined .', 'Then the motion solver ( semi-implicit Euler integrator ) gives the new positions and orientations of the articulated bodies at time .', 'It is advantageous to choose the Exponential Map parametrization [21] for the quaternion integration .', 'In general Eq .', '( 4 ) has to be solved numerically and ODE has two algorithms to do so , one based on the Successive-Over-Relaxation ( SOR ) method [22] and the other based on the Linear Complementary Problem ( LCP ) [23] .', 'LCP time complexity is of order and space complexity ( memory ) of order where is the number of constraint rows [18]; whereas SOR time complexity is of order where is the number of successive-over-relaxation and space complexity of order [18] .', 'Both algorithms have equivalent performances when .', 'But in general LCP is more accurate , although much more time consuming , than SOR .', 'We compared these two algorithms for a chain of length without noticing significant differences in the errors on ( error on the colocalization of joint bearings ) .', 'In order to save computational time , we preferentially run the SOR method with a value of for the relaxation factor and .', 'These values are different from the default values in ODE and work well for a linear chain of rigid bodies connected with ball-and-socket joints .', 'But in some cases , the SOR method does not converge and we then switch to the LCP method , which always converges .', 'However in the case when there are no collisions between the rigid bodies the articulated system is composed of , we were able to derive an exact solution for ( see next subsection ) .', 'Therefore we solve Eq .', '( 4 ) according to the following scheme: For the sake of clarity , let us first consider the example of four rigid cylinders of length each connected with ball-and-socket joints at the extremities with the first one anchored to some fixed point , taken as the origin of the coordinates .', 'The vector is the tangent to the cylinder .', 'The jacobian matrix associated with this system is tridiagonal when there are no collisions , in which case it reads: ( 5 ) For each cylinder , the antisymmetric matrix is associated with the cross product and has the property : ( 6 ) The transpose Jacobian matrix and mass matrix are given by: ( 7 ) Then we get: ( 8 ) and we deduce the final result for the symmetric matrix : ( 9 ) where we define the matrix .', 'We then write in the associated principal axis body frame as : ( 10 ) The vector collects the Lagrange multipliers associated to each joint respectively .', 'The equation gives us a system of coupled equations on .', 'Note that ODE solves in one time all the constraints of this articulated system .', 'This is not the case with the SHAKE algorithm where an unconstrained step is first performed , before correcting the positions and orientations iteratively to get the constraints satisfied eventually .', 'The term from Eq .', '( 4 ) is given by: ( 11 ) where we write and .', 'We can then write the constraint forces and torques ( 12 ) We can now generalise the previous example to the case of a linear chain of rigid cylinders connected with ball-and-socket joints with the first one anchored to the ground .', 'We denote the matrix with the following properties: ( 13 ) ( 14 ) ( 15 ) ( 16 ) Using the following decomposition for the matrix where is a block lower matrix with block identity matrix on the diagonal and where is the block diagonal matrix it is easy to show that for .', 'From these we get the following equations: ( 17 ) ( 18 ) ( 19 ) In order to solve the linear system of equations we define and solve the problem in an iterative way: ( 20 ) ( 21 ) and we get the final solution for by solving the problem : ( 22 ) ( 23 ) The method explained here is the exact solution of the problem where no collisions are present in the system .', 'With this exact resolution the simulation is faster than the SOR algorithm ( and ) with a gain of .', 'To model a DNA molecule , we build a linear chain of rigid cylinders of length each , corresponding to 10 base pairs ( ) , which amounts to the double helix pitch .', 'We connect the cylinders to each other by ball-and-socket joints .', 'The radius of the cylinders is set according to the salt buffer concentration of the experimental data we compare with .', 'Indeed , since the DNA molecule is highly negatively charged , DNA-DNA electrostatic repulsion affects the double helix response in single molecule experiments [9] , [24]–[27] .', 'This effect can be easily and implicitly included in simulations and theoretical models by introducing an effective DNA radius where is the crystallographic radius of the DNA double helix and accounts for the DNA-DNA electrostatic repulsion [28]–[31] .', 'It turns out that may be set equal to the Debye length of the salt buffer solution .', 'As with c the salt concentration given in and in , we set the effective radius to in mmol monovalent salt buffer for comparison with the reference experimental data of Mosconi et al [32] .', 'Alternatively we use an effective radius to fit the experimental data obtained in mmol monovalent salt buffer by Smith et al [1] .', 'We performed all our in silico single molecule experiments with a DNA molecule of contour length .', 'The corresponding number of DNA cylinders in the chain is therefore .', 'The DNA molecule is anchored , at one end , to a planar surface ( mimicking the microscope coverslip ) , and at the other end , to a rotatable bead ( mimicking the magnetic bead ) .', 'We set the bead radius to in order to prevent the DNA from looping around it .', 'At both ends of the DNA chain , the rigid cylinders are tangent to their attachment surface .', 'The final problem that remains to be addressed is how to obtain a correct definition of the bending and twisting behavior of DNA .', 'We have solved this problem by a special choice of the connecting joints and by introducing appropriate restoring torques reacting to the bending and twisting deformations .', 'This has been done based on the bending and twisting energies that are defined according to the usual expressions and respectively .', 'The rigidity constants and are related to the bending and twisting persistence lengths and respectively , through the following equations: ( 24 ) ( 25 ) ( 26 ) where is the Langevin function ( see supplementary Text S1 ) .', 'The bending angle and twisting angle are related to the standard Euler transformation ZXZ and are given by ( 27 ) ( 28 ) ( 29 ) where is the tangent vector of cylinder i , a vector normal to and .', 'These three vectors are the principal axis of cylinder i .', 'We finally get the following expression for the global restoring torque between two connected DNA segments ( see supplementary Text S1 for the complete derivation of this equation ) : ( 30 ) We recall that , for DNA , estimates of the bending persistence length give for salt buffer ( see Refs . [1] , [33]–[35] ) ; whereas estimates of the twisting persistence length give for salt buffer [9] , [32] .', 'According to the size of the unit cylinder we find and .', 'Although well adapted to mechanical simulations , ODE lacks coupling to a thermal bath .', 'As physics engines impose to deal with dynamics equations including inertial terms , in particular for computing constraint forces ( collected in ) , we need to turn to some implementation of stochastic isothermal molecular dynamics in order to thermalize the system: isothermal to simulate the system at constant temperature , stochastic to ensure ergodicity .', 'The corresponding algorithms are all related to Langevin dynamics and can be cast into local and global thermostats .', 'In local thermostats , such as standard Langevin dynamics , a correction force including both a frictional term and a stochastic term is exerted on each particle to drive the system to the canonical distribution at a prescribed temperature .', 'Global schemes of Langevin dynamics are designed to minimize the perturbation introduced by the thermostat on the Hamiltonian trajectory ( so called “disturbance” as defined originally in [36] ) , hence on the dynamical properties , such as autocorrelation functions , and related quantities , such as diffusion coefficients .', 'In these globally applied thermostats the stochastic term of the correction force acting on each particle is proportional to the momentum of that particle .', 'Two main global algorithms have been designed so far:', '( a ) Stochastic Velocity Rescaling methods , most notably the “global thermostat” introduced by Bussi and Parrinello [19] ,', '( b ) the Nosé-Hoover Langevin thermostat [37] .', 'Here we first show how to implement Langevin-Euler equation in the ODE software .', 'Moreover we show that the global thermostat introduced by Bussi and Parrinello is so remarkably adapted to this implementation that it improves quite significantly the sampling efficiency with respect to local Langevin dynamics ( by two orders of magnitude in typical situations ) , while preserving the time-dependent properties such as autocorrelation functions .', 'The sampling efficiency is defined as usual as the number of independent configurations generated during the time necessary to reach thermal equilibrium .', 'To begin with , we add to the “mechanical” forces an additional , thermal contribution containing a frictional term and a random force vector .', 'is the matrix of the coupling frequencies to the thermostat , the matrix of white noise amplitudes and a generalized vector of normalized and independent Wiener processes .', 'and are related through the fluctuation-dissipation theorem , which reads here ( 31 ) where with the temperature of the thermal bath and where the superscript denotes that the matrices and are chosen to be diagonal in the principal axis body frame ( where the matrix is diagonal by definition ) .', 'For simplicity , we choose to fix all the to a common frequency .', 'Note that is the relaxation time of the thermostat , i . e . the autocorrelation time of the kinetic energy ( see supplementary Figure S1 ) .', 'We then improve the sampling efficiency of this Langevin dynamics by extending the “global thermostat” algorithm designed by Bussi and Parinello in 2008 [10] to physics engines .', 'This algorithm allows faster yet correct sampling of the phase space in the canonical ensemble .', 'However , it is designed for the translational degrees of freedom only .', 'In order to apply it to an articulated rigid body system , we therefore have to extend it to the rotational degrees of freedom and adapt it to the ODE software .', 'To this aim , we replace the traditional Langevin-Euler correction force ( local thermostat ) by a corresponding global version , which reads ( 32 ) Eq .', '( 32 ) shows that is proportional to , so that the stochastic force and torque globally associated with the thermostat is in the same direction as .', 'Hence , a free particle , i . e . a particle not connected to any other particle , will move on a straight line between two collisions .', 'Note that , nevertheless , the particle will undergo Brownian motion along this straight trajectory .', 'The global version of the Langevin dynamics minimizes the disturbance induced by the thermostat on the Hamiltonian trajectory ( equal to according to its definition in Ref . [19] , but extended here to the rotational degrees of freedom ) , nevertheless retaining the same thermalization speed as usual Langevin dynamics ( see supplementary Figure S1 ) .', 'When used in the framework of a velocity-based algorithm such as ODE , the global thermostat presents a remarkable advantage .', 'This is because , in this case the global Langevin contribution is decoupled from the constraint forces , in the sense that it cancels out in the equations for .', 'More precisely , with our definition of ( see Eq .', '( 32 ) ) , the contribution to the term in Eq .', '( 4 ) is always zero .', 'In other words , not only minimizes the disturbance of the Hamiltonian trajectory , but also does not disturb the generalized constraint force .', 'Both effects cooperate to achieve a dramatic acceleration of the simulation sampling , that is , in the case of our model , approximately times faster than with the local thermostat ( see below “Preliminary tests of validity and performance of the global thermostat” and Figure 2 ) .', 'Importantly this acceleration is compatible with the correct computation of dynamical properties , such as autocorrelation functions .', 'In all DNA simulations presented in this article , we choose the length of the cylinders as the unit length , the mass of the cylinders as the unit mass and the unit of thermal agitation as the unit of energy , from which we deduce the unit of time .', 'The complete set of parameters of our simulations is given in Table 1 . We also choose to deal collisions with a restitution coefficient equal to without surface friction .', 'Hence , when two rigid bodies collide , the constraint force imposed by the contact joint that temporarily connects them is directed along their common normal at the contact point .', 'We finally choose an error reduction parameter and a constraint force mixing parameter .'], ['We start to validate our methodology by simulating a DNA molecule without any constraint applied on the bead ( neither stretching nor twisting ) .', 'To this aim , we first check the equipartition theorem .', 'When the system is at thermal equilibrium , its temperature is related to the kinetic energy through the equation ( 33 ) where is the number of non-redundant holonomic constraints .', 'This relation is standard since is just the number of degrees of freedom ( dof ) of the system .', 'Moreover the distribution of the kinetic energy of the system at thermal equilibrium follows a Boltzmann law and therefore reads: ( 34 ) with .', 'We checked this relation for a DNA molecule of length µm coupled to the two different Langevin thermostats , local and global respectively .', 'The resulting histograms are shown in Figure 1 , confirming that:', '( i ) the kinetic energy is correctly sampled at thermal equilibrium with both thermostats ,', '( ii ) there is indeed no energy stored in the joints , although these have been softened by effective springs ( see above the error reduction parameter in subsection “Introduction to physics engines” ) .', 'We then quantified the simulation sampling efficiency by means of the autocorrelation function of the end-to-end distance of a DNA molecule of length µm coupled to the two different Langevin thermostats , local and global respectively .', 'Here the average is performed over the time and denotes the lag-time of the autocorrelation function .', 'A demonstration of the performance of the global thermostat in terms of relaxation rapidity is given in Figure 2 . Fitting the exponential decrease of both relaxation curves shows that with use of the global thermostat we reach the saturation value at whereas this value is reached at in the case of the local thermostat , thus resulting in an acceleration factor of about for this system composed of articulated rigid bodies .', 'Note that the dimensionless lag-time is equal to ( see Table 1 ) , with the corresponding number of time steps .', 'Then a typical run using the global thermostat is of the order of tens of millions of time steps , whereas it is of the order of billions of time steps with usual Langevin dynamics .', 'A striking illustration of the sampling acceleration provided by the global thermostat is also given in supplementary Video S1 .', 'We also compute the tangent-tangent correlation function along the polymer with both local and global thermostats .', 'No significant deviations were found between both thermostats .', 'Results obtained with the global thermostat are plotted in supplementary Figure S2 along with the corresponding theoretical curves .', 'A simple calculation shows that the tangent-tangent correlation function decreases as , from which one can calculate the average bending .', 'With the DNA persistence length , this quantity amounts to , to be compared to the result from a fit of the simulation curves , giving .', 'The same comparison can be done for the twist angle ( with ) , for which the simulation average matches the theoretical value .', 'These comparisons show that our simulation results are in very good agreement with the analytical formulae , thus validating', '( i ) our implementation of the bending and twisting rigidities and', '( ii ) the correct sampling of the DNA conformation space by means of the global Langevin thermostat .', 'We then simulate reference force-extension curves , both theoretical and experimental .', 'We thus perform simulations at given stretching force along the z-axis ( normal to the DNA anchor surface ) , and without torsional constraints on the magnetic bead .', 'In order to fit the experimental data obtained by Smith et al [1] in monovalent salt buffer , we set here the DNA radius ( i . e . the radius of the unit cylinders ) to .', 'The resulting force-extension curve is given in Figure 3 where we plot ( red circles ) the dimensionless stretching force as a function of the dimensionless mean relative extension .', 'Here denotes the mean DNA extension , i . e . the mean distance between the bead and the anchor surface , at zero torsional constraint .', 'For comparison , we also plot ( black solid line ) the analytical Worm-Like-Chain ( WLC ) interpolation fitting curve proposed by Bouchiat et al [35] as well as the numerical solution of the WLC model ( black triangles ) obtained by Marko and Siggia with the same persistence length [34] .', 'The simulation reproduces pretty well the WLC behavior , thus validating our implementation of the DNA bending rigidity .', 'Note that , at low forces , the extension saturates at a value greater than zero because of the impenetrable ground and magnetic bead that both confine the DNA molecule .', 'This effect is more pronounced than in the experimental curve [1] because the ratio of the bead radius to the DNA length is higher in our simulations .', 'In Figure 3 , we also show the results obtained in the limit case , for which ( see Eqs .', '( 24–25 ) ) , and when there are no collisions .', 'In this case we expect to observe a Freely-Jointed-Chain response ( with segments of each ) .', 'The analytical force-extension relation for FJC is given by the well-known expression as a Langevin function and it is also reproduced in Figure 3 . Again , the simulation results are in very good agreement with the theoretical formula .', 'More interestingly , magnetic tweezers also allow the application of a torsional strain on a single DNA molecule at constant stretching force .', 'This torsional strain is equal to the number of turns of the magnetic bead around the z-axis due to the rotation of the magnets .', 'The number of turns of the bead is also equal to , the variation of the linking number of the DNA double helix with respect to the intrinsic twist of the DNA double helix with the pitch of the DNA .', 'And we define as usual the DNA relative overtwist as .', 'The method developed here can conveniently describe all physiological situations involving DNA positive supercoiling , which are of main importance for DNA transcription or replication .', 'One limitation of the modeling developed so far is that extensive modifications of the double helix structure are not accounted for , e . g . base pair opening that occurs when negative supercoiling is applied ( DNA denaturation ) , or the S-DNA transition observed at extremely high force , or the P-DNA transition under very positive torque .', 'Nevertheless this methodology gives an unparalleled opportunity to study more complex biological systems , such as protein-DNA complexes: in particular , we are currently addressing the modeling of magnetic tweezers response of chromatin fibres [7] , by associating our in silico DNA model with a rigid body representation of the histone octamer .', 'More recently , in vivo experiments also enable the measurement of the dynamics of chromosome loci in the cell nucleus [11] , [12] , [46] .']]
['Introduction', 'Materials and Methods', 'Results/Discussion']
['We develop a new powerful method to reproduce in silico single-molecule manipulation experiments .', 'We demonstrate that flexible polymers such as DNA can be simulated using rigid body dynamics thanks to an original implementation of Langevin dynamics in an open source library called Open Dynamics Engine .', 'We moreover implement a global thermostat which accelerates the simulation sampling by two orders of magnitude .', 'We reproduce force-extension as well as rotation-extension curves of reference experimental studies .', 'Finally , we extend the model to simulations where the control parameter is no longer the torsional strain but instead the torque , and predict the expected behavior for this case which is particularly challenging theoretically and experimentally .']
['Video game techniques are designed to simulate rigid body dynamics of macroscopic bodies , e . g . characters or vehicles , in a realistic manner .', 'However they are not able to deal with temperature effects , hence they are not able to deal with molecules .', 'In order to extend these powerful techniques to molecular modeling , we implement here Langevin Dynamics in an open source library called Open Dynamics Engine .', 'Moreover we add a “global thermostat” to this Langevin Dynamics , which accelerates the simulation sampling by two orders of magnitude .', 'With these radically new simulation techniques , we prove that we can accurately reproduce single-molecule manipulation experiments in silico , in particular force-extension as well as rotation-extension curves of reference experimental studies .', 'The method developed here represents an unparalleled tool for the study of more complex single molecule manipulation experiments , notably when DNA interacts with proteins .', 'Furthermore the simulation technique that we propose here has all the functionalities required to tackle the nuclear organization of chromosomes at every length scale , from DNA to whole nuclei .']
['physics', 'computer science', 'computer animation', 'biophysic al simulations', 'computing methods', 'biology', 'computational biology', 'biophysics simulations', 'biophysics']