Cooperativity-Dependent Folding of Single-Stranded DNA

The folding of biological macromolecules is a fundamental process of which we lack a full comprehension. Mostly studied in proteins and RNA, single-stranded DNA (ssDNA) also folds, at physiological salt conditions, by forming non-specific secondary structures that are difficult to characterize with biophysical techniques. Here we present a helix-coil model for secondary structure formation, where ssDNA bases are organized in two different types of domains (compact and free). The model contains two parameters: the energy gain per base in a compact domain, $\epsilon$, and the cooperativity related to the interfacial energy between different domains, $\gamma$. We tested the ability of the model to quantify the formation of secondary structure in ssDNA molecules mechanically stretched with optical tweezers. The model reproduces the experimental force-extension curves in ssDNA of different molecular lengths and varying sodium and magnesium concentrations. Salt-correction effects for the energy of compact domains and the interfacial energy are found to be compatible with those of DNA hybridization. The model also predicts the folding free energy and the average size of domains at zero force, finding good agreement with secondary structure predictions by Mfold. We envision the model could be further extended to investigate native folding in RNA and proteins.


Introduction
Single-stranded nucleic acids participate in a myriad of processes [1,2]. Bases in single-stranded DNA (ssDNA) and ssRNA tend to form base pairs of different geometries and stabilities, from the most stable Watson-Crick to the weaker Hoogsteen or wobble base pairs [3]. Base pairs are stabilized by hydrogen bonds and stacking interactions [4,5] that induce high-order tertiary structures [6,7]. Biologically active ssRNAs fold into compact domains, stabilized by base-pairing interactions, connected by flexible unpaired regions facilitating RNA conformational changes [8]. Although DNA is found in a double-stranded form in the cell genome, ssDNA is generated as an essential intermediary during every aspect of DNA metabolism, replication [9,10,11], transcription [12], and repair [13,14,15]. The tendency of ssDNA to form secondary structures interferes with and modulates the action of DNA processing enzymes.
Understanding how molecular interactions drive molecular folding is a fundamental problem in biology, from RNA [16] to proteins [17] but also in naturally occurring [18,19] and synthetic [20] DNA structures. After decades of research, conflicting views on how nucleic acids and proteins fold remain [19,21,22,23,24,25]. Besides the key biological role of ssDNA, understanding how ssDNA folds might shed light on the general mechanisms behind molecular folding in nucleic acids and proteins due to the ubiquitous presence of hydrogen-bonded and stacking interactions.
Secondary structure formation in nucleic acids can be investigated with pulling experiments with optical and magnetic tweezers [26], which allow stretching individual molecules by applying mechanical forces in the pN range and measuring their force-extension curves (FECs) [27,28,29]. These ex-periments show that ssDNA behaves as a semiflexible polymer, with FECs described using ideal elastic polymer models (such as the Worm-Like Chain -WLC-model [30]). However, at low forces ( 15 pN) FECs deviate from the ideal elastic behavior, showing consistent compaction in the form of a force shoulder [31,32]. This compaction indicates that ssDNA folds via secondary structure formation quantified by the reduction in the molecular extension.
We introduce a helix-coil model for ssDNA folding. In the model, bases along the ssDNA chain belong to two types of domains, compact (C) and free (F). C-domains contain bases forming secondary structure, whereas Fdomains contain unfolded free bases. Originally proposed by Zimm-Bragg and Lifson for protein folding [33,34], the helix-coil model has been extended to DNA hybridization [35,36,37,38,39] and RNA folding [40]. Surprisingly though, and to the best of our knowledge, the helix-coil model has never been used before for non-specific secondary structure formation in DNA and RNA.
We performed pulling experiments with optical tweezers on individual ss-DNA molecules over a wide range of experimental conditions, varying molecular length and ionic strength, to test our model. The model reproduces the FECs of individual ssDNA molecules over two-three decades of molecular lengths and sodium and magnesium salt concentrations. In comparison with more complex models [31,41,42], the one presented here is described in terms of two parameters: the energy per base forming C-domains (ǫ) and the interfacial energy between adjacent domains (cooperativity γ). Previous studies addressed the effect of salt on secondary structure formation in ss-DNA [28,32], yet the effect of molecular length has never been addressed. The helix-coil model permits us to determine the finite-size corrections to the free energy of formation of folded ssDNA and the average size of domains at zero force. Despite its simplicity, the model is compatible with Mfold secondary structure predictions, showing its power for predicting cooperativitydependent folding of ssDNA.

2
Materials and Methods

Optical tweezers setup
Experiments were carried out using a miniaturized dual-beam optical tweezers setup described in [43,44]. Briefly, two tightly focused counter-propagating laser beams (P=200 mW, λ=845 nm) created a single optical trap. A DNA molecule was tethered between two micron-sized polystyrene beads, one captured in the optical trap and the other held on top of a glass micropipette [ Fig. 1(a)]. dsDNA handles were labeled with Biotin or Digoxigenin (See Sec. S1 in the Supplemental Material at [URL will be inserted by publisher] for details) to bind selectively to Streptavidin (2.1 µm Kisker Biotech) and anti-Digoxigenin coated beads (3.0-3.4 µm Kisker Biotech). The experiments were performed in a microfluidics chamber. The force exerted on the optically trapped bead was determined by directly measuring the change of light momentum using Position Sensitive Detectors (PSDs). The position of the optical trap was determined by diverting ∼8% of each laser beam to a secondary PSD. The instrument has a resolution of 0.1 pN and 1 nm at a 1 kHz acquisition rate.

DNA substrates
ssDNA easily adsorbs non-specifically on surfaces which makes it difficult to manipulate and to control the number of ssDNA nucleotides that are effectively stretched. In order to avoid this problem, we used DNA hairpins that were unfolded mechanically to generate ssDNA (see Sec. 2.4 for details). In addition, the position at which mechanical unwinding of the hairpin starts can be used as a fiducial marker for alignment of the force-extension curves (FECs, see below). We used eight hairpins, with different sequences and lengths spanning from ∼ 120b to ∼ 14kb, named H N , N being the number of ssDNA bases released by the hairpin unfolding. We performed two different preparations for the DNA substrates (See Fig. S1 in the Supplemental Material at [URL will be inserted by publisher] for details). The shorter hairpins (120b and 204b) were synthesized by annealing and ligating a different set of oligonucleotides (See Sec. S2 in the Supplemental Material at [URL will be inserted by publisher] for their sequences). The longer hairpins (from 700b to 13680b) were synthesized using a long DNA fragment obtained either by a PCR amplification or by a digestion of the linearized λ-phage DNA. All hairpins were annealed to 29bp dsDNA handles that were used as spacers.
Labeling of the handles was achieved by a Digoxigenin/Biotin tailing using a terminal transferase (See Sec. S1 in the Supplemental Material at [URL will be inserted by publisher] for details).

