This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Base pair fluctuations in helical models for nucleic acids

Marco Zoli School of Science and Technology
University of Camerino, I-62032 Camerino, Italy
marco.zoli@unicam.it
Abstract

A statistical method is developed to estimate the maximum amplitude of the base pair fluctuations in a three dimensional mesoscopic model for nucleic acids. The base pair thermal vibrations around the helix diameter are viewed as a Brownian motion for a particle embedded in a stable helical structure. The probability to return to the initial position is computed, as a function of time, by integrating over the particle paths consistent with the physical properties of the model potential. The zero time condition for the first-passage probability defines the constraint to select the integral cutoff for various macroscopic helical conformations, obtained by tuning the twist, the bending and the slide motion between adjacent base pairs along the molecule stack. Applying the method to a short homogeneous chain at room temperature, we obtain meaningful estimates for the maximum fluctuations in the twist conformation with 10.5\sim 10.5 base pairs per helix turn, typical of double stranded DNA helices. Untwisting the double helix, the base pair fluctuations broaden and the integral cutoff grows. The cutoff is found to increase also in the presence of a sliding motion which shortens the helix contour length, a situation peculiar of dsRNA molecules.

pacs:
87.14.gk, 87.15.A-, 87.15.B-, 05.10.-a

I. Introduction

Mesoscopic Hamiltonian models provide a convenient approach for studies of the structural and mechanical properties of nucleic acids as they reduce the number of degrees of freedom required to describe the nucleotide, i.e. the unit comprising the nitrogenous base and the sugar-phosphate group in the molecule backbone pablo07 ; doye11 ; jeon14 ; bish14 ; liang14 . Among the models proposed over the last decades, the Peyrard-Bishop-Dauxois (PBD) model pey93 stands out as it yields a simple and appealing description of the main forces at play in the DNA molecule in terms of a single continuous variable, accounting for the dynamics of the nucleotides on complementary strands through the stretch mode between the base pair mates. Although initially proposed to study the thermally driven DNA denaturation in the thermodynamic limit of an infinite chain, this one-dimensional model with first neighbors interactions has been later applied in different contexts to study DNA thermodynamical and flexibility properties, bubble statistics and dynamics both in short and kilo-base long sequences campa98 ; zdrav06 ; bish09 ; singh09 ; io09 ; io10 ; singh11 ; falo12 ; hando12 ; albu14 ; io14a ; lak19 ; kalos20 . The same model has been also used to estimate the force constants of RNA chains and DNA/RNA hybrids by fitting the melting temperatures of short duplexes weber13 ; weber19 .

While the original PBD model, with onsite Morse potential for hydrogen bonds, predicts base pair lifetimes of open and closed states pey08 much shorter than those inferred from proton–deuterium exchange experiments gueron85 , successive extensions of the model pey09 ; erp18 have added solvent potential terms which i) enhance the threshold for base pair dissociation over the Morse plateau and ii) introduce a re-closing barrier accounting for the hydrogen bonds that the open bases may establish with the solvent, albeit at an higher cost as bases are hydrophobic coll95 . Although such improved model potentials partly reconcile the calculated base pair lifetimes with the experimental estimates, the fact remains that the PBD model is intrinsically 1D whereas a consistent analysis of the thermodynamics and dynamical properties of DNA should not overlook the helical structure of the molecule with its twisting and bending degrees of freedom croq00 ; maiti15 . Furthermore, in a realistic picture for the DNA molecule in solution both strands recombination and large amplitude base pair fluctuations should be allowed but the pair mates should not be free to go infinitely apart. Accordingly, this requires a truncation of the base pair configuration space which is in fact always performed in computational techniques for the PBD models, with or without solvent potentials, e.g. Monte Carlo and molecular dynamics simulations and Transfer Integral methods ares05 ; joy05 . While the truncation removes the divergence of the thermodynamic quantities otherwise encountered for chains with a finite number of base pairs, the cutoff is a measure of the largest separation of a fluctuating base pair identified with the width of the first bound state of the PBD Hamiltonian and generally depends on the Hamiltonian parameters zhang97 . However, even for a specific sequence, the upper bound in the base pair fluctuations is somewhat arbitrary and markedly affects the estimated melting temperature pey95 ; wart85 . Consequently, large discrepancies exist as for the values used in calculations of the partition function and thermodynamical properties of DNA chains with cutoffs ranging from 30Å\sim 30\,\AA up to about ten times the average helix diameter (which is 20Å\sim 20\,\AA) zhang97 ; singh15 . Also, these values hold solely for 1D molecules described by the PBD ladder Hamiltonian whereas the cutoffs may differ significantly in 3D models.

In the finite temperature path integral method for the PBD model, the upper bound in the integration over the base pair paths can be technically determined by the normalization condition for the free particle action io11a ; io14c ; klein04 . Likewise, this condition holds in the path integral for the DNA helical model developed over the last years and applied to compute cyclization, distribution lengths and stretching properties of short chains io16b ; io18c ; io18 ; io18a ; io18b ; io16a . However, for the latter purpose, one may need to take a cutoff larger than the value set by the minimal normalization condition in order to include those large amplitude fluctuations which affect the flexibility of the chain. This points to the importance of defining a rigorous physical criterion which restricts the base pair configuration space selecting a cutoff consistently with the model potential. Here we argue that the cutoff may also vary with the specific helical conformation of the molecule which, in turn, essentially depends on the three variables, i.e. twisting, bending and sliding, describing the motion of any base pair relative to its neighbor in the stack.

Extending to the more realistic and complex three dimensional case an idea previously applied to the PBD ladder model io20 , we notice that the base pair thermal fluctuations are an example of Brownian motion for a particle subjected to the interactions which shape the helix. Following this observation we develop a novel statistical method which sets the integration cutoff independently of the melting temperature of the specific fragment. Taking a single base pair of the chain as the Brownian particle, we compute its time dependent first passage probability making use of the path integral formalism and derive a general benchmark to select the proper cutoff for the base pair fluctuations. As the argument applies in principle to any base pair in the chain, the method permits to truncate the base pairs space configuration and provides a rich relation between the amplitude of the base pair fluctuations and the helical conformation. Further, the method is general enough to be applied to any 3D helical molecule for which a Hamiltonian description is feasible.

The geometrical representation for the helix together with the Hamiltonian model are outlined in Section II while the computational method for the first passage probability is presented in Section III. The results are contained Section IV and some final remarks are made in Section V.

II. Model

