Article Information

• PDF (1915 kb)
• Full Text with Large Figures

PubMed

• Articles by Eran Eyal
• Articles by Ivet Bahar

Related Articles

• Selection of Optimal Variants of Gō-Like Models of Proteins ...

  Selection of Optimal Variants of Gō-Like Models of Proteins through Studies of Stretching Biophysical Journal, Volume 95, Issue 7, 1 October 2008, Pages 3174-3191 Joanna I. Sułkowska and Marek Cieplak Abstract The Gō-like models of proteins are constructed based on the knowledge of the native conformation. However, there are many possible choices of a Hamiltonian for which the ground state coincides with the native state. Here, we propose to use experimental data on protein stretching to determine what choices are most adequate physically. This criterion is motivated by the fact that stretching processes usually start with the native structure, in the vicinity of which the Gō-like models should work the best. Our selection procedure is applied to 62 different versions of the Gō model and is based on 28 proteins. We consider different potentials, contact maps, local stiffness energies, and energy scales—uniform and nonuniform. In the latter case, the strength of the nonuniformity was governed either by specificity or by properties related to positioning of the side groups. Among them is the simplest variant: uniform couplings with no , +2 contacts. This choice also leads to good folding properties in most cases. We elucidate relationship between the local stiffness described by a potential which involves local chirality and the one which involves dihedral and bond angles. The latter stiffness improves folding but there is little difference between them when it comes to stretching. Abstract | Full Text | PDF (286 kb)

• The Dynamics of Signal Triggering in a gp130-Receptor Comple...

  The Dynamics of Signal Triggering in a gp130-Receptor Complex Structure, Volume 15, Issue 4, 1 April 2007, Pages 441-448 Rishi Matadeen, Wai-Ching Hon, John K. Heath, E. Yvonne Jones and Stephen Fuller Summary gp130 is a shared signal-transducing membrane-associated receptor for several hematopoietic cytokines. The 30 Å resolution cryo-electron microscopy (cryo-EM) structure of the Interleukin 11(IL-11)-IL-11 Receptor-gp130 extracellular complex reveals the architecture and dynamics of this gp130-containing signaling complex. Normal-mode analysis reveals a repertoire of conformational changes that could function in signal triggering. This suggests a concerted mechanism of signaling involving all the components of the complex. This could provide a general mechanism of signal transfer for cytokines utilizing the JAK-STAT signaling cascade. Summary | Full Text | PDF (754 kb)

• Stretching to Understand Proteins—A Survey of the Protein Da...

  Stretching to Understand Proteins—A Survey of the Protein Data Bank Biophysical Journal, Volume 94, Issue 1, 1 January 2008, Pages 6-13 Joanna I. Sułkowska and Marek Cieplak Abstract We make a survey of resistance of 7510 proteins to mechanical stretching at constant speed as studied within a coarse-grained molecular dynamics model. We correlate the maximum force of resistance with the native structure, predict proteins which should be especially strong, and identify the nature of their force clamps. Abstract | Full Text | PDF (415 kb)

• …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 94, Issue 9, 3424-3435, 1 May 2008

doi:10.1529/biophysj.107.120733

Biophysical Theory and Modeling

Toward a Molecular Understanding of the Anisotropic Response of Proteins to External Forces: Insights from Elastic Network Models

Eran Eyal and Ivet Bahar^,  

Department of Computational Biology, School of Medicine, University of Pittsburgh, Pittsburgh, Pennsylvania 15213

Address reprint requests to Ivet Bahar, Department of Computational Biology, University of Pittsburgh, W 1040, BST 200 Lothrop St., Pittsburgh 15261. Tel.: 412-648-3332.

Processes in living cells depend on the mechanical properties of biomolecules in addition to their chemistry 1. Biomolecules are often subjected to functionally required mechanical pressures, e.g., within muscle fibers, microtubules, and molecular motors, in addition to their interactions with other molecules in the cell. The evolution of biomolecules presumably led to an optimization of their mechanical behavior to fit their biological function. A recent interesting example is the direct relation between the functional mode of the action of DNA gyrase and the applied mechanical stress 2.

Recent advances in single-molecule atomic force microscopy (AFM) and optical tweezers techniques allow us to examine the response of proteins to uniaxial tensions 3,4,5. The muscle protein titin, for instance, has been extensively investigated using both AFM methods 6,7,8 and optical tweezers 9,10. Experimental studies have also been conducted on proteins such as T4 lysozyme 11, bacteriorhodopsin 12,13,14, a Na⁺/H⁺ antiporter 15, and others 16. Unfolding forces in different pulling directions have been measured for green fluorescent protein (GFP) 17, ubiquitin (Ub) 18, and the lipoyl domain (E2lip3) of acetyl transferase subunit E2p 19. Strikingly, significant differences have been observed between the responses of the same molecule to pulling along different directions, which apparently reflect path-dependent, nonequilibrium events, rather than the passage over an energy barrier (ΔG_unfolding), which is theoretically expected to be independent of the unfolding pathway 20.