ssDNA generation and determination of its molecular extension
Several methods have been developed for obtaining ssDNA [45,46,47,48]. In this work, ssDNA substrate was generated from a DNA hairpin using the blocking oligonucleotide method described in Ref. [32]. Initially, the DNA hairpin was tethered between the two polystyrene beads. An oligonucleotide of ∼ 25 − 30b, complementary to the loop region of the hairpin, was flowed into the chamber (See Table S2 in the Supplemental Material at [URL will be inserted by publisher] for the sequences). Above 15 pN the hairpin was fully unzipped [ Fig. 1(b)], allowing the annealing of the complementary oligonucleotide to the exposed hairpin-loop region. Re-annealing of the hairpin at forces lower than 15 pN was prevented by the annealed oligonucleotide, which worked as a kinetic block. Then, the excess of oligonucleotide was washed out of the chamber to prevent any possible interference with the secondary structure [32]. Force-distance curves were recorded by repeatedly moving the optical trap up and down at ∼ 100 nm/s) between a maximum and a minimum force (∼ 1 pN and ∼ 40 pN), and measuring the force and the distance λ between the trap center and the tip of the micropipette [ Fig. 1(b)]. About 5-10 stretching (force increasing) and releasing (force decreasing) curves were recorded for several tethers (∼ 5-10 in each condition). At low forces (∼ 1 pN) the ssDNA attached non-specifically to the bead surface. However, when low forces were avoided (forces above ∼ 1 pN), no hysteresis was observed between stretching and releasing curves, showing that pulling curves were reversible.
The molecular extension of ssDNA X ssDNA (f ) is related to the extension  tal setup for generating ssDNA. A DNA hairpin, bearing a functionalized dsDNA handles at each strand is tethered between two micron-sized beads; by moving the optical trap the molecule is stretched and the unfolding of the hairpin is observed at f ∼ 15 pN (b) Stretching (red) and releasing (gray) force-distance curves of the hairpin H 7138 in the presence of an oligonucleotide complementary to the loop region. Binding of the oligonucleotide blocks the rezipping of the hairpin, leading to the generation of ssDNA (Sec. 2.4) and secondary structure formation (grey curve). The black line shows the elastic WLC model (Sec. 2.5). (c) ssDNA FECs obtained for the different hairpins (lengths ranging from 120b to ∼14000b, from yellow to dark green) by converting trap-distance (λ) into molecular extension (X ssDNA ) (Sec. 2.4). The x axis is represented in log scale. (d) FECs in (c) normalized by the number of bases (N ). Color code as in panel (c). The secondary structure emerges at forces below ∼ 10 pN, as a shoulder deviating from the WLC [Eq.(2)] (black line). For N 500 FECs collapse into a single curve. The force shoulder for the 204b and 120b molecules is systematically lower, showing finite-length effects. Error bars were obtained by averaging 5-10 releasing cycles of 5-10 different tethers for each hairpin. Experiments were performed at 25°C in 10 mM MgCl2. measured experimentally, λ [ Fig. 1(a)], as follows: where x h (f ) = x h1 (f )+x h2 (f ) is the extension of the two 29bp handles (characterized in [49]) x 0 is an arbitrary shift of the trap position (see App. A.1 for details) and x b (f ) is the displacement of the bead from the center of the optical trap. The latter is given by x b = f /k b with the trap stiffness k b measured with a power spectrum analysis (See Sec. S3 in the Supplemental Material at [URL will be inserted by publisher] for details). The force plotted versus the extension X ssDNA defines the force-extension curve (FEC). The binding of the oligonucleotide to the loop region introduces a correction to Eq. (1) due to the hybridized dsDNA segment. We have verified that this correction is negligible for molecules of 200b or larger (See Sec. S4 in the Supplemental Material at [URL will be inserted by publisher] for details). However, for the shortest 120b hairpin, the correction would be significant. Therefore, we used an alternative approach, the so-called twobranches method [50]. This method implies the use of a longer loop (20b) to slow down the loop formation process, allowing measuring the ssDNA extension between ∼ 3-4 pN and ∼ 15 pN (See App. A.2). The two branches method compares differences in extension between the folded and unfolded branches, without considering corrections due to the optically trapped bead and the handles. Instead, in the blocking oligonucleotide method bead and handles extensions must be subtracted, leading to higher uncertainties.
The FECs of the ssDNA for the eight different molecules studied are shown in Fig. 1(c). The FECs are scaled by dividing the measured extension over the number of nucleotides of each sequence, x ssDNA = X ssDNA /N, as shown in Fig. 1(d).

2.5
Elastic models 2.5.1 Worm-like chain Different polymer models have been proposed to describe the ideal elasticity (i.e. without secondary structure) of nucleic acids [27,51,52]. A widely used model is the Worm-Like Chain (WLC) model with the following interpolation formula between the low-and high-force regimes [51]: where x is the polymer extension at a given applied force f , l is the contour length per base, p is the persistence length and N is the number of monomers of the ideal chain (with total contour length given by L = N · l). Here we used Equation (2) to reproduce the ideal ssDNA elastic behavior at 10 mM MgCl 2 (Fig. 3), with elastic parameters taken from Ref. [32] (p = 0.75 nm, l = 0.69 nm). Previous to fitting the data to the helix-coil model the measured ssDNA extensions were corrected by 5% to match the WLC model at f = 15 pN (See Table S4 in the Supplemental Material at [URL will be inserted by publisher] for the used factors).

