## Abstract

In the present article, we summarize our recent studies of DNA dynamics using the generalized immersed boundary method. Our analysis of the effects of electrostatic repulsion on the dynamics of DNA supercoiling revealed that, after perturbation, a pre-twisted DNA collapses into a compact supercoiled configuration that is sensitive to the initial excess link and ionic strength of the solvent. A stochastic extension of the generalized immersed boundary method shows that DNA in solution subjected to a constant electric field is compressed into a configuration with smaller radius of gyration and smaller ellipticity ratio than those expected for such a molecule in a thermodynamic equilibrium.

- computational model
- DNA dynamics
- immersed boundary method
- polymer
- simulation
- supercoiling

## Introduction

Modelling of many biological systems requires one to take into account the interaction of elastic structures with the surrounding fluid. Among traditional methods for simulation of fluid–structure interaction problems, the IB (immersed boundary) method [1], created to study the hydrodynamics of heart valves, stands out as an efficient numerical procedure that has been used for a variety of such systems. The classical IB method has been extended to the GIB (generalized immersed boundary) method by incorporating the Kirchhoff rod theory to treat long and thin filamentous structures that include both positional and rotational degrees of freedom [2]. The GIB method allows the structure to apply both force and torque to the fluid, and the structure moves at the local fluid velocity and rotates at the local fluid angular velocity.

We have developed an improved GIB method for modelling the dynamics of charged rod-like structures in a viscous fluid [3] and formulated an extension of the GIB method that accounts for thermal fluctuations, with the goal of modelling the dynamics of DNA immersed in a solvent with counterions. DNA is a double-stranded molecule composed of two polynucleotide strands that are bound together by hydrogen bonds between complementary nucleotide bases. The long-range electrostatic repulsion is caused by charged phosphate residues on the DNA backbone. Two of these residues are located every 0.34 nm, each carrying one electronic charge. The main influence of these charges is that the chance of self-contact of DNA, i.e. contact of a point on the surface of DNA with another point on the surface, is greatly reduced compared with electrostatically neutral polymer, which is particularly important during DNA supercoiling.

In our model, the DNA is represented by a homogeneous isotropic, intrinsically straight Kirchhoff elastic rod with charge density distributed uniformly along the axis. In both methods, the self-contact of DNA is treated by including: (i) a long-range electrostatic repulsion, and (ii) a hard-core potential that does not allow two points on the axis of the rod to approach closer than the specified distance (rod diameter) and prevents the passage of the DNA through itself.

## Deterministic model

In the study of DNA rings of length smaller than a few multiples of the persistence length, elastic energy dominates the trajectory of DNA and thermal effects can be neglected. We represent DNA as a homogeneous, naturally straight, electrically charged elastic rod immersed in an incompressible fluid governed by the Navier–Stokes equations. The governing equations are the same as in the model introduced by Lim et al. [3] with electrostatic and steric repulsion forces included in addition to the elastic force. The electrostatic interaction between DNA charges in solution is screened by counterions and hence the strength of any electrostatic effect decreases with increasing ionic strength of the solvent, i.e. the concentration of the dissolved salt. We employ the theory of counterion condensation, due to Manning [4], which concludes that the DNA is surrounded by a concentrated layer of counterions which neutralizes approximately 76% of the DNA charge. The electrostatic energy is then described by the Debye–Hückel term:
(1)
where *q _{i}* is the net charge of the

*i*th base pair, equal to 0.48e,

*d*is the distance between charges

_{ij}*i*and

*j*. The constant ϵ

_{0}is the permittivity of free space, and ϵ

_{w}is the dielectric constant of water at 300 K. The Debye screening parameter, κ, is given by where

*Cs*is the molar salt concentration in mol/litre.

The numerical scheme we employ to solve the system of governing equations uses the standard finite difference method. In particular, we solve the discretized Navier–Stokes equations using fast Fourier transform facilitated by the periodic fluid domain (see [2] for more details). All parameters are the same as in [3].

## Dynamics of supercoiling of DNA minicircles

As an example of results that can be obtained using the GIB method, in the present article, we report on simulations of the supercoiling of a 600 bp DNA minicircle that has initial configuration with a circular axial curve and a uniform initial twist density. Thus the initial configuration is an equilibrium solution of the minicircle, albeit unstable. We perturbed this configuration by varying the twist density along the DNA while leaving the axial curve intact and observed the dynamical behaviour of the minicircle. We chose the excess link Δ*Lk*, i.e. the number of turns that the ends of the DNA are twisted before closure, to be an integer between 1 and 6 and kept *Cs* fixed for the duration of the simulation at 0.01, 0.05 or 0.1 M. Regardless of Δ*Lk* and *Cs*, and the initial perturbation, the DNA trajectories we computed progressed through four characteristic phases: (i) rapid equilibration of the twist density accompanied by a small perturbation of the axial curve, (ii) a collapse of the circle to the point of first ‘contact’ at which DNA comes to close proximity and the motion abruptly slows down, (iii) progression through one or several intermediate configurations with non-increasing number of terminal loops and the creation of additional contacts, and (iv) a slow convergence to a locally stable equilibrium accompanied by writhe relaxation.

