Article

Simulation of DNA Supercoil Relaxation Ikenna D. Ivenso1 and Todd D. Lillian1,* 1

Department of Mechanical Engineering, Texas Tech University, Lubbock, Texas

ABSTRACT Several recent single-molecule experiments observe the response of supercoiled DNA to nicking endonucleases and topoisomerases. Typically in these experiments, indirect measurements of supercoil relaxation are obtained by observing the motion of a large micron-sized bead. The bead, which also serves to manipulate DNA, experiences significant drag and thereby obscures supercoil dynamics. Here we employ our discrete wormlike chain model to bypass experimental limitations and simulate the dynamic response of supercoiled DNA to a single strand nick. From our simulations, we make three major observations. First, extension is a poor dynamic measure of supercoil relaxation; in fact, the linking number relaxes so fast that it cannot have much impact on extension. Second, the rate of linking number relaxation depends upon its initial partitioning into twist and writhe as determined by tension. Third, the extensional response strongly depends upon the initial position of plectonemes.

INTRODUCTION Supercoiling is intimately coupled with a variety of cellular processes including transcription, replication, supercoil regulation, and DNA compaction (1–7). Such processes not only alter the degree of supercoiling in DNA, but also respond to it. During transcription, for example, RNA polymerase is sensitive to the degree of supercoiling (8,9) and induces supercoiling as DNA corkscrews through it (4,8,10,11). Consequently, cells regulate supercoiling with an assortment of vital enzymes known as topoisomerases (2,12–14). Because of their significant cellular role, topoisomerases have become targets for anticancer and antibiotic drug therapies (15–17). While topoisomerases employ a variety of mechanisms, type IB topoisomerases alter DNA topology by nicking one of the two backbones, allowing restrained rotation between the upstream and downstream DNA (18), and re-ligating the nicked backbone. A similar group of enzymes, nicking endonucleases, also alter topology by nicking a single backbone but lack the ability to re-ligate DNA. Characterizing the dynamics of supercoil relaxation, in response to these enzymes, is necessary to advance our fundamental understanding of cellular biology and pave the way for future drug therapies. Single-molecule experimental studies have measured the dynamics of supercoil relaxation in response to nicking by

Submitted October 26, 2015, and accepted for publication March 28, 2016. *Correspondence: [email protected] Editor: Keir Neuman. http://dx.doi.org/10.1016/j.bpj.2016.03.041 Ó 2016 Biophysical Society.

2176 Biophysical Journal 110, 2176–2184, May 24, 2016

type IB topoisomerases (15,18–20) and nicking endonucleases (18,20,21). In such experiments, DNA is anchored at one end to a fixed substrate and at the other to a micronsized paramagnetic bead (see Fig. 1). An external magnetic field is used to manipulate the bead and measure the dependence of extension (bead height) on tension and superhelical density. This calibration is used to calculate a DNA molecule’s instantaneous superhelical density from its extension during relaxation. In such experiments, relaxation by type IB topoisomerases proceeds at ~100 turns/s despite friction between the DNA and enzyme (1,18). A reduced relaxation rate (~3 turns/s) is observed in the presence of the chemotherapeutic drug topotecan (15). Crut et al. (21) estimate that supercoil relaxation could proceed as fast as 10,000 turns/s in response to a nicking enzyme. However, such fast dynamics cannot be resolved accurately by the overdamped motions of a large micron-scale bead (1,21–24). Computational and analytical models have the potential to describe supercoil dynamics with detail inaccessible to experiments. Analytical models efficiently represent the global structure of supercoiled DNA with a handful of parameters (e.g., linking number and plectoneme size). These analytical models provide insight into, for example, the formation and structure of stretched supercoiled DNA (25–29) and the relaxation of supercoils by type IB topoisomerase (25,30). Molecular dynamics simulations offer atomic detail in their representation of supercoiled DNA. Such molecular dynamics simulations reveal, for example, a stress-softening of DNA with bending (31), the formation of kinks and

Simulation of DNA Supercoil Relaxation

number (Lk) and extension after a nicking event. We investigate the influence of the applied tension on Lk relaxation dynamics. Finally we explore the dependence of relaxation dynamics on initial plectoneme position.

MATERIALS AND METHODS Discrete wormlike-chain model

FIGURE 1 Cartoon of recent single-molecule magnetic tweezer experiments in which DNA is attached to a fixed surface below and to a magnetic bead above. The bead is manipulated by a magnetic field. (Inset) We represent supercoiled DNA using a discrete wormlike-chain model. Figure modified and inspired from Crut et al. (21), Lillian et al. (47), and Ivenso and Lillian (60).

bubbles in tightly bent and supercoiled DNA (32–35), and a toroid structure formed by compressed DNA confined to a small cavity (36). Unfortunately, the computational expense of these all-atom representations limit simulations to short length- and timescales. Continuum and coarse-grain models serve as a compromise between the detail and expense of molecular dynamics and the simplicity of analytical models. These models are capable of accounting for details such as DNA sequence-dependent properties (37–41), bending nonlinearities (42–45), and interactions with proteins (36,46–49), for example, while simulating much longer length- and timescales than molecular dynamics. In addition, these models can be paired with molecular dynamics to simulate over multiple scales; see, for example, Taranova et al. (31), Hirsh et al. (36), Villa et al. (46), and Lillian et al. (47). Previously, we developed a multiscale model for supercoil relaxation in response to type IB topoisomerase (47) by representing DNA-protein interactions with a simple potential derived from an all-atom molecular dynamics study (50). Unfortunately, omitting thermal fluctuations in our previous continuum model for DNA limited its application to very short lengths of DNA (approximately one persistence length). In this study, we employ a coarse-grain model (discrete wormlike-chain model) to quantify the dynamics of DNA supercoil relaxation in response to nicking under constant tension. Discrete wormlike-chain (dWLC) models (51–56) are widely used in simulations of supercoiled DNA (see, for example (24,48,56–62)) and many of their predictions agree well with experimental observations. Using our dWLC model, we compare the relaxation rates of linking

