## Abstract

We investigate the effects of electrostatic and steric repulsion on the dynamics of a pre-twisted charged elastic rod immersed in a viscous incompressible fluid. Equations of motion of the rod include the fluid–structure interaction, rod elasticity and a combination of two interactions that prevent self-contact, namely the electrostatic interaction and hard-core repulsion. The governing equations are solved using the generalized immersed-boundary method. We find that after perturbation, a pre-twisted minicircle collapses into a compact supercoiled configuration. The collapse proceeds along a complex trajectory that may pass near several unstable equilibrium configurations, before it settles in a locally stable equilibrium. The dwell time near an unstable equilibrium can be up to several microseconds. Both the final configuration and the transition path are sensitive to the initial excess link, ionic strength of the solvent and the initial perturbation.

## 1. Introduction

Modelling of many biological systems requires one to take into account the interaction of elastic structures with fluid. The immersed-boundary (IB) method (Peskin 2002), created to study the dynamics of heart valves (Peskin 1972), is an efficient mathematical procedure that has been used in the numerical simulation of a variety of such systems. The mathematical formulation of the IB method employs Eulerian independent variables for the description of the fluid and Lagrangian independent variables for the elastic solid. The interaction between the solid and the fluid is mediated by forces applied by the structure on the fluid and vice versa. Dirac delta functions are used to convert forces from one set of independent variables to the others. The numerical scheme for the IB method employs two meshes—a Cartesian mesh for Eulerian variables and a curvilinear mesh for Lagrangian variables, which moves freely through the Cartesian mesh. Smoothed approximations of Dirac delta functions are involved in the interpolation between the two meshes. A distinguishing feature of the IB method is that the velocity of the structure is identical to the velocity of the surrounding fluid owing to no-slip conditions, i.e. the structure is carried by the fluid (Peskin & McQueen 1996; Peskin 2002).

The classical IB method is applicable to the dynamics of two- or three-dimensional structures in a three-dimensional fluid, but cannot be directly applied to a one-dimensional structure (a space curve) because of the presence of a singularity due to delta functions. The generalized IB method (Lim *et al.* 2008) resolves this issue by using a smoothed delta function with a compact support, which has the effect of assigning a cross-sectional radius to the space curve. In addition, the generalized IB method incorporates the special Cosserat rod theory and is distinguished by two new features: (i) the local angular velocity of the fluid at the IB is evaluated, along with the local linear velocity of the fluid, and (ii) the IB applies torque as well as force to the surrounding fluid. It was observed that the generalized IB method prevents the passage of a discretely modelled, infinitesimally thin elastic rod through itself, even without the use of explicit repulsive potentials (Lim *et al.* 2008). This effect appears to be the result of the no-slip condition at the fluid–structure interface, where the interpolated velocity field is continuous, as is also true in the classical IB method. Nonetheless, the prevention of self-penetration did not prove to be reliable in the numerical implementation, since it varied with the grid-spacing and the time step.

Here, we present an improved IB method for modelling the dynamics of elastic rod-like structures in a fluid. In this method, the self-contact is treated naturally by including two types of interactions: (i) 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 (ii) a long-range electrostatic repulsion. Inclusion of these interactions prevents the passage of the rod through itself and enables one to study the effects of the rod diameter and the magnitude of the repulsive force on the dynamics.

The inclusion of electrostatic force was motivated by our application of the model to the dynamics of short DNA molecules in a solvent. DNA is a double-stranded molecule composed of two polynucleotide strands that are bound together by hydrogen bonds between complementary nucleotide bases. In normal conditions, the strands wind around the DNA axis as two right-handed helices with diameters of 2 nm and pitches of 3.57 nm. During biological processes, such as transcription or replication, various mechanical forces act on the DNA and cause its underwinding or overwinding, i.e. an increase or decrease in the helical pitch. The molecule responds by releasing some of the stored twisting energy through a deformation that results in the so-called supercoiled configurations (Bates & Maxwell 1993).

