INTRODUCTION

Top Abstract Introduction Methods Conclusion References

Understanding the role of DNA three-dimensional (3D) structure in its biological function requires a quantitative description of the structure and dynamics of long DNA chains in their various states. Double-helical DNA is a wormlike polyelectrolyte chain, whose conformationneglecting atomic detailcan be described by a simple set of parameters: bending and twisting elasticity, hydrodynamic diameter, and electrostatic interactions between different parts of the same molecule.

The bending elasticity of B-DNA is given by its persistence length of 50 nm and does not strongly depend on monovalent ion concentration between 0.01 and 1 M. In the computations referred to in this paper, we work with a constant persistence length of 50 nm (Hagerman, 1988). For the torsional elasticity we take C = 2.0 · 10¹⁹ erg cm (Schurr et al., 1992). The hydrodynamic properties of DNA can be well described by a hydrodynamic diameter of 2.4 nm (Hagerman and Zimm, 1981). Long-range intramolecular interactions are otherwise due to electrostatic repulsion between the negatively charged phosphate groups on the helix backbone. One of the main points of this paper is to develop an efficient procedure to compute electrostatic interactions in DNA.

A numerical simulation of the time-dependent structure of a long chain molecule like DNA in solution is equivalent to solving the equations of motion for this molecule plus the surrounding solvent. Given the size of the moleculesome 1000 basepairs of DNA, corresponding to ~10⁵ atomsa numerical description at the atomic level on present-day machines is unthinkable over any reasonable time scale. But since one is mainly interested in global features of the structure of the DNA molecule, such as cyclization probabilities, as computed in the accompanying paper, it is possible to view the DNA as a linear chain of segments, each including several tens of basepairs, and to model the solvent as a continuous viscous fluid. The thermal motion of the solvent molecules is included in the equations of motion through a random force.

This approach of describing the dynamics of a system is called Brownian dynamics (BD), a method that has been applied in several cases to modeling DNA dynamics (Allison, 1986; Allison et al., 1989, 1990; Chirico and Langowski, 1994). Up to now, in most cases new BD models were developed for each particular case. However, the power of the BD method and the range of problems in the 3D structure of DNA and other wormlike polymers where it can be successfully applied justifies the implementation of a general-purpose BD package for computing polymer dynamics, with special regard for DNA. In this paper we describe the implementation of such a program. An accompanying paper (Merlitz et al., 1998) shows its application to the problem of DNA loop closure probabilities in the presence of bound proteins and/or bends.

	METHODS

Top Abstract Introduction Methods Conclusion References

The DNA molecule is modeled as an elastic chain with electrostatic interactions. The principal approach taken here is similar for linear and circular DNA; specific differences between the two cases will be mentioned where they apply. The conformation of a chain of N straight segments is specified by the space positions of its vertices, r_i, i = 0, ... , N. The segments are represented by the vectors s_i = r_i+1  r_i, i = 0, ... , N  1. To each segment a system of three orthogonal vectors of unit length is attached (f_i, g_i, e_i), so that the e_i direction coincides with the direction of the ith segment: e_i = s_i/s_i, s_i = |s_i| (Fig. 1). For a circular chain we have the additional restriction that r_N = r₀.

View larger version (10K): [in this window] [in a new window]   FIGURE 1   Chain geometry with the segment vectors si and the segment coordinate systems (fi, gi, ei) that define the relative orientation between segments.

The total energy of a given chain conformation is given as the sum of the stretching, bending, twist, and electrostatic energies.

The stretching energy is defined for each segment i:

(1)

Here k_B is the Boltzmann constant, T the temperature, l₀ the segment equilibrium length, and  the stiffness parameter, so that (l₀)² is approximately equal to the variance of the segment length distribution.

The bending energy is defined for each chain joint. We call a joint unbent if it connects segments that form a straight line at equilibrium, and bent if the angle ^*_i between the segments at equilibrium is nonzero. To each bent joint i we attach an auxiliary unit vector b_i that is fixed in the coordinate system (f_i, g_i, e_i), its polar coordinates being (^*_i, ^*_i). Note that this formalism is different from our first implementation (Chirico and Langowski, 1996) where a heuristic "kink potential" was used that was given in terms of the Euler angles for rotating one segment into the next. Here the bending energy of the ith joint is:

(2)

where _i is the angle between e_i1 and e_i for an unbent joint, or between e_i1 and b_i for a bent joint; _b is the bending rigidity parameter chosen in such a way that the Kuhn length is equal to:

(3)

where The twist energy is defined for each adjacent segment pair:

(4)

where C is the torsional rigidity constant and _i is the twist angle between the (i  1)th and ith segments.

The twist angle _i is calculated by defining a vector p_i = s_i1 × s_i, which is normal both to s_i1 and s_i. Now, we can easily calculate the angles _i between f_i1 and p_i, and _i between p_i and f_i. Then, _i = _i + _i (Fig. 2). During the BD simulation we assume that at time t, _i(t  t)    _i(t)  _i(t  t) + , where _i(t  t) is the twist angle one simulation step ago (t is chosen such that the probability that the twist angle changes by more than ± becomes negligible).

View larger version (8K): [in this window] [in a new window]   FIGURE 2   Definition of the twist angle i = i + i. pi is perpendicular to si1 and si; i is measured between fi1 and pi, i between pi and fi.

The starting point for the electrostatic energy is the expression for the energy of interaction between two uniformly charged nonadjacent segments (i, j) in a 1:1 salt solution in the Debye-Hückel approximation:

(5)

The integration is done along the two segments; _i and _j are the distances from the segment beginnings, r_ij is the distance between the current positions at the segments to which the integration parameters _i and _j correspond;  is the inverse of the Debye length, so that ² = 8e²I/k_BTD, I is the ionic strength, e the proton charge, D the dielectric constant of water, v the linear charge density which for DNA is equal to v_DNA = 2e/, where  = 0.34 nm is the distance between basepairs.

There are two problems to be solved for the Eq. 5: 1) the linear density v should be renormalized from that of DNA to a smaller value in order to ensure the correct excluded volume effects; 2) the integration should be approximated by a more simple procedure to save computation time.

The renormalization of the linear density was done as in Stigter (1977). As pointed out in Schellman and Stigter (1977), the Gouy layer of immobile counterions reduces the effective charge density by a factor of q = 0.73 for NaCl concentrations between 1 and 500 mM. Next, the Debye-Hückel approximation is a linearization of the Poisson-Boltzmann equation and valid only for a very small electric potential:   k_BT/e. We choose the renormalized charge density v* in such a way that the known solution of the Debye-Hückel equation for a straight thin line with charge density v* coincides with the solution of the Poisson-Boltzmann equation for a cylinder of the DNA radius r_ES and the charge density qv_DNA in the regions where   k_BT/e (Fig. 3).

View larger version (11K): [in this window] [in a new window]   FIGURE 3   Electrostatic potential for the interaction between DNA segments. A Debye-Hückel potential (1) was renormalized (3) such as to coincide at large distances with the nonlinear Poisson-Boltzmann equation (2).

In order to save computation time, a tabulation of the double integral (Eq. 5) was used. For simplicity we assume that each segment has the same length l₀. For each segment pair (i, j) we can define a vector l₀R_ij connecting the middle points of the two segments. Then the mutual position of the segments can be defined by the following four dimensionless parameters:

(6a)

(6b)

(6c)

(6d)

Equation 5 can be rewritten in the form:

(7)

where _e = v*²l₀/k_BTD and

(8)

(9)

where v₀, v₁, and v₂ are unit vectors oriented in such a way that v₀v₁ = ₁, v₀v₂ = ₂, and (v₁ × v₀) · (v₂ × v₀)/ |v₁ × v₀v₂ × v₀| = .

In the following we need the partial derivatives of Eq. 8. At first, a four-dimensional table for the values of the integral (Eq. 8) and each of its partial derivatives was constructed numerically. Then, during the simulation, a linear interpolation was used to obtain the values of (f/), (f/₁), (f/₂), (f/s) at particular points of the (, ₁, ₂, ) space. The table steps were chosen in such a way that the values of , arccos ₁, arccos ₂, and arccos  had constant increments. The tabulation range for ₁, ₂, and  is [1,1]. The range of , [_min, _max], and the table size for each argument are parameters of the approximation, which we chose by the following criteria. For the minimum distance, _min, all possible values of the electrostatic energy should be large enough (e.g., >10 k_BT), so that this distance is practically unreachable during the simulations. For distances >_max all possible energy values should be negligible (say, <0.01 k_BT). The mutual displacement of segments corresponding to one -step should not exceed the Debye length, and the same restriction is applied to the displacement of segment ends corresponding to one step in the ₁, ₂, and  dimensions. These criteria are rather "soft," in order to keep the total table size within reasonable limits (the minimum table size for l₀ = 10 nm and I = 1 M is 14 MB of memory, including all partial derivatives of Eq. 8). We should note, however, that even crude hard-core potentials for electrostatic repulsion in DNA can be applied in many cases to predict statistical properties of DNA to good precision (Vologodskii and Cozzarelli, 1995); therefore, the representation of the electrostatic potential according to Eqs. 7 and 8 is probably a good approximation.

