Stretching of long double-stranded DNA and RNA described by the same model

We propose a bead-spring model that accurately reproduces a variety of experimental force-extension curves of long double-stranded DNA and RNA, including torsionally constrained and unconstrained DNA, and negatively supercoiled DNA. A key feature of the model is a specific non-convex energy function of the spring. We provide an algorithm for obtaining five required parameters of the model from experimental force-extension curves. In the plateau region of the force-extension curves, our molecular dynamics simulations show that the polymer separates into a mix of weakly and strongly stretched states without forming macroscopically distinct phases.


Introduction
Double-stranded (ds) DNA and RNA molecules are constantly under different loading conditions in living cells. [1][2][3] Knowledge of the mechanical properties of these polymers is of great importance for understanding many key cellular processes. [4][5][6][7] A widely used experimental technique employed to explore the mechanics of polymers is the single-molecule stretching with optical or magnetic tweezers. 5,[8][9][10] To interpret these experiments under low tension, several variants of the worm-like chain model (WLC) have gained popularity. 8,[11][12][13][14][15][16] The popularity of the models stems, in large part, from their conceptual simplicity, the clear physics behind them, the direct connection with the key mechanical parameters of the polymers, and the facility for extracting these properties from the experimental data. For example, according to the extensible worm-like chain (eWLC) model, 8,17,18 at a given stretching force F, the extension L x = L 0 (1 + F/S − (k B T /4F P ) 1/2 ), where k B is the Boltzmann constant, and T is the absolute temperature.
Thus, the parameters of the polymer, such as the polymer persistence length P, the stretch modulus S, and its contour length L 0 (its maximum end-to-end distance), can be obtained by fitting the model to the experimental force-extension curve. 8,19 However, when it comes to dsDNA and dsRNA under high tension, these popular models have a serious limitation. Specifically, the single-molecule experiments reveal a peculiar plateau region on force-extension curves for dsDNA 8,[20][21][22] and dsRNA. 8,23 Only the stretching behavior of the molecules below the plateau is well described by the eWLC model, 8,17,18 but the plateau itself remains out of reach of this otherwise appealing model (see e.g. Figure 5 of the "Results" section).
To describe full force-extension curves of dsDNA, including the plateau, a variety of theoretical models have been proposed. [24][25][26][27] These models are based on statistical mechanics, and describe a number of possible states with parameters fitted to experimental data. Despite the clear advantage of being more general than the variants of the WLC, these more complex models have not been used nearly as widely in interpreting single-molecule experiments, possibly because the connection to macroscopic parameters that describe properties of the polymers is no longer as intuitive as in the variants of the WLC.
In general, these models rely on empirical potential functions (force-fields) that describe interactions between atoms or beads through a combination of bond stretching, bending, torsional, and electrostatic potentials, to name a few, where each potential has multiple fitted parameters, making the total number of parameters as many a hundred for some general-purpose biomolecular force-fields. 47,48 An advantage of these models, compared to fully phenomenological models such as the variants of the WLC, is that they can provide a microscopic, including atomic-level, interpretation of experiment. The drawback is that direct connection to the macroscopic parameters that describe polymer deformation is harder to establish, and the simulation outcomes may depend on the force-field chosen. These complex force-fields are fine-tuned to work with all of their components present, making it difficult to ask questions that require turning various physical interactions on or off. For the specific case of the DNA and RNA, several of these models have enough complexity built into them to reproduce, to various degrees, the overstretching plateau region, which is a clear advantage over the variants of the WLC. A number of simulation studies of dsDNA and dsRNA stretching at this level of description has been performed, 49-54 but still we are unaware of a mesoscopic or a fully atomistic model that accurately reproduces various regions of experimental force-extension curves of both long dsDNA and dsRNA within the same model.
We suggest a mesoscopic model that may be considered a hybrid between the two approaches, with just enough complexity to faithfully account for the plateau region and its basic physics. We aim to retain the appealing simplicity level of the variants of the WLC with its direct connection to mechanical properties of the polymer, while also retaining the most basic feature of the molecular simulation approach -microscopic granularity. Our longterm motivation is to construct a model that may potentially replace or at least complement the variants of the WLC in the routine interpretation of single-molecule experiments. 3 We propose a mesoscopic model, based on a non-convex potential energy function, that reproduces, closely, the experimental stretching behavior of long dsDNA and dsRNA with one plateau region on their force-extension curves. We use the model to explore the possibility of a phase separation in the plateau region.
The conceptual basis of the proposed model is a recently proposed general framework, 49,55 in which the existence of the plateau region on the force-deformation curve stems from a non-convex profile of the energy of the system. Namely, if the energy as a function of deformation is not convex over some interval, then the system can minimize its energy by separating into two states: weakly and strongly deformed, with the geometric characteristics of the states corresponding to the beginning and the end of the convex hull of the energy function, respectively. 49,55 Our proposed model treats the dsDNA or dsRNA molecule as a chain of N beads connected by non-linear springs, Figure 1A. The stretching behavior of the non-linear spring is governed by the potential, U (r), in the form of a piecewise continuously differentiable function: where S is the stretch modulus of the molecule (in k B T /Å); r, the distance between the centers of neighboring beads -the bond length (inÅ); r 0 , the equilibrium bond length (in A); r a , the bond length at the beginning of the convex hull (inÅ); r b , the bond length at the end of the convex hull (inÅ); ∆, the energy difference between the potential U (in k B T ) and the convex hull at r=(r a +r b )/2. Equivalently, ∆ can be interpreted as the potential energy barrier height between two equally likely states of the spring -weekly and strongly stretched 4 -at an applied load along the r-direction equal to the plateau force of F (r a ) = U (r a ), as depicted in Figure 1B. To account for the bending of the chain, characterized by its persistence length P , we employ the discrete version of the WLC model, see "Methods".
We propose an algorithm for extracting parameters of the new model, P , S, r a , r b , and ∆, from an experimental force-extension curve. The algorithm is based on approximating the force-extension curve by a piece-wise function F L x L 0 , below. This function, the blue curve in Figure 2, describes the minimum energy state of the polymer chain with its fixed ends, governed by the stretching potential U (r).
where F s is the overstretching force at the middle of the plateau (in pN), S is the stretch  We start exploring the model using coarse-grained Molecular Dynamics (MD) simulations, see "Methods", by considering the case of torsionally unconstrained dsDNA in detail, Figure 3. The governing stretching potential U (r) is shown in Figure 3A for two different values of the parameter ∆; the other parameters of the model are obtained according to the algorithm. As illustrated in Figure 3B, the model clearly reproduces all the key features of the experimental force-extension diagram. The parameter ∆ affects the transition of the chain to the plateau region, where the stretching of the chain is bi-modal, Figure 3B. The effect of larger ∆ is a sharper transition to the plateau region in the force-extension diagram, as well as sharper peaks in the probability density ρ of relative bond lengths r/r 0 , as shown in Figure 3B.  Figure 3A; the eWLC fit in the force range of 2-50 pN). Inset: the probability density ρ of relative bond lengths r/r 0 for the model of torsionally unconstrained dsDNA with ∆ = 2 k B T (pink) and ∆ = 8 k B T (blue); L x /L 0 ≈ 1.388.
In the plateau region, the polymer chain is a mix of two states -strongly and weekly stretched, Figure 4. Figure 4: Conformations of the chain at t=0 ns and t=700 ns for torsionally unconstrained dsDNA with ∆ = 8 k B T ; L x /L 0 ≈ 1.388 . Case 1: the simulation begins from the "onephase" conformation, in which the beads are equally spaced with r/r 0 = 1.388. Case 2: the simulation begins from the "two-phase" conformation, in which the first 50 % of the beads are equally spaced with r/r 0 = 1.056 and the other 50 % of the beads are equally spaced with r/r 0 = 1.72. Case 3 (the Debye-Huckel potential is added): the simulation begins from the "mixed state" conformation, in which the bond lengths alternate. The color scale (from red to blue) indicates the relative bond length between the bead of that color and its neighbor to the right. Snapshots are generated using VMD.
Further MD simulations based on the proposed model with the parameters obtained according to the algorithm, Figure 5, demonstrate that the model leads to a good agreement with the experimental force-extension curves of all tested dsDNA and dsRNA molecules, including, most importantly, the plateau region where the eWLC fails. In the plateau regions, the distributions (the probability density ρ) of relative bond lengths r/r 0 demonstrate the co-existence of the two microscopic states: weakly and strongly stretched, Figure 5. In the plateau region of the force-extension curves, the simulations show no formation of distinct (macroscopic) phases in the polymer, as illustrated in Figure 4 (cases 1-2). As shown in Figure 4 (case 3), we also see no phase separation in our model with the adding of long-range electrostatic interactions described by the Debye-Hückel potential.
Our prediction is that long dsDNA and dsRNA molecules do not exhibit (macroscopically) distinct phases in the plateau region at room temperature. We speculate that single-molecule fluorescence resonance energy transfer, with judiciously placed labels, may be used to test the prediction.