Aside from the bending and torsional elasticity, the dynamics of the molecule is influenced by long-range electrostatic interactions between 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 the DNA, i.e. contact of a point on the surface of the DNA with another point on the surface, is greatly reduced as the electrostatic repulsion provides a cushion between any two closely approaching DNA segments. Electrostatic interaction is particularly important for supercoiled DNA dynamics because, upon supercoiling, many parts of the double helix come into close contact. A secondary effect of electrostatics is the so-called electrostatic stiffening of the molecule, i.e. increased resistance to bending (Odijk 1977; Skolnick & Fixman 1977; Baumann *et al.* 1997).

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 dissolved salt. There are several theories describing the extent of such screening. The most commonly employed is the theory of counterion condensation, due to Manning (1978), which concludes that the DNA is surrounded by a concentrated layer of counterions that neutralizes about 76 per cent of the DNA charge. The electrostatic energy is then described by a Debye–Hückel term as a function of the modified charge, molar salt concentration and the distance between the charges.

Equilibrium configurations of supercoiled circular DNA have been studied extensively using elastic-rod models (e.g. Le Bret 1979, 1984; Benham 1989; Yang *et al.* 1993; Tobias *et al.* 1994; Coleman *et al.* 1995, 2000; Westcott *et al.* 1995; Dichmann *et al.* 1996 and also references in Swigon 2009). Schlick *et al.* (1994) incorporated both elastic potential and electrostatic forces in a computational model of DNA that represented the molecule by cubic B-splines. They found that at low salt, the electrostatic term dominates over the bending and twisting terms and the DNA assumes more open configurations, while at high salt, the DNA structure is dominated by the bending term and takes upon compact and bent interwound configurations. Westcott *et al.* (1997) also studied the mechanical equilibria of charged DNA by solving the equations of mechanical equilibrium for elastic DNA at fixed salt concentration. They showed that the extensible DNA yields the same equilibrium configurations as in the inextensible DNA. Recently, Biton *et al.* (2007) and Biton & Coleman (2010) have analysed the influence of electrostatic repulsion on configurations of both open and closed intrinsically curved DNA segments. The influence of electrostatics on DNA configurations in single-molecule manipulation experiments was recently analysed by Clauvelin *et al.* (2009) and Brutzer *et al.* (2010). The departure of a pre-twisted circular elastic rod from the circular equilibrium configuration has been studied using dynamical rod theories that included drag forces to represent the interaction of the rod with fluid (Klapper 1996; Goriely & Tabor 1997; Goyal *et al.* 2005, 2008).

We here use the generalized IB method to study the dynamics of supercoiling of a pre-twisted charged elastic rod and its dependence on the ionic strength of the surrounding solvent, mimicking the supercoiling of DNA minicircles. We model the solvent as a viscous incompressible fluid, assume a no-slip interaction between the solvent and the rod and treat the rod as a homogeneous, transversely isotropic, intrinsically straight and uniformly charged rod obeying special Cosserat theory of rods (Antman 1995), with electrostatic forces governed by the Debye–Hückel theory. The rod is initially circular with a uniform twist density and an imposed excess link corresponding to an integer number of turns of one end with respect to the other before closure. This configuration is perturbed in its twist density (but not the axial curve) and then left to evolve in accord with the governing equations.

## 2. Equations of motion

We consider a closed bent and twisted rod immersed in an incompressible fluid governed by the Navier–Stokes equations. The mathematical model discussed here is based on the model introduced by Lim *et al.* (2008) with electrostatic and steric repulsion forces included in addition to the elastic force. The immersed rod is represented by a three-dimensional closed space curve with an associated orthonormal triad at each point of the curve.

Eulerian coordinates are used in the description of the motion of the fluid, while Lagrangian variables are used to describe the immersed rod and its properties. The interaction between these two types of variables is facilitated by means of integral transformations that involve a smoothed version of the three-dimensional Dirac delta function (Peskin & McQueen 1996).

The coupled system of equations of the rod and the fluid is as follows:2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11and
2.12where **f**^{c}(*s*,*t*) and **f**^{e}_{i}(*t*) are the steric repulsion and electrostatic forces defined below.

Equations (2.1) and (2.2) are the Navier–Stokes equations for the velocity **u**(**x**,*t*) and pressure *p*(**x**,*t*) of an incompressible fluid; here, **x**=(*x*_{1},*x*_{2},*x*_{3}) describes a position in space and *t* is the time. The motion of the fluid is subject to the body force **f**^{b}(**x**,*t*), which here represents the force per unit volume applied to the fluid by the immersed rod. The constant parameters *ρ* and *μ* are the fluid density and the fluid viscosity, respectively.

