Fold Catastrophes and the Dependence of Free-Energy Barriers to

J. Phys. Chem. B 2010, 114, 10821–10825

10821

Fold Catastrophes and the Dependence of Free-Energy Barriers to Conformational Transitions on Applied Force Daniel J. Lacks,* Joshua Willis,† and Michael-Paul Robinson‡ Department of Chemical Engineering, Case Western ReserVe UniVersity, CleVeland, Ohio 44016 ReceiVed: July 14, 2010

Applied mechanical force (f ) can activate conformational change in molecules by reducing the height of a free-energy barrier (∆Gb). In this paper, molecular dynamics simulations are carried out with umbrella sampling and self-consistent histogram methods to determine free-energy profiles for a coarse-grained model of a protein under an applied force. Applied force is shown to cause fold catastrophes, where free-energy minima are destabilized until they disappear. It is well-known that a fold catastrophe at force f ) B implies the scaling ∆Gb ≈ |B - f|3/2 in the limit of ∆Gb f 0, but it is not clear whether this scaling is accurate for physically relevant barrier heights. The simulation results show that the fold catastrophe scaling is in fact accurate in the physically relevant regime and that the two-parameter function ∆Gb ) A(B - f )3/2 is superior to the twoparameter linear function for parametrizing changes in free-energy barriers with applied force. Mechanical force can cause molecules to undergo conformational transitions. For example, proteins under pulling forces unfold from their native states, and these transitions play roles in physiological processes such as muscle function.1 Experimentally, these force-induced conformational transitions are investigated in single molecules using atomic force microscopy2 and optical tweezer techniques.3,4 Mechanical force acts to reduce free-energy barriers,5,6 and this idea forms the basis of quantitative models for conformational transitions.7-13 Basically, the transition rate, r, is approximated as

r ) ν exp(-∆Gb /kT) where ∆Gb is the height of the free-energy barrier, ν is a frequency factor, k is the Boltzmann constant, and T is temperature. The transition is likely within the time interval ∆t when r∆t > ∼1, which occurs when ∆Gb/kT < ∼ln(ν∆t). For experiments where one end of the molecule is pulled at constant velocity V (with the other end fixed), it is more convenient to express this result as ∆Gb/kT < ∼ln(ν∆L/V), where ∆L is an elongation increment and ∆t ) ∆L/V. When the pulling begins, ∆Gb/kT . ln(ν∆L/V), and the molecule stretches without undergoing a conformational transition. However, ∆Gb decreases as the molecule is stretched, and the unfolding transition becomes likely when ∆Gb decreases to the point that ∆Gb/kT ≈ ln(ν∆L/V). For stretching experiments on proteins, ∆L ≈ 1 nm (an upperbound to ∆L is the distance along the stretching coordinate to the barrier for unfolding, which for proteins is usually 1-2 nm21), an estimate23 for ν is ν ≈ 108 s-1, and experiments14 are carried out with V ≈10-104 nm/s; on the basis of these parameter values, unfolding transitions are expected when ∆Gb/kT ≈ 10 at higher V and ∆Gb/kT ≈ 20 at smaller V. To apply these models, it is necessary to have a parametrized function for the change in ∆Gb with applied force f. A linear * To whom correspondence should be addressed. E-mail: daniel.lacks@ case.edu. † E-mail: [email protected]. ‡ E-mail: [email protected].