The data collected using single-molecule force spectrometry boosted the field and made it possible to examine the role of structural features such as secondary structure composition and orientation of hydrogen bonds in determining the unfolding forces 14,18. However, not until recently was it possible to obtain data on the mechanical behavior of proteins in response to different deformation directions (induced by exerting uniaxial tensions at particular pairs of residues) due to the time and labor limitations of these complicated experiments.

Molecular insights on the origins of mechanical responses have been inferred to some extent from theoretical studies. In particular, molecular dynamics (MD) simulations have been advantageously employed for estimating the mechanical unfolding forces and the unfolding paths triggered by forces applied along well-defined pulling directions. Titin, in particular, became a model system for understanding the relation between mechanical stability and biological function by steered 21,22,23,24 and quasiequilibrium 20 MD simulations. MD studies have been performed on other small globular proteins as well, which drew attention to the strong dependence of mechanical stability and unfolding pathways on the linkage through which the force is applied 18,25. For a comprehensive review of simulations to explore the molecular origins of observed mechanical resistance and the possible relation to biological function, the reader is referred to the recent article by Sotomayor and Schulten 5.

Notably, MD simulations highlighted the importance of the native contact topology in determining the stress-induced unfolding behavior of proteins 26. A reasonable agreement, mainly qualitative, with experiments could indeed be achieved even when using coarse-grained (residue-level) models 27 or simple potentials 28 based purely on native contact topologies. An extensive survey of the responses of Protein Data Bank (PDB) structures to uniaxial stretching at their N- and C-termini, predicted by coarse-grained simulations, was recently published by Sulkowska and Cieplak 16, which also provides a comprehensive compilation of experimental data (see http://info.ifpan.edu.pl/BSDB).

Here we demonstrate the use of an analytical methodology, the anisotropic network model (ANM) 29,30, to construct a complete map of the mechanical response of all residue pairs in a given protein to uniaxial deformation. We present results for three proteins (Table 1). The results illustrate how the ANM, recently shown to serve as a useful tool for efficient analysis and visualization of the conformational dynamics and anisotropic fluctuations of PDB structures 31, also captures the experimentally observed anisotropic response of proteins to external stresses. First, we present the use of the ANM to derive the effective force constants associated with uniaxial deformations along any pair of residues. Then, the predictions are compared with available experimental data on GFP from jellyfish, Ub from human, and E2lip3 domain of pyruvate dehydrogenase (PDH) complex from Escherichia coli. We discuss limitations of the approach, specifically the lack of residue specificity and the dependence on equilibrium structure. The former prevents its applicability to proteins where the sequence identity, rather than the overall fold, dominates the behavior; and the latter restricts the theory to deformations near native state.

Most of the results presented here compare the unfolding forces where the same protein is stretched in different directions. The major utility of ANM indeed lies in providing information on the relative mobilities of residues or on the relative responses of different residue pairs under tension, rather than predicting the absolute sizes/strengths of deformation. The comparison of the results for different proteins is complicated by the necessity to calibrate their spring constants and by the varying conditions (e.g., pulling velocity) in different experiments. Yet, some results on different proteins are also presented and discussed.

Two major utilities of the approach are its computational efficiency (it lends itself to deterministic assessment of stress-strain behavior for all residue pairs) and its simplicity, which allows a clear interpretation of the results. The theory essentially delineates the behavior of the protein in the neighborhood of its original (equilibrium state) or early stages of unfolding. Yet a high correlation is observed between the computed resistance to deformation and that which is experimentally observed. This correlation is discussed in light of the shape of the energy landscape near the original energy minimum and the kinetic accessibility of particular deformation directions.

We represent the protein by a structure of N sites, the instantaneous position vectors of which, R_i (1≤i≤N) are identified by the C^α-atoms. Consider two residues, i and j, originally separated by a distance vector:

(1)

where and are the equilibrium position vectors of the two residues (Fig. 1). Suppose an external force (uniaxial tension) F_ij is exerted on them to increase the interresidue distance by a deformation vector d_ij:

(2)

Display large version of this figure

Figure 1

Schematic description of the position vectors at equilibrium, and and their instantaneous deformations induced by mode k, and Stretching of residues i and j (green arrows) gives rise to a deformation vector, d_ij, the direction of which coincidences with that of the equilibrium distance vector, The distance vector induced by mode k is designated

In the absence of external forces, under equilibrium conditions, the structure has 3N−6 internal degrees of freedom and enjoys 3N−6 normal modes of motion. The changes in coordinates driven by these modes will be denoted 1≤k≤3N−6. These modes are conveniently determined by normal mode analysis (NMA), i.e., by the eigenvalue decomposition of the Hessian matrix:

(3)

H is a 3N×3N matrix of the second derivatives of the potential with respect to coordinates, λ_k denotes its kth nonzero eigenvalue, the superscript T designates the transpose, and u^(k) is the kth eigenvector given by

(4)

The ANM lends itself to a simple formulation of H where the respective off-diagonal and diagonal super elements (3×3 matrices) read 29,30,31

(5)

(6)

(7)

The eigenvalue λ_k corresponds to the curvature of the potential along the normal mode k, and hence takes on, in the mode space, the role of the spring constant of the physical space 34,35. Evaluation of absolute λ_k values requires knowledge of force constant γ, whereas their relative values/dispersion and eigenvectors are independent of γ. γ is usually estimated from the experimentally detected ms fluctuations of residues (e.g., B-factors in x-ray structures or deviations in residue positions between NMR models) using the relationships B_i=(8π²/3)〈(ΔR_i)²〉 and 〈(ΔR_i)²〉=trC_ii where trC_ii designates the trace of the ith diagonal superelement of C32,34,35,36,37.

The frequency of the kth mode scales with such that slower modes make larger contributions to observed fluctuations (see Eq. (7)). They also define the directions of deformation that incur the lowest internal resistance. If a pulling direction d_ij coincides with a relaxation mechanism favored by a slow mode, then the work done against these relatively “soft springs” would be expected to be smaller, and this will be reflected by an overall small effective “system-based” spring constant. In this respect, it is of interest to examine the correlation cosine between the direction of the externally applied deformation d_ij and the change in interresidue distance, intrinsically favored by mode k:

(8)

Noting that i), the deformation direction coincides with that of the equilibrium distance vector and ii), scales with the difference between the jth and ith super elements (three-dimensional vectors) of u^(k) (Eq. (4)), Eq. (8) may be rewritten as