Our dWLC model is patterned after (and uses notation largely duplicating) that of Klenin et al. (51) and similar models presented in a variety of published works; see, for example (52–56,60,61). In our dWLC model, DNA is discretized into a chain of N segments with free length l0. The positions of the endpoints (vertices) of the ith segment are denoted by ri ¼ fxi ; yi ; zi gT and riþ1 ¼ fxiþ1 ; yiþ1 ; ziþ1 gT . The indexing of the vertices is shown in Fig. 1 (inset). A set of three orthonormal body fixed vectors attached to each segment is used to track torsional deformations of the DNA. (In the undeformed state, these body-fixed frames are aligned across all segments.) The first of the vectors is denoted ei and is aligned with the ith segment. The other two vectors, denoted f i and gi , are related through the cross product ei ¼ f i  gi . As the molecule’s configuration changes over a time step (Dt), the vector ei remains a unit vector along the segment riþ1  ri . While the vectors f i and gi change with ei , they also depend on changes of the angle fi over a time step Dt, where Dfi ¼ fi ðt þ DtÞ  fi ðtÞ. Physically, Dfi measures the rotation of the ith segment about ei over a time step. (Because ei changes with time, fi is not a physically meaningful angle.) Mathematically,




f i ðt þ DtÞ gi ðt þ DtÞ






 cosðDfi Þ sinðDfi Þ ¼ sinðDfi Þ cosðDfi Þ

 pi , f i ðtÞ pi , gi ðtÞ pi  : pi , gi ðtÞ pi , f i ðtÞ ei ðtÞ  pi

(1) Here, pi is a unit vector along the common normal to ei ðtÞ and ei ðt þ DtÞ, specifically

pi ¼

ei ðtÞ  ei ðt þ DtÞ : jei ðtÞ  ei ðt þ DtÞ j

(2)

Harmonic potentials are used in our dWLC model to represent stretching, bending, and twisting deformations following Klenin et al. (51). The potential representing the stretch of the ith segment is

Esi ¼

kB T 2ðl0 dÞ

2 ðl0

2

 si Þ ;

(3)

where kB is the Boltzmann constant, T is the absolute temperature, d is a parameter representing the stiffness of the ith segment, and si ¼ jriþ1  ri j is the length of the ith segment. The potential representing bending between adjacent segments ((i  1) and i) is given by

Ebi ¼ ab kB Tb2i ;

(4)

where ab is a measure of the bending stiffness of DNA and bi ¼ cos1 ðei1 ,ei Þ is the bend angle between the segments. The elastic energy representing the relative twist between adjacent segments (ði  1Þ and i) is given by

Eti ¼

1C 2 q; 2 l0 i

(5)

Biophysical Journal 110, 2176–2184, May 24, 2016 2177

Ivenso and Lillian where C is the torsional stiffness of DNA and qi is the twist angle between segments. The twist angle ðqi Þ, due to relative rotations between the ði  1Þth and ith segments, is calculated by

tanðqi Þ ¼

f i , gi1  gi , f i1 : f i , f i1 þ gi , gi1

(6)

Electrostatics are approximated by distributing point charges along each segment of the dWLC model and using a Debye-Hu¨ckel potential

Ee ¼

Nl Nl X n2 l20 X l2 Dw i ¼ 1 j ¼ iþNex þ1

   exp krij  H rc  rij ; rij

(7)

3ðN þ 1Þ  3ðN þ 1Þ matrix with entries made up of 33 blocks, Dij; and F(t) is a 3ðN þ 1Þ by 1 vector representing the forces acting on each bead due to the above potentials and corresponding to a state described by rðtÞ, ei ðtÞ, f i ðtÞ, gi ðtÞ, and fðtÞ. In a similar fashion, estimates for the angles fðt þ DtÞ are obtained as

~ þ DtÞ ¼ fðtÞ þ Drd MðtÞ Dt þ F: fðt kB T

Here, fðtÞ ¼ ff1 ðtÞ; f2 ðtÞ; .fN ðtÞgT is an N by 1 vector, and Drd is an ðN þ 1Þ  ðN þ 1Þ diagonal matrix with diagonal entries given by Drd i . Next, the positions and angles at time t þ Dt are obtained as

rðt þ DtÞ ¼ ~rðt þ DtÞ þ see, for example, Vologodskii (59). Here, n is the effective charge density of DNA, l is the number of discrete point charges per segment, Dw is the dielectric constant of water, 1/k is the Debye length, and rij is the distance between point charges i and j. (Here the indices i and j index through point charge numbers rather than segment or vertex numbers.) This approximation excludes contributions from Nex nearest-neighbor point charges. This electrostatic potential is fundamentally different from the previous potentials (stretch, bend, and torsion) because it involves interactions between sites that are distant along the contour length of the DNA. To reduce the computational expense of evaluating negligible contributions to the electrostatic forces from distant point charges, we introduce a cutoff distance, rc, and the Heaviside function, H(), in the potential. When the electrostatic forces are insufficient to prevent the unphysical case of two segments of the model occupying the same volume (and possibly passing through each other), an additional constant repulsive force, Fex, is triggered between point charges; see, for example, Vologodskii (59). To account for hydrodynamics in the dWLC model, a spherical bead is located at each vertex. The bead radius, a, is chosen such that the diffusion coefficient of the chain of beads in a straight configuration approximates that of a cylinder with dimensions similar to DNA (51). Evaluating hydrodynamic interactions in our dWLC model (using the Rotne-Prager tensor, for example) imposes significant computational costs. Consequently, in this study we neglect hydrodynamic interactions (and account for only self-mobility of each bead) in favor of computational speed. This yields the constant isotropic diffusion tensor

