{'fraction_non_alphanumeric': 0.03521744213233575, 'fraction_numerical': 0.010083002104278701, 'mean_word_length': 4.663190996358821, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Twisting short dsDNA with applied tension

The twisting deformation of mechanically stretched DNA molecules is studied by a coarse grained Hamiltonian model incorporating the fundamental interactions that stabilize the double helix and accounting for the radial and angular base pair fluctuations. The latter are all the more important at short length scales in which DNA fragments maintain an intrinsic flexibility. The presented computational method simulates a broad ensemble of possible molecule conformations characterized by a specific average twist and determines the energetically most convenient helical twist by free energy minimization. As this is done for any external load, the method yields the characteristic twist-stretch profile of the molecule and also computes the changes in the macroscopic helix parameters i.e. average diameter and rise distance. It is predicted that short molecules under stretching should first over-twist and then untwist by increasing the external load. Moreover, applying a constant load and simulating a torsional strain which over-twists the helix, it is found that the average helix diameter shrinks while the molecule elongates, in agreement with the experimental trend observed in kilo-base long sequences. The quantitative relation between percent relative elongation and superhelical density at fixed load is derived. The proposed theoretical model and computational method offer a general approach to characterize specific DNA fragments and predict their macroscopic elastic response as a function of the effective potential parameters of the mesoscopic Hamiltonian.


I. Introduction
Understanding the DNA mechanics has important implications as, in cells, DNA is constantly bent, stretched, repaired and processed by proteins which, upon binding, confer to the double helix its biological functions and regulate gene expression [1][2][3][4][5].
The development of optical and magnetic tweezers techniques over the last twenty five years has allowed to gain remarkable insights into the elastic properties of single DNA molecules by studying their response to external forces in the pico-Newton regime [6][7][8][9][10][11].
Such forces are required to oppose the thermal bending fluctuations due to the environment which constantly buffet the molecular bonds thus causing the helix to assume different random walk configurations. In fact, at physiological temperatures, the thermal energy per nano-meter is ∼ 4pN.
Force-extension data of kilo base pairs sequences have been well reproduced by wormlike-chain models that treat DNA as an isotropic rod whose behavior is dominated by entropic elasticity at least up to ∼ 10pN [12]. Instead, at higher external forces, structural changes occur in the intra-strand base pair covalent bonds and the helix is progressively over-stretched to a length larger than its B-form contour length [13]. Later measurements by the rotor bead tracking technique [14] have shown that kilo base pairs DNA molecules over-twist upon stretching up to ∼ 30pN and then untwist above such value.
Importantly, if a torque is applied in order to over-twist the double helix under a constant load, it has also been found that the molecule extends. While analogous results have been obtained by magnetic tweezers experiments, molecular models for DNA in a solvent [15,16] have suggested that a negative inclination of the base pairs towards the minor groove could reduce the helix diameter and elongate the rise distance. Moreover, the molecule extension appears to be a linear function of the applied over-twist in the limited range of those torsional strains which preserve the stable B-form.
Motivated by these findings pointing to a remarkable DNA flexibility together with a rich interplay between its twisting and stretching properties [17,18], we have developed in a previous work [19] a computational method based on a mesoscopic DNA Hamiltonian which treats the helix at the base pair level and retains the fundamental intra-strand and inter-strand base pair interactions, responsible for the helix stability in the presence of a solvent. Mesoscopic models have the capability to predict the thermomechanical behavior of specific sequences through optimization of the potential parameters via direct fitting of experimentally accessible data e.g., the melting profiles [20][21][22][23]. Essentially our method assumes that the single molecule may exist in a broad range of helical conformations, specified through the average number of base pairs per helix turn, and determines the energetically most convenient conformation by free energy minimization.
As the computation is carried out by varying the strength of an external load, one can predict the twisting response of the molecule as a function of the stretching perturbation.
While the method has been applied to short fragments which have been the focus of recent and widespread interest in view of their unexpected flexibility [24][25][26][27][28][29][30][31][32][33][34][35][36][37], the same scheme can be used (compatibly with the available CPU time) for any sequence and length being aware that the latter generally affect the properties of the molecules [38,39].
In this paper, going beyond our previous study, we assume that the helix may be overtwisted (or untwisted) with respect to its equilibrium conformation under a constant load and investigate the ensuing modification on the helical shape. It is emphasized that the load has here the function to align, not that to disrupt [40], the intra-strand stacking bonds.
Accordingly, the external force is tuned within a range of values which do not cause the over-stretching of the molecule backbone. In this way, the method offers a feasible approach to simulate the above described experimental setup. In particular, we derive here the quantitative relation between average helical elongation and superhelical density which, in principle, could be investigated experimentally for sequences of a few tens of base pairs. Furthermore, it is shown that the over-twisting / untwisting transition, observed in kilo-base long sequences as a function of the external force, is predicted by our model and essentially ascribed to a dependence of the helix bending fluctuation on the size of the applied load.
The geometrical representation for the helix is outlined in Section II while the mesoscopic Hamiltonian model is discussed in Section III. The general features of the computational method are given in Section IV and the formulas for the macroscopic helix parameters are defined in Section V. The results are presented in Section VI while some conclusions are drawn in Section VII.

II. Helical Model
In previous Hamiltonian studies of DNA denaturation [41], the double helix has been described by a basic ladder model, see Fig. 1(a), in which the bases are arranged as beads along the complementary strands. The backbone of a molecule with N base pairs is thus a chain of N −1 segments connecting the points O i (i = 1, ..., N). R 0 and d, input parameters of the model, represent the bare helix diameter and rise distance along the molecule stack, respectively, in the absence of fluctuations. Each pair is formed via the hydrogen bond connecting the two mates and only two degrees of freedom per pair, x (1) i and x (2) i , representing the displacements of the pair mates, are included in the model. The in-phase-displacement, i , yields a straightforward harmonic potential energy term in the chain Hamiltonian which can be exactly integrated. Instead the relative distance, r i = x (1) i − x (2) i , measured from the central helical axis, stretches the hydrogen bond and determines the statistical mechanics of the DNA ladder model [42,43]. In fact, r i may even become smaller than R 0 thus compressing the hydrogen bond but the pair mates cannot get too close to each other due to the strands repulsion exerted by the negatively charged phosphate groups [44].
Here we adopt a more realistic picture for the double helix which goes beyond the ladder model assuming that adjacent displacements along the molecule stack, e.g. r i and r i−1 , are allowed to twist and bend as shown in Fig. 1(b). Accordingly, the distance AB between neighbor base pairs is a function of the rotational degrees of freedom with both the torsional angle θ i and the bending angle φ i being integration variables in the calculation of the partition function. Hence twisting and bending fluctuations are incorporated in our model whereas other structural deformations, such as propeller twist that enhances the intra-strand base pair stacking and the presence of grooves relevant to the sequence specific DNA-protein binding [45,46], are not taken into account. While a general description of the base pair degrees of freedom and helical parameters is given e.g., in Ref. [47], the effects of heterogeneous base pair sequence are analyzed e.g. in a comprehensive molecular dynamics simulation for a large set of oligomers [48].
Restricting our study to homogeneous fragments, we further assume that, for short DNA molecules, the most energetically convenient conformations are those for which the helix axis is essentially planar [49,50], the O i 's are pinned to the sheet plane. This is consistent with the fact that, once short molecules close into a ring, the free energy of supercoiling is mostly partitioned into twisting while the writhe contribution to the linking number is negligible [51,52].
As the base pairs are described by the r i 's which depart from the central molecule axis, we assume: 1) to apply the force F ex to such axis, along the direction of the segment O i O i−1 in Fig. 1(b), 2) that the force pulls one end of the chain while the opposite end remains anchored, 3) that, for short molecules, the force acts uniformly on all the chain segments and 4) that the force has the effect to orient the chain segments along a specific direction.
More generally, in a two strand representation and in the experimental setup, one can devise a number of different schemes to pull the chain ends applying the load either to one or both complementary strands [53][54][55]. Moreover, even the pulling direction may become important yielding specific deformations of the helical structure, as shown by molecular force balance experiments on heterogeneous duplexes. This occurs if the forces and loading rates are such to disrupt the intra-strand stacking bonds [56]. These latter cases however do not pertain to our analysis as emphasized in the Introduction. We instead assume that F ex varies in a range of low to moderate values which do not over-stretch the helix and study the response of the homogeneous molecule to the external perturbation.
As detailed in the next Section, the computational technique sums over a large ensemble of molecule configurations and finds by free energy minimization a) the average twist fluctuation i.e., the average number of base pairs per helix turn, b) the average radial fluctuation which provides a measure of the molecule diameter and c) the average intra-strand elongation. This allows to determine the twisting profiles and the molecule structural deformations as a function of F ex .
A major research focus has been the denaturation transition with the associated Watson-Crick base pair openings and formation of bubble profiles which are crucial to the DNA biological functioning [69][70][71][72][73][74][75]. Theoretical analysis of the helix melting transition depart either from Ising-like models describing paired and unpaired complementary bases [76] or from Hamiltonian models treating the hydrogen bonds through a potential energy function of the distance between the pair mates [77,78]. While Ising models for DNA melting have in terms of a continuous variable, i.e. the relative separation between the pair mates on complementary strands. Furthermore, stability properties and flexibility of the helix can be studied in Hamiltonian models by introducing stacking potentials which account for the covalent bonds between neighbor nucleotides along the molecule strands. Certainly, interstrand and intra-strand forces cannot be treated independently as the disruption e.g., of the i-th hydrogen bond moves one (or both) mate(s) of the i-th base pair out of the stack thus affecting also the intra-strand interactions between the adjacent i and i ± 1 bases.
These requirements have been met by the Hamiltonian model used in recent analysis of the flexibility of short DNA sequences which have predicted cyclization values in fair agreement with available experimental data [79]. The same Hamiltonian, consistent with the helical model of Fig. 1(b), is adopted in the present study. In addition, we introduce the effect of a tunable external load which stretches the molecule axis and induces the changes in the optimal twisting configuration as determined by minimization of the free energy.
Then, the Hamiltonian for the helical molecule with N base pairs of reduced mass µ, stretched by a force F ex , is: Note that H a [r 1 ] is treated separately as the first base pair has no preceding neighbor along the molecule backbone.
V 1 [r i ] is one-particle potential accounting for the hydrogen bond between the i-th pair mates and V 2 [r i , r i−1 , φ i , θ i ] is the two-particle stacking term which also depends on the angular degrees of freedom.
Their analytic expressions are: (2) The one-particle potential comprises: a) a hydrogen bond Morse potential which enhances the threshold for base pair dissociation and stabilizes the hydrogen bond through the parameters f s and l s as pointed out in refs. [80,81].
The two-particle potential depends on the distance d i,i−1 between adjacent r i , r i−1 , as marked by the AB segment in Fig. 1(b), which includes the angular variables. The stacking potential contains both an elastic force constant K S and nonlinear terms weighed by the parameters ρ i,i−1 , α i,i−1 which favor cooperative behavior in the formation of local bubbles as it is understood by the following observation. In fact, when the conditions are fulfilled, the effective stacking coupling is ∼ K S · 1 + ρ i,i−1 hence, both the i and i − 1 base pairs are bound. However, thermal fluctuations may violate either-or of the above conditions thus weakening the relative hydrogen bond. In this case, G i,i−1 → 0, the corresponding coupling drops to ∼ K S and also the adjacent base along the stack loosens its hydrogen bond. Accordingly, the fluctuational opening spreads along the strands and untwists the helical molecule. This qualitative picture indicates that the equilibrium twist conformation is sensitive to the interplay between stacking parameters and amplitude of the base pair separations.
For homogeneous sequences, one can shorten, Then, consistently with our previous works [79] we take in the following calculations, For the bare helix diameter and rise distance, we set the values R 0 = 20Å and d = 3.4Å, respectively. This set of model parameters reproduces the experimentally estimated DNA free energies per base pair [62,82]. Further details on the potentials in Eq. (2) are given e.g., in refs. [81,85].
The fourth addendum in the last of Eq. (1) accounts for the presence of the external force field as described in Section II. As F ex has the purpose to straighten the chain thus opposing the coiling effect of the bending fluctuations, we expect that by increasing F ex , the entropy is reduced and, accordingly, the free energy grows.
It is also worth mentioning that our mesoscopic Hamiltonian may be extended to study the thermodynamics and flexibility properties of RNA [83,84] provided that a specific parametrization can be obtained for a model potential with twisting and bending variables.
The Hamiltonian in Eq. (1) basically represents a DNA molecule, in a solvent potential and stretched by an external load, whose equilibrium statistics is obtained by solving the associated partition function. This task is carried out by generating an ensemble of molecule configurations consistent with the model potential and integrating, for each base pair in the sequence, over a distribution of radial and angular degrees of freedom dense enough to stabilize the partition function. Previous studies based on molecular dynamics, Monte Carlo and transfer integral techniques for the Hamiltonian of the simplified ladder model in Fig. 1(a) [58,86,87] have encountered the problem to define an upper bound for the integral over the base pair distances |r i |. The problem stems from the fact that the two-particle potential of the ladder model vanishes for the zero mode (all r i 's equal) whereas the Hamiltonian remains finite for r i → ∞ as the one-particle potential is always bounded. Since V 1 [r i ] is not translationally invariant, the zero mode cannot be removed and the partition function diverges. Such divergence does not occur in our Hamiltonian model as the two-particle potential remains finite due the twist between adjacent base pairs. Nevertheless, a truncation of the configuration space for the base pair amplitudes is required in the computation and the choice of the cutoffs carries some arbitrariness. These issues are handled in the computational method based on the path integral formalism as outlined in Section IV.

IV. Computational Method
In our method, the base pair separations r i are thought of as trajectories r i (τ ) depending on the imaginary time τ = it, with t being the real time for the evolution amplitude in the [88,89]. Hence, τ varies in a range τ b − τ a whose amplitude is set by the inverse temperature β and the partition function for the Hamiltonian in Eq. (1) is obtained by integrating over closed trajectories, ( r i (0) = r i (β) ), defined along the τ -axis, following a route extensively described in previous works [92] .
Accordingly, r i (τ ) can be expanded in Fourier series : whose coefficients generate an ensemble of possible base pair paths and define the associated integration measure Dr i : where Λ 0 T and Λ T are the temperature dependent cutoffs and λ cl is the classical thermal wavelength [90].
As Eq. (4) normalizes the kinetic term in the action, i.e. [91]: our formalism provides a consistent method to derive the cutoffs in the integration over the path configuration space thus avoiding the above mentioned arbitrariness usually encountered in the application of transfer integral techniques to mesoscopic Hamiltonian models.
Furthermore, from Eq. (3), (5), one easily finds that Λ 0 T and Λ T are ∝ √ T hence, the maximum path amplitudes included in the partition function are larger at higher T as expected on physical grounds. Also note that Eq. (5) sets the free energy zero for our system and it is fulfilled for any µ in accordance with the fact that the free energy of a classical system does not depend on µ.
Writing Eq. (1) in terms of the Fourier expansion in Eq. (3) and summing over the bending and twisting degrees of freedom, we finally express the partition function Z N for the helical molecule with N base pairs: where φ M and θ M are the cutoffs on the amplitudes of the bending and twisting fluctuations, respectively.
Then, our computational technique consists of a direct sum over a set of molecule configurations which statistically contribute to the path integral with a Boltzmann weight.
Specifically, the ensemble size, given by the number of trajectories for any base pair in the chain, is enlarged until Z N numerically converges i.e., the state of thermodynamic equilibrium is achieved. This amounts to sum over ∼ 10 8 configurations for each dimer in the chain.
Over such ensemble, we perform the averages to obtain the macroscopic helical parameters as described in the next Section.
While Eq.
From Eq. (6) or Eq. (7), one derives the equilibrium thermodynamics from the free energy, While the results hereafter presented are obtained from Eq. (6), the same qualitative trend is found via computation of Eq. (7).

V. Rise distance and helix radius
As the focus of this work is on the helical torsional response to an applied load, we devise a recursive procedure to express the twist angle in Eq. (6). Precisely, θ i is measured with respect to the ensemble average < θ i−1 > obtained for the preceding base pair along the molecule stack (see Fig. 1(b)) and the latter value is incremented by 2π/h, where h, the number of base pairs per helix turn, is chosen within a physically suitable range. Moreover, for any h in such range, we integrate over a twist fluctuation variable θ f l i around the value < θ i−1 > +2π/h . This idea is formally expressed by: where n is the number of values sampled in the h range and ∆h is the incremental step. Using Eq. (8) the average twists are computed by: while, for the first base pair in the chain, < θ 1 >= 0. Hence, from Eq. (9), we derive the average helical repeat as: Technically, the program sets an increment ∆h and, for any input h, computes a set of values {< h > j , (j = 1, ..., n)} which differ from the initially chosen h-values. For any value < h > j , the corresponding free energy is calculated. As discussed in detail in ref. [19], the procedure is reiterated by taking a finer ∆h until the average twists and associated free energies converge, i.e., they do not further change by increasing n. By minimizing F , one finally selects the equilibrium average twist conformation, denoted in the following by < h > j * , in the presence of a load. By varying F ex , the calculation yields the twist-stretch profile for a specific molecule.
Likewise, one evaluates the average bending angles between adjacent base pair planes, , which also define the average angle between the stacking bonds and are measured with respect to the force direction, as shown in Fig. 1(b). By definition, < φ 1 >= 0.
The numerical convergence of the ensemble integrations should be also tested against the maximum twisting and bending fluctuations. For the twist angles, we find that the appropriate cutoff to stabilize the ensemble averages is, θ M = π/4, which allows for large torsions between adjacent base pairs. As twisting fluctuations may be sizeable also in the presence of F ex , there is no physical reason to introduce a force dependent twisting cutoff.
In the absence of external loads, the bending cutoff φ M = π/2 is large enough to allow for the formation of kinks having the effect to reduce the bending energy between neighbor base pairs [32,[95][96][97]. It is however plausible that, by increasing F ex , the intra-strand bonds should stretch and the amplitude of the bending cutoff should shrink. Although, at this stage, there is not experimental information to quantitatively account for such effect, we have tested some ad hoc functions φ M (F ex ) and chosen, φ M (F ex ) = π[1 − (c · F ex ) z ]/2 with tunable parameters c , z in the calculations hereafter presented.
Analogously to Eq. (9), ensemble averages are carried out to compute the macroscopic helical parameters, i.e., the average distances between adjacent base pairs along the stack, and the average base pair radial fluctuations: Note that the first base pair radial fluctuation, < r 1 >, is calculated with the Boltzmann weight given by the action A a [r 1 ] in Eq. (6).
Before displaying the results, we discuss a possible extension of our model to the over-

VI. Results
The model is applied to a very short homogeneous fragment with N = 10 base pairs.
The calculation is carried out with a fine mesh, i.e. ∆h = 0.015625 in Eq. (8). First, we compute in the absence of external loads, the ensemble averages of the one-particle and two-particle potential energy of the Hamiltonian in Section III. In Fig. 2, the terms > are plotted as a function of < h >. While < V 1 > shows a very weak dependence on < h > ascribable to the interplay between radial and angular fluctuations only in the Boltzmann statistical weight, < V 2 > displays a pronounced minimum versus < h > as the stacking explicitly depends on the twist angle. This points to the fact that some twist conformations are energetically favored.
The precise evaluation of the equilibrium twist is carried out in Fig. 3     Experiments on kilo-base sequences have located the upturn at F ex ∼ 30pN [14]. Analogous experiments on fragments with a few tens of base pairs would help to check whether a similar twisting pattern persists at short length scales and also serve as a criterion to fit the parameters of our model. While in principle the behavior of kilo-base sequences may differ from that of very short fragments, we have repeated the free energy minimization procedure for the same sequence of Fig. 3 varying the parameters of the bending cutoff φ M (F ex ). Some results are displayed in Fig. 5. The computation shows that the occurrence of the over-twisting / untwisting transition is indeed sensitive to the specific dependence of the maximum amplitude of the bending fluctuations on the applied load. By reducing c and (or) increasing z with respect to Fig. 3, we assume that the maximum bending fluctuations decrease more smoothly for large loads and, under these conditions, the over-twisting regime progressively extends up to a few tens of picoNewtons. At last, it is found that for c −1 = 100pN and z = 3, the upturn in the < h > j * plot can be set at F ex ∼ 30pN as in Ref. [14]. For the latter curve, we can estimate the superhelical density as a function of the load, σ(F ex ) = ∆T w/(T w) 0 , where ∆T w is the number of turns added to the helix by increasing the stretching perturbation and (T w) 0 is the unperturbed twist number. For instance, taking the calculated < h > j * at F ex ∼ 8.3pN, we obtain σ(F ex ) = 0.0098.
Our code however, computes the helix structural parameters for the whole set of conformations obtained from Eqs. (8)- (10) and not only for the conformation < h > j * (reported in Figs. 2, 3) which minimizes the free energy. Then we can monitor how the helical shape changes, away from the minimum, assuming an over-twisted conformation (< h > smaller than < h > j * ) or an untwisted conformation (< h > larger than < h > j * ). As this is done at a fixed external load, our reasoning simulates the experimental setups of ref. [14,15] in which small torsional strains are applied to the molecule under constant tension. Fig. 6 shows our findings for the same molecule of Fig. 3 and for two distinct loads. We see that an imposed over-twist causes the average helix diameter to shrink (panel (a)) and the average rise distance to elongate (panel (b)) with respect to the values corresponding to the equilibrium conformation < h > j * . Interestingly, if we slightly untwist the helix, the rise distance shortens (down to the dashed lines in panel (b)) whereas < d > eventually extends if we further untwist the helix (reaching the regime to the right of the dashed lines). This trend holds for any constant load in qualitative agreement with the behavior of kilo-base sequences of Ref. [14].
Likewise, assuming the bending cutoff parameters which set the upturn at F ex ∼ 30pN in Fig. 5, we find that < d > for the fragment with N = 10 base pairs grows by overtwisting the helix with respect to its equilibrium twist conformation, at constant force. The results are shown in Fig. 7(a) for F ex ∼ 9pN as, for such load, some data are available from Ref. [14]. Also the percent relative extension (≡ < d >/< d > j * − 1) is reported in Fig. 7(b) as a function of the superhelical density σ( * ) = (∆T w) j * /(T w) j * . Note that σ( * ) is defined with respect to the equilibrium twist number ((T w) j * = N/ < h > j * ) and (∆T w) j * measures the applied torsional strain at fixed force. Hence, σ( * ) is physically distinct from the previously defined σ(F ex ). The intra-strand distance grows almost linearly above σ( * ) ∼ 0.01 but, at the latter value, we obtain a relative increment which is an order of magnitude lower than that reported in Ref. [14]. It is however remarked that, for any given σ( * ), we are simulating (∆T w) j * values which are too small to allow a strict comparison with data relative to kilo-base sequences.
Finally, one may wonder to which extent the results presented so far are sensitive to the choice of the model potential parameters. In Fig. 8  sequence with base pair dissociation energy D i = 60meV (higher than in Fig. 3) and plot the average helical repeat as function of the external load, for three values of the nonlinear stacking parameter ρ i . All other model parameters are as in Fig. 3. The incremental step, ∆h in Eq. (8), is the same as in Fig. 3   result may suggest that longer sequences have higher intrinsic flexibility (larger ρ i in our model) and therefore they are more resilient to the stretching perturbations. Hence, they begin to untwist only under sizeable loads.
Comparing the green plot (ρ i = 1) with Fig. 3(b), it is found that the < h > j * 's values become slightly larger by increasing D i . This result is apparently surprising as one may expect that higher hydrogen bond energies stabilize the molecule and oppose the helix untwisting. In fact, assuming higher D i 's with all other parameters unchanged, yields more tightly bound complementary strands (< R > shrinks) whereas the intra-strand rise distance < d > grows (not shown here). As a consequence the molecule of the green plot in Fig. 8 appears, with respect to the molecule in Fig. 3, in a narrower albeit more elongated configuration which ultimately favors a slight untwisting and higher < h > j * 's. Note however that the pair dissociation energies are taken within a range consistent with the experimental free energies per base pair [79]. These findings have been qualitatively discussed to point out that, by virtue of the interplay between inter-strand and intra-strand base pair interactions, the model potential parameters are intertwined and should be determined, as a set, by fitting the available data for specific fragments [98].

VII. Conclusions
The mechanical response of DNA molecules to applied loads provides insights into the intrinsic flexibility properties of the helix which may vary with its sequence specificity, length and environmental conditions. While the theory of elastic rods generally accounts for the behavior of kilo-base pair filaments whose contour length largely exceeds their characteristic persistence length, current research is revealing that, at short length scales, all-atom simulations and mesoscopic models should be rather used to capture the elastic behavior of fragments which display strong base pair fluctuational effects. Following our previous investigations on the cyclization properties, end-to-end distance and persistence length of sequences with 100 base pairs or less, we have here studied how the helix macroscopic parameters and its twisting conformation may change under the effect of an external load. The computational method is based on a path integral description for the Hamiltonian model which comprises both inter-strand hydrogen bonds between the pair mates and intra-strand forces between adjacent bases along the molecule stack. Base pair separations are treated as trajectories contributing to the partition function with their specific statistical weight which is essentially determined by the physical constraints of the model potential. Thus, too large contractions of the base pair distance (with respect to the equilibrium helix diameter) are energetically discouraged by the electrostatic repulsions and therefore yield scarce weight to the partition function. On the other hand, the maximum amplitude of the base pair separations encounters in the path integration a temperature dependent cutoff which consistently truncates the configuration space. Crucial to our study are the bending and twisting fluctuations between neighbor base pairs which characterize the form of the stacking potential and also contribute to the partition function. Technically, the density of the base pair paths and angular fluctuations included in the computation is increased until the molecule free energy converges. For any applied force in a suitable pico-Newton range, it is assumed that the molecule may exist in a large set of possible twist conformations each characterized by an average number of base pairs per helix turn. Carrying out the ensemble averages over the base pair degrees of freedom, we eventually obtain the average equilibrium helical twist by minimizing the system free energy. Moreover, we compute the helical twist, away from the equilibrium, associated to over-twisted and un-twisted conformations which can be simulated keeping the external load constant. Thus, our numerical program can predict: a) how the helical twist changes by varying the load and b) how the helical parameters, i.e. average diameter and rise distance, change by over-twisting (untwisting) the helix at constant load.
Applying the method to a very short homogeneous helix, we have found a general pattern similar to that observed in kilo-base long sequences although short fragments may present macroscopic helix parameters which quantitatively differ from the long ones.
Specifically, tuning the cutoff in the bending fluctuations integration, the model can even reproduce the transition between over-twisting and untwisting regime experimentally located at ∼ 30pN in kilo-base sequences, with the caveat that such forces may be excessive to be applied to fragments of only ten base pairs as they could disrupt the stacking bonds.
Furthermore, adding (and subtracting) helix turns to the equilibrium conformation for a fixed stretching perturbation, we have derived the relation between helix elongation and superhelical density which appears essentially linear in the over-twisting regime.
We have also shown how the equilibrium helical repeat, predicted by minimizing the free energy in the presence of a load, could vary with the specific choice of model potential parametrization. Thus, we feel confident that the proposed method makes a valid computational tool to characterize specific fragments for which experiments could provide sufficient information to fit the set of potential parameters.