(9)

The contribution of the kth mode to the deformation reads

(10)

(11)

The macroscopic force is written as a weighted sum over all these contributions such that the effective force constant accompanying the observed deformation becomes

(12)

In the following, 〈κ_ij〉 will be interchangeably referred to as effective spring constant or mechanical resistance. We will examine how the values predicted for three different proteins and different residue pairs in the same protein compare with experimental data.

GFP is an α+β class protein of N=238 residues originally isolated from jellyfish, having a β-barrel structure of 11 strands and a helix that contains a chromophore. Its unique property to fluoresce green light when exposed to blue light makes it an extremely useful probe in biology. Its mechanical properties have been examined in a series of studies 17,38,39. In particular, the anisotropic response of GFP to uniaxial tension along five different deformation directions (Fig. 2) was determined using AFM 17 to observe a strong dependence on the pair of residues on which the forces were applied. Unfolding forces were reported to vary in the range 116≤|F_ij|≤548pN, the lower and upper limits corresponding to the respective residue pairs (i,j)=(Lys-3, Asn-212) and (Asp-117, Tyr-182).

Display large version of this figure

Figure 2

Distribution of the deformations in the distance contributed by each mode k (abscissa) in the presence of extensional forces applied to residues i and j. Each panel corresponds to a particular pair of GFP residues examined in previous experiments 17: (A) 3–132, (B) 3–212, (C) 182–212, (D) 132–212, and (E) 117–182, shown on the right panels. Note the distributions sharply skewed toward the lowest frequency modes in B and D, the same trend to a weaker extent in A, and the relatively uniform distribution in C. The first two cases refer to deformation directions that are more “yielding” as they can be accommodated by low-frequency modes. E, on the other hand, points to the involvement of moderate-to-high frequency modes and thereby the need to apply relatively stronger external forces to induce the same level of deformation. The ordinate scale refers to the force constant γ=0.25 kcal/(molŲ), and the profiles are independent of γ. The cartoon diagrams here and in all figures were generated using Jmol 68.

We calculated the contributions of all modes to the macroscopic force for each pair of residues experimentally examined using Eqs. (10) and (11). The values resulting from all modes are displayed in the five panels of Fig. 2 as a function of mode index. Clearly, the distributions among different modes exhibit a strong dependence on the selected pairs of residues where the external tension is applied. Although in some cases the slower/softer modes make relatively larger contributions (e.g., Lys-3–Asn-212 and Glu-132–Asn-212), others exhibit more uniformly distributed contributions from different modes (e.g., Tyr-182–Asn-212). Notably, a single mode (first mode) is observed in the case of the Lys-3–Asn-212 to induce a deformation of >1Å along the extension direction, whereas the contributions of individual modes generally remain lower than 0.1Å. Note that the absolute size of the modes are based on the normalization of the force constant after the B-factors (see Theory and Methods). Fig. 2 thus provides a detailed description of the deformations along selected directions inherently accessible to GFP as a consequence of its natural vibrational motions in the absence of external forces.