approximation is often used,5,7,8 ∆Gb ) A - Bf, which corresponds to the first-order term in a Taylor series and thus is appropriate for small changes in f. Other investigators have suggested the function ∆Gb ) A(B - f )3/2 because the force is expected to cause a fold catastrophe (where ∆Gb f 0), and this function is theoretically expected in the ∆Gb f 0 limit;15-17 however, the accuracy of this function is not known for physically relevant barrier heights, which are in the range of ∼10 < ∆Gb/kT < ∼20, as discussed above. Previous work addressed this question for potential energy barriers,18,19 but freeenergy barriers differ fundamentally in that they include entropic contributions (which are particularly important for unfolding transitions) and that they are of a much lower dimensionality (potential energy landscapes are in 3N dimensions, where N is the number of particles, while free-energy landscapes may be in one dimension that is described by an order parameter such as the end-to-end distance of a molecule). Molecular simulations are used here to assess the accuracy of these functions for describing free-energy barrier heights in the physically relevant regime. The simulations are carried out for a model of the 27th immunoglobulin domain (Ig27) of the muscle protein titin, which is the benchmark molecule for forceinduced unfolding studies; for example, experiments have been carried out on engineered proteins composed solely of repeated units of Ig27,20 and Ig27 has been addressed in many molecular simulations.21-25 A coarse-grained model of Ig27, shown in Figure 1, is used here so that long simulations can be run to determine free-energy barrier heights to high precision; each amino acid residue is represented by a single site, and interactions between sites are parametrized by fitting to the experimental structure.26,27 Comparison with all-atom models show that this type of coarse-grained model captures the essential physics of the protein.28 A pulling force is applied by including an external electric field, E, oriented along the z-axis and small atomic charges of opposite polarity on the first and last residues (q ) (0.1 e); the end residues are thus pulled in opposite directions along the z-axis with force f ) qE. The atomic charges have a negligible effect on the structure in the absence of the electric field, and the magnitude of the force is varied by changing E.

10.1021/jp106530h  2010 American Chemical Society Published on Web 08/03/2010

10822

J. Phys. Chem. B, Vol. 114, No. 33, 2010

Lacks et al.

Figure 1. Model of the I27 domain. Each site represents an amino acid residue, and the end residues, which are pulled in opposite directions in the simulations, are highlighted in red.

The free energy as a function of the end-to-end-distance of the molecule, G(L), is determined using the method of Li and Makarov.29 Basically, G(L) is related to the probability that the system has a particular end-to-end distance, Prob(L)

G(L) ) -kT ln[Prob(L)] + C The values of Prob(L) are obtained as histograms constructed from a molecular dynamics trajectory. However, this basic approach would be ineffective as the simulation would not sufficiently sample configurations near barriers. Umbrella sampling is used to overcome this problem by adding a harmonic potential, 1/2ku(L - Lu)2, to bias the system to sample configurations near L ) Lu and then removing its effects to obtain the true free-energy profile

1 G(L) ) -kT ln[Prob(L)] - ku(L - Lu)2 + C 2 To obtain G(L) over a range of L, simulations are carried out with various Lu, and the self-consistent histogram method is used to optimally combine these results to obtain the best estimate of G(L). The energy parameters of this coarse-grained model have arbitrary values that have not been scaled to match experimental unfolding temperatures; therefore, free-energy and force results are reported as dimensionless quantities reduced by kT and the initial end-to-end distance, L0 (e.g., dimensionless force f* ) fL0/kT). The temperature is about half of the temperature of the folding transition for this model (due to the arbitrary scale of the energy parameters, this is T ) 50 K). For each value of the applied force, umbrella sampling is carried out using between 8 and 32 simulations with ku ) 2kT and different values of Lu. Each simulations is run for 107 steps with a time step of 0.0005 ps. The simulations are carried out with the Gromacs software package.30

Figure 2. Free-energy landscapes, G(L), for various values of the pulling force. Note that the y-axis offset for each curve is arbitrary and is chosen to allow clear display of the curves. The bottom panels show in more detail the changes near particular barriers.

