Base pair fluctuations in helical models for nucleic acids
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 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.-aI. 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 up to about ten times the average helix diameter (which is ) 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 point-like base pairs, with reduced mass , first proposed in ref.io16b and depicted in Fig. 1(a). Essentially, the coordinate 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 overlaps , the i-th fluctuation vanishes and the pair mates distance equals the average helix diameter . 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 and by the bending between adjacent vectors and . 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 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 angle 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 .


Straightforward geometrical rules permit to derive the distances between adjacent radial displacements, i.e. in Fig. 1(a) and in Fig. 1(b). For instance is given by
(1) |
where 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 is the measure of . 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 | |
---|---|---|
24 | 20 | |
2.8 | 3.3 | |
S | -1.48 | 0.03 |
From Eq. (1), setting to zero both the intrinsic stiffness 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 and barbi99 or a non vanishing twist with zero weber06 .
For the helical linear chain with first neighbors interactions in Fig. 1(a), the Hamiltonian reads:
(2) |
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 and a two-particles intra-strand potential , contains the essential forces which stabilize the helix proho93 ; olson10 ; metz11 ; maher13 ; onuf19 .
is usually expressed by a Morse potential whose hard core represents the repulsive electrostatic interaction between negative phosphates on complementary strands. is the base pair dissociation energy which defines the Morse plateau. sets the potential range through its inverse value. As the bases vibrate, the pair relative distance may become even smaller than . 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 we incorporate in the calculation only fluctuations satisfying the inequality, .
A solvent potential term can be added to the Morse potential in 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 :
(3) |
where is the factor which increases the energy barrier for base pair dissociation and the length defines the spatial range of the solvent.
The two-particle potential 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 and the anharmonic force constants , . Enforcing the condition , the range of the stacking is taken larger than that of the Morse potential. Accordingly, if the inequality holds, the hydrogen bond breaks and the local stacking interaction drops from to 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 , the neighboring fluctuation can take only a restricted range of values so that the energy scale of remains of order 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 value. In this sense, the stacking of the PBD model potential does not sufficiently discourage large base pair fluctuations whereas the twisted stacking 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 and by average values and 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 and , 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, 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 in Eq. (2) is conceived as a trajectory where is the Euclidean time varying in the range and is the inverse temperature. Imposing the closure condition nn1 , we can expand in Fourier series around :
(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 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, . 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 , may exceed or even shrink compatibly with the physical constraint defined by the hard core Morse potential. Accordingly, is defined as the probability that does not return to until whereas is the probability that the path will return for the first time to the initial in the time interval past . The need to introduce two time variables, and , arises from the fact that, for a given , the probabilities are given as a sum over the particle histories in the time lapse maj05 .
For the base pair embedded in the chain and interacting with its first neighbors via the stacking potential in Eq. (2), we write as:
(5) |
where the Heaviside function enforces the condition that has to remain larger than for any . This is implemented in the code by evaluating at any the amplitude of 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 integrals:
(6) |
while the measures of integrations over closed () and open () trajectories are coupled via the two particle potential which connects the and 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), is explicitly defined by
(7) |
where is the classical thermal wavelength and 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 normalizes the kinetic action, i.e.,
(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 , it is numerically found that suffices to satisfy Eq. (8). While this provides a formal criterion to establish a minimum 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 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 , is defined up to which is in fact an open trajectory for any . Hence, a new criterion should be developed to estimate the integral cutoff on the amplitude of the radial fluctuation.
To this purpose we examine Eq. (5) and ask the question: what is the probability that, at the initial time, the fluctuation is larger than ?
From Eq. (4), at , the trajectory is with the Fourier coefficients being integrated on an even domain. Accordingly, the initial probability is expected to be 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 in Eq. (5) by a cutoff , with tunable and then determine the value such that . 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.


IV. Results
Our theory is first tested by setting in Eq. (5) the average bending at 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 , which means an average helical repeat , 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 base pairs is considered which allows for about two turns of the helix whose diameter is .
The potential parameters are those taken in ref.io09 namely, , , , , . 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 and 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 (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 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 values and in Fig. 2(b) which plots the zero time probabilities versus the number of possible paths () for the dimer associated to the base pair. As the time axis partitioning involves 1000 points, the zero time corresponds to .
Fig. 2(b) highlights the dependence of on the path ensemble size and also provides the rule to select the appropriate for a given . In fact, the ’s in Fig. 2(a) are computed by assuming, for each , the values which maximize the corresponding ’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 , the condition is fulfilled provided that the ensemble size for the dimer fluctuations is . Moreover, 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 . Instead, for the terminal base pairs, the first passage probability remains smaller than , e.g. , signalling that the average base pair separation is larger than . This effect is ascribable to the softer stacking force constant which causes looser bonds and fraying at the chain ends berg06 .
A. Twisting


Next we study the interplay between cutoff and twist conformation for the homogeneous chain with . Keeping the same model parameters as above, the molecule unwinding is simulated by increasing the helical repeat and the appropriate is determined for each . The results are displayed in Fig. 3. The ’s in Fig. 3(a) are computed by taking the ’s obtained respectively from the plots in Fig. 3(b). As a main result, markedly decreases for larger . 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 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., . By further increasing the helix unwinds and tends to the ladder representation. Consistently, it is shown in Fig. 3(b), that the selected ’s tend to the cutoff value determined for the PBD model io20 .

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 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 as shown in the inset. For instance, given a chain with i.e. the standard DNA helical repeat at room temperature duguet93 , the obtained value yields a maximum amplitude for the first Fourier component in Eq. (7). This, in turn, corresponds to a reasonable estimate of for the largest breathing fluctuation of the base pair with respect to the average helix diameter in the closed state. This length is of and it is a fair measure for the threshold above which hydrogen bonds are disrupted pey08 ; manghi16 .
B. Sliding and bending

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 is tuned as a fraction of n3 . The probability in Eq. (5) is computed also for this more complex configuration and the obtained cutoff is plotted in Fig. 5 as a function of the ratio , assuming three twist conformations. As in the previous plots, the average bending angle is . It is found that is enhanced in the presence of a finite . This behavior holds for all helical conformations and in particular for the case 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 ’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 is assumed for the twist conformation with . Accordingly the calculated ’s are somewhat reduced with respect to the case .
Altogether these findings are consistent with the interpretation that, at the microscopic level, the contraction of the rise distance 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 and 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 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 , 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 with , 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 . In accordance to the expectation, Ū is reduced with respect to the model without solvent () and this holds both in the case with and without sliding motion. For the intermediate sliding case , we find a cutoff reduction of for a energy barrier increase of i.e., .

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 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.. 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 ’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 ’s included in the computation explains why may get slightly larger than . 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).