The filled circles in Fig. 3 show the correlation between the theoretical 〈κ_ij〉 (see Eq. (12)) and the unfolding forces deduced from experimental measurements 17. A correlation coefficient of 0.94 is observed between the two sets of data. Due to the small number of points, the real correlation might be weaker, yet p-value analysis indicates that it is still significant (p=0.01 to be obtained by chance). Evidently, when a given deformation direction R_ij can be accommodated by moving along relatively softer modes, 〈κ_ij〉 takes on smaller values and vice versa.

Display large version of this figure

Figure 3

Correlation between the theoretical effective force constant 〈κ_ij〉 (ordinate) and experimentally reported spring constants (abscissa) 17,69 for the five studied extensions of GFP. Theoretical spring constants are evaluated using Eq. (12). The theory yields spring constants that are about 10 times softer, attributed to local deformations, rather than those, global, experimentally detected. Data points indicated by open circles refer to experiments 40 performed with two permutations of EYFP (cut at 144/145, higher point (outlier); 173/174, lower point). The gray circle refers to the EYFP wild-type protein.

The same figure displays the results taken from the work of Jimenez et al. 40 on enhanced yellow fluorescent protein (EYFP) using the wild-type structure (shown by a star) and two circular permutations (open circles). EYFP is structurally equivalent to GFP. The two proteins differ by only five substitutions. The two circular permutations each contain a linker of six residues that connect the original C- and N-termini, whereas the peptide bonds at residues 144/145 and 173/174 are broken in the respective structures. These are conceptually hard targets for our ANM-based approach because the forces are applied at residues that are originally covalently bonded but now serve as chain termini. Note that the introduction of a six-residue linker and the dissociation of a peptide bond may have induced some structural changes which would affect the model and predictions. To apply our method to these cases, we have implicitly assumed that the fold does not change with respect to the wild-type GFP structure.

The results show that the theoretical approach overestimates the force in one of the circular permutations. This may be due to some internal relaxation, succeeding the permutation, which the ANM does not incorporate. We also note that the pulling velocity in the study of Jimenez et al. 40 was considerably slower (0.4μm/s vs. 3.6μm/s) than that used by Dietz et al. 17 for GFP. This difference implies that the forces measured with EYFP would be lower than those measured for GFP 17. However, this effect is relatively small, due to the logarithmic dependence of changes in forces on the ratio of pulling velocities. Therefore, the apparent discrepancy between theory and EYFP experiments cannot be explained by the difference in pulling velocities alone.

Although there is a good correlation between the relative sizes of effective force constants and the multidirectional unfolding forces obtained for GFP, the absolute values computed by the ANM (ordinate in Fig. 3) are about one order of magnitude smaller than those inferred from experiments (abscissa). This difference may be attributed to the fact that the former refers to small deformations near the folded state, whereas the latter corresponds to unfolding events. Furthermore, this approach implicitly assumes that the observed mechanical behavior is dominated by the equilibrium state contact topology, and this may not be the case if the transition state between the folded and unfolded state is not close to the native state, or there may even be multiple states/barriers involved in the unfolding event. A typical example is titin I27 which has been shown by Schulten and co-workers to unfold in multiple steps 41,42. The approach here can be used to examine first (closer to native state) transition (dissociation of strands A-B), rather than events occurring at a later stage (e.g., breaking the A′-G patch) in this case. Due to such effects, the difference in the magnitudes of the two sets of force constants is not, therefore, surprising, but their high correlation is, as will be further discussed below.

Using the ANM, we can readily construct a complete map of the mechanical resistance in response to all possible pulling directions (Fig. 4). Efficient assessment of such maps is an advantage of analytical models such as the ANM over numerical approaches such as steered MD or even coarse-grained simulations. The map shows that the residues that belong to secondary structural elements tend to exhibit relatively strong resistance to deformation. The curve under the map displays the results averaged over all pairs for each residue, which provides a profile of the mechanical resistance of individual residues to deformation, in general. Some residues, especially those lying at the loops connecting strands 10–11 and 8–9 are clearly more disposed than others to deformation, as indicated by this profile.

Display large version of this figure

Figure 4