The equilibrium equations (2.3)–(2.8) are employed to relate the force and torque of the immersed rod to the external forces applied on it, where the rod is defined by giving its axial curve **X**(*s*,*t*) and an orthonormal triad embedded in each cross section (**D**^{1}(*s*,*t*),**D**^{2}(*s*,*t*),**D**^{3}(*s*,*t*)). The rod is assumed to have no dynamics of its own other than the one that follows the motion of the fluid. This approximation becomes exact in the limit of a thin elastic rod with density equal to the fluid density and no-slip conditions imposed on the fluid velocity at the rod surface. All variables in equations (2.3)–(2.8) are functions of the linear material coordinate *s* of the rod (not necessarily the arclength) and the time *t*. **F**(*s*,*t*) and **N**(*s*,*t*) are the force and moment (couple) transmitted across the section of the rod at *s*. The expressions −**f**(*s*,*t*) and −**n**(*s*,*t*) are the density of the external force and the torque applied by the rod on the fluid.

The internal force and moment on the perpendicular cross section, **F** and **N**, and also the applied force density **f** and the torque density **n** may be expanded in the basis (**D**^{1},**D**^{2},**D**^{3}). Equations (2.7) and (2.8) are the constitutive relations of the special Cosserat theory of rods (Antman 1995). Here, *a*_{1} and *a*_{2} are the bending moduli of the rod about **D**^{1} and **D**^{2}, respectively, and *a*_{3} is the twisting modulus of the rod. The parameters *b*_{1} and *b*_{2} are the shear-force constants and *b*_{3} is the stretch-force constant. The elastic energy of the rod is given by2.13For simplicity, we assume that *a*_{1}=*a*_{2}≡*a* and *b*_{1}=*b*_{2}≡*b*, which corresponds to the case of a rod with transversely isotropic material properties.