Saleh formula
Previous studies have shown that, at low forces, the measured ssDNA FECs deviate from the WLC predictions [53,54]. These deviations are mainly due to excluded volume effects, and they become specially relevant at low salt concentrations. Saleh and collaborators have proposed different formulae for the regimes of low and high forces in sodium conditions [53,54,55]: at low forces the molecular extension shows a power law dependence with force due to the excluded volume effects; at high forces the elastic response is described by a WLC model with a salt-dependent internal electrostatic tension; finally, a logarithmic dependence interpolates the two force regimes (intermediate force regime). We have combined these results in a single elastic response function covering the three force regimes (low, intermediate, high). The fitting function will be denoted as Saleh-interpolating formula and is given by [53,54,55], if f ≤ f 1 ; with f 1 = a (1000c) γ , a = 0.06 pN, γ = 0.054, l 2 = 0.7 nm, l 1 = l 2 / (2(γ log(f 2 /f 1 ) + 1)) (continuity at f 2 ) and The electrostatic tension f el reads: with an average charge length b = 1.14 nm [53]. l B and κ are the Bjerrum length and the inverse of the Debye length, respectively. Note that the dependencies described in Eq. 3 have been only tested in sodium conditions. However, Eq. 3 can be extrapolated to magnesium by using the 100:1 salt rule for the equivalence of the screening ion effect of monovalent and divalent ions [32,56]. In the App. B we show results of the ssDNA elasticity in presence of glyoxal (which removes secondary structure). At high salts, the WLC and the Saleh formula are compatible with the experimental data, however there are differences at low salts where the Saleh formula reproduces better the data. Consequently, we have used the Saleh formula (Eq. 3) to describe the ideal ssDNA elasticity at different salt conditions (Fig. 4). Comparison of the results for the fitting parameters of the helix-coil model, ǫ and γ, using the WLC and the Saleh formula are shown in Fig. S4 in the Supplemental Material [URL will be inserted by publisher].
The measured ssDNA extensions were corrected by 5% to match the Saleh formula model at f = 15 pN (See Tables S5 and S6 in the Supplemental Material at [URL will be inserted by publisher] for the used factors).

2.6
Helix-coil model for secondary structure The FECs of individual ssDNA molecules deviated from the ideal polymer behavior due to the formation of secondary structures at forces below ∼ 15 pN (Figs. 1B-D). Here, we develop a helix-coil cooperative model that describes the formation of secondary structures and explains the experimental FECs within full force range under study (0-40 pN). We envision the ssDNA as forming a sequence of structured compact domains interspersed by chains of free monomers (Fig. 2). These two domain types are denoted by C (compact) or F (free chain). The ssDNA is then modeled as a chain of N monomers (bases) that can be in two states σ i = ±1: if σ i = −1, the monomer is part of a C-domain; otherwise (σ i = +1), belongs to an F-domain. F-type domains pulled by a force f follow the WLC model or the Saleh formula (Sec. 2.5). C-type domains are stabilized by an average energy per monomer ǫ. The model includes cooperativity through an interfacial energy parameter γ that penalizes (rewards) adjacent monomers in different (equal) domain type. Let N F = N i δ σ i ,1 and N C = N i δ σ i ,−1 be the total number of monomers of F-type and C-type, respectively, fulfilling N = N F + N C . The model Hamiltonian is given by: The first two terms account for the total energy of the compact domains pulled at a force f ; M C (f ) is the number of C-domains, each contributing with an extension d C (f ); the third term is the stretching energy of F-type monomers with the extension per monomer x F (f ), obtained by Eq. (2) or Eq. (3) ( [57]); and the last term stands for the interfacial energy between adjacent monomers. A schematic representation of this model is shown in Fig. 2. At a given force, N F monomers with σ = +1 contribute by x F (f ) to the total extension. The rest of monomers, N C , are distributed into M C compact domains, each one contributing by d C (f ) to the full extension. Hence, the total molecular extension is given by: The Hamiltonian in Eq. (5) is analytically solvable (See Sec. 2.6), leading to the extension per monomer, X ssDNA /N, where the fraction of free monomers, φ F (f ) = N F (f )/N, is given by 1 : with β = 1/k B T . A and B quantities are defined by Eq. (7) can be used to obtain the fraction of bases in C-domains, φ C (f ) = The number of C-domains is given by: .

(11)
The average number of monomers per C-domain reads as: In Eq. (6), x F (f ) is given by the WLC model [Eq. (2)] or the Saleh formula [Eq. (3)], with the parameters described in Sec. 2.5. For the extension d C (f ), we consider each domain as a rigid object of length δ that aligns along the stretching direction (App. C). For example, a hairpin-like C-domain initiating with a stem would have a length equal to the diameter of a B-DNA doublehelix, δ ∼ 2 nm [58]). We checked that this parameter does not affect the predictions of the model significantly (App. C). Therefore, we neglected the domain orientation term (δ = 0) in the analysis and fitted our data to the simpler expression: We fitted the model [Eq. (13)] to the experimental FECs with ǫ and γ as fitting parameters between 2 ≤ f ≤ 15 pN. Lower forces (f < 2 pN) were not considered in the fit to avoid any influence of incorrect tethering (specially in molecules with N 1000b). Importantly, the model reproduces the experimental FECs throughout the whole force range 0 < f ≤ 15 pN in the different conditions tested (different length and salts).

Force-extension curves
For ssDNA measurements a DNA hairpin was fully unzipped by applying a force above 15 pN at the 3' and 5' extremities of the molecule, exposing the loop region. Unzipping of the hairpin was performed in the presence of an oligonucleotide complementary to the loop (Sec. 2.4). Fig. 1(b) shows the pulling curve (red) corresponding to the unzipping of the H 7138 hairpin with the characteristic sawtooth pattern [46]. This is followed by the relaxation curve (grey), after the oligonucleotide hybridizes to the loop region, corresponding to the ssDNA FEC. It differs from the ideal elastic curve (WLC model shown in black) at low forces, where ssDNA compacts and the extension shortens showing the characteristic shoulder of secondary structure formation [32,59]. It is well known that DNA sequence (i.e. GC content) and monovalent and divalent cations (which screen the negatively charged DNA phosphate backbone) affect secondary structure formation [28,32,60]. However, the effect of ssDNA length has never been studied. The nature of the compact structure is an open question due to its highly disordered and dynamic features. A dependence of the FEC with the ssDNA length might be the signature of long-range interactions along the chain, giving information about the nature of the secondary structure. γ(kcal/mol)