Mechanical resistance map for GFP obtained by calculating the effective force constant 〈κ_ij〉 in response to uniaxial extensional forces exerted at each pair of residues. The secondary structure of the protein is shown along the upper abscissa (arrow, β-strand; zigzag line, α-helix). The profile at the lower part of the map displays the mean resistance of each residue, averaged over all values in a given column.

It is also interesting to see which molecular directions are more mechanically resistant to uniaxial tension. Figure 5A displays the distributions of effective force constants 〈κ_ij〉 as a function of the angular deviation of the exerted tension direction away from the cylindrical axis of the β-barrel. For clarity, the histogram for residue pairs that can be most easily pulled (which exhibit the lowest 1% 〈κ_ij〉 values) and that of the most resistant (the top 10% 〈κ_ij〉) pairs are displayed. Only pairs that are separated by more than 25 intervening residues along the sequence and more than a 25Å internode distance and that are part of the “barrel” fold of the protein are included. The protein is observed to be more easily deformable along directions which coincide with its cylindrical axis, whereas radial directions exhibit higher resistances. Examples of residue pairs that exhibit low mechanical resistances are shown in Figure 5B.

Display large version of this figure

Figure 5

Relation between the mechanical resistance (〈κ_ij〉) and the direction of interresidue vector with respect to the cylindrical axis of GFP. (A) Two number distributions are shown for two sets of residue pairs that exhibit opposite behavior: those yielding the lowest 〈κ_ij〉 values in the 1% range (black bars) and those in the upper 10% 〈κ_ij〉 range (gray bars). Clearly, residue pairs which exhibit lower resistance to deformation are oriented along the cylindrical axis (mean at 22°), whereas residue pairs distinguished by their strong resistance are oriented at more perpendicular directions (mean angle of 55° with respect to the cylindrical axis). (B) Illustration of the location of some residue pairs that are predicted to exhibit very low resistance against stretching.

Ub is another (α+β) protein that has been examined by AFM 18 (Figure 6A). This is a compact (76 amino acids) protein, highly conserved in evolution. It has some essential roles, especially in targeting proteins to be degraded but also in other cellular signaling processes 43. The unfolding force of Ub when pulled at the C- and N-termini (Met-1–Gly-76) was measured 18 to be larger than that along the (Lys-48–Gly-76) direction (203 pN vs. 85 pN, respectively). The ANM results for the effective spring constants 〈κ_1,76〉 and 〈κ_48,76〉, on the other hand, are 6.00 N/m and 4.23 N/m, respectively. Calculations yield deformations of the order of 0.15Å (average over all modes), using γ=1.75kcalmol⁻¹Å⁻² deduced from the B-factors reported in the PDB (1ubi 44). The counterpart of Fig. 4 for Ub is presented in Figure 6B.

Display large version of this figure

Figure 6

Mechanical behavior of Ub. (A) Cartoon representation of Ub. Pulling directions which are relevant to the biological function or examined by Carrion-Vasquez and colleagues 18 are indicated, and the pulled residues are labeled. (B) Mechanical resistance map for Ub, equivalent to the one shown in Fig. 4 for GFP. (C) The effective force constant 〈κ_ij〉 for the extensions in the directions shown in (A).

Ub is particularly attractive for analysis since its mechanical resistance is likely to be directly relevant to its biological function. To tag proteins in the Ub-proteosome degradation system, the C-terminal Gly of Ub forms an isopeptide bond with a selected lysine on the target protein. Then, often, additional Ub molecules attach by forming isopeptide bonds between their C-terminal glycines and the surface lysine side chains of the preceding Ub molecules. It was shown that the linkage form of the poly-Ub chain determines the fate of the target molecule in the proteasome. For example, K48-linked chains usually attach to molecules targeted for degradation; K29-linked chains and K11-linked chains do so also 45. On the other hand, K63-linked Ub chains are attached to molecules which are not eventually degraded and serve as a signal for other cellular processes. Several nonspecific processes in the proteasome, before the degradation, involve ATPase-driven mechanical pulling and unfolding of the target protein. The mechanical properties of the Ub chains play a role in these events 46,47,48, and it is suggested that a minimal level of mechanical resistance is necessary for a specific linkage to be functional 18.