The limiting configurations for a 600 bp minicircle can be seen in Figure 1(A). The shapes of limiting configurations at high salt resemble those found for elastic rods without electrostatic repulsion [5]. DNA at low salt has the tendency to limit the total extent in close contact, which results in a preference for toroidal supercoiling. The time of transition from the perturbed circle to a stable equilibrium configuration for the cases shown in Figure 1 varies greatly with Δ*Lk* and *Cs*. The extent of deformation of the DNA ring can be judged from its writhe *Wr*, defined as the average, over all planar projections, of the sum of signed crossings seen in the projection. Alternatively, one can compute writhe as the difference of Δ*Lk* and the total twist in the rod. The graph of *Wr* against time for the configurations in Figure 1(A) is shown in Figure 1(B). Writhe increases monotonically with both increasing Δ*Lk* and increasing *Cs*.

## Stochastic model

Accurate modelling and simulation of long DNA in a fluid requires one to account for intrinsic thermal fluctuations of DNA and Brownian collision forces exerted on DNA by the surrounding solvent. Traditionally, such simulations have used a Brownian dynamics approach [6–11] in which the DNA is represented as a chain of beads endowed with energy accounting for the bending, stretching, electrostatic and/or excluded volume interactions of the DNA. The motion of the beads is computed using a Brownian dynamics algorithm. We have recently developed an alternative method for studying the dynamics of DNA in solution by combining the GIB method [3] with the SIB (stochastic immersed boundary) method of Atzberger et al. [12]. The benefits of this combined SGIB (stochastic generalized immersed boundary) method over the Brownian dynamics simulations is that SGIB accounts for the torsional drag in DNA–fluid interactions, which is essential for accurate study of supercoiling dynamics of DNA, and also that the computing time of SGIB scales with the spatial resolution of the computational domain, not the size of DNA.

In the SGIB method, the non-linear advection term in Navier–Stokes equation is neglected because the Reynolds number, the ratio of inertial force to viscous force, is very small in DNA dynamics. Thus the fluid is treated by the time-dependent Stokes equation for incompressible fluid with velocity **u**(**x**,*t*), pressure *p*(**x**,*t*), density ρ and viscosity μ:
(2)

The force density **f**(**x**,*t*) acting on the fluid consists of two components: the force resulting from mechanical imbalance of the immersed DNA structure and the thermal force, which is the source of both DNA and solvent fluctuations. The mechanical force is computed using balance equations for the forces and moments of DNA structure, represented as an elastic rod, and it is applied to the fluid via generalized delta functions, in the same way as in the GIB method (see [3]). The thermal force is uniform across the entire domain and is represented in the Fourier space by Gaussian white noise derived from standard Brownian motion *B*** _{k}**(

*t*), with

*D*

**, a mode-dependent strength of the forcing (see [12]): (3)**

_{k}The benefit of the Stokes approximation together with Fourier representation can be seen when one transforms the finite-difference discretization of the Stokes equation to a Fourier space: in the standard notation of stochastic differential equations, one obtains
(4)
where α** _{k}** is a mode-dependent decay coefficient and P

_{k}is a projector that guarantees that the solution of eqn (4) is a divergence-free field in the real space. Assuming that the force field is independent of time, this equation can be integrated using Itô calculus to yield a discrete update formula for

**û**

_{k}which is the basis for a numerical algorithm for SGIB simulation. (5)

The last term in eqn (5) represents the projection of a random field consisting of numbers with mean 0 and variance σ_{k}^{2} which is derived from *D*** _{k}** and α

**. (For detailed derivation of equations and the definitions of various constants, see Atzberger et al. [12].) Note that, in the absence of forces or when forces are constant, the formula (5) is exact for any value of the time step Δ**

_{k}*t*. Thus the magnitude of the time step required to maintain accuracy is determined by the dynamics of the immersed structure, not dynamics of the fluid itself. It was shown that the SIB method agrees with Brownian dynamics methods in the case of small particles separated by large distances (far-field hydrodynamics) [13].

The simulation algorithm works as follows. First, the balance equations of the DNA are evaluated, and the resultant forces and moments are computed at the DNA discretization points. Then the forces and moments are converted into forces and moments acting on fluid discretization points using a compact support function that is characteristic of the IB method. The fluid velocity field is incremented according to eqn (5), and then the velocity of the DNA at its discretization points is computed. As the final step, the DNA is propagated along the computed velocity. These steps are repeated until the desired simulation time is reached.

## Dynamics of DNA in electric field

The simulations we report in the present article are inspired by the experiments of Tang et al. [14] in which individual T4 DNA molecules were subjected to moderate electric field and observed to compress into isotropic polymer globules. The compression was observed to occur at lower field intensities when the DNA charge was less neutralized by counterions, and the time scale of compression was found to be proportional to the field intensity. It was also found that alternating current leads to moderate compression at low frequencies and no compression at high frequencies. When the field was turned off, the molecule relaxed back to its natural random coiled state determined by thermodynamic equilibrium.