Equations (2.9)–(2.11) describe the interactions between the fluid and the rod. These interaction equations connect the Lagrangian and Eulerian variables via a three-dimensional smoothed Dirac delta function *δ*_{c}(**x**)=*δ*_{c}(*x*_{1})*δ*_{c}(*x*_{2})*δ*_{c}(*x*_{3}), which acts as a kernel of the integral transformations that appear in the interaction equations. The particular choice of *δ*_{c}(**x**) that we make in this work is the following:2.14where **x**=(*x*_{1},*x*_{2},*x*_{3}) and the function *ϕ* is given by(figure 1). Note that *δ*_{c}(**x**−**X**) is a continuous function of **x** with continuous first derivatives and with support equal to a cube of edge 4*c* centred on **X**. Whenever *c* is an integer multiple of *h*, the function *δ*_{c}(**x**−**X**) satisfies two identities that are of particular importance in this work. Note in particular that these identities hold for all **X**,2.15and2.16where **j** is any vector with integer components, and *h* is the meshwidth of the fluid grid. As mentioned above, these identities hold only if *c*/*h* is a positive integer, and we shall choose *h* so that this is the case. The above identities ensure that force and torque generated by the rod are correctly applied to the fluid by our numerical scheme.

In equation (2.9), the first term describes how the force of the rod acts on the fluid and the second term describes how the torque of the rod acts on the fluid (see Lim *et al.* 2008 for details).

The repulsive force **f**^{c} is included to prevent penetration of the rod surface. It is assumed to be a Hookian force proportional to the amount of penetration, i.e.2.17where2.18for any two material points *s* and *s*′, such that , and (**X**(*s*,*t*)−**X**(*s*′,*t*))·**T**(*s*,*t*)=(**X**(*s*,*t*)−**X**(*s*′,*t*))·**T**(*s*′,*t*)=0. (In practice, for any *s* there are only finitely many *s*′ obeying such conditions, whence the sum in equation (2.17).) The constant *D*=20 Å is the diameter of the DNA, and the constant *g* is chosen sufficiently large so that no significant penetration occurs.

Electrostatics is included in the external force **f**^{b} as a screened Coulombic force , defined by2.19where2.20where **X**_{i} and *q*_{i} are the position and net charge of the *i*th base pair, respectively, *d*_{ij} is the distance between net charges *q*_{i} and *q*_{j} and *n* is the number of base pairs. The constant *ϵ*_{0} is the permittivity of free space and *ϵ*_{w} is the dielectric constant of water at 300 K. For simplicity, we represent each base pair by one charge located on the DNA axis. The condensed monovalent cations around the surface of the DNA reduce the charge of each phosphate group to 0.24e, where e is the elementary charge of an electron (Manning 1978; Williams & Maher 2000). Since there are two such phosphate groups per base pair, *q*_{i} is 0.48e. The Debye screening parameter, *κ*, is given by2.21where *C*_{s} is the molar salt concentration in moles per litre. The screened Coulombic force can be obtained as the first variation of the Debye–Hückel electrostatic potential (Westcott *et al.* 1997),2.22

Equation (2.10) represents the condition that the rod moves at the same velocity as the fluid, after averaging in a manner determined by the smoothed Dirac delta function. Equation (2.11) represents a similar relation between the angular velocity **W**(*s*,*t*) of the rod and a local average of the angular velocity of the fluid, 1/2(∇×**u**). Again, the smoothed delta function is used to determine the appropriate weighted average of the local fluid velocity. The change in orientation of the triad at each point of the rod is given by equation (2.12).

Any deformation of a closed DNA is subject to the constraint of a fixed linking number *Lk*, defined as one half of the number of signed crossings of a DNA strand and the DNA axis in a planar projection (Calladine *et al.* 2004). A related quantity is the excess link Δ*Lk*, which is a topological invariant, just like *Lk*, and is equal to half of the number of signed crossings of the axial curve and the curve that is formed by the endpoints of the vector **D**^{1} that is embedded in the cross section, and is constant in a stress-free configuration. If the axial curve of closed DNA is a simple closed planar curve, then it follows from the theorem of White (1969, 1989) that2.23and hence Δ*Lk* equals the total number of turns of one end of the rod before it is sealed with the other end to produce a ring.

The numerical scheme we employ to solve the system (2.1)–(2.12) and (2.17)–(2.20) is based on the generalized IB method. In particular, since the discretized Navier–Stokes equations are linear in the variables of time level (*n*+1) with constant coefficients that depend on the variables at time level *n*, Fourier transformation decouples the discretized system into *N*^{3} separate 4×4 systems of linear equations that are easily solved (for details, see Lim *et al.* 2008). The parameter values used in the computations are listed in table 1. The elastic and electrostatic constants are based on the physical properties of the DNA (Westcott *et al.* 1997). In particular, the ratio of twisting to bending moduli, *a*_{3}/*a* is here chosen to be 1.4 in accord with the experimental estimates obtained by measurements of the free energy of supercoiling for small DNA rings (Horowitz & Wang 1984) and single DNA manipulation experiments (Strick *et al.* 1996). Each trajectory was continued until the total kinetic energy of the system dropped below 10^{−12}*kT* to make sure that the system converged to a locally stable equilibrium. The computation of one dynamical trajectory takes 2–3 days of CPU time on a common Linux workstation. The order of accuracy of the numerical method was verified by performing calculations on three different grid sizes of 64^{3}, 128^{3} and 256^{3}, comparing the discrete *L*^{2}-norm of the difference in fluid velocities and computing the convergence ratios. For example, for the trajectory of the system with (Δ*Lk*,*C*_{s})=(4,0.1 M) (see below), the convergence ratios for the velocity **u**=(*u*,*v*,*w*) at *t*=0.4 μs were (||*u*_{64}−*u*_{128}||_{2})/(||*u*_{128}−*u*_{256}||_{2})=2.40, (||*v*_{64}−*v*_{128}||_{2})/(||*v*_{128}−*v*_{256}||_{2})=2.56 and (||*w*_{64}−*w*_{128}||_{2})/(||*w*_{128}−*w*_{256}||_{2})=1.98, which implies that the scheme is first-order accurate.

## 3. Results and discussion

In choosing the initial configuration, we start with a configuration that has a circular axial curve of radius *r*_{0} and a uniform twist density Δ*Lk*/(2*πr*_{0}), and hence is an equilibrium solution of equations (2.3)–(2.8) in the absence of external force and moment, i.e. with **f**=0 and **n**=0. We perturb this configuration by varying the twist along the DNA while leaving the axial curve intact. In terms of cylindrical coordinates (*r*,*θ*,*z*) with unit vectors (**r**(*θ*),**θ**(*θ*),**z**), the initial configuration becomes3.1
3.2
3.3
3.4and
3.5whereand3.6Note that integer-valued Δ*Lk* implies that the triad will be continuous at *s*=0, which is the same point as *s*=2*πr*_{0}. The vector **E**(*s*) is a useful reference vector orthogonal to **D**^{3} within the plane spanned by **θ** and **z**. The charges are located on the axis at the discretization points. As each material point of the discretized immersed rod corresponds to a base pair of the circular DNA, the spacing between them is 0.34 nm.