We consider the helical model for a linear DNA chain of NN point-like base pairs, with reduced mass μ\mu, first proposed in ref.io16b and depicted in Fig. 1(a). Essentially, the coordinate rir_{i} is the distance between the bases on complementary strands with respect to the mid helical axis and it measures the base pair radial fluctuations. When BB overlaps OiO_{i}, the i-th fluctuation vanishes and the pair mates distance equals the average helix diameter R0R_{0}. Neighboring base pairs along the molecule stack can be twisted and bent. In our 3D model, the angular degrees of freedom are represented by the twist θi\theta_{i} and by the bending ϕi\phi_{i} between adjacent vectors rir_{i} and ri1r_{i-1}. In the absence of bending, the radial fluctuations occur within the ovals and the model reduces to a fixed-planes representation for the helix io11b . Besides twisting and bending, there is a third relative motion in a dimer i.e., the sliding of the base pairs past each other, whose value depends in general on the specific pyrimidine-purine step and also on the base pair propeller twist, the latter however not accounted for by our point-like model. While the interplay of twisting, bending and sliding generates a manifold of helical conformations at the macroscopic scale calla , the slide motion depicted in Fig. 1(b) has the main effect to shorten the rise along the helical axis and reduce the chain contour length, as found e.g., in dsRNA after comparison with dsDNA molecules having the same sequence gonz13 ; dekker14 . In fact, A-form structures (such as dsRNA) tipically display an average absolute value of slide larger than B-form structures although the average slide is largely dependent on the oligomer sequence. For instance, in dsDNA, a sizeable slide is observed at pyrimidine-purine steps and also at GG/CC steps whereas   SS is essentially zero at AA/TT steps calla . Some average structural parameters found for dsRNA and dsDNA are summarized in Table 1. Note that the bending angle is essentially a measure of the roll an­gle from one base pair to the next, defined in a rigid body parametrization for the local geometry of the base pair steps dick89 . Such parametrization also comprises the tilt angles which however contribute much less to the bending as they are generally smaller than the roll angles olson95 .

Refer to caption
Refer to caption
Figure 1: (Color online) (a) Geometric model for a helical chain with NN point-like base pairs. rir_{i} is the inter-strand distance between the mates of the ithi-th base pair. It is defined with respect to the point OiO_{i} lying along the helix mid-axis. The OiO_{i}’s are separated by the rise distance dd. The angles θi\theta_{i} and ϕi\phi_{i} define the local twist and bending between neighboring base pairs. AB¯\overline{AB} is the distance between the radial displacements rir_{i}, ri1r_{i-1}. (b) Neighboring base pairs with relative sliding perpendicular to the helix mid-axis. The upper pair (ith)(i-th) may go either to the left or to the right with respect to the i1i-1 pair. Note that, due to the slide SS, the rise along the helix axis (HSH_{S}) is smaller than dd.

Straightforward geometrical rules permit to derive the distances between adjacent radial displacements, i.e. AB¯\overline{AB} in Fig. 1(a) and AB¯\overline{A^{\prime}B^{\prime}} in Fig. 1(b). For instance AB¯\overline{AB} is given by

di,i1¯=[(d+risinϕi)2+ri12+(ricosϕi)22ri1ricosϕicosθi]1/2,\displaystyle\overline{d_{i,i-1}}=\,\bigl{[}(d+r_{i}\sin\phi_{i})^{2}+r_{i-1}^{2}+(r_{i}\cos\phi_{i})^{2}-2r_{i-1}\cdot r_{i}\cos\phi_{i}\cos\theta_{i}\bigr{]}^{1/2}\,, (1)

where dd is the bare rise distance in the absence of fluctuations, i.e. the bond length for beads arranged along a linear chain in the freely jointed model kovac73 ; busta92 . Likewise di,i1(S)¯\overline{d_{i,i-1}(S)} is the measure of AB¯\overline{A^{\prime}B^{\prime}}. Note that also the twist and bending angles in Fig. 1(b) may differ from those in Fig. 1(a) due to the sliding motion.

dsRNA dsDNA
R0R_{0} 24 20
HSH_{S} 2.8 3.3
S -1.48 0.03
θ¯\bar{\theta} 32o32^{o} 34.6o34.6^{o}
ϕ¯\bar{\phi} 8.6o8.6^{o} 3.5o3.5^{o}
Table 1: Average structural parameters reported for A-form dsRNA and B-form dsDNA helices dekker14 ; olson03 ; wang79 . R0R_{0}, in units Å, is the average helix diameter. HSH_{S} and SS, in units Å, are respectively the rise distance per base pair and the slide as depicted in Fig. 1(b) n3 . θ¯\bar{\theta} and ϕ¯\bar{\phi} are respectively the average twist and bending angles between adjacent base pairs. The rise, slide, twist and bending angle may vary significantly along the specific sequence according to the dinucleotide step oroz10 ; lavery14 .

From Eq. (1), setting to zero both the intrinsic stiffness dd and the angular variables, one recovers the original PBD model. Later works, which have gone beyond the PBD ladder model introducing helicity without bending of the base pair planes, have assumed either finite values for dd and θi\theta_{i} barbi99 or a non vanishing twist θi\theta_{i} with zero dd weber06 .

For the helical linear chain with first neighbors interactions in Fig. 1(a), the Hamiltonian reads:

H=Ha[r1]+i=2NHb[ri,ri1,ϕi,θi],\displaystyle H=\,H_{a}[r_{1}]+\sum_{i=2}^{N}H_{b}[r_{i},r_{i-1},\phi_{i},\theta_{i}]\,,
Ha[r1]=μ2r˙12+V1[r1],\displaystyle H_{a}[r_{1}]=\,\frac{\mu}{2}\dot{r}_{1}^{2}+V_{1}[r_{1}]\,,
Hb[ri,ri1,ϕi,θi]=μ2r˙i2+V1[ri]+V2[ri,ri1,ϕi,θi],\displaystyle H_{b}[r_{i},r_{i-1},\phi_{i},\theta_{i}]=\,\frac{\mu}{2}\dot{r}_{i}^{2}+V_{1}[r_{i}]+V_{2}[r_{i},r_{i-1},\phi_{i},\theta_{i}]\,\,,
V1[ri]=Di[exp(bi(|ri|R0))1]2,\displaystyle V_{1}[r_{i}]=\,D_{i}\bigl{[}\exp(-b_{i}(|r_{i}|-R_{0}))-1\bigr{]}^{2}\,,
V2[ri,ri1,ϕi,θi]=Ki,i1(1+Gi,i1)di,i1¯2,\displaystyle V_{2}[r_{i},r_{i-1},\phi_{i},\theta_{i}]=\,K_{i,i-1}\cdot\bigl{(}1+G_{i,i-1}\bigr{)}\cdot\overline{d_{i,i-1}}^{2}\,,
Gi,i1=ρi,i1exp[αi,i1(|ri|+|ri1|2R0)].\displaystyle G_{i,i-1}=\,\rho_{i,i-1}\exp\bigl{[}-\alpha_{i,i-1}(|r_{i}|+|r_{i-1}|-2R_{0})\bigr{]}\,. (2)

For the chain in Fig. 1(b), di,i1¯\overline{d_{i,i-1}} is replaced by di,i1(S)¯\overline{d_{i,i-1}(S)} in the fifth of Eqs. (2).