Effects of molecular length on secondary structure
We measured the FECs of eight different ssDNA molecules [ Fig. 1(c)] with lengths spanning from 120 to 13680b in 10 mM MgCl 2 (Sec. 2.3) and similar GC content, ∼ 50% 2 . FECs are shown in Fig. 1(d) with the extension rescaled by the number of ss bases of each molecule. Molecules with N 500b approximately collapse into a single curve, deviating from the ideal elastic behavior at forces f 10 pN and showing a remarkable shoulder. Small differences are due to different GC content (from 44% to 53%) which correlates with the height of the force shoulder observed in the FECs (see  below and App. E). In contrast, for N 500b FECs do not collapse into the same curve and the height of the shoulder decreases with N, indicating finitelength effects. This might be due to the depletion of secondary structure arrangements in the regime N 500b.
As shown in Fig. 3(a), the FECS for the studied molecules (yellow to dark green circles for increasing lengths) are well reproduced by Eq. (13) in the proposed model (green curves). The best-fitting values for the average energy per base of C-domains, ǫ, and the cooperativity parameter, γ, as a function of the molecular length are shown in Figs. 3(b)-(c). ǫ is found to follow a phenomenological 1/N-correction, with ǫ 0 = 0.18 (1)  Equation (14) quantifies corrections to the stabilization energy of C-domains due to the free ends of the polymer. In generic statistical physics models, corrections to the equilibrium free energy are determined by the surface to volume ratios, being on the order of 1/N for the one-dimensional helix-coil model. For shorter molecules (N 500b), ǫ starts deviating from the saturation value ǫ 0 . We stress that the FECs in Fig.3(a) cannot be fitted to the model with γ = 0, as shown in Fig. S5 in the Supplemental Material at [URL will be inserted by publisher]. Therefore, cooperativity is key to reproduce the experimental results, suggesting the suitability of the helix-coil model.