Figure 6A shows the five biologically relevant pulling directions between the C-terminal glycine (Gly-76) and four lysines (11,29,48,63) and the N-terminus (Met-1). The theoretically predicted 〈κ_ij〉 values for these extension directions are presented in Figure 6C. The pair G76–K48 exhibits a moderate mechanical resistance, which is apparently strong enough to maintain the poly-Ub chain intact, whereas the target protein is processed by the ATPases in the 19S unit of the proteasome. The G76–K11 direction, on the other hand, shows relatively low resistance and its role in the process is indeed questionable 45. The pairs (Gly-76, Met-1) and (Gly-76, Lys-63) pairs exhibit relatively high resistances to extension. Interestingly, no N-terminal linked and K63-linked chains have been characterized to date among the poly-Ub chains that target proteins for degradation. This suggests that a strong resistance to uniaxial tension might not present a suitable setup for the degradation machinery. As for the K29 linkage that is predicted to be equally resistant to deformation, we note that although K29-linked chains exist and are involved in proteasome targeting, they are proposed to switch to K48 linkage during extension 45. The exact biological functional role of K11-, K63-, and K29-linked as well as N-terminal linked Ub chains remains to be further explored. The analysis here suggests that overall an intermediate flexibility, lending itself to conformational adaptability while maintaining stability, may provide an optimal framework for ubuiqitination reactions. Yet, many factors other than mechanical stability may affect the selection of the appropriate linkage.

Finally, we note that atomic simulations performed for Ub 18,25 show that the two pulling directions (N-C) and (48-C) yield effective (maximal) forces of 2000 and 1200pN, respectively. These two relative values are in accord with the ratio (6.0vs. 4.2N/m) of the force constants predicted here (Figure 6C).

Yet another protein that has been investigated with respect to its mechanical resistance in different pulling directions is E2lip3, the lipoyl domain (E2p) of the PDH. This domain forms the structural core of PDH and is responsible for the transfer of an acetyl group in the complex. In a recent study 19, its response to pulling in two different directions (Figure 7A) was studied. The mechanical resistance along the Met-1–Lys-41 direction was shown to be significantly stronger than that along the Met-1–Ala-80 direction (177pNvs. <25pN). The unfolding force along the latter was so weak that it was hardly detectable 19.

Display large version of this figure

Figure 7

Mechanical resistance of E2lip3. (A) Cartoon representation of E2lip3. Pulling directions examined by Brockwell and colleagues 19 are indicated, and the pulled residues are labeled. (B) The complete resistance map for E2lip3. The secondary structure of the protein is shown along the upper abscissa and right ordinate (arrow, β-strand; zigzag line, α-helix). The mean resistance of each residue is shown in the profile at the lower part of the map. The inset in the top right corner shows the equivalent resistance map obtained with a Gō potential in coarse-grained MD simulations 27. The color code of the matrix is as in Fig. 4, and in the inset, yellow to blue colors indicate high to weak mechanical resistance, respectively 27. The portion of the map corresponding to the diagram in the inset is indicated by the black triangle.

Our results are in qualitative agreement with these observations. The difference between the two directions is not, however, as large as suggested by experiments, the respective 〈κ_ij〉 values being in the ratio of 1.25:1.

The orientation of hydrogen bonds with respect to the pulling direction was suggested in previous work 19 to explain the anisotropic mechanical response. The departure of ANM results from experimental data may be attributed in this respect to the effect of hydrogen bonds, which are not explicitly taken into consideration in the ANM. On the other hand, a steered MD study 19,49 yielded results similar to those obtained by ANM, showing that the 1–41 direction is more mechanically resistant to deformation but not by one order of magnitude, as implied by AFM measurements.

To investigate the sensitivity of the results to the initial coordinates, we repeated the calculations with various NMR models of E2lip3. The first five models in the coordinates file (PDB code 1qjo) yielded the same qualitative results with force constant ratios with respect to the original model being in the range 1.45±0.25. This shows that at least in this case the results are insensitive to the resolution of the structure and suggests that the method can be applied to low-resolution models.

Figure 7B displays the mechanical resistance map for E2lip3. The inset shows previous results obtained using a Gō model (Fig. 7 in West et al. 27). The comparison with the results here (the lower triangular portion of the map indicated by the black frame) reveals similarities at particular regions (e.g., relatively strong resistance to deformation for residue pairs belonging to the respective strands 4 and 8), whereas other regions (e.g., N-terminal) show different behavior.

In this study, we focused mainly on the relative responses of a given protein to deformations along different pulling directions. However, such data exist for a few systems only. Indeed, the large majority of experiments reported in this area refer to stretching proteins at their C- to N-termini 16. The comparison of our predictions with the unfolding forces reported for different proteins is complicated by several factors. First, to make a theoretical assessment, we need to know the generic spring constant γ corresponding to interresidue interactions in each protein. In principle, this is an adjustable parameter that uniformly rescales the absolute magnitude of residue motions without affecting the shape of the normal modes or the distributions of residue fluctuations and their cross correlations.