Dij ¼ D0 Idij :

kB T 2

4phðr rd Þ l0

:

(9)

A second-order Brownian dynamics algorithm following Iniesta and de la Torre (63) is used to simulate the dynamic response of the dWLC model. In this two-part algorithm, first an estimate (denoted by a tilde) of the vertex positions, rðt þ DtÞ, at time t þ Dt is obtained as

~rðt þ DtÞ ¼ rðtÞ þ

 Dt  ~ D Fðt þ DtÞ  FðtÞ ; (12) 2kB T

and

~ ~ þ DtÞ þ Drd Mðt þ DtÞ  MðtÞ Dt: (13) fðt þ DtÞ ¼ fðt 2kB T Here, MðtÞ is an ðN þ 1Þ by 1 vector representing the torque on each segment due to the above potentials when the dWLC model is in the state ~ þ DtÞ and described by rðtÞ, ei ðtÞ, f i ðtÞ, gi ðtÞ, and fðtÞ. The force Fðt ~ þ DtÞ vectors are due to the above potentials when the moment Mðt dWLC model is in the state described by ~rðt þ DtÞ, ~ei ðt þ DtÞ, ~f i ðt þ DtÞ, ~ þ DtÞ. RðtÞ and FðtÞ represent random thermally ~gi ðt þ DtÞ, and fðt induced displacements and rotations sampled from distributions with the following properties:

hRi ¼ 0;

(14)

hR5Ri ¼ 2DDt;

(15)

hFi ¼ 0;

(16)

hF5Fi ¼ 2Drd Dt:

(17)

(8)

Here, D0 ¼ kB T=6pha is the diffusion constant of a single bead, h is the solvent viscosity, I is the 33 identity matrix, and dij is the Kronecker delta function. (In the Supporting Material, we provide evidence that the conclusions of this study would be qualitatively unchanged with the inclusion of hydrodynamic interactions.) The hydrodynamics of rotating a dWLC segment about its length is represented by a cylinder of radius rrd with rotational diffusion constant

Drd i ¼

(11)

Dt DFðtÞ þ R: kB T

(10)

Here, rðtÞ ¼ fx1 ðtÞ; y1 ðtÞ; z1 ðtÞ; .zNþ1 ðtÞgT is a 3ðN þ 1Þ by 1 vector containing the positions, ri , of all (Nþ1) vertices at time t; D is a

2178 Biophysical Journal 110, 2176–2184, May 24, 2016

Simulation In this article, our dWLC simulations are intended to parallel recent singlemolecule experiments such as Crut et al. (21) and van Loenhout et al. (64). Specifically, we represent a L ¼ 21 kilo-basepair (7 mm) length of DNA with N ¼ 700 segments. The ends of the dWLC model are connected to two parallel surfaces: a stationary surface representing the fixed substrate and a movable surface representing a large paramagnetic bead. The mobile surface is constrained such that it can only translate in the direction of the applied tension, which is along the vertical axis (z). To represent impenetrability of these surfaces, we incorporate potentials that repel any vertices (excluding the two end vertices at which the DNA is attached to the surfaces) that come within a specified distance. Specifically, the potentials are

Eup

8 > > > >
Lb > > > : Fb Lb þ Fb ðzi  zNþ1 Þ;

6Lb < ðzNþ1  zi Þ; 0 < ðzNþ1  zi Þ%6Lb ; ðzNþ1  zi Þ%0; (18)

Simulation of DNA Supercoil Relaxation

Elow

8 > > > >
Lb > > > : Fb Lb þ Fb ðz1  zi Þ; ðzi  z1 Þ%0; (19)

where z1 and zNþ1 are, respectively, the vertical components of the positions of the first and last vertices (corresponding to the positions of the surfaces); zi is the vertical component of the position of the ith vertex; Lb is the length scale over which these surface forces develop; Fb is the full repulsive force from the surfaces; and Eup and Elow are the potentials corresponding to the upper and lower surfaces. The two end segments (1 and N) are constrained to be normal to both surfaces and incapable of rotating. A constant tension, Fs, is applied to vertex (N þ 1) associated with the movable surface. No attempt is made to represent the hydrodynamics associated with the two surfaces. Our simulations are intended to represent experiments in which stretched and supercoiled DNA is at thermal equilibrium before a nicking event. Consequently, we use an ensemble of Ns ¼ 100 distinct initial conditions that are obtained from Ns ¼ 100 equilibration simulations. In each equilibration simulation, an arbitrary plectonemic state is allowed to anneal for >1 s (corresponding to over 900 h of computer time per simulation using a single computational core) under constant tension, Fs ¼ 1.6 pN, and linking difference, DLk ¼ 86. Linking difference, DLk, is a measure of DNA supercoiling and is the sum of excess twist, DTw, and writhe, Wr. (Specifically, DLk and DTw measure the change of DNA linking number and twist, respectively, from an undeformed double helical state.) For our simulations, DTw is calculated as

DTw ¼

N 1 X qi ; 2p i ¼ 2

(20)

