Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 11756.
Published online 2017 September 18. doi:  10.1038/s41598-017-10849-2
PMCID: PMC5603602

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.


G-Quadruplexes (G4) are non-canonical conformations of DNA or RNA sequences rich in guanine nucleobases. Unlike the classic double helix1, the basic structural unit in the G4 motif is the G-tetrad, a planar arrangement of four guanines (G) nucleobases held by Hoogsten hydrogen bonds2. The presence of monovalent cations like K+ and Na+ stabilises the highly electronegative central channel along the axis of the G4 stem. G4s have been observed in some biological sequences within telomeres and in promoter regions3,4, being their high mechanical stability one relevant property in their potential role in processes such as gene expression and chromosome maintenance5,6.

The mechanical stability of G4s have been studied at the single molecule level by means of AFM and optical tweezers (OT)712. Mechanical stability in this framework is defined in terms of the unfolding force or unfolding free energy. Typical values of the unfolding forces for G4 range between 20–30 pN and depend on both the G4 topology and the type of ions present in the G4 channel7. Single-molecule characterisation is at present limited by the temporal and spatial resolutions of the experimental techniques. Complementary information as the structural changes during the unfolding is obtained by molecular simulations.

The dynamics of the G4 has been modelled at different temporal and length scales, from quantum calculations13,14 and molecular dynamics simulations1520 to mesoscopic approaches21,22. Molecular dynamics allows an atomistic description of G4 properties. In this regard, we previously studied the mechanical unfolding of a fragment of the human telomeric sequence that can be folded in different geometries by using Steered Molecular Dynamics. We showed that the unfolding pattern in the force-extension curves is correlated with the loss of coordination of the central ions in the G4 and that its stability is significantly decreased if the ions are removed20. These results cannot be compared directly with the experimental results due to the high pulling velocity used in this molecular dynamics simulation (around 6 orders of magnitude higher than in the experiments), which is known to affect the unfolding forces.

Larger temporal and spatial scales can be achieved by means of mesoscopic models. Unlike for double-stranded (ds) DNA or single stranded (ss) structures, like DNA hairpins, only few mesoscopic models have been developed for G42123 and none, to our knowledge, to characterise the mechanical unfolding. Margaret et al.21 developed a 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 approach22. 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.


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 Ustr = ∑i½ks(lil0)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 G424.
  • A Morse potential between consecutive beads of the same plane to simulate the Hoogsten hydrogen bonds between the guanines: UGG = ∑pgDG(eαG(dGGgr0)−1)2, 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 formed12,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 d GGg are the four distances in the same plane and δ sets the length scale for the decay of D G. Thus D G varies from D G = 3D 0 when the G4 is folded and the side of the plane is r 0, to D 0 when any of the distances from the equilibrium position (d GGg  r 0) increases beyond the length scale δ −1, where D 0 is the Morse potential depth of a single guanine couple.
  • An interaction potential between each central ion and its neighbour guanines UIG=ig(A/dig12Qig/dig), where the sums go over the two ions i = 1, 2 and the eight neighbour guanines g = 1,..., 8, A/dig12 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 2(r0/2)2+(p0/2)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 UI1I2QI1I2/dI1I2 is also included. We set Qig = 2.5QI1I2. 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 interaction25. 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 Uben=ikb/2[1cos(θiθ0)], 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.
  • A repulsive Lennard Jones interaction VLJ = ∑ij>i4ε[(σ/rij)12 − (σ/rij)6] if r ij < 1.122σ between all the beads of the G4. This term accounts for the excluded volume effect.

Figure 1
Mesoscopic model for G4. Left: Scheme of a parallel G4 assembly, taken as exemplary structure, 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 ...

The dynamics of the system is obtained from the overdamped equations of motion of the j-th bead of the G4 (j = 1,..., 21)


and for the i-th bead representing the ions (i = 1, 2)


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 107 time steps, from which, the first 106 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.

Figure 2
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 ...

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 plasmid11. 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 ceiling26. 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 conservation27. 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.


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 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.

Figure 3
Melting curves in terms of different magnitudes as a function of the dimensionless temperature T* (thermal energy): (a) Average coordination number of each ion Co I. I 1 and I 2 stand for each of the ions, whereas δ = 0 indicates ...

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 G420.

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 Tm. In each trajectory, we use the mean of the two energies TI1 and TI2 at which the coordination Co I of ions 1 and 2 respectively drop below the value 0.5, namely TI1/I2. To account for the stochastic effects in the unfolding, we repeated the melting simulations N s = 100 times from which we set Tm=12(TI1+TI2)=TI1/I2, 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 COI: 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 Tm=0.4875, represented by a solid line in each panel of Fig. 4.

