Mesoscopic model for DNA G-quadruplex unfolding

Genomes contain rare guanine-rich sequences capable of assembling into four-stranded helical structures, termed G-quadruplexes, with potential roles in gene regulation and chromosome stability. Their mechanical unfolding has only been reported to date by all-atom simulations, which cannot dissect the major physical interactions responsible for their cohesion. Here, we propose a mesoscopic model to describe both the mechanical and thermal stability of DNA G-quadruplexes, where each nucleotide of the structure, as well as each central cation located at the inner channel, is mapped onto a single bead. In this framework we are able to simulate loading rates similar to the experimental ones, which are not reachable in simulations with atomistic resolution. In this regard, we present single-molecule force-induced unfolding experiments by a high-resolution optical tweezers on a DNA telomeric sequence capable of adopting a G-quadruplex conformation. Fitting the parameters of the model to the experiments we find a correct prediction of the rupture-force kinetics and a good agreement with previous near equilibrium measurements. Since G-quadruplex unfolding dynamics is halfway in complexity between secondary nucleic acids and tertiary protein structures, our model entails a nanoscale paradigm for non-equilibrium processes in the cell.

bead-spring model for dsDNA with three beads per nucleotide. They studied thermal properties of dsDNA and also those of ssDNA, specifically the melting of a DNA hairpin and the folding of the thrombin aptamer (a G4 with two G-tetrad planes). Rebic et al. developed a mesoscopic model following a bottom up approach 22 . Their model presents three different beads: one bead for the guanines, another for the nucleotides in the loop, and the last one for the ions, which interact by means of tabulated potentials.
In this work we propose a mesoscopic model for a fragment of the human telomeric DNA capable of forming a single G4 with a resolution of one bead per base and investigate both its thermal and its mechanical stability. The simplified representation of the model allows, on the one hand, the exploration of the mechanical unfolding at lower velocities than permitted in atomistic simulations and, on the other hand, the direct comparison with a high resolution mechanical unfolding experiment we herein present. Methods DNA G4 model. Each nucleotide and the two central ions are represented by a single bead as depicted in Fig. 1 where l i is the distance between two consecutive beads in the chain; d GG the distance between two guanines of the same plane; d IG the distance between each central ion and its neighbour guanines (this distance is defined for the eight closest guanines to each ion) and θ i the angle between three consecutive guanines. The Hamiltonian of such system is composed by the following terms.
• A harmonic interaction between two consecutive beads = ∑ − U k l l ( ) str i s i 1 2 0 2 , where l i = |r i+1 − r i |, r i is the position vector of the i-th particle, l 0 is the equilibrium separation between beads and k s the elastic constant of the chain. We take l 0 = 0.65 nm, which is approximately the distance between consecutive phosphor atoms of the backbone according to the crystal structure of the human telomeric parallel G4 24 .
• A Morse potential between consecutive beads of the same plane to simulate the Hoogsten hydrogen bonds between the guanines: , where the sums go over the number of planes p = 1,..., 3 and the number of guanines in each plane g = 1,..., 4, d GGg is the distance between two consecutive guanines in a G-quartet plane and r 0 the equilibrium length of the side of each plane. The strength of the hydrogen bond interaction highly depends on their environment. In the case of the G4, it has been shown that the hydrogen bonds between guanines are stronger when all the bonds of the same plane are formed 12,13 . To take into account this cooperative effect, we take the depth of the Morse potential D G to be dependent on the distances between the guanine of the same plane as follows:  where each nucleotide is represented by a single bead (not to scale). For simplicity, the twist between successive G-tetrad planes is not represented. Such rotation makes the distance between two consecutive planes (p 0 = 0.5l 0 , see the text for details) lower than the equilibrium distance between two consecutive beads in the chain (l 0 = 0.65 nm). Right: Potential between guanines of the same plane U GG and between one ion and its neighbour guanines U IG . The red lines mark the threshold distances used to define the coordination numbers between guanines and between the central ions and the neighbour guanines. • An interaction potential between each central ion and its neighbour guanines where the sums go over the two ions i = 1, 2 and the eight neighbour guanines g = 1,..., 8, A d / ig 12 is a repulsive term that accounts for the excluded volume effect and Q ig describes the strength of the effective attraction between the ion and the guanine. The constant A is selected in such a way that the minimum of U IG is located at + r p 2( /2) ( /2) 0 2 0 2 (p 0 is the distance between the centres of two consecutive planes). The Coulomb contribution Q ig /d ig is shifted linearly to zero in the interval a < b indicated in Fig. 1 with the two black vertical lines. A purely repulsive interaction between the two ions = U Q d /   I I  I I  I I   1 2  1 2 1 2 is also included. We set = . Q Q 2 5 ig I I 1 2 . It is worth noting that even if the interactions between the two ions and between each ion and the guanine have an electrostatic origin, their features are different. In fact, the interaction between the guanine base and the ion is due to the metal-ion coordination between the ion and the oxygen O6 of the guanine, and then is more sensitive to their distance separation than if it were a pure electrostatic interaction 25 . For this reason we introduce, as a typical procedure adopted in these cases, a cutoff distance for the interaction between the ion and the guanines.
• A bending energy interaction between the three consecutive beads that belong to each side of the G4 stem where k b is a bending elastic constant, θ 0 ≈ 150° due to the twist between planes (not represented in Fig. 1) and θ i is the angle between vectors l i+1 and l i . This term accounts for the stacking interactions between the consecutive guanines and confers stability to the G4 stem. For the beads of the loops the bending can be neglected.
122σ between all the beads of the G4. This term accounts for the excluded volume effect.
The dynamics of the system is obtained from the overdamped equations of motion of the j-th bead of the G4 (j = 1,..., 21) The last terms in eqs 2 and 3 represent the thermal contribution as a Gaussian uncorrelated noise. The damping is taken implicit into the time units. The mass of the ion m i and of the nucleotide m j are taken as m j = 3.85m i . We use the following dimensionless units: length is given in units of l 0 = 0.65nm and energy in units of D 0 . The energy and time units are derived in the next sections in order to match the experimental values of the G4 melting temperature and unfolding force with the simulations. The Lennard-Jones parameters are σ = 0.9 and ε = 0.5, while the rest of the parameters are D 0 = 1, k s = 100, k b = 2, α G = 10 and δ = 0.3, all in dimensionless units. The value of the latter parameter δ is chosen such that the cooperativity takes effect -by reducing the potential depth of every bond in the plane -when even a single guanines bond stays at a distance lying in the plateau of the Morse potential, i.e. when the bond is completely open. The first parameter D 0 has been set to 3, in the same order of magnitude as the value reported in ref. 12 , where the quantum contributions have been taken into account, but without thermal fluctuations.
The system of equations 2 and 3 are integrated with the stochastic Euler algorithm with dt = 10 −4 . For the melting simulations, were the dynamics is studied at different temperatures, the simulations are started from the lowest temperature and a folded conformation of the G4. The final positions and velocities at each temperature are used as the initial conditions for the next temperature. Each simulation lasts for 10 7 time steps, from which, the first 10 6 steps are for thermalisation. Pulling simulations are conducted at a temperature lower than the melting until the G4 is unfolded.