The use of a circular initial configuration is standard in studies of DNA supercoiling, and corresponds to the situation in which a sudden change in the stress-free value of twist density pushes a stable circular minicircle above the instability threshold. This can be accomplished experimentally, for example, by adding intercalating agents to the solution.

We study the dynamics of initially circular DNA with excess link Δ*Lk* fixed at a value between 1 and 6 and the molar salt concentration *C*_{s} fixed at 0.01, 0.05 or 0.1 M. Regardless of Δ*Lk*, *C*_{s} and the initial perturbation, the majority of DNA trajectories we computed progressed through four characteristic phases described in figure 2.

### (a) Dependence of limiting configurations on Δ*Lk* and *C*_{s}

We first explore the dependence of trajectories and limiting configurations on Δ*Lk* and *C*_{s} while the initial perturbation is held fixed. For Δ*Lk*=1, the circular configuration is stable even in the absence of electrostatics (Coleman & Swigon 2000). A familiar result in the theory of inextensible, unshearable elastic rods states that the circular equilibrium configuration is unstable if Δ*Lk* exceeds a critical value (Michell 1889; Goriely 2006). The equivalent critical value for a charged rod with *a*_{3}/*a*=0.67 has been estimated by Schlick *et al.* (1994) as a function of the ionic strength: Δ*Lk**(0.1 M)=2.27, Δ*Lk**(0.05 M)=2.38, Δ*Lk**(0.01 M)=2.81, Δ*Lk**(0.005 M)=3.25. We find that for our choice of *a*_{3}/*a*=1.4, the circular configuration with Δ*Lk*=2 is unstable for . Furthermore, we find that and .

When the circular state is unstable, the DNA deforms into a supercoiled configuration. Each column in figure 3 depicts limiting steady-state equilibrium configurations at specified values of excess link and salt concentration when the initial perturbation is chosen in accord with equation (3.6) with *ϵ*_{1}=1, and *ϵ*_{m}=0 for *m*≥2. (Although each of these configurations is locally stable, they should not be viewed as the minimum energy configurations for there may be other configurations at the same value of Δ*Lk* with even lower energy.)

The shapes of equilibrium configurations, supercoiled plasmids can be classified according to their symmetry (Coleman & Swigon 2000). In the absence of electrostatics, the minimum elastic-energy configuration for any value of excess link has the *D*_{2} symmetry group (dihedral group of order 2) with three axes of rotational symmetry. Additional stable configurations with higher symmetries (*D*_{3}, *D*_{4}, and so on) can be found for larger values of excess link. There are two types of supercoiled configurations: plectonemically supercoiled configurations have a small number of terminal loops and large two-ply plectonemic regions in which two helical DNA segments run in parallel. Toroidally supercoiled configurations have multiple terminal loops and DNA segments at the points of self-contact run transversally. The distinction is meaningless for Δ*Lk*<3. In figure 3, among the limiting configurations with Δ*Lk*≥3 those with (Δ*Lk*,*C*_{s})=(3,0.01 M), (5,0.1 M), (6,0.05*M*) and (6,0.1 M) are plectonemically supercoiled, while the remaining ones are toroidally supercoiled.

The shapes of limiting configurations at high salt resemble those found for elastic rods without electrostatic repulsion (Coleman & Swigon 2000). The DNA at low salt has the tendency to limit the total extent in close contact, which results in a preference for toroidal supercoiling.

In figure 4, we compare the bending, twisting and electrostatic energies of configurations in figure 3. (The electrostatic energies are shown as a difference from the electrostatic energy of the circle at the same *C*_{s}.) Bending energy increases with Δ*Lk*, but does not vary much with salt concentration. On the other hand, the twisting energy is independent of Δ*Lk* at medium and high salt concentration, but depends strongly on Δ*Lk* at low salt concentration. The differences in electrostatic energy among the configurations are much smaller than the differences in either bending or twisting energy—there is only a 6.43*kT* increase in electrostatic energy at *C*_{s}=0.1 M as Δ*Lk* varies from 1 to 6, although the configurations vary substantially; the bending energy, for example, changes by 34.52*kT*. This reflects the charge screening—charges that are located more than several Debye lengths apart (which, in this case, equals *κ*^{−1}=10 Å) do not interact with each other. Writhe increases monotonically with both increasing Δ*Lk* and increasing *C*_{s}.

The stretching energy (not shown) plays a minor role, in accord with observations of Westcott *et al.* (1997), who showed that supercoiling dynamics of DNA is not much affected by DNA extensibility.

To assess the compaction of the supercoiled molecule, we measured the distance of closest approach, i.e. the closest distance between the centre of a given base pair and the centre of another base pair that has at least 20 base pairs separation along the axis of the DNA (figure 5). The distance of closest approach decreases with increasing molar salt concentration and increasing Δ*Lk*, as a result of increased supercoiling and decreased electrostatic repulsion. The distance of closest approach can be compared in figure 5 with experimental estimates of the effective DNA radius, based on measurements of knotting probability in the absence of supercoiling (Rybenkov *et al.* 1993). In those experiments, no supercoiling was imposed and hence we plot their estimated values at an arbitrary value Δ*Lk*=1.6. The values of closest approach we calculate are smaller than the effective DNA diameter at the same ionic strength because supercoiling induces forces that push the DNA into closer proximity. The values of plectonemic radius versus ionic strength reported by Ubbink & Odjik (1999) and Rybenkov *et al.* (1997) are about twice as large as those we find, most likely because our values are for equilibrium configurations and do not account for thermal fluctuations that would tend to increase the plectonemic radius.

Numerical results also show that for a fixed salt concentration, the distance between two adjacent base pairs (i.e. the stretching of the DNA) is essentially constant, regardless of different initial twist number Δ*Lk* (in accord with Westcott *et al.* 1997) and it slightly increases as the salt concentration decreases. Figure 6 shows five different projections of the equilibrium configurations for (Δ*Lk*,*C*_{s})=(5,0.01 M) and (5,0.05 M), to demonstrate that there are no actual self-contact points anywhere along the configuration as the distance of closest approach is larger than *D*.

In the theory of elastic rods with self-contact, multiple locally stable equilibrium configurations exist for Δ*Lk*>1 (Coleman & Swigon 2000; Coleman *et al.* 2000). In the presence of electrostatics, two types of locally stable equilibria (lobed and interwound) were found by Schlick *et al.* (1994) for several values of Δ*Lk*. It is not known, however, how many local equilibria exist or whether they all can be attained by a path from a perturbed circular configuration, i.e. whether the domain of attraction of each local equilibrium extends up to a neighbourhood of the circular equilibrium. We have examined this question by sampling over different perturbations of circular configurations with excess link Δ*Lk*=6 and the molar salt concentration *C*_{s}=0.1 M, and observed which configurations appeared as the limits of corresponding trajectories. Our search led to a total of five distinct limiting configurations shown in figure 7. Among those are the 2-loop plectoneme, which minimizes the total energy, two types of 3-loop configurations, one 4-loop configuration and one 5-loop configuration. The 4-loop configuration was the most commonly occurring limiting configuration.

### (b) Dynamics of transitions

The time of transition from the perturbed circle to a stable equilibrium configuration for the cases shown in figure 3 varies greatly with Δ*Lk* and *C*_{s}. The dynamics of each trajectory that produced a configuration in figure 3 can be followed on the plots of total energy, kinetic energy and writhe versus time in figures 8⇓–10.

As expected, for every transition path, the total energy decreases with time due to the strong dissipation in the system. The plateaus in the graphs of total energy coincide with local minima in the kinetic energy and correspond to transient states near unstable equilibrium configurations. The transient states resemble the unstable equilibria of elastic rods with self-contact in the absence of electrostatics, which can be found using explicit solutions of the governing equations of Kirchhoff’s theory of rods (Coleman & Swigon 2000; Coleman *et al.* 2000).

The time of convergence to the limiting configuration can be estimated from the graphs of writhe versus *t* (figure 10) as the largest *t* for which , where *δ* is some small positive number, say 0.01. The convergence time for each trajectory that produced a configuration in figure 3 can be found in table 2. The convergence time is longest for configurations with Δ*Lk*=2 regardless of ionic strength. For Δ*Lk*>3, the convergence times differ depending on whether the limiting configuration shows plectonemic or toroidal supercoiling. Toroidally supercoiled configurations are attained much faster at around 10.1 μs on average, independent of Δ*Lk* and *C*_{s}. Plectonemic configurations are attained on average in 47.3 μs.

The most prominent transient state, found as a local minimum in the plot of the kinetic energy in figure 9, can be seen in the trajectory for (Δ*Lk*,*C*_{s})=(3,0.01 M). Other transient states are in (Δ*Lk*,*C*_{s})=(5,0.05 M), (5,0.1 M) and (6,0.05 M) (figure 11).

The fastest transition is observed for Δ*Lk*=1 as that trajectory involves only equilibration of the twist along the DNA and no deformation of the axial curve. In this case, the relaxation of kinetic energy obeys a power law, i.e. *E*_{K}∼*t*^{−m}, with *m*=1.8±0.1, which is independent of *C*_{s}.

For Δ*Lk*=2, the transition paths are also independent of *C*_{s}. They start with an equilibration of twist and approach to the unstable circular configuration, followed by a buckling transition along the branch of configurations with *D*_{2} symmetry and terminating at a figure 8 configuration.

For Δ*Lk*=3, the transition paths depend on *C*_{s}. After the initial twist equilibration, all three trajectories proceed towards the unstable equilibrium similar to the configuration shown in figure 11*b*i. Then, for *C*_{s}=0.01 M, the trajectory continues towards a plectonemic configuration in figure 11*c*i. For *C*_{s}=0.05 M and 0.1 M, however, the trajectory continues towards the configuration with *D*_{3} symmetry, seen in figure 3. This transition is slower than the one for *C*_{s}=0.01 M, but because it does not dwell in a transient state, the final writhe is attained at about the same time (table 2).

For Δ*Lk*=4, the transition paths are again essentially independent of the ionic strength and similar to the path for (Δ*Lk*,*C*_{s})=(3,0.05 M), with no significant dwell time near any locally stable equilibrium.

For Δ*Lk*=5, the transition paths depend on *C*_{s}. After the initial twist equilibration, the trajectory for *C*_{s}=0.01 M proceeds directly towards a configuration with four terminal loops. For *C*_{s}=0.05 M and 0.1 M, the trajectories proceed to a transient shown in figure 11*b*ii and 11*b*iii. Then, for *C*_{s}=0.05 M, the trajectory proceeds towards a second transient state figure 11*c*ii and continues towards a limiting configuration with four lobes. For *C*_{s}=0.1 M, however, the trajectory departs towards a transient state figure 11*c*iii with three terminal loops and continues with a slow equilibration of the writhe towards the final configuration with two plectonemic loops.

Finally, for Δ*Lk*=6, the transition paths again depend on *C*_{s}. After the initial twist equilibration, the trajectory for *C*_{s}=0.01 M proceeds directly towards a configuration with four terminal loops, similarly as for Δ*Lk*=4. For *C*_{s}=0.05 M and 0.1 M, the transition proceeds along a path similar to (or exactly like) the one in figure 2 with one transient state at *t*=2.6 μs (for *C*_{s}=0.1 M). The final equilibration of the writhe is again a slow process because of the presence of a long plectonemic loop.

Figure 11*a*ii–*d*ii and 11*a*iii–*d*iii illustrates the remarkable sensitivity of the supercoiling process to ionic strength, whereby a small difference in salt concentration causes the system to depart the same transient configuration along markedly different paths.

With the exception of the initial twist relaxation, the approach of a configuration **X**(*t*) to a local equilibrium **X**_{eq} is close to exponential, and hence we can define the relaxation time of a local equilibrium, *τ*, as ‖**X**(*t*)−**X**_{eq}‖≈*e*^{−t/τ}. The relaxation times can be estimated from the graph of the kinetic energy, which scales as *e*^{−2t/τ}, when equilibrium is approached. Relaxation times for the limiting configurations with Δ*Lk*>1, shown in figure 3, are given in table 3. They are computed for the transition phase after which no further contacts between DNA segments are formed and no further rearrangement of contacts occurs. With the exception of Δ*Lk*=4, for which we observe the fastest relaxation times, for any other given Δ*Lk*, the relaxation time gets larger as salt concentration increases. The toroidally supercoiled configurations have shorter relaxation times than the plectonemically supercoiled ones.

## 4. Summary and conclusion

In this paper, we present an improved generalized IB method for simulation of the dynamics of elastic rods in a fluid. By including electrostatic repulsion and steric (hard-core) interaction between rod segments, we completely eliminate any self-penetration of the modelled rod during simulations and assure the conservation of linking number. The model can be applied to macroscopic structures, such as underwater cables, ropes and fishing lines (charged or not), to cellular components, such as bacterial filaments and microtubules, and to short semiflexible polymers, such as DNA.

We here focus on supercoiling of twisted DNA rings, and explore the dependence of the dynamical trajectories and limiting configurations on the excess link, ionic strength and perturbation. Our goal is to obtain a better idea about the kinetics of supercoiling processes in short molecules, motivated by the possibility that loop-formation kinetics, not just equilibrium, may be important in biological processes in which binding and unbinding of proteins happens on fast time scales. We confirm the earlier findings (Schlick *et al.* 1994; Westcott *et al.* 1995) that the type of limiting configuration is influenced by ionic strength in which higher salt gives rise to plectonemic supercoiling, while low salt gives rise to toroidal supercoiling. In all the limiting structures we found, the repulsion is sufficiently strong that the DNA surface remains free of direct contact.

We find that the dynamics of the ring 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 total energy. Similar plateaus were observed in the dynamics of elastic rods with self-contact (Goyal *et al.* 2008). 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 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.

We confirm the existence of multiple locally stable configurations for the given ionic strength and excess link, and find that each such configuration can be attained with an appropriate initial perturbation of the circle. Our search recovered all locally stable equilibria known from the theory of elastic rods with self-contact (Coleman & Swigon 2000), and hence, we conjecture that any locally stable equilibrium configuration of a closed rod can be attained from the circular configuration with the same excess link by a dynamical trajectory that starts with an appropriate initial perturbation in twist density.

Our findings on the departure from the circular configuration are in accord with the results from perturbation analyses of looping and ring collapse transitions (Goriely & Tabor 1997) and numerical analyses of supercoiling of DNA plasmids (Klapper 1996), in the sense that if the excess link is large enough, the departure from the unstable twisted circular configuration is fastest along the mode with four loops. It is not the only unstable mode and other modes can be excited with different perturbations—these modes generally lead to other locally stable configurations.

At the length scales corresponding to the size of our minicircles, random fluctuations, not accounted by our method, have an additional effect on DNA supercoiling dynamics, primarily by enabling transitions between various locally stable configurations and perturbing the configurations in such a way that the probability of occurrence of any configuration is given by the Boltzmann distribution. Stochastic dynamics of 1000–1500 bp DNA plasmids has been studied using Brownian dynamics with excluded volume (Chirico & Langowski 1994) and electrostatics (Jian *et al.* 1998). The average relaxation times for writhe reported in Chirico & Langowski (1994) are 6 μs for Δ*Lk*=4 and 3.5 μs for Δ*Lk*=6, those reported in Jian *et al.* (1998) are 20 μs for Δ*Lk*=4 and 10 μs for Δ*Lk*=5.5, both at *C*_{s}=0.01 M. Our results indicate that relaxation time depends on the limiting configuration, but is in the range of 0.8–8 μs. Thermal fluctuations can be accounted for in the framework of the generalized IB method by adding a random force with a Boltzmann distribution acting on the fluid (Atzberger *et al.* 2007). We plan to combine electrostatic repulsion and thermal fluctuations in a future implementation of the generalized IB method.

The theory of counterion condensation, employed in this paper, is based on approximations that are valid only in the limit of low salt concentrations and for an infinitely long, straight polyion. Although the theory has been shown to provide good estimates of DNA repulsive forces up to the values of ionic strength considered here, a more accurate and potentially interesting extension of the current work would be to account explicitly for the dynamics of ion distribution and its interaction with the DNA. It is possible that the ions would localize near the plectonemic regions of DNA, which would make configurations showing such regions energetically favourable.

## Acknowledgements

We are grateful to Charles Peskin, David McQueen and Estarose Wolfson for the use of their visualization software. Y.K. was supported by the National Research Foundation of Korea, grant 2010-0006165. D.S. acknowledges support by an Alfred P. Sloan Fellowship and grant RGP0051/2009 by Human Frontiers in Science Project.

- Received March 31, 2010.
- Accepted July 19, 2010.

- © 2010 The Royal Society