3.3
Salt dependence of secondary structure Values (kcal/mol) parameters ǫ, γ and the FECs is shown in Figs. S4 and S6 in the Supplemental Material at [URL will be inserted by publisher]. Salt-induced stabilization of secondary structure is captured by the model as an increase in ǫ with NaCl (magenta) and MgCl 2 (blue) concentrations. Note that, for ǫ = 0 and zero force, F-and C-domains are equally probable meaning that half of the bases are of C-and F-type. ǫ values can be compared with the Watson-Crick (WC) base-pair energies of the nearest-neighbour (NN) model. Therefore, we define ǫ bp = 2ǫ as the equivalent of the NN base-pair energy for DNA hybridization. Similarly to the energies in the NN model [46,64], we find that ǫ bp depends logarithmically with salt concentration [ Fig. 4(c)]: with c being the salt concentration (in M units), ǫ 0 bp the value at 1M and m ǫ the salt correction parameter for secondary structure. Equation (15) is a phenomenological expression that follows from applying standard physical chemistry theories of thermodynamic activity to diluted ionic solutions. Fits of ǫ bp to Eq. (15) are shown in Fig. 4(c) as dashed lines. The results from the fits are ǫ 0 bp = 0.34(3) kcal/mol, m ǫ = 0.11(2) kcal/mol in sodium and ǫ 0 bp = 0.75(9) kcal/mol, m ǫ = 0.09(2) kcal/mol in magnesium. The values for ǫ 0 bp are ∼ 5 times smaller than the average NN base-pair energy (∼ 1.6 kcal/mol) and almost 3 times smaller than the most unstable NN bp (AT/TA∼ 0.85 kcal/mol at 1 M NaCl and 10 mM MgCl 2 ) [46,64]. Such a lower stabilization is expected for random sequences that lack full complementarity regions. In fact, C-domains might consist of aggregates of base pairs, but also large loops and mismatches, which decrease base-pairing stability [65,66]. In agreement with this, unzipping forces in dsDNA are ∼ 15 pN [ Fig. 1(b)] (see also [67]), much larger than the force at the shoulder in the FEC (∼ 5-10 pN). Interestingly, the slope m ǫ =0.11 (2) kcal/mol [Eq. (15)] in sodium is compatible with the homogeneous salt correction parameter of the unified oligonucleotide dataset [68] and with that derived from unzipping experiments [46] (m = 0.114 kcal/mol and m = 0.104 kcal/mol, respectively). For the MgCl 2 case, m ǫ = 0.09(2) kcal/mol is about twice the Mfold value [69] (m = 0.055 kcal/mol) but closer to the unzipping value [70] [m = 0.07(2) kcal/mol]. Note that Mfold assumes that the salt correction for sodium is exactly twice that of magnesium (m NaCl Finally, for the cooperativity term γ, we have also assumed the phenomenological logarithmic salt dependence of Eq. (14) [ Fig. 4(d)]: Notice that the switching of a base (from F to C or the opposite) involves an interfacial energy cost/gain of ±2γ [Eq. (5)]. Note the opposite sign of the salt corrections for ǫ and γ. In both cases, salt screens the electrostatic interactions by weakening phosphates repulsion. This leads to an increased base pairing stabilization and a (much weaker) decreased cooperativity.
A compendium of all results in this section is shown in Table 1 and the fitting parameters obtained with the WLC model are shown in Table S7 in the Supplemental Material at [URL will be inserted by publisher].

C-and F-domains organization
We use the proposed helix-coil cooperative model (Sec. 2.6) to investigate how bases are distributed between F-and C-domains. To this end, we compute the fraction of bases belonging to C-domains at each force f , φ C (f ) [Eq.   7) and (10) using a salt and length dependent ǫ and γ parameters as given in Table 1. (b) Average number of bases or size per C-domain, n C [Eq. (12), with salt dependent ǫ and γ parameters given in Table 1], as a function of force for 10, 100, 250 and 1000mM NaCl. The schematic depictions show a highly compacted form of the ssDNA at low forces that is stretched at higher forces, decreasing the size of the C-domains. The arrow indicates the threshold force f * . (c) n C dependence on salt concentration (magenta for NaCl and blue for MgCl2) at 0 and 5 pN. Symbols correspond to Eq. (12) using the best fitting values for ǫ and γ  Table 1. Shadowed areas show the uncertainty region of n C , from the errors in ǫ and γ. The schematic depictions show the negatively charged ssDNA (black line) screened by the cations (grey dots) in the solvent stabilizing C-domains.  Fig. 5(a) show φ C for each condition at two different forces, f = 0 and 5 pN. φ C (f ) presents its maximum value at 0 pN, φ C (0) = 0.5-0.9, depending on salt and length. φ C (f ) decreases until a threshold force is reached, f * . Above f * the secondary structure does not form (φ C (f * ) ∼ 0) and the experimental FEC matches the ideal elastic model prediction Eqs. (2)(3). f * depends on molecular length and salt concentration: for long molecules (N ≥ 500b) and high salts ([NaCl]=1000 mM, [MgCl 2 ]=10 mM), f * ∼ 10 pN, whereas for either short molecules or low salts ([NaCl]=10 − 50 mM, [MgCl 2 ]=0.5 mM), f * ∼ 5-7 pN. Note that φ C (0) is non-negligible, even at the lowest NaCl concentrations (10-50 mM NaCl, dashed lines in the insets of Fig. 5(a), middle). In spite of this, a shoulder in the FEC is not observed [ Fig. 4(a)], because ǫ is too small. Fig. 5(b) shows n C (f ) versus force obtained by fitting the model to the experimental data, and a similar trend to φ C (f ) is observed: it is maximum at zero force, n C (0), and decreases as force increases, until the threshold force is reached, f * , for which n C → 1 (the minimum value allowed by the model). The values of f * for φ C (f ) and n C (f ) are the same, indicating that secondary structure is disrupted above f * . As shown in Fig. 5(c), n C increases with sodium and magnesium for all forces. The values for n C range from ∼ 10b (low salt) to ∼ 30b (high salt). Overall, the model predicts that the applied force destabilizes the secondary structure, by reducing both the size and number of C-domains [sketch in Fig. 5 . For a given ssDNA sequence, Mfold computes folded structures of free-energy ∆G (with respect to the random coil) in a given range above the absolute minimum, defining the energy spectrum of that sequence. We analyzed each folded structure in the spectrum [ Fig. 6(a)], and assigned nucleotides as being of F-or C-type depending on whether they belong to a folded hairpin-like structure (C-domain) or they are part of a free linker (Fdomain). This permits us to calculate φ C , n C , by counting the number of bases forming C-domains (φ C ) and the number of bases per C-domain (n C ).
To compare our model predictions with Mfold, we derive an analytical expression for the free energy of the folded ssDNA at zero force from the Eq. (17) Fig. 6(c) corresponds to a 1004b random sequence with the same GC content as in H 13680 , since the Mfold analysis cannot be performed with very long sequences (larger than ∼5000b, [72]). The same results are obtained for different sequence lengths whenever the GC content is conserved. It is important to stress the strong dependence of the free energy ∆G with the GC content (see Fig. S7 in the Supplemental Material at [URL will be inserted by publisher]), which must be taken into account for a proper comparison. The match between the Mfold prediction and Eq. (17) works better when using the Saleh formula [Eq. (3)] as compared to the WLC [Eq. (2)] for the ideal ssDNA elasticity (see Fig. S7 in the Supplemental Material at [URL will be inserted by publisher]). The equivalent of Fig. 6(c) in magnesium is not predictive because Mfold uses a salt correction that is known to be inaccurate [70].
In addition, for short ssDNA molecules (N < 1000b), the helix-coil model and Mfold also predict similar values for φ C (0) and n C (0)  Fig. 6(a), right structure]. To test this, we have computed low-lying energy structures with Mfold by constraining the maximum distance between base pairs. For 150b distance, the values for φ C (0) and n C (0) predicted by Mfold for long molecules (N > 1000b) agree with our model [red circles in Figs. 6(d)-(e)], supporting this interpretation. As shown in Fig. 6(f)-(g), the helix-coil model also reproduces the Mfold predictions for φ C (0) and n C (0) as a function of the salt concentration for short enough molecules [Mfold predictions are obtained for the same 1004b sequence as in Fig. 6(c)].
It is remarkable that, despite its simplicity, the helix-coil model compares so well with Mfold. Mfold is a predictive tool that uses the NN model with energy parameters derived from large datasets obtained from DNA melting experiments. The excellent agreement between our model and Mfold reflects the fact that both build upon the energetics of hybridization. The relevant motifs in Mfold (base pairs, loops, single-stranded segments, etc.) are also the ultimate building blocks of the compact disordered structures. Interestingly, the combined effects of these blocks are effectively described by only two parameters (ǫ and γ) in the helix-coil model.

Discussion and Conclusion
We have developed a cooperativity-dependent folding model that reproduces experimental FECs of ssDNA over different molecular lengths and salt conditions. The model assumes that ssDNA folds forming compact domains (C) interspersed by free bases (F). While free bases are described by the ideal elastic model, C-domains are assumed to have a negligible extension. Compared to more complex theoretical models [31,71], the model presented here includes only two parameters: ǫ, the average energy per base in a C-domain, and γ the interfacial energy between adjacent domains (cooperativity term). With only two parameters, the model reproduces the experimental FECs of ssDNA spanning two decades of contour length (100 b N 14 kb) and three decades of sodium and magnesium concentrations (10 ≤ [NaCl] ≤ 1000 mM and 0.5 ≤ [MgCl 2 ] ≤ 10 mM). Remarkably, γ = 0 is unable to fit the data suggesting that cooperativity is essential for ssDNA folding, as shown in Fig. S5 in the Supplemental Material at [URL will be inserted by publisher]. The simplest interpretation of the cooperativity term is the nature of the stacking interaction itself. In nucleic acids two distinct forces stabilize a duplex: hydrogen bonding and stacking. Hydrogen bonding acts transversely to the phosphate backbone bringing bases close to each other. In contrast, stacking acts longitudinally along the backbone, which is the nat-ural direction of the single stranded polymer. It is this geometric alignment between the stacking force and the polymer direction that generates cooperativity. This explains why ssDNA and ssRNA show a cooperative plateau in the FECs [71,48] and many RNAs form tertiary structures with coaxial stacking between different helices (making RNA folding a cooperative process [72,73]).
To further check the applicability of the model, we carried out pulling experiments to predict the FECs on another ssDNA (∼ 8 kb, ∼ 50% GC content) obtained by unpeeling an 8kbp dsDNA [45,47] in a different optical tweezers instrument [27]. In App. F, we show the predicted FEC from the model using the ǫ, γ parameters derived in our study [ Fig. 3(b)-(c)] without any fitting procedure. The good match with the experimental data further validates our model's predictive power and the derived energy parameters.
In a previous work [32], we found a phenomenological formula for the FEC at different salt conditions. This formula included two parameters, a critical force (f † ) and a force-width parameter (δ), akin to ǫ and γ in the helix-coil model. The phenomenological formula gives comparable fits to the helix-coil model; however, the former has a limited scope at the level of physical interpretation, especially at low forces where their behaviors differ (as shown in Sec. S10 in the Supplemental Material at [URL will be inserted by publisher]).
A remarkable feature of the ssDNA FECs is the absence of force rips, which is due to the low values of ǫ ( 0.3kcal/mol) compared to k B T (∼0.6kcal/mol at T = 298K). For a random sequence, Brownian forces smear out any force rips appearing during the formation of compact domains. This gives rise to the smooth and reversible FECs observed in the experiments. For ǫ k B T we should observe FECs with specific sequence-dependent force rips as it is the case in unzipping experiments [ Fig. 1(b)] where base-pair energies are ∼ 2 − 4k B T . Still, sequence effects are observed in the height of the shoulder of the FEC that increases with GC content. Concomitantly, the value of ǫ increases too (App.D).
Regarding salt dependence effects, we observe that ǫ and γ have logarithmic salt corrections. The salt-correction parameter for ǫ in sodium [m ǫ = 0.11(2) kcal/mol] is compatible with that of dsDNA hybridization [m = 0.11(1) kcal/mol]. This indicates that monovalent cations play a similar role in secondary structure formation and dsDNA hybridization. The stabilization effect of magnesium appears to be slightly larger for ssDNA folding [m ǫ = 0.09(2) kcal/mol] as compared to duplex hybridization [m = 0.06 (2) kcal/mol]. This demonstrates the distinct role of divalent cations in bringing together distant nucleotides to fold ssDNA. Such a long-range effect should be less prominent for duplex formation where the stem always grows in the proximity of the hybridization junction.
Apparently, the interfacial energy, γ, decreases much weakly with the ionic strength, |m γ | ≪ m ǫ raising the question whether the logarithmic dependence of Eq. (16) is the right phenomenological description. It would be very interesting to extend thermodynamic activity theories to describe salt corrections to the cooperativity term. Using the less-accurate WLC model to fit the data, a systematic salt dependence in both sodium and magnesium that is compatible with Eq. (16) is obtained (see Fig. S6 in the Supplemental Material at [URL will be inserted by publisher] for details). This would be an indication that salt similarly screens the interaction between adjacent bases in secondary structure and dsDNA hybridization.
The model predicts the free energy of the folded structure and its nature at zero force (φ C (0) and n C (0)). Free energies are in agreement with Mfold predictions for all investigated lengths, while φ C (0) and n C (0) are not well captured for long molecules (N 1000b) due to the underestimation of large C-domains (∼ 100b) by our model. Interestingly, the model also predicts that ssDNA folds at zero force even for very short lengths N ∼ 10-100b at physiological conditions (10 mM MgCl 2 ) [i.e. φ C (0) ∼ 0.5 in Figs. 5(a) (top), 6(d)]. This prediction is of biological relevance because short ssDNA fragments (15 − 30b) play important roles during DNA metabolic processes.
The agreement between the helix-coil model and Mfold is remarkable. This suggests that Mfold includes cooperativity as an essential ingredient for secondary-structure prediction. Indeed cooperativity is implicit in the energy of the stem-loop structures formed by stacking adjacent base pairs in the NN model and other secondary motifs. Mfold does not account for base-pair mismatches and tertiary motifs in general. This limitation of Mfold excludes the possibility of predicting disordered structures made of a large number of small folded regions of low thermodynamic stability but entropically favorable because of their large number. On the other hand, our model does not include heterogeneity in the sequence while Mfold does, predicting large secondary domains. These two features make the most important differences between the helix-coil model and Mfold.
Thermal molecular folding (at zero force) cannot be directly measured with force spectroscopy. However techniques such as hydrogen exchange combined with NMR and mass spectrometry might grant access to naturally folded domains and folding kinetics [74,75]. The kind of structures expected in ssDNA should differ from those reported in native RNAs and proteins [19,25]. Still, the cooperativity mechanisms might be similar and ultimately related to the fundamental nature of electrostatic interactions [76,77].
It might be interesting to apply the current model to investigate ssRNA folding, where base stacking interactions are stronger. Future work should also extend this study to incorporate sequence disorder while keeping the simplicity of a few energy parameters. This might be useful for models of folding in native RNAs and proteins showing cooperativity-dependent sequential folding.  1), where x 0 is a shift of the trap position, λ. x 0 results from imposing λ = 0, at zero force by extrapolating the force-distance curve (FDC) in the initial slope (i.e. before the first rip) to zero force (left scheme in Fig. 7). For this procedure to work, tethers must be aligned along the pulling axis. To verify the tether alignment, we compare the pattern of the experimental unzipping curve with the theoretically predicted FDC in equilibrium [46], as shown in Fig 7. If both curves match, the alignment is correct.  Figure 7: Checking tether alignment. Experimental force-distance unzipping curves obtained for an aligned tether (left, blue) and a misaligned tether (right, orange). The agreement between the theoretical curve (purple) and the experimental one (blue)is a signature of correct tether alignment (left panel). Misaligned tethers usually lead to the lack of some of the initial characteristic unzipping rips (right panel). A schematic depiction of aligned and misaligned tethers are shown on top. For obtaining the theoretical curves, we use the NN energies for the base pair energies [46], the elastic parameters from [32] and a bead stiffness of k b ∼ 0.07 pN/nm.