Force-induced unfolding experiments.
To adjust and validate the mesoscopic model, we performed constant velocity pulling experiments in which a DNA telomeric sequence that yields a G4 is unfolded by means of a high-resolution OT device, as depicted in Fig. 2. The central hexanucleotide-repeat sequence is flanked by dsDNA handles for its manipulation in the optical setup (see below). The trap stiffness is k = 0.135 ± 0.004pN/nm. The micropipette is moved relative to the optically trapped bead with a pulling velocity v exp = 11.8 ± 1.4 nm/s near the rupture event. The elastic constant acting on the G4 due to the dsDNA handles is estimated from the slope of their force-extension curve in the enthalpic elasticity regime before the rupture, which gives k DNA ≈ 0.4 ± 0.1 pN/nm. Synthesis of telomeric DNA molecules. The molecular construction for OT experiments consists of a telomeric G4-forming sequence (5′-TATA (GGGTTA) 5 TAGT-3′) flanked by two dsDNA handles, one of 650 bp at the 5′ end (handle A) and the other of 401 bp at the 3′ end (handle B), for their attachment to beads through digoxigenin-antidigoxigenin and biotin-streptavidin labelling, respectively (Fig. 2). The extra repetition and the flanking ssDNA nucleotides provide configurational flexibility for the formation of the G4 in the presence of the dsDNA handles. The three DNA fragments were obtained by PCR amplifications of a conveniently modified pUC18 plasmid 11 . Handle B was labeled during its PCR amplification using a 5′-biotinylated primer (Integrated DNA Technologies, IDT, Coralville, IA) and handle A was modified afterwards adding a very short tail of digoxigenin-dUTP at its 3′ end (DIG Oligonucleotide Tailing Kit, 2nd Generation (Roche); 37 °C-15 min). Both DNA templates and labeled handles were purified using Wizard SV Gel and PCR Clean-Up System (Promega). Finally, equimolar amounts of the telomeric DNA and the two DNA labeled handles were mixed and annealed in the presence of 100 mM KCl in a PCR apparatus using the annealing procedure described in ref. 11 .
Optical setup. Measurements have been performed in a dual-beam optical trap in which two diode lasers (250 mW at maximum power, 808 nm wavelength) and associated optics are compacted into a miniaturised instrument suspended from the ceiling 26 . Each beam is delivered by an optical fibre and its position in the plane perpendicular to the optical axis is controlled by bending the optical fibres using piezoelectric crystals. The two laser beams are counter-propagating and brought to the same focus with orthogonal polarisations, which allow their optical paths to be separated using polarising beam splitters. Each beam is passed through a pellicle beam-splitter that redirects about 5% of the intensity, which is used to measure the position of the trap. The remaining light is focused through water-immersion objectives (NA = 1.20) to form the optical trap in a microfluidics chamber, which also contains a micropipette. The light exiting the trap from each beam is collected by the opposite objective lens, which redirects it to position-sensitive detectors that monitor the three force components acting on the trapped bead. Force is measured using the light momentum conservation 27 . This setup design reduces the mechanical drift and allows a large measurement stability over time thus enhancing the discernibility of the rupture events associated with the unfolding of the DNA structure. It allows an approximate force resolution of 0.1 pN and a distance resolution of 0.5 nm.