Preliminary simulations of a 6000 bp linear DNA molecule in a fluid in the presence of counterions and electric field show that, in accordance with experiments, the simulated molecule has a tendency to undergo compression. As in Tang et al. [14], we monitor the compression and the molecule shape by computing the radius of gyration *R*_{g} and the average ratio of between the major and minor axes of the molecule, *R*_{M}/*R*_{m}, when approximated by an ellipsoid. Figure 2 shows the results of one such simulation for a DNA molecule that was initially in a shape consisting of four turns of a helix with radius 81 nm and pitch 25 nm. The uniform electric field is applied and generates a force acting on each charged group of the molecule. Figure 2(B) shows the changes in the DNA configuration induced by the combined effect of the electric field and the fluid. Since charges are constrained in position, this force has a constant linear density and is parallel to the direction of the field. As a result, the DNA accelerates in the direction of the field; this acceleration is opposed by the drag force of the fluid. The fluid flows around the molecule and also passes though the molecule with somewhat decreased velocity. This difference in velocities results in a moment that is exerted on the outermost parts of the molecule. This moment results in a rotation of the molecule about a circular axial curve that is perpendicular to the direction of the flow; the overall motion resembles a toroidal vortex ring. The ring-shaped configuration collapses further to form a tight circular bundle due to hydrodynamics. The whole process is also reflected in Figure 2(A) in the graphs of the radius of gyration, which decreases steadily to approximately 50 nm, and in the graph of *R*_{M}/*R*_{m} which oscillations describe the rotations of the molecule about the centre circle of the torus.

## Conclusions

Combination of the elastic rod model for DNA with electrostatic repulsion and the IB method provides a novel tool for the study of the dynamics of DNA that complements existing techniques of Brownian dynamics simulations or lattice Boltzmann methods. We find that, in the absence of thermal fluctuations, supercoiling dynamics is overdamped and that the time of transition from the circle to a supercoiled configuration can vary greatly depending on excess link and ionic strength and that the transition slows down at various times as the configuration approaches unstable equilibria, with characteristic plateaus in the writhe. The dwell time near such equilibria is of the order of 1–10 μs, which is enough time for a protein to bind and connect DNA segments that are in close proximity. The transition paths are sensitive to ionic strength; small changes in ionic strength may cause divergence of the trajectories in the neighbourhood of transient states.

The technique developed by Atzberger et al. [12] allows us to extend the IB studies of DNA dynamics to include thermal fluctuations. In preliminary simulations, we found that the compression of DNA in electric field is caused primarily by the interaction of the DNA with a fluid through which the molecule is moving. The difference in fluid speed around the molecule and through the molecule gives rise to moments that fold the molecule into a compact shape. Electrostatic forcing is essential for this process because it causes a uniform force density to be applied to every segment of the molecule, and hence it does not, by itself, give rise to a force differential that would tend to extend the molecule. Such force differential is observed, for example, when the molecule is moved through the fluid by attaching it at one end to an object, such as a magnetic or glass bead. The attached end of the molecule then experiences a force in one direction while the remaining part of the molecule experiences a drag force in the opposite direction generated by the fluid. We plan to investigate this effect in detail in further studies and describe the effect of changes in ionic strength and the electric field intensity on the timescale and amount of compression.

The model of DNA elasticity used in this research is based on an elastic model for DNA structure in which the bases remain bonded in base pairs during all deformations of the structure. The method as presented is therefore not suitable for DNA that is strongly negative supercoiled, sharply kinked or overstretched, since those deformations are known to result in DNA denaturation or disruption of base-pairing. The elastic model could be potentially generalized to account for denaturation and then paired with the SGIB method. The SGIB method described in the present study differs significantly from the method presented by Zhang et al. [15] (also named the ‘immersed boundary method’): our DNA is fine-grained (the separation between DNA discretization points is 4 bp, i.e. 1.36 nm), our energy includes bending, twisting and electrostatic energies, and we update the position of the DNA by using explicit integral of the equation of motion (5) instead of Brownian dynamics algorithm based on a mobility matrix.

## Funding

Y.K. was supported by the National Research Foundation of Korea [grant number 2010-0006165]. D.S. acknowledges support by the Human Frontiers in Science Project [grant number RGP0051/2009].

## Footnotes

Topological Aspects of DNA Function and Protein Folding: An Independent Meeting held at the Isaac Newton Institute for Mathematical Sciences, Cambridge, U.K., 3–7 September 2012, as part of the Isaac Newton Institute Programme Topological Dynamics in the Physical and Biological Sciences (16 July–21 December 2012). Organized and Edited by Andrew Bates (University of Liverpool, U.K.), Dorothy Buck (Imperial College London, U.K.), Sarah Harris (University of Leeds, U.K.), Andrzej Stasiak (University of Lausanne, Switzerland) and De Witt Sumners (Florida State University, U.S.A.).

**Abbreviations:**
GIB, generalized immersed boundary;
IB, immersed boundary;
SGIB, stochastic generalized immersed boundary;
SIB, stochastic immersed boundary

- © The Authors Journal compilation © 2013 Biochemical Society