Figure 4
Distribution of temperature Tmelt at the melting point calculated from different magnitudes over 100 realisations. (a) Temperature at which the coordination of the ions Co I goes below 0.5; (b) Temperature at which the coordination of the planes ...

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: Eu=D0kBTm/Tm=2.33kBTroom, (T room = 298K). The unit of force is F u = E u/l u = 14.7 pN. 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 ions20. 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 direction is calculated as F = k A(x A(t)  x 1(t)), 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.

Figure 5
Scheme of the pulling simulations. The first and last nucleotides in the G4 chain are represented by black beads and labeled as 1 and 21, respectively. Bead 21 is fixed to the origin of coordinates, whereas bead 1 is attached to a harmonic spring k A ...

The unfolding force value F u obtained in the all-atom simulations is in the order of 102 pN20, one order of magnitude higher than previous experiments of G4 mechanical unfolding8,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 6d and e show three examples of experimental curves, and three examples of simulation curves, respectively.

Figure 6
(a) Force as a function of the extension x 1 with the mesoscopic model at different pulling velocities. (b) Force as a function of the position of the pulling spring x A with the mesoscopic model at different pulling velocities. The velocity is given ...

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 7
Diagrams of the structural changes in the G4 as a function of the position of the pulling spring (vt + x 0) at different velocities. (a) Distance between beads 1 and 21. (b) Radius of gyration. (c,d) Average coordination number of the ...

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 sampling29 and the weighted histogram analysis method30. 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.

Figure 8
Potential of Mean Force (PMF) calculated from umbrella sampling and the weighted histogram analysis method.

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.0788lu/tu  0.051 nm/tu. Equaling this value to the experimental one v exp = 11.8 nm/s, we get the time unit to t u = 0.0043s.

Figure 9
Distribution of the unfolding force over 100 realisations at different velocities.
Figure 10
Mean value of the unfolding force as a function of r = k A v and fitting to the kinetic models (note that the three lowest values of the loading rates are excluded from the DHS  fitting because in this region some rebinding ...

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 model33. 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: feq=2kAG+. 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:



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 E1(z)=zessds is the exponential integral which is interpolated by e z E 1(z)  ln(1 + e γ/z). ku(f)=k0exp(β(fxu0.5kAxu2)) is the force dependent unfolding rate from the Bell model. In our simulations the pulling is performed with a harmonic spring and then r=dF(t)/dt=ddtkA(vt+xA(0)x1(t))kAv.

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 model32. 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.

Table 1
Fit parameters obtained from the mean unfolding force vs the loading rate.

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 × 105. According to the behaviour of d 1−21 as a function of time we can split the folded conformations in two main elongation ranges: (i) 0.8 < d 1−21 < 2.2 nm (t < 0.05 × 105, visible in the right inset of the figure), and (ii) 2.2 < d 1−21 < 2.9 nm (0.05 × 105 < t < 9 × 105 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 × 105 < t < 9 × 105). 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.

Figure 11
Pulling simulation at E = 0. Top: End to end distance d 1−21 and force F as a function of the time. The dashed green lines mark the equilibrium value of d 1−21 at v = 0. The solid black lines delimit the ...

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.


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 G420, on the melting temperatures of DNA G428, 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 experiments7.

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 repeats11, like those existing in human telomeric sequences.


Work supported by the Spanish Ministry of Economy and Competitiveness (MINECO), grant No. FIS2014-55867, co-financed by FEDER funds. We also thank the support of the Aragón Government and Fondo Social Europeo to FENOL group. Work in J.R.A.-G. laboratory was supported by a grant from MINECO, No. MAT2015-71806-R).

Author Contributions

Author Contributions

A.E.B.-P. conducted the simulations, A.F. and F.F. designed the project, I.G. conducted the experiment, J.R.A.-G. designed the experiment, A.F., F.F., A.E.B.-P., and J.R.A.-G. analysed the results. All authors reviewed the manuscript.


Competing Interests

The authors declare that they have no competing interests.


Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