Ha[r1]H_{a}[r_{1}] is taken out of the sum, as the first base pair is coupled only to the successive base pair along the chain. Being the sum of a one-particle inter-strands potential V1[ri]V_{1}[r_{i}] and a two-particles intra-strand potential V2[ri,ri1,ϕi,θi]V_{2}[r_{i},r_{i-1},\phi_{i},\theta_{i}], HH contains the essential forces which stabilize the helix proho93 ; olson10 ; metz11 ; maher13 ; onuf19 .

V1[ri]V_{1}[r_{i}] is usually expressed by a Morse potential whose hard core represents the repulsive electrostatic interaction between negative phosphates on complementary strands. DiD_{i} is the base pair dissociation energy which defines the Morse plateau. bib_{i} sets the potential range through its inverse value. As the bases vibrate, the pair relative distance may become even smaller than R0R_{0}. However, radial fluctuations have a lower bound which follows from the fact that the potential hard core reduces the statistical weight of too negative base pair contractions. Hence, by imposing the condition   V1[ri]DiV_{1}[r_{i}]\leq D_{i}   we incorporate in the calculation only fluctuations satisfying the inequality, |ri|R0ln2/bi|r_{i}|-R_{0}\geq-\ln 2/b_{i}.

A solvent potential term can be added to the Morse potential in V1[ri]V_{1}[r_{i}] to enhance the height of the energy threshold below which the base pair is closed and introduce a hump over the Morse plateau thus accounting for the re-closing barrier and possible strand recombination with the solvent. While several analytical choices have been proposed in the literature pey08 , we take in the following a widely used expression which well accounts for the mentioned physical effects coll95 ; druk01 ; io11b :

VSol[ri]=Difs(tanh((|ri|R0)/ls)1)\displaystyle V_{Sol}[r_{i}]=\,-D_{i}f_{s}\bigl{(}\tanh((|r_{i}|-R_{0})/l_{s})-1\bigr{)} (3)

where fsf_{s} is the factor which increases the energy barrier for base pair dissociation and the length lsl_{s} defines the spatial range of the solvent.

The two-particle potential V2[ri,ri1,ϕi,θi]V_{2}[r_{i},r_{i-1},\phi_{i},\theta_{i}] extends to the 3D helical model the peculiar stacking term chosen in the 1D PBD model. The choice was motivated to account for the cooperativity effects in the melting transition of DNA chains, in the thermodynamic limit pey93 . The potential displays the stacking dependence on the angular variables and contains three parameters per dimer   i.e., the elastic Ki,i1K_{i,i-1} and the anharmonic force constants ρi,i1\rho_{i,i-1}, αi,i1\alpha_{i,i-1}. Enforcing the condition   αi,i1<bi\alpha_{i,i-1}<b_{i}, the range of the stacking is taken larger than that of the Morse potential. Accordingly, if the inequality   |ri|R0αi,i11|r_{i}|-R_{0}\gg\alpha_{i,i-1}^{-1}   holds, the hydrogen bond breaks and the local stacking interaction drops from Ki,i1(1+ρi,i1)\sim K_{i,i-1}(1+\rho_{i,i-1}) to Ki,i1\sim K_{i,i-1} thus favoring the breaking also of the adjacent base pair. This may trigger cooperatively the formation of local bubbles.

The Hamiltonian in Eq. (2) differs substantially from the original PBD. First, the inclusion of the bending degree of freedom permits to address the cyclization and flexibility properties of open ends chains which are beyond the range of a ladder model. Second, the model with bending can be adapted to deal with circular molecules whereas a ladder model cannot io14b . Third, the presence of the twist angle provides a restoring force which stabilizes the stacking potential io12 . Precisely, if a base pair undergoes a large fluctuation rir_{i}, the neighboring fluctuation ri1r_{i-1} can take only a restricted range of values so that the energy scale of V2V_{2} remains of order DiD_{i} or below. These fluctuations are those which contribute most to the partition function. Instead, in the absence of a twist, both adjacent base pair may fluctuate independently and still yield a low V2V_{2} value. In this sense, the stacking of the PBD model potential does not sufficiently discourage large base pair fluctuations whereas the twisted stacking V2V_{2} in Eq. (2) confers stability to the double stranded molecule. Importantly, it also follows that a twisted model avoids the well known divergence of the partition function found in the 1D ladder model zhang97 . Such divergence can be tackled in 1D only by truncating the real space available to base pair fluctuations.

Then, Eq. (2) models the essential interactions through five parameters which can be tweaked to account for sequence heterogeneity in chains of any length. For short linear chains, finite size effects are relevant and should be implemented by taking open boundary conditions joy07 and parameter values at the chain ends consistent with looser hydrogen bonds and/or weaker stacking, although the strength of the latter may vary with the specific dimer zgarb14 .

Computation of the thermodynamics and mechanical properties for the model in Eq. (2) requires integration over both the radial and the angular variables with suitable choice of the respective integration cutoffs. In general, to perform such calculations, one has first to minimize the free energy over a broad range of possible twist angles and then select the equilibrium twist conformation, as explained in detail in ref. io17 . As a result, this procedure can be significantly time consuming even for short fragments io19 .

However, for the task of establishing a consistent relation between model parameters and radial cutoff, the details of the base pair fluctuations over the twist and bending angles make a minor contribution. Accordingly we chose to replace ϕi\phi_{i} and θi\theta_{i} by average values ϕ¯\bar{\phi} and θ¯\bar{\theta} which are taken as input parameters. This permits to retain the 3D nature of the model whereas the computational time is markedly reduced. Further, by tuning ϕ¯\bar{\phi} and θ¯\bar{\theta}, one can study: i) the relation between radial cutoff and macroscopic helical conformation, ii) to which extent such relation is affected by the sliding motion depicted in Fig. 1(b). To this purpose we refer hereafter to the helical repeat i.e., the average number of base pair per helix turn defined by, h= 2π/θ¯h=\,2\pi/\bar{\theta} as a macroscopic measure of the twist conformation of the molecule.

III. First-passage probability for a base pair

While DNA breathing can be tuned by binding enzymes that enhance the lifetime of the open states bonnet03 ; marcus13 and change the helical repeat stasiak97 , it is recognized that nucleic acids molecules display an intrinsic dynamics due to the conformational fluctuations around the average helix diameter which expose the bases to the solvent and allow regulatory proteins to access the code kalos11 ; wang12 ; marko15 ; lee19 .

The idea underlying the computational method is that the base pair thermal fluctuations are trajectories, providing an example of Brownian motion albeit subjected to the constraint that the base pairs are organized, at room temperature, in a stable helical structure. Accordingly, in the finite temperature path integration fehi , the fluctuation rir_{i} in Eq. (2) is conceived as a trajectory ri(τ)r_{i}(\tau) where τ\tau is the Euclidean time varying in the range [0,β][0,\beta] and β\beta is the inverse temperature. Imposing the closure condition ri(0)=ri(β)r_{i}(0)=\,r_{i}(\beta) nn1 , we can expand ri(τ)r_{i}(\tau) in Fourier series around R0R_{0}:

ri(τ)=R0+m=1[(am)icos(2mπβτ)+(bm)isin(2mπβτ)],\displaystyle r_{i}(\tau)=\,R_{0}+\sum_{m=1}^{\infty}\Bigl{[}(a_{m})_{i}\cos\bigl{(}\frac{2m\pi}{\beta}\tau\bigr{)}+(b_{m})_{i}\sin\bigl{(}\frac{2m\pi}{\beta}\tau\bigr{)}\Bigr{]}\,, (4)

whereby a set of coefficients defines a base pair state and measures the fluctuational amplitude. The path integral method applied to DNA has been discussed in a number of papers, e.g. refs.io11b ; io14c , to which we refer for details. While the expansion in Eq. (4) holds for any base pair in the chain, the focus is now on the motion of the jthj-th base pair identified, for instance, with the mid-chain base pair in Fig. 1(a). For the latter, the assumed initial condition for the average pair mates separation is, <rj>=R0<r_{j}>=\,R_{0}. It is remarked that, for the argument proposed hereafter, any other base pair in the chain could have been picked besides the mid-chain one.

Due to thermal fluctuations, at any successive time tt, rjr_{j} may exceed R0R_{0} or even shrink compatibly with the physical constraint defined by the hard core Morse potential. Accordingly, Pj(R0,t)P_{j}(R_{0},\,t) is defined as the probability that rjr_{j} does not return to R0R_{0} until tt whereas Fj(R0,t)=dPj(R0,t)dtF_{j}(R_{0},\,t)=\,-dP_{j}(R_{0},\,t)dt is the probability that the path will return for the first time to the initial R0R_{0} in the time interval dtdt past tt. The need to introduce two time variables, tt and τ\tau, arises from the fact that, for a given tt, the probabilities are given as a sum over the particle histories rj(τ)r_{j}(\tau) in the time lapse [0,t][0,t] maj05 .

For the jthj-th base pair embedded in the chain and interacting with its first neighbors via the stacking potential in Eq. (2), we write Pj(R0,t)P_{j}(R_{0},\,t) as:

Pj(R0,t)=Dr1exp[Aa[r1]]i=2,ijNDriexp[Ab[ri,ri1,ϕ¯,θ¯]]\displaystyle P_{j}(R_{0},\,t)=\,\oint Dr_{1}\exp\bigl{[}-A_{a}[r_{1}]\bigr{]}\cdot\prod_{i=2,\,i\neq j}^{N}\oint Dr_{i}\exp\bigl{[}-A_{b}[r_{i},r_{i-1},\bar{\phi},\bar{\theta}]\bigr{]}\cdot\,
rj(0)rj(t)Drjexp[Ab[rj,rj1,ϕ¯,θ¯]]τ= 0tΘ[rj(τ)R0],\displaystyle\int_{r_{j}(0)}^{r_{j}(t)}Dr_{j}\exp\bigl{[}-A_{b}[r_{j},r_{j-1},\bar{\phi},\bar{\theta}]\bigr{]}\cdot\prod_{\tau=\,0}^{t}\Theta\bigl{[}r_{j}(\tau)-R_{0}\bigr{]}\,,
(5)

where the Heaviside function Θ[..]\Theta[..] enforces the condition that rj(τ)r_{j}(\tau) has to remain larger than R0R_{0} for any τ[0,t]\tau\in[0,t]. This is implemented in the code by evaluating at any   τ\tau the amplitude of rj(τ)r_{j}(\tau) in Eq. (4) and discarding those sets of coefficients which don’t comply with such condition.

The action functionals in Eq. (5) are obtained by the following dτd\tau integrals:

Aa[r1]=0β𝑑τHa[r1(τ)],\displaystyle A_{a}[r_{1}]=\,\int_{0}^{\beta}d\tau H_{a}[r_{1}(\tau)]\,,
Ab[ri,ri1,ϕ¯,θ¯]=0β𝑑τHb[ri(τ),ri1(τ),ϕ¯,θ¯],\displaystyle A_{b}[r_{i},r_{i-1},\bar{\phi},\bar{\theta}]=\,\int_{0}^{\beta}d\tau H_{b}[r_{i}(\tau),r_{i-1}(\tau),\bar{\phi},\bar{\theta}]\,,
Ab[rj,rj1,ϕ¯,θ¯]=0t𝑑τHb[rj(τ),rj1(τ),ϕ¯,θ¯],\displaystyle A_{b}[r_{j},r_{j-1},\bar{\phi},\bar{\theta}]=\,\int_{0}^{t}d\tau H_{b}[r_{j}(\tau),r_{j-1}(\tau),\bar{\phi},\bar{\theta}]\,, (6)

while the measures of integrations over closed (Dri\oint{D}r_{i}) and open (Drj\int Dr_{j}) trajectories are coupled via the two particle potential which connects the ithi-th and jthj-th base pairs along the stack. For the two measures however, different criteria should be applied to set the integration cutoffs on the radial fluctuations.

Consistently with Eq. (4), Dri\oint{D}r_{i} is explicitly defined by

Drim=1(mπλcl)2×Λi(T)Λi(T)d(am)iΛi(T)Λi(T)d(bm)i,\displaystyle\oint{D}r_{i}\equiv\prod_{m=1}^{\infty}\Bigl{(}\frac{m\pi}{\lambda_{cl}}\Bigr{)}^{2}\times\int_{-\Lambda_{i}(T)}^{\Lambda_{i}(T)}d(a_{m})_{i}\int_{-\Lambda_{i}(T)}^{\Lambda_{i}(T)}d(b_{m})_{i}\,,\, (7)

where λcl\lambda_{cl} is the classical thermal wavelength and Λi(T)\Lambda_{i}(T) is the temperature dependent cutoff which sets the maximum amplitude for the base pair Fourier coefficients. In turn, the cutoff can be determined by using the property that the path integral measure Dri\oint{D}r_{i} normalizes the kinetic action, i.e.,

Driexp[0β𝑑τμ2r˙i(τ)2]= 1.\displaystyle\oint{D}r_{i}\exp\Bigl{[}-\int_{0}^{\beta}d\tau{\mu\over 2}\dot{r}_{i}(\tau)^{2}\Bigr{]}=\,1\,.\, (8)

In fact, using Eqs. (4), (7), the l.h.s. of Eq. (8) yields a product over the Fourier components of decoupled Gaußian integrals. Setting Λi(T)Uiλcl/mπ3/2\Lambda_{i}(T)\equiv\,{{U_{i}\lambda_{cl}}/{m\pi^{3/2}}}, it is numerically found that Ui= 2U_{i}=\,2 suffices to satisfy Eq. (8). While this provides a formal criterion to establish a minimum UiU_{i} for the closed trajectories in the chain, it is understood that larger cutoffs may be also assumed for practical tasks e.g., for the computation of flexibility properties which essentially depend on large amplitude base pair fluctuations io11a .

As for the jthj-th base pair, while the path expansion in Eq. (4) can be performed, one cannot apply the normalization condition in Eq. (8). This follows from the observation that, for any τ\tau, rj(τ)r_{j}(\tau) is defined up to rj(t)r_{j}(t) which is in fact an open trajectory for any t<βt<\beta. Hence, a new criterion should be developed to estimate the integral cutoff on the amplitude of the jthj-th radial fluctuation.