A.2 Two-branches method
The two-branches method is used for obtaining the ssDNA FEC of the shortest molecule (H 120 ). As shown in Fig. 8, the molecule is unzipped and rezipped with pulling velocities of ∼ 100 nm/s. Two branches are defined depending on whether the hairpin is folded (F) or unfolded (U). The optical trap position in the F ( λ F ) and U (λ U ) branches, reads as: where x d is the extension of the folded hairpin and λ 0,F , λ 0,U are drift corrections.
x ssDNA (f ) can be obtained along the force range where the two branches  Figure 8: Two-branches method. Force-distance curves for the H 120 , with several stretching/releasing cycles. Experimental configurations corresponding to F and U branches are schematically shown in the top. The folded (λ F ) and unfolded (λ U ) branches are shown in red and blue, respectively, with unfolding and refolding events shown in black.
coexist (∼4-17 pN) from Eqs. 19 and 18: where λ 0 (t) stands for the drift correction (See Sec. S11 in the Supplemental Material at [URL will be inserted by publisher] for details).

B Experiments with glyoxal
Experiments with chemically treated ssDNA have been performed in order to check the characterization of the free ssDNA elasticity. Following references [28,53,54] we have carried out experiments with glyoxal using the following protocol. We first unzipped dsDNA in the absence of glyoxal to generate ssDNA. Since glyoxal prevents base pairing the blocking loop oligonucleotide was not necessary to prevent hairpin rezipping. Once unzipped we kept the force constant above the unzipping force (∼20 pN at 1M NaCl and 10mM MgCl2, ∼15pN at 10mM NaCl and 0.5mM MgCl2) and flowed glyoxal at 1M concentration.
We waited for 15 minutes for glyoxal to covalently attach to the exposed ssDNA bases. Glyoxal changes the index of refraction of the solvent, modifying force calibration. Therefore, once glyoxal has coated ssDNA we washed the chamber by flowing the original buffer to remove glyoxal and measure the glyoxal-ssDNA FECs. Figure 9 shows the results we obtained for different molecules at the lowest and highest salt conditions in sodium and magnesium. For comparison, we also show the results obtained without glyoxal. The differences between the FECs obtained in the two conditions (with and without glyoxal) is due to the presence or absence of secondary structure.  In Figure 10 we show the FECs averaged over the molecules at each condition together with three elastic models: an extensible-FJC model (with parameters from [32]), the WLC model in Eq. (2) and the Saleh formula in Eq. (3). The results show that at high salts the inextensible-WLC and the Saleh formula are both compatible with the data, however there are differences at low salts where the Saleh formula works better. This confirms the results shown in Fig. S4 in the Supplemental Material at [URL will be inserted by publisher] where the largest differences among the elastic models appear at the lowest salt concentrations.  Figure 10: Comparison between elastic models and the experimental FEC with glyoxal. Green points show the averaged re-scaled FECs of the H 7138 ssDNA molecule for different salt concentrations (from 2-3 different tethers). The continuous lines correspond to the elastic curve obtained from the WLC (blue, elastic parameters from [32]), extensible Freely Jointed Chain (magenta, elastic parameters from Ref. [32]) and orange (Saleh formula, described in Sec. 2.5.2).