Therefore the relative mobilities of residues are uniquely computed for a given protein, but their absolute sizes depend on the γ. The usual approach for a quantitative comparison with experiments is to evaluate γ based on the average B-factors reported for a given protein, and we adopt the same approach here. Note that the B-factors may be sensitive to experimental conditions (crystallization temperature, crystal form/symmetry group) and refinement errors. We also note that the effective forces increase with the pulling velocity, which varies between experiments. This dependence is, however, logarithmic for a given extension 50, which may be comparable to the range of uncertainty in the predictions.

In view of these problems, we restricted our analysis to a subset of proteins that have i), available x-ray structures, and ii), recorded pulling velocities in the relatively narrow range of 0.3–0.6μm/s. The results are shown in Fig. 8. Each point refers to a different protein (see the caption) except for fibronectin, three different domains/repeats of which are displayed: the labels c, f, and g, correspond to the respective domains ¹⁰FNIII, ¹²FNIII, and ¹³FNIII. We have not shown the results for Ub as the corresponding data (6.0N/mvs. 203pN) lie outside the range of the figure. Calculations performed using the high-resolution structure of ¹⁰FNIII available in the PDB yield good agreement with the measurements of Oberdorfer et al. 51 (c), more or less consistent with the relative unfolding forces observed for other proteins (titin, spectrin, EYFP, and protein L). On the other hand, the theoretical forces appear to be weaker than those observed in experiments for the other two repeats, ¹²FNIII and ¹³FNIII. In these two cases, we used the coordinates from the relatively higher (2.8Å) resolution crystal structure of heparin and integrin binding segment of human fibronectin structure 52,53.

Display large version of this figure

Figure 8

Comparison of theoretically predicted force constants (ordinate) and experimentally measured unfolding forces (abscissa) for a series of proteins resolved by x-ray crystallography and subjected to pulling experiments at their N- and C-termini, with a velocity of 0.3–0.6 μm/s. Results are presented for (a) Spectrin 70,71, (b) EYFP 40,72, (c) Fibronectin repeat/domain 10 (¹⁰FNIII) 51, (d) Titin (I1) 73, (e) Protein L 72, (f) ¹³FNIII 12, and (g) ¹²FNIII 12. See Table 1 for the PDB structures used in ANM calculations. The coordinates for ¹²FNIII and ¹³FNIII are taken from the same crystal structure (1fnh), whereas those of ¹⁰FNIII refer to a different PDB structure (1fnf) with considerably smaller temperature factors 52,53.

As the structures of the two repeats are very similar (root mean-square deviation of 1.4Å), the predicted unfolding forces are very similar (Fig. 8, open circles). The experimental unfolding forces differ, however (89vs. 124pN), presumably due to their difference in sequence. Despite the strong structural similarity, the two domains indeed share only 25% sequence identity. This example illustrates a case where the theory here fails. In fact, the ANM exclusively depends on the fold (or contact topology); and as such, structurally homologous structures that exhibit different unfolding forces (due to their different sequences) cannot be explained by the ANM-based model. We note that the different values in this case for ¹⁰FNIII, ¹²FNIII, and ¹³FNIII originate partly from their different fluctuation amplitudes indicated by their B-factors. The B-factors experimentally measured for ¹⁰FNIII (1fnf) average out to 36Ų as opposed to the mean values of 55Ų for the other domains (1fnh).

The analysis here is based on the basic premise that the intrinsic dynamics of a protein near equilibrium conditions has an impact on its anisotropic mechanical behavior. To elucidate the intrinsic, structure-encoded dynamics, we resorted to an NMA with a coarse-grained (anisotropic network) model. NMA yields an ensemble of modes of motion that the protein enjoys under equilibrium conditions, some more easily accessible than others (as implied by their low eigenvalues/force constants). The externally observed (macroscopic) force to initiate a deformation would be expected to be small, to the extent that the external stress complies with those “easier” movements. So, the problem reduced to the examination of the correlation between the direction of the externally applied force and that of the modes inherently induced/favored by the individual structures to make an assessment of the kinetic accessibility of extensions along particular directions.

Clearly, NMA describes the motions in the neighborhood of an energy minimum approximated by a harmonic well. Those motions, which require a transition over an energy barrier or involve any nonlinear effects, cannot be theoretically accounted for by NMA. Therefore, one might not expect to see a correlation between predictions made with a simple NMA (ANM) and the unfolding forces observed for proteins. Yet, a correlation of 0.94 is observed between the two sets in Fig. 3, and qualitative features observed for other proteins are confirmed by the NMA predictions.

Why do results from NMA exhibit a correlation with the data upon “unfolding”? Does the preferential dynamics of the protein near its original (folded) state also affect, if not dominate, the evolution of motions beyond the early stages of deformation? Does a small curvature in the energy landscape along a particular mode direction also entail a lower barrier to be surmounted in many cases? We are not in a position, yet, to answer all these questions; but we present results from simple calculations that suggest that the resistance to deformation experienced at early stages (along a given direction away from the original energy minimum) affects to a large extent the behavior at longer times or larger deformations, at least in the examined cases. The slow modes appear to constitute “nucleation seeds” for unfolding and their stiffness presumably correlates to some extent with the effective forces experimentally measured.