To this purpose we examine Eq. (5) and ask the question: what is the probability that, at the initial time, the jthj-th fluctuation is larger than R0R_{0} ?

From Eq. (4), at t= 0t=\,0, the jthj-th trajectory is   rj(0)=R0+m=1(am)jr_{j}(0)=\,R_{0}+\sum_{m=1}^{\infty}(a_{m})_{j} with the Fourier coefficients being integrated on an even domain. Accordingly, the initial probability Pj(R0, 0)P_{j}(R_{0},\,0) is expected to be 1/2\sim 1/2 n2 . This is the benchmark to be met in the computation of the first-passage probability as a function of time. In order to implement this criterion, we first truncate the Fourier integration rj(0)rj(t)Drj\int_{r_{j}(0)}^{r_{j}(t)}Dr_{j} in Eq. (5) by a cutoff   Λj(T)=Ujλcl/mπ3/2\Lambda_{j}(T)=\,{{U_{j}\lambda_{cl}}/{m\pi^{3/2}}}, with tunable UjU_{j} and then determine the value   UjU_{j} such that Pj(R0, 0)1/2P_{j}(R_{0},\,0)\sim 1/2. In this way one selects the cutoff on the base of a physical constraint for a specific set of model parameters and for a given helical conformations.

Refer to caption
Refer to caption
Figure 2: (Color online) (a) First-passage probability versus time for the mid-chain base pair in equilibrium with the N1N-1 base pairs at room temperature. A homogeneous short chain with average twist and bending angles is considered. Various values for the integral cutoff (UjUU_{j}\equiv{U}) over the radial base pair fluctuations are assumed. (b) For any integral cutoff, the zero time probability is computed as a function of the number of trajectories for the chain dimer containing the jthj-th base pair.

IV. Results

Our theory is first tested by setting in Eq. (5) the average bending at ϕ¯= 6o\bar{\phi}=\,6^{o} consistently with the indications of Fluorescence Resonance Energy Transfer studies probing the DNA bending elasticity at short length scales kim14 . The average twist angle is initially taken as θ¯= 36o\bar{\theta}=\,36^{o}, which means an average helical repeat   h= 10h=\,10, close to the usual experimental value for kilo-base long DNA chains at room temperature wang79 . Although the model potential contains only first neighbors radial interactions, the base pairs are correlated along the stack due to the helical conformation. In the following calculations, a short homogeneous chain of N= 21N=\,21 base pairs is considered which allows for about two turns of the helix whose diameter is R0= 20ÅR_{0}=\,20\AA.

The potential parameters are those taken in ref.io09 namely, Di= 30meVD_{i}=\,30\,meV, bi= 4.2Å1b_{i}=\,4.2\,\AA^{-1}, KiKi,i1= 60meVÅ2K_{i}\equiv K_{i,i-1}=\,60\,meV\cdot\AA^{-2}, ρiρi,i1= 1\rho_{i}\equiv\rho_{i,i-1}=\,1, αiαi,i1= 0.35Å1\alpha_{i}\equiv\alpha_{i,i-1}=\,0.35\,\AA^{-1}. This set is in line with the parameters used by other groups in investigations of the PBD model zhang97 ; campa98 and derived by fitting the experimental melting temperatures, though some discrepancies persist mostly as for the stacking force constants eijck11 ; io20 . For the dimers containing the terminal base pairs, subjected to end fraying effects weber15 , the stacking parameters   K2,1K_{2,1} and KN,N1K_{N,N-1} are taken one half of the value assumed for the internal dimers. While this choice is arbitrary, it permits to weigh the impact of chain end effects on the base pair cutoff.

We proceed by making, at first, the simplifying assumption to take a vanishing dd (like in the PBD ladder model) whose consequence will be pointed out below.

Given the set of parameters, the output of the calculation still essentially depends on a) the cutoff UjU_{j} and b) the number of base pair trajectories included in the integrals of Eq. (5). This is shown in Fig. 2(a) which plots the probabilities as a function of time for different UjU_{j} values and in Fig. 2(b) which plots the zero time probabilities versus the number of possible paths (NpN_{p}) for the dimer associated to the jthj-th base pair. As the time axis partitioning involves 1000 points, the zero time corresponds to t/β= 103t/\beta=\,10^{-3}.

Fig. 2(b) highlights the dependence of Pj(R0, 0)P_{j}(R_{0},\,0) on the path ensemble size and also provides the rule to select the appropriate NpN_{p} for a given UjU_{j}. In fact, the Pj(R0,t)P_{j}(R_{0},\,t)’s in Fig. 2(a) are computed by assuming, for each UjU_{j}, the NpN_{p} values which maximize the corresponding Pj(R0, 0)P_{j}(R_{0},\,0)’s thus ensuring that the calculation is sampling the largest set of base pair fluctuations, i.e. of molecule configurations, consistent with the model parameters.

It is found that, for   Uj 7.3U_{j}\sim\,7.3, the condition   Pj(R0, 0)1/2P_{j}(R_{0},\,0)\sim 1/2 is fulfilled provided that the ensemble size for the dimer fluctuations is   Np3.2106N_{p}\sim 3.2\cdot 10^{6}. Moreover, NpN_{p} is constant for all dimers in the stack but the terminal ones as all base pairs are subjected to the same physical constraints imposed by the model potential. Since the chain is homogeneous, the selected value is a good cutoff for all internal base pairs, thereby identified as UjU¯U_{j}\equiv\bar{U}. Instead, for the terminal base pairs, the first passage probability remains smaller than 1/21/2, e.g. Pj=N(R0, 0)0.018P_{j=N}(R_{0},\,0)\sim 0.018, signalling that the average base pair separation is larger than R0R_{0}. This effect is ascribable to the softer stacking force constant which causes looser bonds and fraying at the chain ends berg06 .

A. Twisting

Refer to caption
Refer to caption
Figure 3: (Color online) (a) First-passage probability versus time for the mid-chain base pair in equilibrium with the N1N-1 base pairs at room temperature. The chain is taken as in Fig. 2 except the average twist for which four hh values are considered. The respective average twist angles θ¯\bar{\theta} are (from top to bottom in the legend)   40o, 36o, 32.7o, 30o40^{o},\,36^{o},\,32.7^{o},\,30^{o}. For each twist conformation, the probability is computed assuming the respective cutoff Ū (in the inset) that fulfills the initial time condition. (b) Zero time probabilities versus integral cutoff for three twist conformations. The UU values for which the plots intersect the dashed line, are the U¯\bar{U}’s reported in the inset in (a). \large\bullet   marks the cutoff found for the PBD ladder model.