Conclusions
In this work, we have suggested a model that could replace the variants of the WLC in interpreting single-molecule stretching experiment. The model, built on a non-convex potential energy function, is able to reproduce well the force-extension curves of the long DNA and RNA with a single force plateau under low and high tension. The model requires five parameters that can be easily derived from a single force-extension curve using the suggested algorithm. The algorithm is based on an approximation of a force-extension curve by a function describing the minimum-energy state of the polymer chain with its fixed ends.
Although the proposed mesoscopic model has its limitations, including no explicit account for sequence dependence or torsional rigidity, we showed that it (i) accurately reproduces the The total potential energy of the chain, U total , has two terms: The first term is summed over all adjacent pairs of beads, and r i is the distance between the centers of bead i and bead (i + 1). The second term is the harmonic bending energy of the chain (the discrete version of the WLC model), where P is its persistence length (inÅ); φ i , the angle between the vectors from the centers of bead i to bead (i − 1) and bead (i + 1), Figure 1A; the equilibrium bond length r 0 is 11h, where the distance between base pairs h is given in Table 1 below.
In Figure 4 (case 3 only), the Debye-Hückel potential, U DH (r) = l b k B T q 2 exp(−κr)/re 2 , is added as a perturbation to U total to model long-range interactions between charged beads in salt solution for r < r c , where r is the distance between the centers of beads, e is the elementary charge, the Bjerrum length l b is 10Å, the charge of each bead q is −2|e|, the inverse Debye length κ is 0.01Å −1 , and the force cutoff r c is 10r 0 . At t = 0 ns, the contribution of the Debye-Hückel potential to the total potential energy of the chain does not exceed 4 %.
The force in the simulations is measured with a model of the optical trap, represented as a harmonic spring potential, E trap (r 1 ) = (K trap /2)(r 1 − r 2 ) 2 , where the equilibrium length r 2 is 5r 0 , r 1 is the distance between the right fixed end and the center of the rightmost bead of the chain (Figure 1A), and K trap is the stiffness of the optical trap, Table 1. We make the following assumptions: (i) the average mass of a base pair is 650 Da, (ii) the contour length of dsDNA and dsRNA is fixed at L 0 = (N − 1)r 0 .
Each simulation is 700 ns long, and begins from the starting conformation of the chain, case 1, 2 or 3 in Figure 4, in which the force in the optical trap is zero, Figure 1A. We find that 700ns is enough to ensure convergence to the same final state regardless of the initial conformation of the chain, Figure 4. To obtain each point on the force-extension curves shown in Figure 3B and 5D-5F, we select the "one-phase" conformation (case 1 in Figure 4) with r/r 0 = 0.8 as the starting conformation of the chain. The duration of each simulation is composed of the 50 ns equilibration run followed by the production run. During the 650 ns production run, the values of the force in the optical trap and the relative extension L x /L 0 are recorded at each integration time step, and then averaged over 650 ns. To obtain the next point on the force-extension curves, the above protocol is repeated starting from the same starting conformation, modified by a 0.025r 0 increment to each bond length. In the insets of Figure 3B and 5D-5F, the relative bond lengths are recorded every 0.1 ns during the 650 ns production run. Table 1 lists other simulation details.