Coarse-grained NMA (e.g., ANM) has been observed in a number of recent applications to sample structural changes beyond those that would be confined to an energy minimum in the full atomic description of the protein, e.g., passages between substates separated by relatively low energy barriers within a given global energy minimum. Classical examples are the transitions between the T and R states of allosteric proteins such as hemoglobin 54,55 and aspartate transcarbamylase 56,57, and many other examples can be found in the literature 58,59,60,61. Interestingly, the intrinsic flexibilities of amino acids predicted by the Gaussian network model were also shown to correlate with the hydrogen-deuterium exchange free energy costs observed in a series of proteins, although experiments were conducted under denaturing conditions 62. There is also a large body of literature in which the native contact topology or associated potentials are utilized to examine folding/unfolding kinetics and mechanisms. We note that in previous studies, such as those of Kleiner and Shakhnovich 28 and Cieplak 63,64,65, realistic results were obtained for unfolding pathways and associated forces, despite the assumption of the ground state of the protein located at the native structure and the simple -type potentials adopted in simulations.

These studies, and observations here, point to the utility of examining the dynamics of proteins intrinsically favored by their fold using coarse-grained models. The models that lend themselves to a unique analytical solution such as the ANM also provide the advantage of a thorough sampling of the energy landscape near the original energy minimum, on a low resolution but broad scale. The results here also reveal that the overall topology of the protein plays a major role in determining the anisotropic mechanical resistance. This also implies that evolutionary related proteins are expected to show similar qualitative mechanical behavior in general.

We should, however, bear in mind that the theoretical data obtained with such coarse-grained models essentially provide a qualitative assessment of the relative behavior of different residues or residue pairs. Effective force constants or displacements were reported not to make a quantitative comparison with their counterparts derived from unfolding experiments but to give an estimate of the stress-strain dependence at early stages of deformation (or uniaxial tension). Besides, we note that the experimentally detected forces increase with pulling (constant) velocity applied in AFM measurements. And, MD simulations, which necessarily deform the protein within timescales much shorter than those occurring in experiments (due to the timescale of simulations), usually yield effective forces that are one order of magnitude larger than those experimentally measured 16. The discrepancy with experimental values decreases only if simulations are conducted under quasiequilibrium conditions 20. In view of these uncertainties in absolute values, we have focused mainly on relative forces (or force constants) associated with different pairs in a given protein.

A major practical advantage of the approach here is its computational efficiency. In molecular simulations, even those which utilize coarse-grained models and simple potentials 16, the estimation of the unfolding force in a given direction takes hours at least. In the approach here, on the other hand, we can estimate the effective force constant, 〈κ_ij〉 for all pairs of residues within seconds. This type of information can be readily tested by AFM and optical tweezers, which may permit us to improve our model and gain a better understanding of the molecular origins of observed behavior. The fast computation time makes it feasible to easily estimate the mechanical resistance in proteins for which there are no experimental data, as well as to generate the complete mechanical resistance maps for proteins of interest. One utility of constructing such maps would be the possibility of carefully designing AFM pulling experiments, i.e., selecting the residues where the external forces should be applied, based on the computational data that can be readily evaluated for any PDB structure.

Despite the relative success of the ANM-based approach here, it is important to stress its limitations. First, rather than predicting an absolute force for deforming a protein, the theory provides a reasonable description of the relative (anisotropic) responses to deformations along different directions in a given protein. The comparison of the unfolding forces of different proteins, on the other hand, is complicated by the differences in the intrinsic (generic) force constants that best reproduce the mechanical behavior of each protein, and by the differences in experimental setups, with the pulling velocities playing a role. Second, the theory is based purely on the contact topology, irrespective of the identity of amino acids.

Although we maintain the view that the overall topology plays a crucial role in mechanical response, as proposed in earlier studies 50,66, it should also be recognized that there exists cases where sequence details become important, and even dominant, such as the fibronectin domains (Fig. 8), the immunoglobulin domains, and recombination of immunoglobulin fragments 67. The theory here cannot explain the differences in the mechanical behavior of proteins that have the same fold but exhibit mechanical responses to deformation. Finally, because the models utilize native state coordinates, the predictions would be expected to agree with experiments to the extent that the transition state is close to the original equilibrium state. Events far from the original state, or transitions that involve multiple barriers and/or intermediates (such as those observed by Schulten and collaborators for titin 41,42), are beyond the applicability range of the theory.