Ana Elisa
Bergues-Pupo‡
,
Reinhard
Lipowsky
and
Ana
Vila Verde§
*
Max Planck Institute of Colloids and Interfaces, Department of Theory & Bio-Systems, Am Mühlenberg 1, 14476 Potsdam, Germany. E-mail: ana.vilaverde@mpikg.mpg.de; Fax: +49 (331) 567 9602; Tel: +49 (331) 567 9601
First published on 30th September 2020
Single alpha helices (SAHs) stable in isolated form are often found in motor proteins where they bridge functional domains. Understanding the mechanical response of SAHs is thus critical to understand their function. The quasi-static force–extension relation of a small number of SAHs is known from single-molecule experiments. Unknown, or still controversial, are the molecular scale details behind those observations. We show that the deformation mechanism of SAHs pulled from the termini at pull speeds approaching the quasi-static limit differs from that of typical helices found in proteins, which are stable only when interacting with other protein domains. Using molecular dynamics simulations with atomistic resolution at low pull speeds previously inaccessible to simulation, we show that SAHs start unfolding from the termini at all pull speeds we investigated. Unfolding proceeds residue-by-residue and hydrogen bond breaking is not the main event determining the barrier to unfolding. We use the molecular simulation data to test the cooperative sticky chain model. This model yields excellent fits of the force–extension curves and quantifies the distance, xE = 0.13 nm, to the transition state, the natural frequency of bond vibration, ν0 = 0.82 ns−1, and the height, V0 = 2.9 kcal mol−1, of the free energy barrier associated with the deformation of single residues. Our results demonstrate that the sticky chain model could advantageously be used to analyze experimental force–extension curves of SAHs and other biopolymers.
Natural SAHs typically contain a large fraction of positively and negatively charged amino acids, and (i, i + 3) or (i, i + 4) salt bridges (where i denotes the position of an amino acid in the linear amino acid sequence) are thought to be particularly important for their thermodynamic stability.7,8 The pattern of positively and negatively charged amino acids is not indispensable for SAH formation, however: polyglutamic acid9 and polylysine10,11 also form SAHs at pH levels where the side chain residues are uncharged.
The mechanical response to tensile forces has been investigated using single molecule force spectroscopy (SMFS) or optical tweezers experiments for only a few SAHs.9–13 Current evidence suggests that both the magnitude of the unfolding force and the qualitative features of the force–extension curve depend on the amino acid sequence, for reasons that are not yet understood. Single molecule experiments show that the stretching of the six turn (≈22 amino acids) helical linker connecting the My12 and My13 domains of muscle protein myomesin occurs around 40 pN in a reversible way.13 A comparable unfolding force (F < 30 pN) was also observed for the much longer SAH domain (97 residues) from myosin 10.12 Substantially more intense forces, however, have been reported from atomic force microscopy for unfolding of homopolypeptide SAHs: 150 < F/pN < 200 for polylysine (in the form Cys3–Lys30–Cys)10,11 and up to 1000 pN for polyglutamic acid (Glu80–Cys).9 For the natural SAHs, for polylysine and for coiled coils under tension, experimental evidence indicates that unfolding takes place at constant force, making these structures important force absorbers.10–14 In contrast, for polyglutamic acid the force increases monotonically with extension.9
Steered molecular dynamics simulations can provide great mechanistic insight about the deformation of single helices under tension. To date, most simulations have focused on metastable helices (either homopolypeptides or natural sequences);15–18 whether SAHs respond as metastable helices remains an open question. A constant force plateau was observed in atomistic, explicit solvent simulations of a 52-residue helix belonging to the vimentin AH1 domain simulated isolated in solution (i.e., in a metastable state).18 Hydrogen bond breaking was proposed to be the key barrier to unfolding, in a velocity-dependent manner: at low pull speeds helices appeared to unfold turn-by-turn via the simultaneous breaking of 3-to-4 hydrogen bonds, whereas at large pull speeds (>1 nm ns−1) hydrogen bonds appeared to break individually and helices unfold residue-by-residue. This velocity-dependent mechanism was thought to be behind the two regimes of dependence of the unfolding force on pulling speed that were observed in the simulations. Unfolding started randomly throughout the helix and not just at the pulled ends, i.e., cooperativity was low.17,18 To date we are aware of only a single simulation study of a SAH, with 97 residues.12 That study, with atomistic resolution and implicit solvent representation, indicated this helix unfolds at constant force, consistent with experiments on the same system and with simulations of metastable helices.17,18 In contrast with simulations of metastable helices, unfolding started predominantly from the ends of the SAH suggesting that the cost of creating an interface between helical and unfolded domains is high, i.e., cooperativity in SAHs is high. Hydrogen bond breaking was identified as the key barrier to unfolding, but in a velocity-independent manner. The plateau force observed in simulations at the lowest simulated pull speed (0.01 nm ns−1) was ≈50 pN, comparable to the experimentally detected force (<30 pN) for the same SAH12 despite the large difference in pull speed between experiment and simulation.
There is a clear need to clarify the mechanism of unfolding of SAHs, the extent to which it is cooperative, its dependency on the pull speed and the dominant interactions behind the barrier to unfolding, among other open questions; here we address some of these issues using all-atom molecular simulation studies of a short SAH. Reliable and complete understanding, however, can only be obtained with a strong connection between experiment and simulation. But, how to facilitate this connection given that pull speeds rarely overlap? Both simulation and experimental studies suggest that, for helices as well as for many other polymers under tension, rupture events are well described by a 1-dimensional free energy landscape consisting of two stability wells separated by a barrier.19–21 Coarse insight into this free energy landscape is often extracted, for experimental and for simulation data, via dynamic force spectroscopy models: most often the Bell–Evans,22 Yoreo23 and Dudko24 models. These models yield the distance, xE, from the first stable state to the transition state, and the zero-force off rate, koff, of the rupture event; the Yoreo model in addition includes the equilibrium unfolding force. One way, albeit indirect, of connecting experiment and simulation is by comparing these parameter values. In addition, for a given estimate of the natural vibration frequency, the height, V0, of the free energy barrier can be calculated from koff. This procedure results in estimates of V0 with substantial uncertainty, though. This uncertainty is problematic because our ability to predict the mechanical response of helices – and thus directly connect experiments and simulations – depends critically on this parameter. Bad estimates of V0 lead to predictions of the mechanical response that are not even qualitatively correct.25–27 Another issue with using these models to treat data is that both the Bell–Evans and the Dudko models assume refolding does not take place, an assumption that is not met at low pulling speeds. A final issue is that the models do not consider the case of cooperative helix unfolding.
The sticky chain model proposed by Jäger28 can resolve these issues and thus provide a better connection between SMFS experiments and molecular simulations. This model, based on a 1-dimensional free energy landscape as the Bell and Yoreo models, can in principle be used to fit the full force–extension curve of a polymer under tension and yields, among other parameters, reliable values of V0. It is a kinetic model that allows for unfolding and refolding and is for that reason valid also at low pull speeds, where the Bell–Evans model may fail.20 It comes in two flavors – cooperative or not – enabling insight into this point. It is simple to implement, and thus advantageous relative to Monte Carlo models that could also be used for the same purpose.19 The greater detail of the sticky chain model relative to Bell–Evans-like models thus offers a pathway to connect experimental and simulation force–extension curves. To our knowledge, however, this model has only been used to fit force–extension curves of bacterial pili;29 its suitability for SAHs has not yet been tested.
In this work we combine atomistic, implicit solvent simulations of a SAH with the sticky chain model, to gain insight into the deformation mechanism of SAHs under tension and to test the potential of the sticky chain model to fit force–extension curves of protein helices. We use steered molecular dynamics and umbrella sampling simulations to systematically characterize the mechanical response and deformation mechanism of a short SAH in a range of pull speeds, starting from the high-speed, non-equilibrium unfolding regime all the way to quasi-static unfolding; our lowest pull speed (10−3 nm ns−1) is one order of magnitude lower than reported all-atom simulation studies. We then use the simulation data to test the extent to and conditions under which the sticky chain model28 can be used to fit force–extension curves of SAHs. We evaluate how the parameters of the 1-dimensional free energy landscape associated with residue unfolding obtained from the fit compare with expectations, and with values obtained from fits of the Bell–Evans and the Yoreo models.
Pulling simulations were conducted at constant speed, v ∈ {10−3, 10−2, 10−1, 1} nm ns−1. For the three highest pull speeds, thirty independent runs were started from the same initial configuration, with velocity reassignment at t = 0. For the lowest pull speed, only four independent runs were performed because of their high computational cost. Force–extension curves averaged over 4 runs already overlap with the 30-run average for each of the two intermediate pull speeds (comparison not shown). The 4 runs done for the lowest pull speed should thus be sufficient. Pulling was done via a virtual spring of constant ks connected to the Cα carbon of glycine at the C-terminus of the helix; the Cα of the glycine at the N-terminus was restrained via a virtual spring of constant K, see Fig. 1. The virtual springs had values K = 10ks = 1000 kJ mol−1 nm−2. Configurations were saved every vt = 0.01 nm unless otherwise noted, resulting in 1000 saved configurations per simulation run. For the analysis of the hydrogen bond breaking events, simulations at the lowest pull speed were repeated to save configurations every 10 ps (i.e., the same as for the highest pull speed) but using a limited strain range, 0.5 < ε = vt/L0 < 0.8, to avoid excessively large trajectory files.
Umbrella sampling simulations were performed to build the potential of mean force as a function of the extension. Starting conformations were obtained from a pulling trajectory at v = 10−2 nm ns−1, by selecting configurations at equally spaced displacements (vt = 0.2 nm), yielding 30 separate umbrella simulations. The umbrella restraint was applied to the Cα of the terminal amino acids and the restraint constant was 100 kJ mol−1 nm−2. Each umbrella was simulated for 50 ns.
![]() | ||
Fig. 2 Double well potential of the chain model under a constant force F.28 The folded and unfolded states have equilibrium distances b0 and b0 + ΔL, respectively and elastic constants k1 and k2. The position of the unfolding barrier is xE and is force-independent. |
The transitions between the two states are thermally activated and are described by Bell-like kinetics, namely the transition rate depends exponentially on the external force. The equations describing the model are for the force–strain curve:
![]() | (1) |
![]() | (2) |
ε = vt/L0, | (3) |
A particularity of this double-well potential model is that at zero force, the energy of the unfolding barrier V0 coincides with the energy of the unfolded state, i.e., the barrier to refolding at zero force is zero. This choice has negligible impact on the dependence of the unfolding force on the pulling speed20 but has the advantage of allowing a direct determination of the attempt frequency of opening of each link:
koff = ν0 × exp−βV0. | (4) |
![]() | (5) |
For the conditions of our simulations, only 3 parameters of the two state model (k1, k2 and ν0) are obtained from fitting the force–extension curves obtained in the atomistic simulations, because the remaining parameters can be inferred from the atomistic model as explained below.
![]() | ||
Fig. 3 Atomistic simulations: (A) average force–strain curves at the indicated pull speeds, obtained from the molecular dynamics simulations. For the lowest pull speed, the average is over 4 force traces; for the other pull speeds, it is over 30 force traces. The strain is defined as ε = vt/L0. Average values of the plateau forces are shown in Table 1. (B) Potential of mean force from umbrella sampling (black), work (W–v = 10−3; cyan) realized at the lowest pull speed, and potential energy for the 4 pull velocities (same colors as in (A)), as a function of the chain strain (ε = x/L0), obtained from the molecular dynamics simulations. The different definition of strain used here is necessary to directly compare the PMF with the quantities extracted from the pull simulations. Every (E,ε) point of the potential energy curves is obtained by binning E(t) values of all simulation runs performed for each v in the corresponding strain bin and averaging over the number of E-points per bin (this number is not constant); the shadow show the standard error of the mean. An analogous procedure is used to obtain the W–v = 10−3 curve. |
Analytical and numerical models have demonstrated that force plateaus naturally arise from a two-state free energy landscape associated with extension of links in a chain; the extent of cooperativity determines whether the force remains truly constant or increases slowly with extension in the plateau or quasi-plateau region.21,28 Force plateaus have also been seen in simulations and experiments of different molecules under tension and in a wide range of pulling velocities, suggesting two-state free energy landscapes associated with stretching are common to many (bio)polymers. For example, force plateaus have been observed in simulations of helices in tension, even at higher pulling velocities16,18,35 or under adiabatic conditions;34,36,37 in single molecule measurements of polylysine SAHs,10 of helical segments connected by joints,37,38 of coiled coils,14,39 of dextran19 and of double stranded DNA;40,41 in steered molecular dynamics simulations of duplex DNA,42 and simulations of coiled coils in tension39 and shear.43,44
The average plateau force (Fp; see Table 1), increases substantially with increasing pull speed, varying between at the lowest pull speed up to ≈200 pN at the highest pull speed. The plateaus at the two lowest pull speeds, v ∈ {10−3, 10−2} nm ns−1, almost overlap, suggesting that unfolding at the lowest pull speed may take place under quasi-static conditions. To test this possibility, we calculated the work realized on the helix at the lowest pull speed and compare it with the potential of mean force (PMF) as a function of strain calculated using the umbrella sampling simulations. Both curves (light blue and black) are shown in Fig. 3(B) and clearly do not overlap, indicating that unfolding at the lowest velocity in fact does not proceed under quasi-static conditions.
v (nm ns−1) | ||
---|---|---|
10−1 | 10−2 | 10−3 |
126 ± 0.5 | 97 ± 0.5 | 87 ± 1 |
The difference between the PMF and the work realized at v = 10−3 nm ns−1 (Fig. 3(B)) together with the fact that the plateau forces at the two lowest pulling velocities are very similar, suggest that there is a timescale separation between the processes playing a significant role at equilibrium and those important for v ≥ 10−3 nm ns−1. We attempted to gain insight into these processes by plotting the secondary structure as a function of time for two fixed values of strain of the umbrella simulations. The results, shown in ESI,† Fig. S8, indicate that the secondary structure substantially changes as a function of time in some parts of the peptide and that some of the changes occur only once within the simulated 50 ns. We speculate that these infrequent but substantial structural changes increase the chain entropy in the PMF calculations and thus soften the helix. These rare events are not well-sampled even at the slowest pull speed used in the pull simulations, which would explain why the work to extend the helix realized at this pull speed is larger than the PMF.
A plateau force <30 pN was reported from experiment for a SAH rich in negative and positive charged amino acids;12 quasi-static experiments done on a SAH of comparable composition reported a similar unfolding force.13 Those forces are significantly lower than . This difference between experiment and simulation may be an artifact of the simulations, such as the limited timescales accessed as discussed above, or shortcomings of the force field. Alternatively, it may reflect true differences between the helix simulated here and those studied experimentally. Clarifying this point requires closer connection between experiment and simulation setups and is beyond the scope of this work. Investigating the non-equilibrium mechanical response of SAHs to applied tension, which will likely require high-speed SMFS setups,45 will be an important step towards this aim.
The large contribution of entropy to the mechanical response of helices under tension may have substantial implications: if unfolding is cooperative – i.e., the price to create multiple interfaces between folded and unfolded sections throughout the helix is high, so the helix unfolds residue-by-residue starting from the termini – then the entropic contribution arises at the level of individual amino acids. This contribution is by definition included in 1-dimensional models of the free energy landscape of amino acid stretching, and so these models should be sufficient to model the mechanical response of helices under tension. On the contrary, if unfolding under tension is non-cooperative – i.e., multiple unfolding/refolding events occur randomly throughout the length of the helix – then chain entropy plays a role, and models of the mechanical response of helices under tension must include it. We note that the random opening version of the sticky chain model indeed allows for unfolding/refolding throughout the helix but does not correctly account for chain entropy.28
![]() | ||
Fig. 4 Secondary structure changes as a function of strain at the highest (top) and lowest (bottom) velocities for four different realizations. The y-axis shows the secondary structure of each amino acid in the chain, determined using STRIDE:48 turn (turquoise); α-helix (pink); 310-helix (dark blue) and coil (white). The light blue circle indicates the restrained end of the helix (amino terminus); the acid terminus, connected to the traveling virtual spring, is at y = 0. |
Fig. 4 gives the impression that unfolding occurs in clusters of 3 or 4 consecutive amino acids. This impression is misleading and results from the manner in which STRIDE assigns secondary structures. For example, if amino acids i + 1, i + 2, i + 3 and i + 4 initially meet all the criteria used to define a α-helix, but the hydrogen bond between i and i + 4 subsequently breaks, all four amino acids will lose their α-helix status. To ascertain whether unfolding occurs in clusters or residue-by-residue, we examined the backbone ϕ and ψ angles of each amino acid as a function of time. Results for selected residues are shown in ESI,† Section S3. The results confirm that at both speeds, unfolding occurs predominantly residue-by-residue.
![]() | (6) |
v (nm ns−1) | I | II | III |
---|---|---|---|
1 | 0.10 ± 0.02 | 0.26 ± 0.02 | 0.64 ± 0.03 |
10−3 | 0.13 ± 0.02 | 0.28 ± 0.03 | 0.59 ± 0.04 |
α-Helical structure is determined by backbone angles ϕ and ψ. To assess whether the frequent hydrogen bond breaking events seen in the HB-break time series in Fig. 5 are accompanied by loss of helical structure, we plotted these angles as Ramachandran plots for all residues forming the loop stabilized by the backbone hydrogen bonds (i, i + 4) indicated in that figure. These results are shown in ESI,† Fig. S4–S7. The Ramachandran plots confirm that transient hydrogen bond breaking is not accompanied loss of α-helical structure. Permanent destruction of backbone (i, i + 4) hydrogen bonds only occurs after at least one of the relevant ϕ and ψ transitions out of α-helical configuration.
These results demonstrate that hydrogen bonds, while clearly contributing to the thermodynamic stability of helices, are not the main barrier to helix stretching. We (and others)52 propose that this conclusion is general because experiments do not suggest a substantial dependence between the strength of intrahelical hydrogen bonds and amino acid identity,53 but polylysine, polyglutamate and SAHs have very different mechanical response under tension as discussed in the introduction. Which interaction(s) contribute the most to the free energy barrier to unfolding of helices under tension remains for now an open question and is being investigated by our group.
![]() | (7) |
The cross-covariance functions are shown in the bottom plot of each panel in Fig. 5. A large peak at t = 0 is never observed, which demonstrates that consecutive hydrogen bonds along the helical backbone do not preferentially break simultaneously at any pull speed. Instead, all cross-covariance functions show distributions of time lags of several nanoseconds. These distributions are, as expected, narrower for the low pull speed because strain values leading to substantial loss of helical structure and thus permanent breakage of hydrogen bonds (see also discussion in ESI,† Section S2) are reached faster.
Parameter | Value |
---|---|
N | 27 |
b 0 | 0.15 nm |
ΔL | 0.23 nm |
V 0 | 20 pN nm (4.8kBT) |
The fitting parameters obtained in each case are shown in Table 4, and the corresponding curves are shown in Fig. 6. Both fitting attempts yield comparable fitting parameters and reasonable force–strain curves for the 4 pull speeds, confirming that the model and fitting procedure are robust. These results indicate that reliable insight into the 1-dimensional free energy landscape associated with helix unfolding under tension may be obtained already with one quasi-static and one non-equilibrium force–strain curve. Fits using curves at more than one pull speed are nevertheless expected to yield better estimates of the fit parameters, so the following discussion focuses on the results of the two-curve fit.
Parameter | v (nm ns−1) | |
---|---|---|
10−1 | {10−1, 10−13} | |
k 1 (pN nm−1) | 1321 | 2118 |
k 2 (pN nm−1) | 5188 | 4792 |
ν 0 (ns−1) | 0.23 | 0.82 |
k off (s−1) | 0.18× 107 | 0.63× 107 |
x E (nm) | 0.17 | 0.13 |
![]() | ||
Fig. 6 Top: Force as a function of strain obtained with the sticky chain model (lines) using the two set of parameters in Table 3; to facilitate comparisons, the force–strain curves from all-atom simulations, already presented in Fig. 3, are repeated here (points). (A) Fitting one curve, v = 10−1 nm ns−1. (B) Fitting 2 curves, at v ∈ {10−1, 10−3} nm ns−1. Bottom: Number of folded links as a function of strain obtained from the sticky chain model. (C) Fitting one curve. (D) Fitting 2 curves. |
The bond constant, k2 = 4792 pN nm−1, characterizing the free energy basin of the unfolded state is substantially larger than k1 = 2118 pN nm−1, which characterizes the free energy basin of the folded state. This difference is consistent with the fact that stretching an extended amino acid chain involves deforming angle potentials (bonds are constrained in our model), which are substantially stiffer than any stabilizing contribution to the folded state (dihedrals, hydrophobic effect, hydrogen bonds and salt bridges).
The natural vibration frequency, ν0 = 0.82 ns−1, obtained from the fit is substantially lower than the typical value of 104 ns−1 used as estimate of this quantity.55 This difference demonstrates that the commonly used estimate of ν0 is in fact far from suitable in this case, where each link has much higher mass and much lower bond stiffness (k1) than covalent bonds between 2 atoms, for which the estimate was originally made. A reliable value of ν0 is indispensable to calculate V0 from values of koff (eqn (4)), which can be easily obtained from Bell–Evans fits of rupture force vs. velocity data. This was the procedure followed by Buehler et al.18 to analyze force–strain curves from simulation; from a kBuehleroff = 0.807 × 105 s−1 they obtain VBuehler0 = 11.1 kcal mol−1 for a metastable helix and the same pull speed range investigated here. Using our estimate of ν0 would result in VBuehler0 = 5.5 kcal mol−1. This value is still substantially higher than our own estimate, V0 = 4.8kBT = 2.9 kcal mol−1, because kBuehleroff is 2 orders of magnitude lower than what we obtain (see Table 4). Whether these differences reflect true differences between the SAH investigated here and the AH1 helix from the vimentin intermediate filament investigated by Buehler et al., or are instead artifacts arising from the different force fields used, remains unclear.
![]() | (8) |
The Yoreo model23 relates the same two observables,
![]() | (9) |
In the top panel of Fig. 7 we show force–strain curves obtained with the sticky chain model and the parameters in Tables 3 and 4 (for fits done at v ∈ {10−1, 10−3} nm ns−1) for a wide range of pull velocities. Values of the plateau force at each pull speed, indicated by the black markers, are re-plotted as a function of the loading rate in the bottom panel of the same figure. We fit the Bell–Evans model to the 3 points with higher loading rate, for which re-binding is expected to be minimal, and fit the Yoreo model to the entire curve. The Bell–Evans model is linearized for the fit (F = p1 × x + p2), with Δx and koff obtained as Δx = kBT/p1 and koff = (p1exp(−p2/p1))−1.
![]() | ||
Fig. 7 (A) Force–strain curves generated from the sticky chain model using the parameters in Tables 3 and 4 (for fits using the curves at v ∈ {10−1, 10−3} nm ns−1); the black crosses indicate the values taken as the plateau force shown in panel (B). The pull speed is indicated as a multiple of the quasi-static velocity, veq = 10−3 nm ns−1. (B) Plateau force vs. loading rate, ksv, from the sticky chain model (points) and fits using the Bell–Evans and Friddle–Yoreo model (lines). |
The parameter values obtained from the fits are shown in Table 5. Parameters extracted with both models are comparable, but not identical, to those used as input (Table 4) to the sticky chain model, indicating that both the Yoreo and the Bell–Evans models are reasonable choices for initial analyses, but yield only coarse estimates of Δx and koff. These results illustrate that, for helical systems where a (quasi-)force plateau can be observed as a function of strain, fits using the sticky chain model are preferable.
Bell–Evans | Yoreo | |
---|---|---|
Δx (nm) | 0.14 ± 0.06 | 0.11 ± 0.01 |
k off (s−1) | (0.09 ± 0.2) × 108 | (0.19 ± 0.04) × 108 |
Footnotes |
† Electronic supplementary information (ESI) available: Individual force–strain traces at the lowest pull speed; hydrogen bonds and salt bridges vs. strain; Ramachandran plots as a function of time; secondary structure vs. time from umbrella sampling simulations. See DOI: 10.1039/d0sm01166e |
‡ Present address: Max Delbrück Center for Molecular Medicine, Hannoversche Str. 28, 10115 Berlin, Germany. |
§ Present address: University of Duisburg-Essen, Faculty of Physics, Lotharstr. 1, 47057 Duisburg, Germany. |
This journal is © The Royal Society of Chemistry 2020 |