where the summation is over all interior vertices except that representing nicked DNA. Wr for the simulated conformations may be calculated by using method 1a of Klenin and Langowski (65), for example, on a virtually closed chain similar to that used in Vologodskii and Marko (66). We restrict our simulations to one moderate value of DLk because equilibration is computationally expensive. Because our current model lacks chirality, symmetry dictates that similar results would be obtained for DLk ¼ þ86. In our simulations, nicking endonuclease activity is represented by eliminating the torsional potential at one randomly selected vertex of the dWLC for the duration of the simulation. Experiments suggest that base-stacking interactions at a nick support a small degree of torsional stress (see, for example, Shore and Baldwin (67), Zhang and Crothers (68), and Laurens et al. (69)), while permitting full rotation about the nicked strand. Unfortunately, the torsional mechanics at a nick are poorly characterized. Consequently, in our model we neglect any torsional potential associated with base stacking to facilitate simulation. All other potentials remain unchanged at the nicked vertex in our model because nicked DNA retains much of its bending stiffness; see, for example, Mills et al. (70), Demurtas et al. (71), Qu et al. (72), and Cong et al. (73). The choice of discretization length scale, or segment length l0, largely determines the resolution of our simulations and significantly impacts the computational cost. Only DNA structures with length scales larger than ~l0 can be represented accurately with our model. Unfortunately, decreasing l0 to improve accuracy requires an increase in computational effort with more vertices needed to represent a fixed length of DNA. In addition, reducing l0 increases stiffness at the segment/vertex level requiring more time steps to account for the corresponding faster dynamics. For example, the segment stretching stiffness depends on l0 and largely dictates the maximum permitted time-step size because DNA is stiff in extension. To balance tradeoffs between accuracy and computational efficiency,

we take l0 ¼ 10 nm and d ¼ 0.1 (consistent with other simulation studies (51,55)). While this value for d corresponds to a softer stretching stiffness than physically expected for DNA, it allows for larger time steps with little expected impact on supercoil dynamics; see, for example, Klenin et al. (51) and Jian et al. (55). Using l0 ¼ 10 nm, we expect our model to reasonably predict relaxation timescales for the system despite having limitations on representing the mechanics of sharply bent DNA primarily localized to plectoneme end loops. Other simulation parameters are summarized in Table 1.

RESULTS AND DISCUSSION Extension is a poor dynamic measure of DLk relaxation In Fig. 2, we show extension, zNþ1 ðtÞ=L, and DLkðtÞ as functions of time for each of Ns ¼ 100 simulations following a nicking event at time t ¼ 0 ms. In addition, we show the average of these quantities over all simulations, hzNþ1 ðtÞ=Li and hDLkðtÞi. These simulations were performed with Fs ¼ 1.6 pN. The extensional response is a poor means of characterizing DLk relaxation dynamics, as evidenced by the disparate timescales in Fig. 2. Interestingly, it appears that DLk relaxation does not require a significant change in extension as it relaxes much faster than extension. To quantify the difference in relaxation rates, we fit an exponential function (of the form c1 et=t þ c2 , where t is a time constant) to the average extension and average DLk (black curves in Fig. 2). The time constants associated with extensional response (t hz=Li ¼ 7:8 ms) and DLk relaxation (t hLki ¼ 0:78 ms) differ by approximately an order of magnitude. Single-molecule experiments often rely upon the extensional response (as measured by the displacement of a micron scale magnetic bead) to quantify DLk relaxation; TABLE 1

Table of Simulation Parameters

Parameter DNA contour length Absolute temperature Segment equilibrium length Stretching stiffness parameter Bending stiffness parameter Torsional stiffness Monovalent salt concentration Dielectric constant Effective charge density (75) Inverse Debye length Point charges per segment Electrostatic cutoff distance Excluded volume force Viscosity Bead hydrodynamic radius Rotational diffusion radius Surface force length scale Surface force Simulation timestep size Number of simulations Applied tension

Symbol

Value

L T l0 d ab C

21 kilo-basepairs (7 mm) 298 K 10 nm 0.1 2.41 3  1028 N-m2 150 mM 8.9  109 F/m 8.01 e/nm 1.261 nm1 5 8.72 nm 5 pN 0.001 kg/s-m 2.24 nm 1.2 nm 1 nm 5 pN 0.5 ns 100 1.2–3.2 pN

Dw n k l rc Fex h a rrd Lb Fb Dt Ns Fs

Biophysical Journal 110, 2176–2184, May 24, 2016 2179

Ivenso and Lillian

where N above is the number of vertices corresponding to the DNA above the largest plectoneme and amag is the radius of the magnetic bead. The spring stiffness in our lumped model ðklump Þ is taken as k lump ¼

Fs ; Dz

(22)

where Dz is the average change in extension between the initial state and the equilibrium nicked state. Consequently the time constant for this lumped system is t lump z=L ¼ FIGURE 2 (Top) Extensional response as a function of time, after a nicking event at t ¼ 0 ms. All Ns ¼ 100 simulations are depicted in gray. The average extension is depicted in black. The average initial end-to-end extension of the supercoiled DNA is ~80% of L. (Bottom) DLk as a function of time for all simulations in gray and average DLk in black.

see, for example, Seol et al. (19) and Crut et al. (21). Our results suggest that relying on the extensional response as a means of estimating DLk relaxation rates could introduce significant errors. Furthermore, the difference in timescales between extension and DLk suggests that the location at which a nicking endonuclease acts on the DNA has very little effect on the extensional response timescales ðt z=L Þ. Specifically, DLk relaxes so fast that it cannot have much impact on extension. From the experimental data presented in Crut et al. (21), we estimate the relaxation time for extension to be ~ t exp z=L z50 ms for an applied tension of 1.4 pN and initial DLk ¼ 100. (This timescale was estimated by reading the time it takes for the change in extension to reach ~63% of its total change to equilibrium in Fig. 5b of Crut et al. (21).) This is much longer than t hz=Li ¼ 7:8 ms, obtained from our simulations carried out at an applied tension of 1.6 pN and initial DLk ¼ 86. This difference can be explained, in part, by the lack of a large bead in our simulations. Incorporating the hydrodynamic drag of a large magnetic bead in our simulations must slow extension. (However, relaxation of DLk is expected to be insensitive to the inclusion of this drag.) To estimate the change in t hz=Li due to a large magnetic bead, we present a lumped (single degree of freedom, bead-spring) model for the system illustrated in Fig. 1. The hydrodynamic drag of the bead in our lumped model ðglump Þ depends upon the drag of both the magnetic bead and the DNA. The contribution from DNA is based on the observation that DNA above the largest initial plectoneme displaces the most in the vertical direction upon nicking. Specifically, we write    glump ¼ 6ph N above  1 a þ amag ; (21)