Free-energy landscapes are shown in Figure 2 for various values of f*. Increasing f* destabilizes the more compact structures and changes the free-energy landscape as follows: (1) For f* < 3.50, the landscape has a local minimum corresponding to the native state (N), and a local minimum corresponding to the unfolded state (U) at large L. State N becomes increasingly destabilized with respect to state U as f* increases. While GN < GU at f* ) 0, the value of (GU - GN) decreases with increasing f*, such that GU < GN by the time f* ) 3.50; we did not attempt to find the precise value where GU ) GN because it is unimportant as the very high barrier between states N and U (∆Gb,NU > 50kT for f* < 3.50) prevents possible transitions to the state U. (2) At f* ) 3.50, a third local minimum appears between states N and U, corresponding to an unfolding intermediate state (I). States N and I are separated by a barrier with height ∆Gb,NI in the forward direction (increasing L) and height ∆Gb,IN in the reverse direction. Similarly, states I and U are separated by a barrier with height ∆Gb,IU in the forward direction and height ∆Gb,UI in the reverse direction. The values of ∆Gb,NI, ∆Gb,IN, and ∆Gb,IU are shown in Figure 3 (the values of ∆Gb,UI were not determined as they are not necessary for our analysis). (3) State N becomes increasingly destabilized with respect to state I as f* increases. While GI > GN at f* ) 3.50, the value of (GI - GN) decreases with increasing f*, and GI < GN for f* > 4.71; this result is evident in Figure 3, noting that (GI - GN) ) (∆Gb,NI - ∆Gb,IN). Also, ∆Gb,NI decreases with increasing f*, as shown in Figure 3; at f*

Free-Energy Barriers

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10823

Figure 3. Heights of free-energy barriers as a function of pulling force. Squares represent ∆Gb,NI, triangles represent ∆Gb,IN, and circles represent ∆Gb,IU.

) 5.38, the value of ∆Gb,NI decreases to 0, and state N disappears, corresponding to a fold catastrophe. (4) Similarly, state I becomes increasingly destabilized with respect to state U as f* increases. When state I first appears, ∆Gb,IU is very large, but it decreases with increasing f*, as shown in Figure 3, and is in the physically relevant range (10 < ∼∆Gb/kT < ∼20) for 5.9 < f* < 6.5. A fold catastrophe occurs at f* ) 7.71, where the value of ∆Gb,IU decreases to 0 and state I disappears. The response of this protein to the pulling force under experimentally relevant pulling rates can be inferred from these free-energy landscapes. For f* < 4.71, the molecule remains in state N since this has the lowest free energy of an accessible state (even though state U has a lower free energy, the large barrier, >35kT, prevents the transition); in this regime, the most probable end-to-end distance, corresponding to the position of the minimum of the G(L) curve, increases with increasing f*. At f* ≈ 4.71, the molecule undergoes a transition to state I as the free energy of state I becomes lower than that of state N; since the barrier ∆Gb,NI is small at this point ( 5 (modeling smaller barriers is not as important as the system quickly surmounts these barriers). It is clear that the fold-catastrophe-motivated form gives a better fit to the data than the linear form with the same number of fitting parameters; the mean-squared error in the fits is 4 times larger for the linear function than that for the fold-catastrophe-motivated function. Why is the function motivated by fold catastrophe scaling, ∆Gb ) A(B - f )3/2, accurate outside of the ∆Gb f 0 regime for which it is rigorously valid, and is this accuracy expected to be general? Previous studies of potential energy barriers, ∆Eb (rather than free-energy barriers, ∆Gb), give insight into these questions.18,19 First, the energy barrier scaling is expected to be

10824

J. Phys. Chem. B, Vol. 114, No. 33, 2010

Lacks et al. for yield stresses of metallic glasses follow a (T/Tg)2/3 temperature dependence (where Tg is the glass transition temperature), which they attribute to a ∆Gb ) A(B f )3/2 dependence of barrier heights (here, f represents the stress).31 In their work, it was simply assumed that this scaling is reasonable at the large barrier heights relevant to yielding; the present results show direct support of this assumption. Acknowledgment. This material is based upon work supported by the National Science Foundation under Grant No. DMR-0705191. References and Notes

Figure 5. Comparison of two-parameter fits to free-energy barrier heights for the barrier ∆Gb,IU. The dashed line represents ∆Gb ) A + Bf, and the solid line represents ∆Gb ) A(B - f )3/2. Only data with ∆Gb/kT > 5 are used for the fitting.