Next we study the interplay between cutoff and twist conformation for the homogeneous chain with d= 0d=\,0. Keeping the same model parameters as above, the molecule unwinding is simulated by increasing the helical repeat and the appropriate U¯\bar{U} is determined for each hh. The results are displayed in Fig. 3. The Pj(R0,t)P_{j}(R_{0},\,t)’s in Fig. 3(a) are computed by taking the U¯\bar{U}’s obtained respectively from the plots in Fig. 3(b). As a main result, U¯\bar{U} markedly decreases for larger hh. This is understood by observing that our helical chains are considered to be stable at room temperature also in the untwisted conformations although, in the latter, large amplitude base pair fluctuations would easily disrupt the hydrogen bonds and unstack the helix. Accordingly, helical molecules in a large hh conformation sustain only short scale fluctuations in order to preserve the overall stability. This result follows from the choice of a model which has been assumed to have no intrinsic stiffness i.e., d= 0d=\,0. By further increasing hh the helix unwinds and tends to the ladder representation. Consistently, it is shown in Fig. 3(b), that the selected U¯\bar{U}’s tend to the cutoff value determined for the PBD model io20 .

Refer to caption
Figure 4: (Color online) As in Fig. 3(a) but with a finite rise distance dd. Five twist conformations are considered. The respective average twist angles θ¯\bar{\theta} are (from top to bottom in the legend)   40o, 36o, 34.3o, 32.7o, 30o40^{o},\,36^{o},\,34.3^{o},\,32.7^{o},\,30^{o}. For each hh value, the first passage probability is computed assuming the respective cutoff Ū (in the inset) that fulfills the initial time condition.

This trend is however reversed once we introduce a finite rise distance which provides a more realistic model for the helical molecule and accounts for the intrinsic stiffness of the chain. The results are shown in Fig. 4 where the time dependent probability is computed as a function of time by varying the average twist angle. Now the finite dd confers stability to the helix which can sustain large amplitude base pair fluctuations also in the untwisted conformations. Accordingly the cutoff value, for which the zero time probability condition is fulfilled, grows versus hh as shown in the inset. For instance, given a chain with h= 10.5h=\,10.5 i.e. the standard DNA helical repeat at room temperature duguet93 , the obtained value U¯= 5.2\bar{U}=\,5.2 yields a maximum amplitude Λj(T)= 1.08Å\Lambda_{j}(T)=\,1.08\AA for the first Fourier component in Eq. (7). This, in turn, corresponds to a reasonable estimate of 2.2Å\sim 2.2\AA for the largest breathing fluctuation of the jthj-th base pair with respect to the average helix diameter in the closed state. This length is 10%\sim 10\% of R0R_{0} and it is a fair measure for the threshold above which hydrogen bonds are disrupted pey08 ; manghi16 .

B. Sliding and bending

Refer to caption
Figure 5: (Color online) The cutoff Ū, which enforces the constraint for the initial time probability, is plotted as a function of the ratio between absolute value of slide SS and bare rise distance dd. The schematic of the molecule with sliding motion is given in Fig. 1(b). Three twist conformations are considered while the average bending angle is ϕ¯=6o\bar{\phi}=6^{o}. The effect of a broader bending is shown in the inset for the conformation with h= 10h=\,10.

Now we focus on the geometrical model of Fig. 1(b) which contains a further parameter, i.e. the sliding motion of neighboring base pairs having the overall effect to shorten the helix. The model has a finite rise distance and the sliding SS is tuned as a fraction of dd n3 . The probability in Eq. (5) is computed also for this more complex configuration and the obtained cutoff U¯\bar{U} is plotted in Fig. 5 as a function of the ratio |S|/d|S|/d, assuming three twist conformations. As in the previous plots, the average bending angle is ϕ¯= 6o\bar{\phi}=\,6^{o}. It is found that U¯\bar{U} is enhanced in the presence of a finite SS. This behavior holds for all helical conformations and in particular for the case   h= 11h=\,11 with a sizeable sliding. The fact that a double helix with a pronounced base pair sliding should sustain broader radial fluctuations seems in accordance with the observation that dsRNA is wider than dsDNA dekker14 .

If however the short chain has a larger average bending flexibility, for instance due to melting bubbles formation or large angle deformations at specific dinucleotide steps widom04 ; archer06 ; vafa12 , then lower U¯\bar{U}’s should suffice in order to to fulfill the first-passage probability constraint. This is in fact the situation emphasized in the inset, whereby a chain with ϕ¯= 9o\bar{\phi}=\,9^{o} is assumed for the twist conformation with h= 10h=\,10. Accordingly the calculated U¯\bar{U}’s are somewhat reduced with respect to the case ϕ¯= 6o\bar{\phi}=\,6^{o}.

Altogether these findings are consistent with the interpretation that, at the microscopic level, the contraction of the rise distance HSH_{S} due to the slide may confer an enhanced flexibility to the molecule which is then capable to sustain large scale fluctuations while maintaining the stable helical conformations. It is emphasized that our calculation takes SS and hh as independent variables whereas the structural parameters are intertwined in nucleic acids that is, a larger slide is generally observed in helical molecules with a smaller average twist angle. To meet these experimental observation, one may tune SS and set up a self-consistent procedure to compute, by free energy minimization, both radial cutoff and ensemble averaged helical repeat.

Finally, we re-consider the Hamiltonian in Eq. (2) and add the solvent term in Eq. (3) to the one-particle potential. By increasing fsf_{s}, one can simulate the effect of a higher salt concentration in the solvent which screens the electrostatic repulsion between complementary strands thereby enhancing the base pair dissociation energy hence, an intrinsically more stable molecule configuration albu14 ; singh15 . Assuming the twist conformation h= 10h=\,10 with ϕ¯=6o\bar{\phi}=6^{o}, we have re-run the program in order to test how the cutoff Ū varies in the presence of the solvent factor. The results are shown in Fig. 6 where Ū is plotted as a function of |S|/d|S|/d. In accordance to the expectation, Ū is reduced with respect to the model without solvent (fs= 0f_{s}=\,0) and this holds both in the case with and without sliding motion. For the intermediate sliding case   |S|/d= 0.33|S|/d=\,0.33, we find a cutoff reduction of 18%\sim 18\% for a energy barrier increase of 10%\sim 10\% i.e., fs= 0.1f_{s}=\,0.1.

Refer to caption
Figure 6: (Color online) The cutoff Ū is plotted as a function of the ratio between slide SS and bare rise distance dd. The effect of the solvent potential in Eq. (3) is shown by tuning fsf_{s} . The twist conformation with h= 10h=\,10 is considered while the average bending angle is ϕ¯=6o\bar{\phi}=6^{o}.

V. Conclusions