1. Arias-Gonzalez JR. Single-molecule portrait of DNA and RNA double helices. Integr. Biol. 2014;6:904. doi: 10.1039/C4IB00163J. [PubMed] [Cross Ref]
2. Burge S, Parkinson GN, Hazel P, Todd AK, Neidle S. Quadruplex DNA: sequence, topology and structure. Nucleic Acids Res. 2006;34:5402. doi: 10.1093/nar/gkl655. [PMC free article] [PubMed] [Cross Ref]
3. Lam EY, Beraldi D, Tannahill D, Balasubramanian S. G-quadruplex structures are stable and detectable in human genomic DNA. Nat. Commun. 2013;4:1796. doi: 10.1038/ncomms2792. [PMC free article] [PubMed] [Cross Ref]
4. Siddiqui-Jain A, Grand CL, Bearss DJ, Hurley LH. Direct evidence for a G-quadruplex in a promoter region and its targeting with a small molecule to repress c-MYC transcription. Proc. Natl. Acad. Sci. USA. 2002;99:11593. doi: 10.1073/pnas.182256799. [PubMed] [Cross Ref]
5. Endoh T, Sugimoto N. Mechanical insights into ribosomal progression overcoming RNA G-quadruplex from periodical translation suppression in cells. Sci. Rep. 2016;6:1. doi: 10.1038/srep22719. [PMC free article] [PubMed] [Cross Ref]
6. Hänsel-Hertsch R, Di Antonio M, Balasubramanian S. DNA G-quadruplexes in the human genome: detection, functions and therapeutic potential. Nat. Rev. Mol. Cell Biol. 2017;18:279. doi: 10.1038/nrm.2017.3. [PubMed] [Cross Ref]
7. de Messieres M, Chang JC, Brawn-Cinani B, La Porta A. Single-molecule study of G-quadruplex disruption using dynamic force spectroscopy. Phys. Rev. Lett. 2012;109:058101. doi: 10.1103/PhysRevLett.109.058101. [PubMed] [Cross Ref]
8. Koirala D, et al. A single-molecule platform for investigation of interactions between G-quadruplexes and small-molecule ligands. Nat. Chem. 2011;3:782. doi: 10.1038/nchem.1126. [PMC free article] [PubMed] [Cross Ref]
9. Long X, et al. Mechanical unfolding of human telomere G-quadruplex DNA probed by integrated fluorescence and magnetic tweezers spectroscopy. Nucleic Acids Res. 2013;41:2746. doi: 10.1093/nar/gks1341. [PMC free article] [PubMed] [Cross Ref]
10. Ghimire C, et al. Direct Quantification of Loop Interaction and pi-pi Stacking for G-Quadruplex Stability at the Submolecular Level. J. Am. Chem. Soc. 2014;136:15544. doi: 10.1021/ja503585h. [PubMed] [Cross Ref]
11. Garavís M, et al. Mechanical Unfolding of Long Human Telomeric RNA (TERRA) Chem. Commun. 2013;49:6397. doi: 10.1039/c3cc42981d. [PubMed] [Cross Ref]
12. Fonseca Guerra C, Zijlstra H, Paragi G, Bickelhaupt FM. Telomere structure and stability: covalency in hydrogen bonds, not resonance assistance, causes cooperativity in guanine quartets. Chemistry-A European Journal. 2011;17:12612. doi: 10.1002/chem.201102234. [PubMed] [Cross Ref]
13. Yurenko YP, Novotn J, Sklen V, Marek R. Exploring non-covalent interactions in guanine-and xanthine-based model DNA quadruplex structures: a comprehensive quantum chemical approach. Phys. Chem. Chem. Phys. 2014;16:2072. doi: 10.1039/C3CP53875C. [PubMed] [Cross Ref]
14. Poudel L, et al. Implication of the solvent effect, metal ions and topology in the electronic structure and hydrogen bonding of human telomeric G-quadruplex DNA. Phys. Chem. Chem. Phys. 2016;18:21573. doi: 10.1039/C6CP04357G. [PubMed] [Cross Ref]
15. Li MH, Luo Q, Xue XG, Li ZS. Toward a full structural characterization of G-quadruplex DNA in aqueous solution: Molecular dynamics simulations of four G-quadruplex molecules. J. Mol. Struct-Theochem. 2010;952:96. doi: 10.1016/j.theochem.2010.04.035. [Cross Ref]
16. Islam B, et al. Conformational dynamics of the human propeller telomeric DNA quadruplex on a microsecond time scale. Nucleic Acids Res. 2013;41:2723. doi: 10.1093/nar/gks1331. [PMC free article] [PubMed] [Cross Ref]
17. Stadlbauer P, Krepl M, Cheatham TE, Koca J, Sponer J. Structural dynamics of possible late-stage intermediates in folding of quadruplex DNA studied by molecular simulations. Nucleic Acids Res. 2013;41:7128. doi: 10.1093/nar/gkt412. [PMC free article] [PubMed] [Cross Ref]
18. Li H, Cao E, Gisler T. Force-induced unfolding of human telomeric G-quadruplex: a steered molecular dynamics simulation study. Biochem. Bioph. Res. Co. 2009;379:70. doi: 10.1016/j.bbrc.2008.12.006. [PubMed] [Cross Ref]
19. Yang C, Jang S, Pak Y. Multiple stepwise pattern for potential of mean force in unfolding the thrombin binding aptamer in complex with Sr2+ J. Chem. Phys. 2011;135:225104. doi: 10.1063/1.3669424. [PubMed] [Cross Ref]
20. Bergues-Pupo AE, Arias-Gonzalez JR, Morón MC, Fiasconaro A, Falo F. Role of the central cations in the mechanical unfolding of DNA and RNA G-quadruplexes. Nucleic Acids Res. 2015;43:7638. doi: 10.1093/nar/gkv690. [PMC free article] [PubMed] [Cross Ref]
21. Linak MC, Tourdot R, Dorfman KD. Moving beyond Watson-Crick models of coarse grained DNA dynamics. J. Chem Phys. 2011;135:205102. doi: 10.1063/1.3662137. [PubMed] [Cross Ref]
22. Rebi M, Mocci F, Laaksonen A, Ulin J. Multiscale simulations of human telomeric G-quadruplex DNA. J. Phys. Chem. B. 2014;119:105. doi: 10.1021/jp5103274. [PubMed] [Cross Ref]
23. Stadlbauer P, et al. Coarse-Grained Simulations Complemented by Atomistic Molecular Dynamics Provide New Insights into Folding and Unfolding of Human Telomeric G-Quadruplexes. J. Chem. Theory Comput. 2016;12:6077. doi: 10.1021/acs.jctc.6b00667. [PubMed] [Cross Ref]
24. Parkinson GN, Lee MP, Neidle S. Crystal structure of parallel quadruplexes from human telomeric DNA. Nature. 2002;417:876. doi: 10.1038/nature755. [PubMed] [Cross Ref]
25. Bhattacharya D, Arachchilageand GM, Basu S. Metal Cations in G-Quadruplex Folding and Stability. Frontiers in Chemistry. 2016;4:38. doi: 10.3389/fchem.2016.00038. [PMC free article] [PubMed] [Cross Ref]
26. de Lorenzo S, Ribezzi-Crivellari M, Arias-Gonzalez JR, Smith SB, Ritort F. A Temperature-Jump Optical Trap for Single-Molecule Manipulation. Biophys. J. 2015;108:2854. doi: 10.1016/j.bpj.2015.05.017. [PubMed] [Cross Ref]
27. Smith SB, Cui Y, Bustamante C. Optical-trap force transducer that operates by direct measurement of light momentum. Methods Enzymol. 2003;361:134. doi: 10.1016/S0076-6879(03)61009-8. [PubMed] [Cross Ref]
28. Mergny JL, Phan AT, Lacroix L. Following G-quartet formation by UV-spectroscopy. FEBS letters. 1998;435:74. doi: 10.1016/S0014-5793(98)01043-6. [PubMed] [Cross Ref]
29. Torrie GM, Valleau JP. Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling. J. Comput. Phys. 1977;23:187. doi: 10.1016/0021-9991(77)90121-8. [Cross Ref]
30. Kumar S, Bouzida D, Swendsen RH, Kollman PA, Rosenberg JM. The weighted histogram analysis method for free-energy calculations on biomolecules I. The method. J. Comput. Chem. 1992;13:1011. doi: 10.1002/jcc.540130812. [Cross Ref]
31. Evans E, Ritchie K. Dynamic strength of molecular adhesion bonds. Biophys. J. 1997;72:1541. doi: 10.1016/S0006-3495(97)78802-7. [PubMed] [Cross Ref]
32. Dudko OK, Hummer G, Szabo A. Intrinsic rates and activation free energies from single-molecule pulling experiments. Phys. Rev. Lett. 2006;96:108101. doi: 10.1103/PhysRevLett.96.108101. [PubMed] [Cross Ref]
33. Friddle RW, Noy A, De Yoreo JJ. Interpreting the widespread nonlinear force spectra of intermolecular bonds. Proc. Natl. Acad. Sci. 2012;109:13573. doi: 10.1073/pnas.1202946109. [PubMed] [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group