Results
Melting simulations. In this section we study the thermal stability of the G4 with our model. To this end, different magnitudes are calculated as a function of the temperature T*: the average coordination number of each ion Co I , the average coordination number of each plane Co P , the radius of gyration of the molecule R G and the heat capacity C v . The coordination number of each ion is 1 if the distance between it and its neighbour guanine is lower than 1.63nm (2.5l 0 ) and 0 otherwise. In each plane the coordination number between two consecutive guanines is put equal to 1 if the distance between them is lower than 1.75 nm (a distance around the plateau of the Morse potential is reached) and 0 otherwise. Co I and Co P are defined as the average over the total possible contacts for each ion-guanine (8) and in each plane (4), respectively and over the simulation time at each temperature. Figure 3 shows the melting curves of the G4 described by the above mentioned magnitudes. For comparison, the results of a simulation without the central ions are also included. In both panels a) (Co I ) and b) (Co P ) we observe that the curves decrease in a stepwise fashion approximately at the same temperature for the two measures. Analogously, in the presence of the ions, R G (panel c), and C V (panel d) show a similar stepwise change at the same temperature. The narrow temperature interval in which Co P goes to zero, as well as the sharp transition of R G , are a consequence of the cooperativity term included in the definition of the Morse potential (Eq. 1). Conversely, in absence of the cooperativity term (δ = 0), the measures undergo their changes at higher temperatures (the breaking of the planes Co P is not shown for simplicity): in the case of the gyration radius the transition is quite smooth and the dissociation of the three planes Co P takes also place in a wide interval instead of appearing as sharp transitions. The transition temperatures are then very different from each other and do not permit any Experimental configuration for the mechanical unfolding of the G4 at the single-molecule level. The two ends of the G4 are tethered by two dsDNA handles, which are in turn attached to micronsized beads via biotin-streptavidin and digoxigenin-antidigoxigenin bonds. One of the beads is kept fixed by air suction to a micropipete while the other one is trapped in the laser-beam focus forming the OT. "B" stands for biotin and "D" for digoxigenin.
definition of a melting temperature. In other words, the contribution of the cooperative term appears important in determining the steepness and so the uniqueness of the melting temperature in the thermal unfolding.
Moreover, the abrupt changes in Co P , R G and C v , characteristic of the melting, is highly influenced by the presence of central ions. We observe in Fig. 3 (panels b,c and d) that the transitions of all the measures occur at a lower temperature if the ions are not present ('no I' labels), indicating that their coordination increases the thermal stability of the G4 20 .
It is important to note that the sharp changes in the behaviour of all the curves of the magnitudes described above occur at the same temperature. That gives the possibility to use indistinctly any of those magnitudes to quantify the thermal unfolding and define the melting temperature ⁎ T m . In each trajectory, we use the mean of the two energies ⁎ T I1 and ⁎ T I2 at which the coordination Co I of ions 1 and 2 respectively drop below the value 0.5, namely ⁎ T I I 1/ 2 . To account for the stochastic effects in the unfolding, we repeated the melting simulations N s = 100 times from which we 1/ 2 , where ⋅ denotes the ensemble average. The distribution of melting temperatures is presented in Fig. 4 (panel a). The other panels of the figure present analogue distributions by using magnitudes different than C O I : the loss of the plane coordination Co p (panel b), the increase of the gyration radius (panel c), and the position of the peak in the heat capacity max{C v } (panel d). By using the ion coordination Co I , the resulting melting temperature is = . We can now adjust the energy unit of the model. Taking the experimental value T m = 65 °C = 328K of the melting temperature of a telomeric sequence in [K + ] solution reported in ref. 28 we get: = = .
In the next section, the pulling simulations are conducted at T* = 0.4269 which corresponds to T = 296 K in real units. This temperature value is represented with a dashed line in Fig. 4. Note that this value lies inside the range of distribution of melting temperatures. Thus, when doing the pulling at this temperature, there is a non-negligible probability that the unfolding occurs due to thermal effects.
Mechanical unfolding at physiological conditions. In a previous work we studied the mechanical unfolding of different G4 conformations at the atomistic level by means of steered molecular dynamics and showed that the force measured during the unfolding was correlated with the loss of coordination of the central ions 20 . In that work, a harmonic spring k A is attached to one atom of the extreme of the G4 and displaced at constant velocity v, while another atom at the opposite end is fixed. The component of the force along the pulling , where x A (t) and x 1 (t) are the components of the distance along the pulling direction of the spring end (point A) and the pulled bead 1, respectively, as depicted in Fig. 5. The bead 21 is fixed to the origin of coordinates.
The unfolding force value F u obtained in the all-atom simulations is in the order of 10 2 pN 20 , one order of magnitude higher than previous experiments of G4 mechanical unfolding 8,11 and the experiment we present here where F u ≈ 20.4 ± 6.9 pN (mean ± s.d., N = 163 experiments). This difference is due to the high values of both the parameters velocity v = 1 nm/ns and the elastic constant k A = 1650 pN/nm necessary for all-atom simulations in order to reach a reasonable simulation time. With our mesoscopic model, we are able to decrease both values and to obtain unfolding forces comparable with the experiments. The elastic constant used in the mesoscopic simulations is set to k A = 0.4 pN/nm, in accordance to our experimental value, as explained in the previous section. Figure 6a and b show the force F as a function of the extension x 1 and the position of the spring x A for the different pulling velocities, respectively. All the curves exhibit a clear jump that coincides with the abrupt increase of the G4 extension (Fig. 6a), so revealing the unfolding of the G4 structure. In those conditions, the unfolding patterns reveal a unique unfolding force F u that we define as the maximum force measured in the spring before the jump. Different to the atomistic simulations, we find that F u is in the same order of magnitude as the experimental values and that decreases when lowering the velocity. This behaviour is due to the presence of thermal fluctuations, which facilitate the unfolding at lower pulling velocity. The nearly saturating behaviour at low velocities indicates that the unfolding occurs in the near equilibrium regime where F u is independent of the velocity. To express the velocity in real units, we need to specify the time units, which will be defined later when considering the mean value of the unfolding force as a function of the velocity. Figure 6c show the comparison between one simulations at v = 0.08 and the experiment, showing a clear agreement on the values of the unfolding force, and, in the inset, the distribution of unfolding forces obtained from the simulations, compared with the Gaussian drawn by using the experimental mean unfolding force and the corresponding standard deviation (drawn as a solid red line).  Figure 7 shows different magnitudes that characterise the mechanical unfolding: the distance between beads 1 and 21 d 1−21 (panel a), the gyration radius R G (panel b) and the average coordination number between the two central ions and their eight neighbour guanines (panels c and d). We note that the mechanical unfolding pathway presents a similar behaviour as the thermal one for all the simulated velocities. Specifically, the gyration radius R G increases abruptly almost at the same time as the ions coordination number goes to zero, showing the equivalence of the different measures. Moreover, the cooperativity term in the Morse potential also presents a similar effect in the mechanical unfolding compared with the thermal unfolding: if the cooperativity is removed the unfolding occurs at higher values of the force and a multi-peak structure is observed in the force extension curves, corresponding to the consecutive rupture of the different planes. Figure 8 shows the Potential of Mean Force (PMF). The PMF is the free energy along the extension x 1 and gives an equilibrium measure of the mechanical stability. It is calculated by using umbrella sampling 29 and the weighted histogram analysis method 30 . The initial conformations for calculating the PMF are taking from a pulling simulation at v = 0.01. The PMF exhibits a change of the convexity around x 1 = 2 nm which is in correspondence with the extension at which the force jumps during the pulling simulations.
To account for the influence of the stochastic effects during the unfolding, N R = 100 simulations at each velocity are performed. The distributions of the unfolding forces F u at each velocity are shown in Fig. 9. In agreement with the single realisation results, the distributions displace towards lower values of the forces as the velocity is decreased. This behaviour is better observed from the mean value of the unfolding force as a function of the velocity, or equivalently, as a function of the pulling rate r = k A v as shown in Fig. 10. In this figure we have included two   experimental force values: the unfolding force F u = 20.4pN obtained in our experiment at v exp = 11.8nm/s and the equilibrium force F u = 2.5 pN obtained in the constant force experiment of Long et al. 9 . The loading rate in the dimensionless units of the mesoscopic model corresponding to this first force value is r = 0.0014, so the velocity is v = r/k A = 0.0788l u /t u ≈ 0.051 nm/t u . Equaling this value to the experimental one v exp = 11.8 nm/s, we get the time unit to t u = 0.0043s. The dependence of F u vs r is similar to the observed in force dynamic spectroscopy experiments, where the strength of molecular bonds or the mechanical resistance upon unfolding is studied at different loading rates. Several analytical theories based on the Kramer's one-dimensional theory of diffusive barrier crossing in the presence of force have been proposed in order to get the kinetic parameters of a simplified energy landscape of the molecule at zero force from the F u vs r data: specifically, the height of the barrier separating the bounded/folded and unbounded/folded states G + , the position of the barrier x u and the unfolding rate constant in the absence of force k 0 . Two of the most widely used models in this analysis are the Bell-Evans-Ritchie (Bell) 31 and the Dudko-Hummer-Szabo (DHS) 32 . The Bell model predicts a linear behaviour for both the mean and the most probable rupture force as a function of the logarithm of the loading rate. In our data this linear behaviour is not observed and we use, more consistently, the DHS model, which predicts a non linear behaviour for high pulling velocity values. However, this model is not valid at low velocities where rebinding/refolding events can occur. For this reason we exclude the three lowest velocities from the fitting analysis with this model. Another kinetic model that takes into account refolding events in the dynamics, and then is valid in the region of low velocities, is the Yoreo model 33 . From this model the equilibrium unfolding force f eq , which is independent of the pulling velocity, is also obtained. f eq is the force at which the force dependent folding and unfolding rates are equal and depends on the elastic constant k A of the pulling spring: We will use this model to fit the simulated mean unfolding forces in the whole interval of loading rates.
The mean unfolding force as a function of the pulling rate with the DHS and Yoreo model reads, respectively: u eq B u In the DHS model (Eq. 4), r is the rate of variation of the applied pulling force, γ ≈ 0.577 is the Euler-Marchesoni constant and ν is a parameter that sets the shape of the free energy potential to a cusp potential (if ν = 1/2) or linear-cubic potential (if ν = 2/3). β = 1/k B T. In the Yoreo model (Eq. 5), f eq is the velocity independent equilibrium force and 2 is the force dependent unfolding rate from the Bell model. In our simulations the pulling is performed with a harmonic spring and then = = The values of the parameters obtained from the fitting are summarised in Table 1. Both variants of the DHS model, with ν = 1/2 and ν = 2/3, give similar results in terms of the model parameters and the goodness of the fit. Differently, with the Yoreo model the estimated distance x u is lower than DHS, while the transition rate k 0 is higher. A similar trend has been observed in the unfolding of both Titin and RNA when these parameters are obtained using the Bell model, where lower values of x u and higher values of k 0 are obtained with respect to the DHS model 32 . This behaviour agrees with the fact that both Yoreo and Bell models have the same dependence of the unbinding rate k u (F) as a function of the force. Force spectroscopy experiments performed in other G4 systems validate the order of magnitude of the parameters obtained with our model. For instance, Messieres et al. 7 obtained G + = 5.3k B T, x u = 0.9 nm and k 0 = 0.004 s −1 for a parallel G4 with four guanine tetrads by simultaneously fitting the unfolding force distributions at different loading rates r = 2, 7, 24 pN/s with the probability distribution function of the DHS model.