2180 Biophysical Journal 110, 2176–2184, May 24, 2016

glump : klump

(23)

In our simulations, the average number of beads above the largest plectoneme is ~300 ðN above z300Þ and the average change in extension is ~1 mm (Dzz1 mm). When the magnetic bead is omitted (i.e., replacing amag with a), the time constant for this lumped model, t lump z=L ¼ 8 ms, closely matches the relaxation time observed in simulations, t hz=Li ¼ 7:8 ms. For a bead of radius amag ¼ 1:5 mm, the relaxation time increases substantially to t lump z=L ¼ 26 ms and becomes comparable to the experimental observation, t exp z=L z50 ms. Now, we extend our simple lumped model to estimate the influence of DNA-DNA and DNA-magnetic bead hydrodynamic interactions on t hz=Li . (We remind the reader that our simulations neglect hydrodynamic interactions due to their significant computational expense.) To approximate hydrodynamic interactions in our lumped model, we assume a fraction of the DNA moves together with the magnetic bead, as a rigid body, and remains in a nearly undeformed state. Therefore, we replace Eq. 21 with  1 glump ¼ kB TET Dlump E: (24) Here, ET ¼ f0; 0; 1; .; 0; 0; 1g and Dlump is the RotnePrager-Yamakawa diffusion matrix (see, for example, Zuk et al. (74)) for a system of N above beads arranged along the z axis. The first ðN above  1Þ beads represent DNA with a radius a while the last represents the magnetic bead with radius amag . Accounting for hydrodynamic interaction in this lumped model yields t lump z=L ¼ 2 ms without a magnetic bead ðamag ¼ aÞ and t lump z=L ¼ 18 ms with a magnetic bead (amag ¼ 1:5 mm). This calculation suggests that accounting for hydrodynamic interactions in our simulations would yield dramatically reduced relaxation times for extension. However, the timescale for extension is expected to remain substantially larger than that for linking number relaxation. (In the Supporting Material, we simulate a much shorter length of DNA, for which we can reasonably compute hydrodynamic interactions, and find results in agreement with the lumped model presented here.)

Simulation of DNA Supercoil Relaxation

The value tLk approaches a limit with Increasing Fs Here we investigate the effect of changing the magnitude of Fs on t hLki . As before, we obtain t hLki by fitting an exponential function to hDLkðtÞi for a range of Fs. Results are presented in Fig. 3. In Fig. 3, we observe that t hLki decreases and appears to approach a limit of ~0.5 ms with increasing tension (within the range reported). Here, t hLki is insensitive to changes in tension for Fs R2:4 pN. In this high tension regime, plectonemic supercoils are very rare, and therefore DLk is stored almost entirely as DTw. A simple model in which DNA is assumed a straight torsionally elastic cylinder of length L, radius rrd, and stiffness C can explain the high tension limit for DLk relaxation, t limit hLki . In this simple model, nicking at a random site along the length yields two segments of DNA: the average length of the longer of the two is 3L=4 while its average squared length is 7L2 =12. Consequently, the relaxation time for the first mode of this segment can be expressed as t limit hLki ¼

16hðr rd Þ Cp

27 2 L 12

:

(25)

This evaluates to t limit hLki ¼ 0:48 ms using the simulation parameters. Within the lower tension regime, DLk is distributed between DTw and Wr. Because our current model maintains bending stiffness at nicked sites, all DLk stored in Wr must be converted to DTw before relaxing. This extra step in the relaxation process requires t hLki to be larger in the lower tension regime. We believe the conversion rate of Wr to DTw, and consequently t hLki, depend upon the tension. For this section, equilibration simulations at forces other than Fs ¼ 1.6 pN are substantially shorter to alleviate the

FIGURE 3 The dependence of DLk relaxation rates on applied tension, Fs. For Fs R 2.4 pN, DLk is stored almost entirely as DTw (extended), while for Fs % 2.4 pN, DLk is distributed between DTw and Wr (plectonemic).

computation burden. Note that the position and structure of plectonemes evolve relatively slowly during an equilibration simulation. While this may have a quantitative impact on the results of Fig. 3, we believe these results to be qualitatively insensitive to the initial structure and position of plectonemes. In contrast, we will show next that the timescale for extension is sensitive to the initial plectonemic state. Extensional response strongly depends upon initial plectoneme position Here we investigate the dependence of relaxation time on the initial position of the largest plectoneme. To do so, we divide runs into three groups based on the position, along the contour length, of the largest plectoneme at time t ¼ 0 ms. (One simulation began with two plectonemes of equal size, for which we use the midpoint between the plectonemes as the position of the largest plectoneme.) The first group includes runs in which the largest plectoneme is in the lower-third while the second and third groups correspond to runs in which the largest plectonemes are in the middle- and upper-thirds. Fig. 4 shows hzNþ1 ðtÞ=Lii and hDLkðtÞii as functions of time for each of these groups; here we use the subscript i to index through the three groups. While plectoneme position has very little observable effect on the t hLki , it has a very marked effect on t hz=Li . As shown in Fig. 4 (top), the nearer the plectoneme to the upper-end of the DNA, the faster the extensional response. Our explanation for this is that the position of the dominant plectoneme determines the length of the DNA contour that

FIGURE 4 The effect of initial plectoneme position on relaxation is observed by dividing simulation runs into three groups based on the position, along the contour length, of the largest plectoneme at time t ¼ 0 ms. The first group (i ¼ 1) includes runs in which the largest plectoneme is in the lower-third (24 simulations) while the second and third groups correspond to runs in which the largest plectonemes are in the middleand upper-thirds (46 and 30 simulations respectively). (Top) Average extension as a function of time for each group, hzNþ1 ðtÞ=Lii . (Bottom) Average DLk as a function of time for each group, hDLkðtÞii .