So far, we neglected the fact that the segment length is only approximately equal to l₀. That means that the "charged" segment does not coincide exactly with the "geometrical" segment. This seems to be a good approximation for the excluded volume effects, since the chain is supposed to be stiff with respect to stretching (  1). One has to be careful, however, to avoid that during the simulation parts of the chain cross each other through discontinuities between adjacent charged segments, if the length of their phantom geometrical counterparts is greater than l₀. In order to exclude this possibility we define the R_ij vector for two segments (i, j) of arbitrary length in the following way:

(10)

where r_i(j)^(m) is a "shifted" middle point of the segment:

(11)

and r_j^(m) = (r_j + r_j+1)/2 is the actual middle point of the segment. This means that the geometrical segment that is used to calculate the excluded volume interaction is shifted with its end toward the segment joint that is closest to the other segment in the interaction. Thereby, any gaps that might appear at the joint due to stretching are automatically closed and the chains cannot cross.

Since forces and torques are the partial derivatives of the energy over the system coordinates, the latter should be specified formally. For each segment i we chose the following four coordinates: the three space coordinates r_i of the segment beginning and the angle _i of rotation of the local vector basis (f_i, g_i, e_i) around the e_i axis. For the angle coordinate the zero position is not defined. Therefore, in addition, we need to specify how to keep the _i coordinates unchanged while a displacement of the r_i coordinates takes place. The following rule was used to derive forces and torques and to perform moves in the simulation procedure. If, as a result of r_i displacement, the e_i vector has a new value e'_i and A_i is a rotation matrix, so that e'_i = A_ie_i, then the new value of the f_i vector is f'_i = A_if_i.

In a linear chain three additional degrees of freedom r_N are required for the last segment. Here are the expressions for the forces and torques for the ith segment obtained by differentiating Eqs. 1, 2, 4, and 7:

Stretching force

The stretching force acting on the ith vertex from the (i + 1)th one is parallel to the ith segment:

(12)

The total stretching force for the ith vertex is therefore

(13)

Bending force

The contribution of the energy stored in the bending angle _i to the bending force acting on the ith vertex from the (i + 1)th one is perpendicular to the ith segment and lies in the bend plane:

(14a)

for an unbent joint, and

(14b)

for a bent joint. Here _i = p_i/p_i; p_i = s_i1 × s_i, so that _i = e_i1 × e_i/sin ^*_i, where ^*_i is the angle between e_i1 and e_i (for an unbent joint, ^*_i = _i).

The analogous contribution to the bending force acting on the ith vertex from the (i  1)th one is perpendicular to the (i  1)th segment and also lies in the bend plane:

(15a)

for an unbent joint, and

(15b)

for a bent joint.

Note that Eqs. 14a and 15a for an unbent joint are symmetric with respect to renumbering the vertices in opposite order (one force can be obtained from the other by substituting s_i  s_i1, e_i  e_i1, _{i  i1}). Obviously, there is no such symmetry for a bent joint because the b_i vectors are not defined in a symmetrical way. We note also that the expressions for a bent joint (14b, 15b) become those for the unbent joint upon substituting b_i  e_i.

The total bending force for the ith vertex is

(16)

Bending torque

The torque on segment i induced by bending is for a bent joint:

(17)

For an unbent joint, this torque is equal to zero.

Twisting force

The force on the ith vertex from the (i + 1)th one induced by mutual twisting of the (i  1)th and ith segments by an angle _i is perpendicular to the plane of the ith bend:

(18)

The symmetric expression for the twisting force acting on the ith vertex from the (i  1)th one is

(19)

The total twisting force for the ith vertex is then

(20)

Twisting torque

The torque on the ith segment induced by twisting the (i + 1)th segment by an angle _i+1 with respect to the ith one is

(21)

The total twisting torque for the ith segment is then

(22)

Electrostatic force

The contribution of the electrostatic interaction between segments i and j to the force acting on the ith vertex is where f_ij = f(_ij, _ij, _ji, _ij). An analogous relationship is valid for the force F_ij,2^(e) acting on the (i + 1)th vertex from segment j. We assume that F_ij,1^(e) = F_ij,2^(e) = 0 when i = j or i = j ± 1. The partial derivatives of f (Eq. 8) are tabulated as described above. The expressions for the derivatives of _ij, _ij, _ji, and _ij with respect to r_i and r_i+1 are given below, where the following auxiliary notations are used: _ij = R_ij/_ij is a unit vector in the R_ij direction; _i(j) and _i(j) are the coefficients in the expression for the "shifted" middle-point of a segment: so that (see Eq. 11): The upper variant corresponds to the condition (a) in Eq. 11, the lower variant corresponds to the condition (b).