Mechanical unfolding at T = 0.
To better understand the meaning of the parameters obtained from the fitting in the context of our model we look at one pulling simulation without thermal noise (T* = 0). Figure 11 shows the behaviour of d 1−21 (distance between beads 1 and 21) and the force F as a function of the time during a pulling simulation. We can identify different folded conformations before the unfolding occurrence at t ≈ 9 × 10 5 . According to the behaviour of d 1−21 as a function of time we can split the folded conformations in Dudko

Exp.
Messiers v const. 7 5.3 0.9 0.004 Long F const. 9 (feq = 2.5pN) 0.6 0.24 Table 1. Fit parameters obtained from the mean unfolding force vs the loading rate. R 2 is the coefficient of determination that accounts for the goodness of the fitting. k A = 0.4 pN/nm. two main elongation ranges: (i) 0.8 < d 1−21 < 2.2 nm (t < 0.05 × 10 5 , visible in the right inset of the figure), and (ii) 2.2 < d 1−21 < 2.9 nm (0.05 × 10 5 < t < 9 × 10 5 in the main figure). When looking at the folded configurations corresponding to these two groups, we note that the increase in the extension for the first one is mainly related to a rearrangement of the distances and relative orientation between the planes, almost without difference in the distances of the very G4 planes: the planes have rotated along their axes. In the second interval, the increase in the extension is mainly related with the stretching of the sides of the planes. The possible extension of a guanine still bonded to the others in the G4 plane is ruled by the width of the Morse potential α −1 = 0.065 nm. Looking at the d 1−21 -values for different velocities at room temperature in Fig. 7a we realise that the extension of the molecule before the unfolding depends on the velocity: we observe that for low velocities, the unfolding is more likely to occur when the G4 lies in configurations belonging to the first elongation subset 1 < d 1−21 < 2.2 nm and for higher velocities when belonging to the second elongation subset 2.2 < d 1−21 < 2.9 nm (0.05 × 10 5 < t < 9 × 10 5 ). In fact, at low velocities (v < 0.01) the system does not stabilise at lengths of the second d 1−21 subset, while, conversely, for v > 0.01 the G4 does stabilise in the second subset before the G4 unfolds (see Fig. 7a). This means that the unfolding pathway depends on the velocity, and particularly at low velocities, the dynamics is strongly assisted by the thermal fluctuations, with also the presence of refolding events. In these conditions the DHS theory does not work, and for this reason the points at the lowest velocities have been excluded from our fitting analysis. The parameters we got from the fitting are in agreement with the transition from the second folded subset of lengths to the unfolded state. The values we calculated through the DHS model fit appear to be in agreement with the following conclusions: x u ≈ 0.61 nm for ν = 2/3 and x u ≈ 0.84 nm for ν = 1/2 are both in the order of the expected extension length (x u = 0.7 nm) of the second subset, which, in fact, is our estimation of the distance between the barrier and the state of the G4 before the unfolding. In addition, though the abrupt denaturation the G4 undergoes, the fit with the parameter ν = 2/3 reproduces a little better the expected barrier position, so suggesting that a smooth potential appears a better description than a cusp potential for this unfolding mechanism. The values obtained for G + with the DHS model are lower than the free energy value obtained in the PMF at x 1 = 2.2 nm which appears to be the unfolding point, while the value obtained with the Yoreo model is closer to the value of the PMF in this point. These results corroborate the fact that the DHS model is describing the transition from the second subset length that requires less energy than the transition from the first substate, which is better described by the energy estimated from the Yoreo model. The equilibrium force f eq calculated from this model is however larger than the experimental value reported by Long et al. 9 which can be due to the fact that f eq depends on the elastic constant used for the pulling. Finally we note that the unfolding force obtained without the thermal contribution (T* = 0) is F u ≈ 500 pN.