C Hairpin orientation
A DNA hairpin has a helix diameter of d = 2 nm, corresponding to that of dsDNA (B-form) [78]. Under an externally applied force, the hairpin extension can be modelled as a single bond of length d that orients along the stretching direction [67,49,79]: where k B is the Boltzmann constant and T is the temperature. We investigated the effect of varying δ in our model. Starting from δ = 0 and upon increasing its value we found that the fitting energy parameters ǫ and γ remain almost unchanged (Fig. 11): ǫ and γ increase 3% when 0 ≤ δ ≤ 0.6 nm. We conclude that the effect of δ is very small in the obtained fitting parameters. For the sake of simplicity, we have used δ = 0nm throughout the analysis presented in this work.
Notice that our model, as presented in Sec.2.6, is only applicable for δ ≤ l ∼ 0.7 nm, otherwise for δ > l the model predicts, at sufficiently high forces, a structure made of alternating C and F-domains of only one base (φ C = 1/2, n C = 1). This is clearly incorrect as secondary structure is fully disrupted at sufficiently high forces. This inaccuracy of the helix-coil model is consequence of its simplicity, and might be corrected by imposing a minimum number of bases, n min , to form a C-domain such that n C > n min . With this restriction, the helix-coil model becomes too complicated and is not analytically solvable anymore. This fact and the negligible effect of δ on ǫ, γ (Fig. 11) explain why we took n min = 1 and δ = 0, keeping the simplicity of our model.

E GC content effect on secondary structure
Sequences with high GC content present the shoulder in the FEC at larger forces [28]. The molecules studied in this work (Sec. 3.2) present small differences in the GC content (44-53%). For long molecules (N 500b, squares), ǫ as a function of the GC content increases (Fig. 12), compatible with the larger stability of secondary structures with high GC content, as shown in Fig. 12. For shorter molecules (N 500b, yellow circles), ǫ decreases strongly, which we interpret as finite size effects (see main text).

F Model validation
To test the model, we used experimental FEC data of a 8022 bases ssDNA molecule obtained by unpeeling a dsDNA molecule using a different optical tweezers setup [27] (Fig. 13). In this method, one strand of the dsDNA molecule is attached by its 5' and 3'ends to the beads. Pulling the molecule above 80-100 pN (∼ 200 nm/s) induces the mechanical denaturalization of the double helix and promotes the releasing of the free strand (the one not attached to the beads) [45,47]. We use our model without any fitting parameters to predict the FEC of this molecule at the experimental conditions (10 mM MgCl 2 and 20 mM NaCl). The good match between the model prediction and the experimental data (Fig. 13) further validates the model. Best fitting values for ǫ as a function of the GC content of the sequence, for the 8 molecules studied. Long molecules (N 500b) are shown in squares, whereas short molecules (N 500b), deviating from the collapsed FEC trend observed in Fig. 1(d), are shown in circles. The linear fit (dashed line) performed with long molecules (700b ≤ N ≤ 14kb), shows an increase of ǫ with GC content, as expected from the higher stability of GC bps as compared to AT bps. Color code as in Fig 11.
where D is a constant, and therefore will not be considered in the calculations below. The partition function reads as: where tr is the trace and V is the transfer matrix given by: By diagonalizing V , we find the eigenvalues: The partition function [Eq. (27)] can then be written as: which, in the thermodynamic limit N → ∞ (λ − < λ + ), leads to:

G.1 Free-energy computation
From the partition function [Eq. (32)], we can compute the free-energy per monomer (base) as: with A, B and C given in Eqs. (23)(24)(25). The free-energy per monomer at zero force, G 0 N , is obtained by imposing f = 0 in Eq. (33): For a given DNA sequence, Mfold gives the free energy of formation of secondary structures with respect the random coil (i.e. configuration with all bases free). In order to compare the free-energy of the helix-coil model with Mfold predictions, we subtract to Eq. (34) the energy per base of the random coil at zero force [Eq. (5) with f = 0, σ i = 1, ∀i], ǫ/2 − γ: The fraction of bases in F-domains can be written as: where stands for the ensemble average. Using the transfer matrix definition in Eq. (28), the average σ for the first monomer reads as: where and A, B and C are given in Eqs. (23)(24)(25). Considering Pauli matrices, V ′ = σ z V , we can write: Using the cyclic property of the trace, Eq. (39) can be extended to any site, n: By diagonalizing V , we can write: where |v ± are the eigenvectors: with the angle θ given by: In the thermodynamic limit Z = λ N + and Eq. (41) reduces to: .
Finally, using Eq. (45) in Eq. (36) leads to: The fraction of bases in C-domains φ C (f ) = (1 − N F (f )) /N, then reads as: On the other hand, the number of compact domains, M C , is given by: Using the transfer matrix formalism, we can write σ n σ n+r as: The number of compact domains [Eq. (48)], then reads as: where we have used Eq. (45).