The partial derivatives are thus:

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

In Eqs. 29 and 30 the positive value of the square root is taken when R_ij · (s_i × s_j) > 0, and its negative value otherwise.

The total electrostatic force acting on the ith vertex is

(31)

The boundary conditions for the force expressions, Eqs. 12-31, are different for linear and circular chains. For a linear chain all parameters with indexes out of range do not exist and can be formally set to zero. The allowed ranges are 0  i  N  1 for F_i,next^(s), T_i,next^(t); 1  i  N  1 for F_i,next^(b), F_i,prev^(b), T_i^(b), F_i,next^(t), F_i,prev^(t); 0  i  N  1, 0  j  N  1 for F_ij,1^(e), F_ij,2^(e). For a closed circular chain circular boundary conditions are in effect, such that all indices have to be taken modulo N.

Hydrodynamic interactions

In order to model hydrodynamic interactions defined between spherical objects, a bead with radius a was attached to each chain vertex. With N' = N + 1 beads for a linear and N' = N beads for a circular chain, we describe the hydrodynamic interaction between beads i and j by a 3 × 3 Rotne-Prager tensor:

(32a)

(32b)

(32c)

where r_ij = r_j  r_i, r_ij = |r_ij|, D₀ = k_BT/6a,  is the water viscosity, I is a unit 3 × 3 matrix, and r  r denotes a matrix with the components (rr); ,  = x, y, z.

The bead radius a is connected with the DNA hydrodynamic radius r_HD in the following way. Let us consider two geometrical objects, both modeling a piece of DNA of Kuhn length B. The first object is a straight line with n beads attached to it so that the distance between neighboring beads is l₀. The number of the beads is equal to the closest integer to B/L₀, and the bead radius is a. The second object is a cylinder of the radius r_HD and the length nl₀. We choose the bead radius a in such a way that the diffusion coefficients of the both objects have the same value. The diffusion coefficient for the string of beads is [see, for example, Hagerman and Zimm (1981)]: where the tensors S_ij are related to the diffusion tensors D_ij for the bead row in the following way: _j=1ⁿ S_ijD_jk = _ikI, _ik is the Kronecker delta symbol.

According to Tirado and Garcia de la Torre (1979, 1980), the diffusion coefficient of a cylinder is:

(33a)

with

(33b)

(33c)

(33d)

Brownian dynamics algorithm

The simulations were performed using a second-order Brownian dynamics algorithm. A chain displacement from a given conformation {r_i(t), _i(t)} at time t was done in two half-steps. A tentative first-order displacement is:

(34a)

(34b)

where D_ij(t), F_j(t), T_i(t) are calculated for the given conformation as outlined above; D^rot = k_BT/4r_RD²l₀ is the rotational diffusion coefficient assumed to be equal to that of a cylinder of the DNA radius r_RD and the length l₀ (without end-effect corrections), t is the time step, and the random displacements R_i, _i have the following covariances:

(35a)

(35b)

The R_i values were obtained as R_i = _i=0^N A_ikX_k, where X_k are uncorrelated Gaussian-distributed random vectors with zero mean and unit variance for each component. The A_ik tensors are related to the D_ij tensors as: _k=0^N A_ikA_jk^T = D_ij. A_ik was obtained through a Cholesky factorization.

The final half-step is:

(36a)

(36b)

Here F'_j(t + t) and T'_i(t + t) are the forces and the torques calculated for the conformation {r'_i(t + t), '_i(t + t)}. Since the values of the diffusion tensor D_ij depend on the chain conformation much weaker than forces and torques, they are kept the same as in the first half-step.

The initial chain conformation is of no importance for a linear chain. However, in modeling a closed circular DNA, the initial conformation determines the linking number difference (White, 1989)

(37)

that, being a topological invariant, does not change during the simulation procedure. In Eq. 37 Tw = (1/2) _i=0^N1 _i is the deviation of the DNA twist from its equilibrium value in a linear unstressed chain and Wr is the writhe of the chain.

As an initial conformation we use a flat regular N-sided polygon with Wr = 0 and _i = 2Lk/N for all segments. Lk does not necessarily have to be an integer; in that case, ₀ = (₀ + ₀ + 2 Lk) mod 2  . The invariance of