Conclusions
We have developed a mesoscopic model that captures the main features of DNA G4 thermal and mechanical stability. It is characterised by single beads representing each nucleotide of the ssDNA chain and the monovalent cations located at the central channel of the G4 stem. Model parameters are obtained based on our previous atomistic study of telomeric DNA G4 20 , on the melting temperatures of DNA G4 28 , and on the mechanical unfolding experiment conducted in this work. Among the many potential terms in the model, the cooperativity term between the guanines in the bond of the plane interactions has a key role. In fact it allows the modification of the shape of the thermal transition from sharp -when the cooperativity is strong -to smooth -when the cooperativity term is either low or removed. This phenomenological contribution can be adjusted to describe other conformations of the G4, such as the antiparallel or mixed arrangements of the strands.
The model correctly evidences the importance of the ions in the stabilisation of the G4 structure, whose rupture force are increased with their presence. More importantly, the model gives a very good description of the system under mechanical stretching. In this context, it is able to reproduce unfolding forces in the same order of magnitude as in experimental studies, forces that are impossible to be reached with atomistic calculations. Related to that, we are able to explore wide time scales, and study the unfolding pathways, by using loading rates up to five orders of magnitude lower than those allowed in microscopic simulations. The evaluated mean rupture force as a function of the loading rate nicely reproduces the nonlinear increasing behaviour observed in dynamic force spectroscopy experiments at high velocities and the almost force independent regime at low velocities, behaviours that are well fitted by the DHS and Yoreo models, respectively. The values of the parameters of the one dimensional free energy landscape assumed in both DHS and Yoreo models appears to be in very good agreement with both the umbrella sampling simulations (which permit to determine the energy barrier G + ), and to the pulling simulations at zero temperature (that allowed the estimation of the position of the potential barrier x u ). Moreover they are also in the order of magnitude of the corresponding parameters obtained in some other experiments 7 .
The proposed model, eventually with some extension of it, can be used as a tool for performing systematic studies on mechanical stability of different G4 conformations. In this regard it may be applied to understand different G4 geometries according to the loop orientation (parallel vs antiparallel) or allow the mechano-chemical analysis of G4 tandem repeats 11 , like those existing in human telomeric sequences.