Biophysical Journal 110, 2176–2184, May 24, 2016 2181

Ivenso and Lillian

Although in this article we simulate the response of supercoiled DNA to a single strand nick induced by a nicking enzyme, in the future we plan to simulate supercoil relaxation in response to various topoisomerases. Unlike nicking endonucleases, the type IB topoisomerases, for example, continue to interact with the swiveling DNA after a nicking event. These interactions include re-ligation that may occur before full relaxation, and friction that slows the relaxation process. In addition, we hope that incorporating new algorithms and faster computational resources will make it feasible to represent hydrodynamic interactions for the lengths of DNA we consider here. SUPPORTING MATERIAL FIGURE 5 Dependence of extensional timescale thz=Li on initial plectoneme position. A linear fit to the data is provided.

is displaced significantly (in the vertical direction) following a nicking event. Specifically, all of the DNA above the plectoneme will be displaced upward by a distance roughly equal to the length of DNA within the plectoneme, while the DNA below will experience very little displacement. Consequently, t hz=Li decreases with increasing distance of the plectoneme from the lower (stationary) boundary. Even the fastest t hz=Li3 (where t hz=Li3 ¼ 4.6 ms) is significantly longer than t hLki (where t hLki ¼ 0.78 ms). Fig. 5 presents a scatterplot of t z=L versus position of the largest plectoneme for all Ns ¼ 100 simulations. Here, the values t z=L for each simulation were obtained by extracting the first time at which the extension makes 63% of its total change to thermal equilibrium. A linear fit to the data suggests that t z=L decreases from ~15 ms near the stationary boundary to 2 ms at the movable boundary: a difference of a factor of almost 8. We are unaware of any single molecule experiments that have detected this trend, perhaps because plectoneme position is undetectable in many experimental systems and this effect may be obscured by the slow dynamics of a large magnetic bead.

Supporting Materials and Methods and one figure are available at http:// www.biophysj.org/biophysj/supplemental/S0006-3495(16)30180-1.

AUTHOR CONTRIBUTIONS I.D.I. and T.D.L. designed the research, performed the research, analyzed the data, and wrote the article.

ACKNOWLEDGMENTS The authors thank Pawan Selokar, David Bell, and Justin Polk for their helpful contributions to this research. The authors acknowledge the High Performance Computing Center at Texas Tech University at Lubbock for providing high performance computing resources that have contributed to the research results reported within this article. This research was funded, in part, from startup funds provided to T.D.L. by Texas Tech University.

REFERENCES 1. Koster, D. A., A. Crut, ., N. H. Dekker. 2010. Cellular strategies for regulating DNA supercoiling: a single-molecule perspective. Cell. 142:519–530. 2. Champoux, J. J. 2001. DNA topoisomerases: structure, function, and mechanism. Annu. Rev. Biochem. 70:369–413. 3. Wang, J. C. 1987. Recent studies of DNA topoisomerases. Biochim. Biophys. Acta. 909:1–9.

CONCLUSIONS In this article, we employ a discrete wormlike chain model to represent the relaxation dynamics of supercoiled DNA in recent single-molecule experiments. Our simulations resolve supercoil dynamics following a single strand nick with detail at levels unattainable in recent experiments. We find that the timescale of the extensional relaxation is much larger than that of DLk. Consequently, extension is a poor dynamic measure of DLk, although it is routinely used as such in experiments. Our simulations show that the rate of relaxation of DLk is sensitive to tension and approaches a limit as tension increases. Finally, we show that the relaxation rate for extension depends strongly on the initial position of the largest plectoneme.

2182 Biophysical Journal 110, 2176–2184, May 24, 2016

4. Liu, L. F., and J. C. Wang. 1987. Supercoiling of the DNA template during transcription. Proc. Natl. Acad. Sci. USA. 84:7024–7027. 5. ten Heggeler-Bordier, B., W. Wahli, ., J. Dubochet. 1992. The apical localization of transcribing RNA polymerases on supercoiled DNA prevents their rotation around the template. EMBO J. 11:667–672. 6. Teves, S. S., and S. Henikoff. 2014. DNA torsion as a feedback mediator of transcription and chromatin dynamics. Nucleus. 5:211–218. 7. Baranello, L., D. Levens, ., F. Kouzine. 2012. The importance of being supercoiled: how DNA mechanics regulate dynamic processes. Biochim. Biophys. Acta. 1819:632–638. 8. Pruss, G. J., and K. Drlica. 1989. DNA supercoiling and prokaryotic transcription. Cell. 56:521–523. 9. Ma, J., L. Bai, and M. D. Wang. 2013. Transcription under torsion. Science. 340:1580–1583. 10. Tsao, Y.-P., H.-Y. Wu, and L. F. Liu. 1989. Transcription-driven supercoiling of DNA: direct biochemical evidence from in vitro studies. Cell. 56:111–118.

Simulation of DNA Supercoil Relaxation 11. Samul, R., and F. Leng. 2007. Transcription-coupled hypernegative supercoiling of plasmid DNA by T7 RNA polymerase in Escherichia coli topoisomerase I-deficient strains. J. Mol. Biol. 374:925–935.

34. Curuksu, J., M. Zacharias, ., K. Zakrzewska. 2009. Local and global effects of strong DNA bending induced during molecular dynamics simulations. Nucleic Acids Res. 37:3766–3773.

12. Zechiedrich, E. L., A. B. Khodursky, ., N. R. Cozzarelli. 2000. Roles of topoisomerases in maintaining steady-state DNA supercoiling in Escherichia coli. J. Biol. Chem. 275:8103–8113.

35. Harris, S. A., C. A. Laughton, and T. B. Liverpool. 2008. Mapping the phase diagram of the writhe of DNA nanocircles using atomistic molecular dynamics simulations. Nucleic Acids Res. 36:21–29.