more accurate away from the ∆Eb f 0 regime than other fold catastrophe scaling relations (i.e., those for ∆x and the curvature at the minimum) because the energy barrier is obtained from integrations of the other quantities, and integrations act to prolong earlier behavior.18 Second, a fitting procedure in which the parameters A and B are fit to the entire set of results increases the accuracy of the scaling relation far from ∆Eb f 0 limit (even though the fitted parameters may differ from those obtained for the ∆Eb f 0 regime where the scaling is rigorously valid).19 Third, ∆Eb ) A(B - f )3/2 was found to be accurate outside of the ∆Eb f 0 regime in several diverse systems, including models of proteins and metallic glasses,18 as well as simple model potentials.19 Thus, the ∆Gb ) A(B - f )3/2 functional form, which is exact in the ∆Gb f 0 limit, is also accurate for barrier heights that are physically relevant. This result has the implications that can be observed experimentally (1) The applied force required to activate a conformational change, Fmax, is known to depend on the pulling speed V. The dependence of Fmax on V is governed by the functionality of ∆Gb(f ); Dudko et al. show that the functional form ∆Gb ) A + Bf implies the relationship Fmax ) a + b* ln(V), while the functional form ∆Gb ) A(B f )3/2 implies the relationship Fmax) a + b* ln(V2/3).16 The present results show that free-energy barrier heights are best represented by the form ∆Gb ) A(B - f )3/2, thus suggesting that the relationship Fmax ) a + b* ln(V2/3) should be most accurate. (2) To extract kinetic parameters from single-molecule pulling experiments, a model must be assumed for ∆Gb(f ). Dudko et al. show that values obtained for kinetic parameters strongly depend on the assumed model (varying by more than an order of magnitude between models).10 The present results suggest that the model with ∆Gb ) A(B - f )3/2 would give the most reliable kinetic parameters. (3) Glassy materials undergo plastic deformation when subjected to a stress exceeding the yield stress of the material. This plastic deformation occurs through molecular rearrangements that are limited by free-energy barriers; these barriers are reduced by applied stress, and the yield stress corresponds to the stress that is required to lower the barrier to the point that the rearrangements can occur (the yield stress is analogous to Fmax). Johnson and Samwer show that experimentally determined values