Analysis of the thermodynamics and flexibility of nucleic acids at the mesoscopic scale requires computation of a partition function obtained as an integral over the space configuration available for the base pair fluctuations. Focusing on a three dimensional Hamiltonian model for double stranded helical molecules, we have proposed a statistical method to select the integration cutoff on the amplitude of the base pairs separation, a crucial parameter for quantitative predictions of the physical properties. The method conceives the base pair thermal fluctuations as a Brownian motion for a particle subjected to the interactions which render the molecule stable hence, the base pair vibrates around the equilibrium helix diameter. Setting a zero time condition for the particle trajectory, we have calculated the time dependent probability for the base pair to return to the initial position by summing over those particle histories which are physically consistent with the model potential. In turn, the first-passage probability obeys a zero time condition which sets the constraint to establish the proper cutoff on the base pair fluctuations. Applying the method to a short homogeneous molecule, we find that the cutoff depends significantly on the macroscopic helical conformation defined by the three variables i.e. twist, bending and sliding for the relative motion of two base pairs in any dimer of the chain. In particular, for a realistic Hamiltonian model incorporating an intrinsic stiffness in the stacking potential, the room temperature base pair fluctuations are estimated to be 10%\sim 10\% of the helix diameter. It is also found that the cutoff gets larger upon increasing the average twist angle in line with the expectation that the base pair fluctuations should be broader in untwisted helices. Likewise, the radial cutoff grows in molecules with a significant slide which shortens the rise per base pair along the molecule axis. This seems consistent with the observation that A-form helices are indeed wider and shorter than B-helices in which the slide is small. These main conclusions hold in general and are not restricted to the chain here considered, although quantitative estimates of the cutoff may vary according to the set of parameters chosen to model hydrogen bonds, stacking forces and solvent environment for a helical molecule. Thus, the first-passage probability method provides a sound statistical benchmark to select potential parameters and maximum base pair fluctuations for specific sequences and conformations of nucleic acids and, in particular, may foster mesoscopic studies of RNA properties so far less investigated than DNA.