H DNA hairpin synthesis
H 120 and H 204 were synthesized following the protocol in Ref. [79]. The oligonucleotide forming the 3' end of the hairpin was labelled with a digoxigenin tailing. After a purification step (QIA Nucleotide removal kit), all the oligonucleotides forming the hairpin (

J Drift correction
The drift is a low frequency noise due to macroscopic effects, such as temperature changes or mechanical vibrations. Correcting the drift is important to precisely determine the molecular extension. For long molecules (∼ 1000b), a constant shift λ 0 is imposed in each folding-unfolding cycle to ensure the collapse between consecutive cycles, as described in [32]. However, for shorter molecules (e.g. H 120 ), a time dependent shift, λ 0 (t), needs to be considered to obtain an accurate estimation of the molecular extension. To do so, we first take a reference unfolding curve [e.g. magenta curve in Fig. 15(a)] and find the position of the trap, λ ref , at a certain aligning force where the hairpin is always unfolded (∼ 20 pN, black points in the magenta curve). Next, for all cycles of the tether, the position of the trap at the aligning force, λ align , is computed, and the shift calculated as λ 0 = λ align − λ ref [ Fig. 15(a)]. The continuous function λ 0 (t) is build by performing two spline interpolations (for unfolding and refolding curves) to λ 0 as a function of time data, [ Fig. 15(b)], as in Ref. [46]. This allows us to perform sub-nanometric corrections, as shown in Fig.15(c).

K Elastic parameters and scaling factors
At high forces (around ∼ 15 pN or above) ssDNA is expected to behave as an ideal polymer. The ideal regime can be described by the WLC model [main text, Eq. (2)]. Since the scope of the work is not the characterization of this regime, we use the WLC elastic parameters from Ref. [32] (Sec. E in main text) for the bases forming F-domains (i.e. ideal regime) in the different salt conditions. In this way, we can reduce the number of free parameters of the model to just two: ǫ and γ. However, this restriction leads to small differences between the experimental FECs and the theoretical WLC curves at high forces. To correct these differences, we impose that the ssDNA extension at f = 15 pN coincides with that of the WLC model by multiplying the measured extension by a scaling factor, comprised between 0.95 and 1.05. The elastic parameters (from Ref. [32]) and scaling factors used for each molecule and salt condition are given in Tables 5, 6 and 7.   The binding of the blocking oligonucleotide in the loop region slightly changes the elastic response of the molecule. However, the change in the extension per base induced by the hybridized oligonucleotide decreases with the length of the molecule and it is negligible for long enough molecules [32]. In order to investigate the effect of the blocking oligonucleotide, we measured the FEC (following the procedure described in App. A, main text) with and without the blocking oligonucleotide for the shortest molecule studied with this method, H 204 . As it is shown in Fig. 16, both FECs as well as the elastic fits (persistence, p, and contour length, l) are compatible within error bars. This indicates that the dsDNA hybridized region has a minimal effect on the elasticity of the molecule, for molecules with N 204.

M Obtaining trap stiffness
For forces of the order of a few tens of pN (i.e. small bead displacements in the trap), the optical trap can be approximated as an harmonic potential, with a trap stiffness k b . This allows us to relate the displacement of the bead in the trap, x b , with the applied force as f = k b x b . The bead in the optical trap can be modeled as a Brownian particle in a harmonic potential, and the power spectrum S f (ν) of the measured force follows a Lorenzian distribution: where γ is the drag coefficient of the particle, ν is the frequency and ν c is the so-called corner frequency: ν c = k b /(2πγ). From fitting Eq. (51) to the zero-force power spectrum, γ and k b are determined. N Comparison with the phenomenological formula from [32] In our helix-coil model, the extension of C-domains in the stretching direction can be neglected (App. C, main text), and therefore the ssDNA extension per monomer Eq. (12) (main text) can be written in terms of the fraction of bases in F-domains and its extension: where x F (f ) is the extension of a single monomer given by the WLC model [inverting main text Eq. (2), as described in Ref. [57]] and φ F (f ) can be written as [Eq. (6), main text]: with β = 1/k B T and A and B are given by: Note that φ F corresponds to the fraction of unpaired bases φ defined in Ref. [32], and described by the following phenomenological formula: In order to compare our model with this phenomenological description, we fitted Eq. (56) (f † and δ) to φ F curves [Eq. (53)] with ǫ and γ obtained from the fits of the helix-coil model to the experimental FECs [ Fig. 3(a) and Fig. 4 agreement between the phenomenological description (dashed lines) and the helix-coil model (continuous lines) is observed . Overall, the phenomenological formula reproduces successfully the experimental data, although some differences appear in the predictions for the fraction of free bases at low forces, especially at high salts. The force drops faster from ∼ 2 − 3 pN to zero in the predicted FECs from the helix-coil model than these obtained for the phenomenological formula, as observed in Fig. 19. This is caused by the difference in the leading term in the extension [Eq. (52)]: at the limit f → 0, x F (f ) ∼ f (Hookean limit), which corresponds to a leading term in the extension of ∼ f 2 for the phenomenological formula [Eq. (56)], and of ∼ f 3 for the helix-coil model [Eq. (53)].
The values for f † and δ obtained from fitting the φ F predicted from the helix-coil model are reported in Table 9. These values are similar to those obtained from directly fitting the phenomenological formula Eq. (56) to experimental FECs [32], also shown in Table 9. Note that the φ F fitting procedure allows to estimate f † and δ at low salt concentrations (i.e. 10, 25 and 50mM NaCl) where the direct fitting to Eq. (56) cannot be used (because secondary structure is formed at low forces and FECs do not present any clear force shoulder).
As shown in Fig. 20(a) f † is found to depend logarithmically with salt concentration:      Table 8 recompiles the formulae and values used for obtaining f † and δ for varying length and salt concentration. Table 9: Parameters (ǫ and γ) obtained from the fits of Fig. 3(a) and Fig. 4(a)-(b) (main text). The values of f † and δ have been obtained by fitting Eq. (56) to the obtained φ F [Eq. (53)] using the values of ǫ and γ for each condition. The values of f † and δ given in Ref. [32] are also shown.