13. Wang, J. C. 1996. DNA topoisomerases. Annu. Rev. Biochem. 65:635–692.

36. Hirsh, A. D., M. Taranova, ., N. C. Perkins. 2013. Structural ensemble and dynamics of toroidal-like DNA shapes in bacteriophage f29 exit cavity. Biophys. J. 104:2058–2067.

14. Schoeffler, A. J., and J. M. Berger. 2008. DNA topoisomerases: harnessing and constraining energy to govern chromosome topology. Q. Rev. Biophys. 41:41–101. 15. Koster, D. A., K. Palle, ., N. H. Dekker. 2007. Antitumour drugs impede DNA uncoiling by topoisomerase I. Nature. 448:213–217. 16. Koster, D. A., F. Czerwinski, ., N. H. Dekker. 2008. Single-molecule observations of topotecan-mediated TopIB activity at a unique DNA sequence. Nucleic Acids Res. 36:2301–2310. 17. Pommier, Y., E. Leo, ., C. Marchand. 2010. DNA topoisomerases and their poisoning by anticancer and antibacterial drugs. Chem. Biol. 17:421–433. 18. Koster, D. A., V. Croquette, ., N. H. Dekker. 2005. Friction and torque govern the relaxation of DNA supercoils by eukaryotic topoisomerase IB. Nature. 434:671–674. 19. Seol, Y., H. Zhang, ., K. C. Neuman. 2012. A kinetic clutch governs religation by type IB topoisomerases and determines camptothecin sensitivity. Proc. Natl. Acad. Sci. USA. 109:16125–16130. 20. Bai, H., J. E. Kath, ., J. F. Marko. 2012. Remote control of DNAacting enzymes by varying the Brownian dynamics of a distant DNA end. Proc. Natl. Acad. Sci. USA. 109:16546–16551. 21. Crut, A., D. A. Koster, ., N. H. Dekker. 2007. Fast dynamics of supercoiled DNA revealed by single-molecule experiments. Proc. Natl. Acad. Sci. USA. 104:11957–11962. 22. Goshen, E., W. Z. Zhao, ., M. Feingold. 2005. Relaxation dynamics of a single DNA molecule. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71:061920. 23. Daniels, B. C., and J. P. Sethna. 2011. Nucleation at the DNA supercoiling transition. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 83:041924. 24. Vologodskii, A. 2009. Determining protein-induced DNA bending in force-extension experiments: theoretical analysis. Biophys. J. 96:3591–3599. 25. Marko, J. F. 2007. Torque and dynamics of linking number relaxation in stretched supercoiled DNA. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76:021926. 26. Daniels, B. C., S. Forth, ., J. P. Sethna. 2009. Discontinuities at the DNA supercoiling transition. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 80:040901. 27. Purohit, P. K. 2008. Plectoneme formation in twisted fluctuating rods. J. Mech. Phys. Solids. 56:1715–1729. 28. Emanuel, M., G. Lanzani, and H. Schiessel. 2013. Multiplectoneme phase of double-stranded DNA under torsion. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 88:022706. 29. Moroz, J. D., and P. Nelson. 1997. Torsional directed walks, entropic elasticity, and DNA twist stiffness. Proc. Natl. Acad. Sci. USA. 94:14418–14422. 30. Schurr, J. M., and B. S. Fujimoto. 2013. Supercoiled pseudocircular domains in single-twisted DNAs under tension: elastic constants and unwinding dynamics in complexes with Topo I. Biopolymers. 99:1046–1069.

37. Goyal, S., T. Lillian, ., N. C. Perkins. 2007. Intrinsic curvature of DNA influences LacR-mediated looping. Biophys. J. 93:4342–4359. 38. Wedemann, G., C. Mu¨nkel, ., J. Langowski. 1998. Kinetics of structural changes in superhelical DNA. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics. 58:3537–3546. 39. Klenin, K. V., and J. Langowski. 2001. Diffusion-controlled intrachain reactions of supercoiled DNA: Brownian dynamics simulations. Biophys. J. 80:69–74. 40. Swigon, D., B. D. Coleman, and W. K. Olson. 2006. Modeling the Lac repressor-operator assembly: the influence of DNA looping on Lac repressor conformation. Proc. Natl. Acad. Sci. USA. 103:9879–9884. 41. Yang, Y., I. Tobias, and W. K. Olson. 1993. Finite-element analysis of DNA supercoiling. J. Chem. Phys. 98:1673–1686. 42. Sivak, D. A., and P. L. Geissler. 2012. Consequences of local interstrand dehybridization for large-amplitude bending fluctuations of double-stranded DNA. J. Chem. Phys. 136:045102. 43. Du, Q., C. Smith, ., A. Vologodskii. 2005. Cyclization of short DNA fragments and bending fluctuations of the double helix. Proc. Natl. Acad. Sci. USA. 102:5397–5402. 44. Zheng, X., and A. Vologodskii. 2009. Theoretical analysis of disruptions in DNA minicircles. Biophys. J. 96:1341–1349. 45. Scho¨pflin, R., H. Brutzer, ., G. Wedemann. 2012. Probing the elasticity of DNA on short length scales by modeling supercoiling under tension. Biophys. J. 103:323–330. 46. Villa, E., A. Balaeff, ., K. Schulten. 2004. Multiscale method for simulating protein-DNA complexes. Multiscale Model. Simul. 2:527–553. 47. Lillian, T. D., M. Taranova, ., N. C. Perkins. 2011. A multiscale dynamic model of DNA supercoil relaxation by topoisomerase IB. Biophys. J. 100:2016–2023. 48. Mielke, S. P., W. H. Fink, ., C. J. Benham. 2004. Transcription-driven twin supercoiling of a DNA loop: a Brownian dynamics study. J. Chem. Phys. 121:8104–8112. 49. Swigon, D., and W. K. Olson. 2008. Mesoscale modeling of multi-protein-DNA assemblies: the role of the catabolic activator protein in Lacrepressor-mediated looping. Int. J. Non-linear Mech. 43:1082–1093. 50. Wereszczynski, J., and I. Andricioaei. 2010. Free energy calculations reveal rotating-ratchet mechanism for DNA supercoil relaxation by topoisomerase IB and its inhibition. Biophys. J. 99:869–878. 51. Klenin, K., H. Merlitz, and J. Langowski. 1998. A Brownian dynamics program for the simulation of linear and circular DNA and other wormlike chain polyelectrolytes. Biophys. J. 74:780–788. 52. Chirico, G., and J. Langowski. 1996. Brownian dynamics simulations of supercoiled DNA with bent sequences. Biophys. J. 71:955–971. 53. Chirico, G., and J. Langowski. 1994. Kinetics of DNA supercoiling studied by Brownian dynamics simulation. Biopolymers. 34:415–433.