DATA AVAILABLITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

  • (1) T.A. Knotts, N. Rathore, D.C. Schwartz, J.J. de Pablo, J. Chem. Phys. 126, 084901 (2007).
  • (2) T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).
  • (3) J.-H. Jeon, W. Sung, J. Biol. Phys. 40, 1-14 (2014).
  • (4) C. Nisoli, A.R. Bishop, J. Chem. Phys. 141, 115101 (2014).
  • (5) H. Li, Z. Wang, N. Li, X. He, H. Liang, J. Chem. Phys. 141, 044911 (2014).
  • (6) T. Dauxois, M. Peyrard, A.R. Bishop, Phys. Rev. E 47, R44-R47 (1993).
  • (7) A. Campa, A. Giansanti, Phys. Rev. E 58, 3585-3588 (1998).
  • (8) S. Zdravković, M.V. Satarić, Phys. Rev. E 73, 021905 (2006).
  • (9) B.S. Alexandrov, V. Gelev, Y. Monisova, L.B. Alexandrov, A.R. Bishop, K.\varnothing. Rasmussen, A. Usheva, Nucleic Acids Res. 37, 2405-2410 (2009).
  • (10) S. Srivastava, Y. Singh EPL 85, 38001 (2009).
  • (11) M. Zoli, Phys.Rev. E 79, 041927 (2009).
  • (12) M. Zoli, Phys.Rev. E 81, 051910 (2010).
  • (13) S. Srivastava, N. Singh, J. Chem. Phys. 134, 115102 (2011).
  • (14) R. Tapia-Rojo, D. Prada-Gracia, J.J. Mazo, F. Falo, Phys.Rev. E 86, 021908 (2012).
  • (15) A. Sulaiman, F.P. Zen, H. Alatas, L.T. Handoko, Phys. Scr. 86, 015802 (2012).
  • (16) D.X. Macedo, I. Guedes, E.L. Albuquerque, Physica A 404, 234-241 (2014).
  • (17) M. Zoli, Soft Matter 10, 4304-4311 (2014).
  • (18) I.V. Likhachev, V.D. Lakhno, Chem. Phys. Lett. 727, 55-58 (2019).
  • (19) M. Hillebrand, G. Kalosakas, Ch. Skokos, A. R. Bishop, Phys. Rev. E 102, 062114 (2020).
  • (20) G. Weber, Nucleic Acids Res. 41, e30 (2013).
  • (21) E. de Oliveira Martins, V.B. Barbosa, G. Weber, Chem. Phys. Lett. 715, 14-19 (2019).
  • (22) M. Peyrard, S. Cuesta-López, G. James, Nonlinearity 21, T91 (2008).
  • (23) M. Guéron, J.-L. Leroy, in Methods in Enzymology 261, Academic Press, San Diego (1985).
  • (24) M. Peyrard, S. Cuesta-López, D. Angelov, J. Phys.: Condens. Matter 21, 034103 (2009).
  • (25) M. Marty-Roda, O. Dahlen, T.S. van Erp, S. Cuesta-López, Phys. Biol. 15, 066001 (2018).
  • (26) F. Zhang, M.A. Collins, Phys. Rev. E 52, 4217 (1995).
  • (27) T. Strick, J.-F. Allemand, V. Croquette, D. Bensimon, Prog. Biophys. Mol. Biol. 74, 115–140 (2000).
  • (28) A. Garai, S. Saurabh, Y. Lansac, P.K. Maiti, J. Phys. Chem. B 119, 11146-11156 (2015).
  • (29) M. Joyeux, S. Buyukdagli, Phys. Rev. E 72, 051902 (2005).
  • (30) S. Ares, N.K. Voulgarakis, K.Ø. Rasmussen, A.R. Bishop, Phys. Rev. Lett. 94, 035504 (2005).
  • (31) Y.L. Zhang, W.M. Zheng, J.X. Liu, Y.Z. Chen, Phys. Rev. E 56, 7100-7115 (1997).
  • (32) T. Dauxois, M. Peyrard, Phys. Rev. E 51, 4027 (1995).
  • (33) R.M. Wartell, A.S. Benight, Phys. Rep. 126, 67-107 (1985).
  • (34) A. Singh, N. Singh, Physica A 419, 328-334 (2015).
  • (35) M. Zoli, Eur. Phys. J. E 34, 68 (2011).
  • (36) M. Zoli, J. Theor. Biol. 354, 95-104 (2014).
  • (37) H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physycs and Financial Markets, World Scientific Publishing, Singapore (2004).
  • (38) M. Zoli, J. Chem. Phys. 144, 214104 (2016).
  • (39) M. Zoli, J. Chem. Phys. 148, 214902 (2018).
  • (40) M. Zoli, Physica A 492, 903-915 (2018).
  • (41) M. Zoli, EPL 123, 68003 (2018).
  • (42) M. Zoli, EPL 130, 28002 (2020).
  • (43) M. Zoli, Phys. Chem. Chem. Phys. 18, 17666 (2016).
  • (44) M. Zoli, Phys. Chem. Chem. Phys. 22, 26901 (2020).
  • (45) M. Zoli, J. Chem. Phys. 135, 115101 (2011).
  • (46) C.R. Calladine, H.R. Drew, Understanding DNA, Academic Press, San Diego, USA, 1992.
  • (47) E. Herrero-Galàn, M.E. Fuentes-Perez, C. Carrasco, J.M. Valpuesta, J.L. Carrascosa, F. Moreno-Herrero, J.R. Arias-Gonzalez, J. Am. Chem. Soc. 135, 122-131 (2013).
  • (48) J. Lipfert, G.M. Skinner, J.M. Keegstra, T. Hensgens, T. Jager, D. Dulin, M. Köber, Z. Yu, S.P. Donkers, F.-C. Chou, R. Das, and N.H. Dekker, Proc. Natl. Acad. Sci. U.S.A. 111, 15408-15413 (2014).
  • (49) X.-J. Lu, W.K. Olson, Nucleic Acids Res. 31, 5108-5121 (2003).
  • (50) J.C. Wang, Proc. Natl. Acad. Sci. USA 76, 200-203 (1979).
  • (51) The slide is conventionally defined negative if the upper base pair in the dimer shifts to the right with respect to the lower base pair calla .
  • (52) I. Faustino, A. Pérez, M. Orozco, Biophys. J. 99, 1876-1885 (2010).
  • (53) M. Pasi, J.H. Maddocks, D. Beveridge, T.C. Bishop, D.A. Case, T. Cheatham, P.D. Dans, B. Jayaram, F. Lankas, C. Laughton, J. Mitchell, R. Osman, M. Orozco, A. Pérez, D. Petkevičiūtė, N. Spackova, J. Sponer, K. Zakrzewska and R. Lavery, Nucleic Acids Res. 42, 12272-12283 (2014).
  • (54) R.E. Dickerson, Nucleic Acids Res. 17, 1797-1803 (1989).
  • (55) A.A. Gorin, V.B. Zhurkin, W.K. Olson, J. Mol. Biol. 247, 34-48 (1995).
  • (56) M. Fixman, J.J. Kovac, J. Chem. Phys. 58, 1564 (1973).
  • (57) S. Smith, L. Finzi, C. Bustamante, Science 258, 1122-1126 (1992).
  • (58) M. Barbi, S. Cocco, M. Peyrard, Phys. Lett. A 253, 358 (1999).
  • (59) G. Weber, Europhys. Lett. 73, 806 (2006).
  • (60) Y.Z. Chen, E.W. Prohofsky, Phys. Rev. E 47, 2100 (1993).
  • (61) G. Zheng, L. Czapla, A.R. Srinivasan, W.K. Olson, Phys. Chem. Chem. Phys. 12, 1399-1406 (2010).
  • (62) S. Talukder, P. Chaudhury, R. Metzler, and S.K. Banik, J. Chem. Phys. 135, 165103 (2011).
  • (63) J.P. Peters, S.P. Yelgaonkar, S.G. Srivatsan, Y. Tor, L.J. Maher III, Nucleic Acids Res. 41, 10593-10604 (2013).
  • (64) A.V. Drozdetski, A. Mukhopadhyay, A.V. Onufriev, Front. Phys. 7, 195 (2019).
  • (65) K. Drukker, G. Wu, G.C. Schatz, J. Chem. Phys. 114, 579-590 (2001).
  • (66) M. Zoli, J. Chem. Phys. 141, 174112 (2014).
  • (67) M. Zoli, J. Phys.: Condens. Matter 24, 195103 (2012).
  • (68) S. Buyukdagli, M. Joyeux, Phys. Rev. E 76, 021917 (2007).
  • (69) M. Zgarbová, M. Otyepka, J. Sponer, F. Lankas, P. Jurecka, J. Chem. Theory Comput. 10, 3177-3189 (2014).
  • (70) M. Zoli, J. Phys.: Condens. Matter 29, 225101 (2017).
  • (71) M. Zoli, Phys. Chem. Chem. Phys. 21, 12566 (2019).
  • (72) G. Altan-Bonnet, A. Libchaber, O. Krichevsky, Phys. Rev. Lett. 90, 138101 (2003).
  • (73) C. Phelps, W. Lee, D. Jose, P.H. von Hippel, A.H. Marcus, Proc. Natl. Acad. Sci. U.S.A. 110, 17320-17325 (2013).
  • (74) K. Kiianitsa, A. Stasiak, Proc. Natl. Acad. Sci. U.S.A. 94, 7837-7840, (2013).
  • (75) A. Apostolaki, G. Kalosakas, Phys. Biol. 8, 026006 (2011).
  • (76) J.L. Killian, M. Li, M.Y. Sheinin, M.D. Wang, Curr. Opin. Struct. Biol. 22, 80–87 (2012).
  • (77) J.F. Marko, Physica A 418, 126-153 (2015).
  • (78) S.-R. Choi, N.-H. Kim, H.-S. Jin, Y.-J. Seo, J. Lee, J.-H. Lee, Comput. Struct. Biotech. J. 17, 797-804 (2019).
  • (79) R.P. Feynman, A.R. Hibbs, Quantum Mechanics and Path Integrals, (Mc Graw-Hill, New York, 1965).
  • (80) In Transfer Integral methods applied to the PBD model, the partition function is customarily obtained by applying periodic boundary conditions which amount to close the linear chain into a loop joy05 ; kalos20 . This procedure may be questionable in short chains due to the relevance of finite size effects. This drawback is solved by the Path Integral method which imposes the closure condition on the time axis whereas the chain maintains the open ends in real space. Furthermore, the PI method can be also extended to deal with circular molecules io13 .
  • (81) M. Zoli, J. Chem. Phys. 138, 205103 (2013).
  • (82) S.N. Majumdar, Curr. Sci. 89, 2076 (2005).
  • (83) While the Fourier coefficients are integrated on an even domain, too negative ama_{m}’s are discarded due to the physical condition associated to the hard core of the one particle potential, as discussed after Eq. (2). Such condition depends on the parameter which regulates the range of the Morse potential. This asymmetry in the choice of ama_{m}’s included in the computation explains why Pj(R0, 0)P_{j}(R_{0},\,0) may get slightly larger than 1/21/2. Hence the approximation sign used in the text.
  • (84) T.T. Le, H.D. Kim, Nucleic Acids Res. 42, 10786-10794 ( 2014).
  • (85) L. van Eijck, F. Merzel, S. Rols, J. Ollivier, V.T. Forsyth, M.R. Johnson, Phys. Rev. Lett. 107, 088102 (2011).
  • (86) I. Ferreira, T.D. Amarante, G. Weber, J. Chem. Phys. 143, 175101 (2015).
  • (87) D. Andreatta, S. Sen, J.L. Pérez Lustres, S.A. Kovalenko, N.P. Ernsting, C.J. Murphy, R.S. Coleman, M.A. Berg, J. Am. Chem. Soc. 128, 6885-6892 (2006).
  • (88) M. Duguet, Nucleic Acids Res. 21, 463-468 (1993).
  • (89) M. Manghi, N. Destainville, Phys. Rep. 631, 1 (2016).
  • (90) T.E. Cloutier, J. Widom, Mol. Cell 14, 355-362 (2004).
  • (91) C. Yuan, E. Rhoades, X.W. Lou, L.A. Archer, Nucleic Acids Res. 34, 4554-4560 (2006).
  • (92) R. Vafabakhsh, T. Ha, Science 337, 1097-1101 (2012).