(1) Erickson, H. P. Reversible Unfolding of Fibronectin Type-III and Immunoglobulin Domains Provides the Structural Basis for Stretch and Elasticity of Titin and Fibronectin. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 10114–10118. (2) Rief, M.; Gautel, M.; Oesterhelt, F.; Fernandez, J. M.; Gaub, H. E. Reversible Unfolding of Individual Titin Immunoglobulin Domains by AFM. Science 1997, 276, 1109–1112. (3) Kellermayer, M. S. Z.; Smith, S. B.; Granzier, H. L.; Bustamante, C. Folding-Unfolding Transitions in Single Titin Molecules Characterized with Laser Tweezers. Science 1997, 276, 1112–1116. (4) Tskhovrebova, L.; Trinick, J.; Sleep, J. A.; Simmons, R. M. Elasticity and Unfolding of Single Molecules of the Giant Muscle Protein Titin. Nature 1997, 387–308. (5) Eyring, H. Viscosity, Plasticity, and Diffusion as Examples of Absolute Reaction Rates. J. Chem. Phys. 1936, 4, 283–291. (6) Kramers, H. A. Brownian Motion in the Field of Force and the Diffusion Model of Chemical Reactions. Physica 1940, 7, 284–304. (7) Bell, G. I. Models for the Specific Adhesion of Cells to Cells. Science 1978, 200, 618–627. (8) Evans, E.; Ritchie, K. Dynamic Strength of Molecular Adhesion Bonds. Biophys. J. 1997, 72, 1541–1555. (9) Merkel, R.; Nassoy, P.; Leung, A.; Ritchie, K.; Evans, E. Energy Landscapes of Receptor-Ligand Bonds Explored with Dynamic Force Spectroscopy. Nature 1999, 397, 50–53. (10) Dudko, O. K.; Hummer, G.; Szabo, A. Intrinsic Rates and Activation Free Energies from Single-Molecule Pulling Experiments. Phys. ReV. Lett. 2006, 96, 108101. (11) Hummer, G.; Szabo, A. Kinetics from Nonequilibrium SingleMolecule Pulling Experiments. Biophys. J. 2003, 85, 5–15. (12) Dudko, O. K. Single-Molecule Mechanics: New Insights from the Escape-over-a-Barrier Problem. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 8795–8796. (13) Freund, L. B. Characterizing the Resistance Generated by a Molecular Bond as it is Forcibly Stretched. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 8818–8823. (14) Best, R. B.; Fowler, S. B.; Toca Herrera, J. L.; Steward, A.; Paci, E.; Clarke, J. Mechanical Unfolding of a Titin Ig Domain: Structure of Transition State Revealed by Combining Atomic Force Microscopy, Protein Engineering and Molecular Dynamics Simulations. J. Mol. Biol. 2003, 330, 867–877. (15) Garg, A. Escape-Field Distribution for Escape from a Metastable Potential Well Subject to a Steadily Increasing Bias Field. Phys. ReV. B 1995, 51, 15592–15595. (16) Dudko, O. K.; Filippov, A. E.; Klafter, J.; Urbakh, M. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 11378–11381. (17) Dias, C. L.; Dube, M.; Oliveira, F. A.; Grant, M. Scaling in Force Spectroscopy of Macromolecules. Phys. ReV. E 2005, 72, 011918. (18) Maloney, C. E.; Lacks, D. J. Energy Barrier Scalings in Driven Systems. Phys. ReV. E 2006, 73, 061106. (19) Tshiprut, Z.; Urbakh, M. Exploring Hysteresis and Energy Dissipation in Single-Molecule Force Spectroscopy. J. Chem. Phys. 2009, 130, 084703. (20) Marszalek, P. E.; Lu, H.; Li, H.; Carrion-Vazquez, M.; Oberhauser, A. F.; Schulten, K.; Fernandez, J. M. Mechanical Unfolding Intermediates in Titinmodules. Nature 1999, 402, 100–103. (21) Lu, H.; Schulten, K. The Key Event in Force-Induced Unfolding of Titin’s Immunoglobulin Domains. Biophys. J. 2000, 79, 51–65. (22) Paci, E.; Karplus, M. Unfolding Proteins by External Forces and Temperature: The Importance of Topology and Energetics. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 6521–6526. (23) Li, P.-C.; Makarov, D. E. Theoretical Studies of the Mechanical Unfolding of the Muscle Protein Titin: Bridging the Time-Scale Gap between Simulation and Experiment. J. Chem. Phys. 2003, 119, 9260–9268.

Free-Energy Barriers (24) Cieplak, M.; Hoang, T. X.; Robbins, M. O. Thermal Effects in Stretching of Go-Like Models of Titin and Secondary Structures. Proteins 2004, 56, 285–297. (25) Duff, N.; Duong, N.-H.; Lacks, D. J. Stretching the Immunoglobulin 27 Domain of the Titin Protein: The Dynamic Energy Landscape. Biophys. J. 2006, 91, 3446–3455. (26) Clementi, C.; Nymeyer, H.; Onuchic, J. N. Topological and Energetic Factors: What Determines the Structural Details of the Transition State Ensemble and En-Route Intermediates for Protein Folding? An Investigation for Small Globular Proteins. J. Mol. Biol. 2000, 298, 937–953. (27) SMOG@ctbp: Structure-based MOdels in Gromacs. http://smog. ucsd.edu/. (28) Whitford, P. C.; Noel, J. K.; Gosavi, S.; Schug, A.; Sanbonmatsu, K. Y.; Onuchic, J. N. An All-Atom Structure-Based Potential for Proteins:

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10825 Bridging Minimal Models with All-Atom Empirical Forcefields. Proteins 2009, 75, 430–441. (29) Li, P.-C.; Makarov, D. E. Ubiquitin-like Protein Domains Show High Resistance to Mechanical Unfolding Similar to That of the I27 Domain in Titin: Evidence from Simulations. J. Phys. Chem. B 2004, 108, 745– 749. (30) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. (31) Johnson, W. L.; Samwer, K. A Universal Criterion for Plastic Yielding of Metallic Glasses with a (T/Tg) 2/3 Temperature Dependence. Phys. ReV. Lett. 2005, 95, 195501.

JP106530H