31. Taranova, M., A. D. Hirsh, ., I. Andricioaei. 2014. Role of microscopic flexibility in tightly curved DNA. J. Phys. Chem. B. 118:11028–11036.

54. Allison, S. A. 1986. Brownian dynamics simulation of wormlike chains. Fluorescence depolarization and depolarized light scattering. Macromolecules. 19:118–124.

32. Lankas, F., R. Lavery, and J. H. Maddocks. 2006. Kinking occurs during molecular dynamics simulations of small DNA minicircles. Structure. 14:1527–1534.

55. Jian, H. M., A. V. Vologodskii, and T. Schlick. 1997. Combined wormlike-chain and bead model for dynamic simulations of long linear DNA. J. Comput. Phys. 136:168–179.

33. Mitchell, J. S., C. A. Laughton, and S. A. Harris. 2011. Atomistic simulations reveal bubbles, kinks and wrinkles in supercoiled DNA. Nucleic Acids Res. 39:3928–3938.

56. Jian, H., T. Schlick, and A. Vologodskii. 1998. Internal motion of supercoiled DNA: Brownian dynamics simulations of site juxtaposition. J. Mol. Biol. 284:287–296.

Biophysical Journal 110, 2176–2184, May 24, 2016 2183

Ivenso and Lillian 57. Wocjan, T., K. Klenin, and J. Langowski. 2009. Brownian dynamics simulation of DNA unrolling from the nucleosome. J. Phys. Chem. B. 113:2639–2646.

67. Shore, D., and R. L. Baldwin. 1983. Energetics of DNA twisting. I. Relation between twist and cyclization probability. J. Mol. Biol. 170:957–981.

58. Huang, J., T. Schlick, and A. Vologodskii. 2001. Dynamics of site juxtaposition in supercoiled DNA. Proc. Natl. Acad. Sci. USA. 98:968–973. 59. Vologodskii, A. 2006. Brownian dynamics simulation of knot diffusion along a stretched DNA molecule. Biophys. J. 90:1594–1597. 60. Ivenso, I., and T. D. Lillian. 2014. Brownian dynamics simulation of the dynamics of stretched DNA. American Society of Mechanical Engineers, IDETC/CIE 2014, Buffalo, NY. DETC2014–35487.

68. Zhang, Y., and D. M. Crothers. 2003. High-throughput approach for detection of DNA bending and flexibility based on cyclization. Proc. Natl. Acad. Sci. USA. 100:3161–3166.

61. Ikenna, I., and T. D. Lillian. 2015. The dynamics of DNA supercoiling: a Brownian dynamics study. American Society of Mechanical Engineers, IDETC/CIE 2015, Boston, MA. DETC2015–47444. 62. Lillian, T. D., I. Ivenso, ., J. Polk. 2015. Supercoil diffusion along stretched DNA by Brownian dynamics. Biophys. J. 108:231a. 63. Iniesta, A., and J. G. de la Torre. 1990. A second-order algorithm for the simulation of the Brownian dynamics of macromolecular models. J. Chem. Phys. 92:2015–2018. 64. van Loenhout, M. T. J., M. V. de Grunt, and C. Dekker. 2012. Dynamics of DNA supercoils. Science. 338:94–97.

69. Laurens, N., D. A. Rusling, ., G. J. L. Wuite. 2012. DNA looping by FokI: the impact of twisting and bending rigidity on protein-induced looping dynamics. Nucleic Acids Res. 40:4988–4997. 70. Mills, J. B., J. P. Cooper, and P. J. Hagerman. 1994. Electrophoretic evidence that single-stranded regions of one or more nucleotides dramatically increase the flexibility of DNA. Biochemistry. 33:1797–1803. 71. Demurtas, D., A. Amzallag, ., A. Stasiak. 2009. Bending modes of DNA directly addressed by cryo-electron microscopy of DNA minicircles. Nucleic Acids Res. 37:2882–2893. 72. Qu, H., C. Y. Tseng, ., G. Zocchi. 2010. The elastic energy of sharply bent nicked DNA. Europhys. Lett. 90:18003. 73. Cong, P., L. Dai, ., J. Yan. 2015. Revisiting the anomalous bending elasticity of sharply bent DNA. Biophys. J. 109:2338–2351.

65. Klenin, K., and J. Langowski. 2000. Computation of writhe in modeling of supercoiled DNA. Biopolymers. 54:307–317.

74. Zuk, P. J., E. Wajnryb, ., P. Szymczak. 2014. Rotne-Prager-Yamakawa approximation for different-sized particles in application to macromolecular bead models. J. Fluid Mech. 741:R5.

66. Vologodskii, A. V., and J. F. Marko. 1997. Extension of torsionally stressed DNA by external force. Biophys. J. 73:123–132.

75. Vologodskii, A., and N. Cozzarelli. 1995. Modeling of long-range electrostatic interactions in DNA. Biopolymers. 35:289–296.

2184 Biophysical Journal 110, 2176–2184, May 24, 2016