INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The giant protein titin, also known as connectin, is a roughly 30,000-amino acid-long filament that spans half of the muscle sarcomere and plays a number of important roles in muscle contraction and elasticity (Labeit et al., 1997; Maruyama, 1997; Kellermayer and Granzier, 1996; Wang et al., 1993), as well as controlling chromosome shape in the cell nucleus (Machado et al., 1998). During muscle contraction, titin, which is anchored at the Z-disk and at the M-line, exerts a passive force that keeps sarcomere components uniformly organized. The passive force developed in titin during muscle stretching restores sarcomere length when the muscle is relaxed. Titin is composed of ~300 repeats of two types of domains, immunoglobulin (Ig) domains and fibronectin type III (FnIII) domains, and a PEVK (70% proline, glutamic acid, valine, and lysine residues region (Labeit and Kolmerer, 1995). The FnIII domains are located only in the A-band of the molecule, the PEVK region is located in the I-band, while the Ig domains are distributed along the whole length of titin.

The region of titin located in the sarcomere I-band is believed to be responsible for titin extensibility and passive elasticity (Erickson, 1994; Granzier et al., 1996; Greaser et al., 1996). The I-band region of titin consists mainly of two tandem regions of Ig domains, separated by the PEVK region. The Ig domains each form -sandwich structures, while the PEVK region maintains a random coil conformation due to the charges on its glutamic acid and lysine residues. When titin is stretched, the PEVK region unfolds and elongates. Under extreme conditions, the Ig domains in the titin I-band unfold to provide the necessary extension. When forces are released, the unfolded Ig domains refold quickly (Rief et al., 1997; Carrion-Vazquez et al., 1999b).

A recent immunoelectron microscopy study showed that a unique N2B sequence in I-band of rat cardiac titin, located between tandem Ig and PEVK region (Labeit and Kolmerer, 1995), is responsible for muscle extension under physiological conditions ( 0.5 µm per scarcomere) (Helmes et al., 1999), and that the tandem Ig domains in the same titin are kept in the folded state. However, one immunoelectron micrograph (Fig. 7 in the above paper) shows the tandem-Ig domains in different titins of the same sarcomere to be stretched unevenly, some titins being stretched 300 Å longer than the average extension. This observation suggests that Ig domain unfolding occurs in a small percentage of titins during sarcomere stretching in the physiological range and, hence, in part of the normal function of muscle.

Titin Ig and FnIII domains have been observed, using single molecule techniques such as atomic force microscopy (AFM) and optical tweezers experiments, to be protected against strain-induced domain unfolding (Rief et al., 1997, 1998; Carrion-Vazquez et al., 1999b; Kellermayer et al., 1997; Tskhovrebova et al., 1997). AFM experiments have demonstrated that rather strong forces, of the order of 100 pN, must be exerted before Ig and FnIII domains rupture and unfold on a millisecond time scale (Rief et al., 1997; Oberhauser et al., 1998). Other proteins that do not encounter mechanical strain under physiological conditions have been found to exhibit much less resistance against strain-induced unfolding as demonstrated, for example, through AFM experiments on the helical protein spectrin (Rief et al., 1999) and T4 lysozyme (Yang et al., 2000).

AFM experiments on protein domain stretching do not resolve atomic-level detail of the domains' conformational changes during the enforced unfolding. We have used a computer simulation technique, steered molecular dynamics (SMD) (Grubmüller et al., 1996; Izrailev et al., 1997, 1998; Isralewitz et al., 1997; Kosztin et al., 1999; Wriggers and Schulten, 1999; Gullingsrud et al., 1999; Bryant et al., 2000), to complement the AFM observations by providing a detailed atomic picture of stretching and unfolding of individual protein domains (Lu et al., 1998; Krammer et al., 1999; Lu and Schulten, 1999a, b; Marszalek et al., 1999).

The previous SMD studies of force-induced unfolding of titin immunoglobulin domain I27, the sequence and structure of which is shown in Fig. 1, a and b, qualitatively reproduced the two main features observed in AFM force-extension profiles: the number of the force peaks and the distances between the force peaks (Lu et al., 1998). The SMD studies accomplished quantitative agreement with observations in two respects. First, AFM and chemical denaturation experiments showed that a potential barrier exists between the folded and unfolded states, the barrier height estimated to be 22 kcal/mol (Carrion-Vazquez et al., 1999b); simulations that sampled up to 18 SMD runs carried out under identical stretching conditions matched the height of this potential barrier (Lu and Schulten, 1999b). Second, a mechanical unfolding intermediate that elongates titin's Ig domain by 6 to 7 Å and that has been resolved as a `hump` in AFM's force-extension profile, has also been revealed through SMD simulation and reported together with the observations (Marszalek et al., 1999).

View larger version (41K): [in this window] [in a new window]   FIGURE 1   (a) Sequence and secondary structure of titin I27. (b-d) Representative snap shots of the three-phase unfolding scenario of titin immunoglobulin domain 127: (b) native structure, (c) at extension of 10 Å (unfolding intermediate), and (d) at extension of 25 Å. In (b) all backbone H-bonds between -strands A and B and between -strands A' and G are formed. In (c) H-bonds between -strand A and B are broken, but those between -strands A' and G are formed. In (d) all backbone H-bonds between -strands A and B and between -strands A' and G are broken. The two -sheets are colored in green (sheet ABED) and orange (sheet A'GFC), except the residues involved in H-bonds between -strands A and B are colored in blue, and residues involved in H-bonds between -strands A' and G are colored in red. The same color code is used in Figs. 3, 7, 8 and 14. H-bonds are represented as thick dotted lines and broken H-bonds as thin dashed lines. [Protein structures shown in this figure and in the following figures were generated by the program VMD (Humphrey et al., 1996)].

From SMD results and related AFM data we can describe titin I27 unfolding as a three-phase process (Fig. 1, b-d). Phase I is the pre-burst extension. In this phase the H-bonds between -strands A and B are broken and the domain is extended by 6 to 12 Å, depending on the forces applied. Phase I results in a relatively stable mechanical unfolding intermediate (depicted in Fig. 1 C). Phase II is the burst event in which the domain must overcome the dominant potential barrier along the unfolding pathway. The barrier is crossed when all six H-bonds between -strands A' and G are broken (from Fig. 1 c to Fig. 1 d) at an extension only 3 Å longer than that reached in phase I. Phase III is the post-burst extension where the domain can be extended with less resistance.

Recently, several theoretical groups have published studies on forced unfolding of protein domains using approaches and methods different from those applied in the present study. Rohs et al. (1999) studied the stretching of -helix and -hairpin systems using molecular mechanics; they estimated the magnitude of forces involved in the unfolding of these secondary structures of proteins. Socci et al. (1999) and Klimov and Thirumalai (1999) studied the relation between force and the reaction coordinate by stretching proteins described by a lattice model. Evans and Ritchie (1999) modeled the Ig domain unfolding as a single bond-breaking event and approached the problem using Kramers-Smoluchowski theory, which had also been used in analyzing SMD results (Lu and Schulten, 1999b).

Paci and Karplus (1999) studied FnIII domain unfolding by means of so-called biased molecular dynamics. The authors used an implicit solvent model to reduce computational effort and suggested that the two potential barriers resulting from their study, one more than AFM experiments have revealed, are due to van der Waals (vdW) interactions and not due to hydrogen bonds (H-bonds). These authors questioned the simulation in Lu et al., 1998 on the ground that forces needed in simulated Ig domain stretching were 10 times larger than those in AFM experiments. The authors suggested that at the high speed of simulated stretching, e.g., at the speed used in Lu et al., 1998, hydrogen bond forces are overestimated and the role of water molecules breaking interstrand hydrogen bonds spontaneously is ignored, because the respective fluctuations in water-protein interactions occur too rarely to be detected in SMD simulations.

However, SMD simulation at reduced stretching speed (Lu and Schulten, 1999b) reproduced AFM stretching forces within a factor of two when scaled, revealing in all cases the same hydrogen bond-breaking scenario as in the earlier investigation (Lu et al., 1998). The new results argue strongly for hydrogen bond protection as an explanation of the force peak in AFM experiments and of the key molecular mechanism of titin elasticity. In contrast, SMD simulations on the stretching of helical proteins (Lu and Schulten, 1999a) revealed that these proteins, when stretched, can unfold without the requirement to break backbone hydrogen bonds concurrently, exhibiting little resistance against stretching, which is in agreement with observation (Rief et al., 1999; Yang et al., 2000).

In this study we will analyze the H-bond breaking as the key event initiating force-induced domain unfolding. We will show that relevant water-protein interactions fluctuate on the time scale of SMD simulations and demonstrate that water-mediated interstrand (A'-G) hydrogen bond-breaking events precede force-induced Ig domain unfolding.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Stretching forces were applied to titin Ig domain I27 with two SMD protocols: constant-velocity moving restraints and constant force restraints.

SMD using constant-velocity moving restraints simulates the action of a moving AFM cantilever on a protein. One atom of the protein, the C-atom of the C-terminus residue, is restrained to a point in space (restraint point) by an external, e.g., harmonic, potential. The restraint point is then shifted in a chosen direction at a predetermined constant velocity, forcing the restrained atom to move from its initial position in the protein.

SMD simulations of constant force stretching were implemented by fixing the C-atom of the N-terminus residue of the domain I27 and by applying a constant force to the C-atom of the C-terminus residue along the direction connecting the initial positions of the N-terminus to the C-terminus.

Solvation of I27 was modeled as described in previous studies (Lu et al., 1998; Lu and Schulten, 1999b). The resulting structure of the protein-water system with 12,500 atoms was then equilibrated for 1 ns before SMD simulations were performed.

The molecular dynamics simulations were carried out using the programs NAMD (Nelson et al., 1996) and X-PLOR (Brünger, 1992) with the CHARMM22 force field (MacKerell, Jr. et al., 1998), using an integration time step of 1 fs and a uniform dielectric constant of 1. The non-bonded Coulomb and vdW interactions were calculated with a cutoff using a switching function starting at a distance of 10 Å and reaching zero at 13 Å. The TIP3P water model was used for the solvent (Jorgensen et al., 1983). The atomic coordinates of the entire system were recorded every picosecond. The simulations presented in this work will be referred to as SMD-(velocity value), e.g., SMD-(1.0 Å/ps), for constant velocity simulations and as SMD-(force value), e.g., SMD-(1000 pN), for constant force simulations. Multiple runs under the same condition are differentiated by a subscript, e.g., SMD-(750 pN)₂.

The data analysis, including analysis of distances and interaction energies between the hydrogen bond partners, was performed with the program XPLOR. During the SMD simulation, the hydrogen bond energy is implicitly included in the CHARMM22 force field by appropriate parametrization of the partial charges and the vdW parameters (MacKerell, Jr. et al., 1998). An explicit hydrogen-bonding energy term, however, was used in the trajectory analysis, with the parameters adopted from param11.pro in XPLOR. This energy term depends on the distance between the hydrogen bond partners and the donor-hydrogen-acceptor angle. The minimum energy of the hydrogen bond nitrogen-hydrogen-oxygen (NHO) is 3.5 kcal/mol. This energy is reached when the OH distance is 1.9 Å and the NHO angle is 180°. When the OH distance is larger than 3 Å, or when the NHO angle reaches 90°, an H-bond is considered broken.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Titin I27 has been simulated with 1000 ps free dynamics at 300 K, two constant velocity SMD simulations with pulling speeds 0.1 and 0.5 Å/ps, and four constant force SMD simulations with pulling force 1200 pN, 1000 pN, and 750 pN (twice).

Equilibration

During the 1000 ps equilibration run at 300 K the overall structure remained stable, as reflected by a root-mean-square-deviation (RMSD) from the starting structure of ~1.5 Å for heavy atom positions.

The oxygen hydrogen (OH) distances of the interstrand backbone H-bonds between -strands A and B and between -strands A' and G are presented in Fig. 2. There are three H-bonds between A and B and six H-bonds between A' and G. The diagrams show that all the hydrogen bonds, except K85(O)V13(H) and V13(O)K87(H), break and reform more than once during the equilibration. H-bonds E3(O)S26(H) and K87(O)V15(H), the bonds nearest to the domain termini, are of marginal stability during the last 500 ps of the equilibration. For the six H-bonds between A' and G, the middle four bonds are very stable, but the two H-bonds at both ends of this -sheet, i.e., H-bonds Y9(O)N83(H) and KST(O)V15(H), appear to be more flexible. Fig. 3 shows a breaking event for H-bond Y9(O)N83(H). At 650 ps of the equilibration run, this backbone H-bond is still formed; the closest water molecule is at least 3 Å away (Fig. 3 a) from either of the H-bond partners. At 680 ps, a pair of water molecules approach and form H-bonds with Y9(O) and N83(H) separately (Fig. 3 b), and at this time the backbone H-bond between Y9(O) and N83(H) is broken. At 750 ps, these water molecules have left and the backbone H-bond is reformed (Fig. 3 c).

View larger version (58K): [in this window] [in a new window]   FIGURE 2   O-H distance of the backbone H-bonds versus time during the 1000 ps equilibration. The data for H-bonds between -strands A and B are placed in the top row; the data for H-bonds between -strands A' and G are placed in the two bottom rows. The H-bond partners are labeled in each graph.

View larger version (19K): [in this window] [in a new window]   FIGURE 3   Snapshots of backbone H-bonds between -strands A' and G during the equilibration. (a) At 650 ps, all backbone H-bonds are formed: water molecules are at least 3 Å away from the H-bond partners. (b) At 700 ps, the backbone H-bond between Y9(O) and N83(H) is broken; two water molecules form H-bonds with Y9(O) and N83(H). (c) At 730 ps, backbone H-bond Y9(O)N83(H) reforms. In all figures only the backbone atoms and the hydrogen atoms involved in backbone H-bonds are shown; oxygen is drawn in purple and hydrogen in black; the remaining backbone atoms are drawn in red, matching the color code in Fig. 1; the residues involved in backbone H-bonds are marked in (a).

Constant velocity stretching

The SMD simulations using a constant velocity stretching protocol closely mimic the force application that takes place in an AFM experiment. The force extension profiles (Fig. 4) from the simulations presented here are very similar to those reported in Lu et al. (1998) and Lu and Schulten (1999b), both in shape and maximum force values. The force peak values recorded in SMD-(0.5 Å/ps) and SMD-(0.1 Å/ps) runs measure 2100 pN and 1350 pN, respectively. The values are ~10% higher than the peak force values reported in Lu et al. (1998) and Lu and Schulten (1999b). This difference may be due to fact that the newly equilibrated structure resulting from a 1000-ps simulation is more stable than the previously equilibrated one resulting from a 100-ps simulation. In both cases the force peaks occur at extensions between 11 and 14 Å. Fig. 5 shows the RMSD of coordinates of heavy atoms in I27 during simulation SMD-(0.1 Å/ps). The protein is close to the native structure at the force peak region, as indicated by an RMSD value of 2.5 Å from the native structure for the whole protein. If one excludes from the above RMSD calculation the residues in -strands A/A' and G, then the resulting RMSD values are between 1.5 and 1.7 Å at the force peak extension. This number is very close to the RMSD value of the equilibration run, which shows that for its most part I27 is not perturbed at the force peak region. The RMSD value of I27 without -strands A/A' and G increases to above 2 Å only at extensions greater than the force peak region.

View larger version (30K): [in this window] [in a new window]   FIGURE 4   Force extension profiles from SMD simulations with constant velocity stretching. In both cases, the dominant force peak is located around extensions of 13 to 15 Å. Simulation SMD-(0.5 Å/ps) had been stopped after 320 ps, when the extension reached 160 Å. Simulation SMD-(0.1 Å/ps) had been stopped after 700 ps, when the domain extension reached 70 Å. The gray area highlights the region where force peak occurs.

View larger version (34K): [in this window] [in a new window]   FIGURE 5   The RMSD value of heavy atom coordinates of I27 during constant velocity stretching SMD-(0.1 Å/ps) for the first 400 ps of simulation. Two alignment results are presented, one with the whole protein aligned, the other with part of the protein (residues 15-77), i.e., the residues not belonging to the terminal -strands A/A' and G, aligned. The gray area highlights the period corresponding to the force peak region shown in Fig. 4.

Monitoring the individual H-bond distances along the unfolding trajectory, we found concurrent breaking of the H-bonds between -strands A and B and between A' and G at extension of 14 Å. A set of distance versus extension plots showing similar results from a previous simulation has been presented elsewhere (Fig. 6 in Lu et al., 1998).

View larger version (30K): [in this window] [in a new window]   FIGURE 6   Extension profiles from simulations SMD-(1000 pN), SMD-(750 pN)1, and SMD-(750 pN)2. Three phases can be distinguished: I (pre-burst) at extensions from 0 to 12 Å; II (burst) at extensions around 12-15 Å, corresponding to the mechanical unfolding intermediate; III (post-burst) at extensions larger than 15 Å. The lengths of phase II (lifetime of the intermediate) are 200 ps, 900 ps, and 1000 ps for SMD-(1000 pN), SMD-(750 pN)1, and SMD-(750 pN)2, respectively. Snapshots at the times indicated by arrows are shown in Fig. 8.

View larger version (14K): [in this window] [in a new window]   FIGURE 7   Representative snapshots of H-bond breaking events in the -sheet formed by strands A and B during simulation SMD-(750 pN)1. (a) At 50 ps, all three H-bonds are formed. (b) At 100 ps, H-bond E3(O)S26(H) is broken and E3 forms an H-bond with a water molecule. (c) At 300 ps, all three H-bonds in -sheet A-B are broken; these bonds never reform during the remainder of the simulation. In all cases only the backbone atoms and the hydrogen atoms involved in backbone H-bonds are shown; oxygen is drawn in purple and hydrogen in black; the remaining backbone atoms are drawn in blue, matching the color code in Fig. 1.

Constant force stretching

In another set of SMD runs, I27 was subjected to a constant stretching force during the entire simulation. Depending on the magnitude of the applied force, I27 elongates at different rates and to different extensions. The extension profiles from constant force stretching simulations, SMD-(1000 pN), SMD-(750 pN)₁, and SMD-(750 pN)₂, are shown in Fig. 6 and demonstrate a clear three-phase process that corresponds to the three-phase unfolding scenario described in the Introduction. In phase I (pre-burst extension) the domain extends rapidly from 0 to 12 Å, within 30 ps for SMD-(1000 pN) and within 200 ps for SMD-(750 pN)_1,2. This phase, clearly discernable in Fig. 6, involves the breaking of H-bonds between -strands A and B, as shown in Fig. 1 c. In phase II (burst) the domain extension fluctuates between 12 and 15 Å for a time much longer than is spent in phase I, namely 75 ps, 200 ps, 800 ps, and 900 ps for SMD-(1200 pN), SMD-(1000 pN), SMD-(750 pN)₁, and SMD-(750 pN)₂, respectively. This phase, also readily recognized in Fig. 6, corresponds to the lifetime of the mechanical unfolding intermediate and ends with the breaking of six H-bonds between -strands A' and G, as shown in Fig. 1 d. In phase III (post-burst extension) the domain resumes fast extension after the six H-bonds between -strands A' and G are already broken.

The plateau region in phase II of the extension profile, i.e., the lifetime of the mechanical unfolding intermediate, has been modeled using the theory of mean first passage times for a barrier-crossing event (Schulten et al., 1981; Lu and Schulten, 1999b). For the intermediate, characterized through an extension of 12 to 15 Å, I27 experiences a major potential barrier counteracting the stretching force. Constant force simulation trajectories SMD-(750pN)₁ and SMD-(750pN)₂ have been selected for a detailed analysis of the H-bond breaking in the plateau region because I27 remains in this extension for ~1000 ps. The constant-velocity stretching trajectories are not suited for a corresponding analysis because in these simulations I27 passes the extension of 12 to 15 Å (the force peak extension in constant velocity SMD) quickly, i.e., within 60 ps for SMD-(0.1 Å/ps)).

Hydrogen bond breaking

As shown in Fig. 7, the breaking of the three H-bonds between -strands A and B occurs at the end of phase I in SMD-(750 pN)₁. H-bond E3(O)S26(H) breaks at 100 ps; H-bonds K6(O)E24(H) and E24(O)K6(H) break at 200 ps. These H-bonds remain broken for the entire plateau region (phase II) and beyond, implying that they do not contribute to the dominant potential barrier protecting I27 from strong external forces.

Fig. 8, a and b present snapshots of -strands A' and G at the initial stage of the plateau region for SMD-(750 pN)_1,2. At 230 ps in SMD-(750pN)₁, H-bond V11(O)K85(H) is broken while the other five H-bonds connecting the two strands are maintained. One of the broken H-bond partners, V11(O), forms an H-bond with a nearby water molecule, as seen in Fig. 8 a. At 110 ps in SMD-(750pN)₂, two H-bonds, Y9(O)N83(H) and V13(O)K85(H), are broken, while the other four H-bonds are still formed. Three of the four partners of the broken H-bonds form H-bonds with nearby water molecules, as shown in Fig. 8 b. In both cases the -strands A' and G are still linked together through the remaining backbone H-bonds, and the domain does not elongate. Several picoseconds after the bonds are broken, water molecules leave and the mentioned interstrand H-bonds reform; the H-bonds become only permanently broken at the end of the plateau region when all six H-bonds break and the -strands A' and G separate.

View larger version (27K): [in this window] [in a new window]   FIGURE 8   Representative snapshots of H-bond breaking events in the -sheet formed by strands A' and G during simulations SMD-(750 pN)1 and SMD-(750 pN)2. The panel letters correspond to the arrows in Fig. 6. (a) At 230 ps in SMD-(750 pN)1, H-bond V11(O)K85(H) is broken and VII(O) forms an H-bond with a nearby water molecule. (b) At 110 ps of SMD-(750 pN)2, H-bonds 9O83H and V13(O)K87(H) are broken; three backbone H-bond partners, Y9(O), V13(O), and N83(H), form H-bonds with nearby water molecules. (c) At 1000 ps of SMD-(750 pN)1, all six backbone H-bonds are broken. (d) At 1080 ps of SMD(750 pN)2, all six backbone H-bonds are broken. In all cases only the backbone atoms and the hydrogen atoms involved in backbone H-bonds are shown; oxygen is drawn in purple and hydrogen in black; the remaining backbone atoms are drawn in red, matching the color code in Fig. 1. The residues involved in backbone H-bonds are marked.

Fig. 8, c and d present snapshots of -strands A' and G at the end of the plateau region for SMD-(750 pN)₁ and SMD-(750 pN)₂, respectively. At 1000 ps in SMD-(750 pN)₁, all six interstrand H-bonds become broken (Fig. 8 c); similarly, at 1100 ps in SMD-(750 pN)₂, all six interstrand H-bonds become broken (Fig. 8 d). In both cases, water molecules form H-bonds with at least one of the H-bond partners for all six broken H-bonds. Six or seven water molecules diffuse in between strands A' and G. At this time, without the force from interstrand H-bonds to hold I27 together, the -strands A' and G are pulled apart by external forces, after which the backbone H-bonds never reform.

Fluctuation of backbone H-bonds

Examining the SMD trajectories, one readily recognizes fluctuating water molecules that interact with the atoms involved in interstrand H-bonds. During the equilibration, as shown in Fig. 2, backbone H-bonds fluctuate between formed and broken states due to `attacks' of the water molecules, but remain intact for most of the time. To investigate the dynamics of the interstrand H-bonds, we monitored the respective H-bond energies as described in Methods.

Fig. 9 presents the energy versus time plots for H-bonds between -strands A and B for SMD-(750 pN)₁ and SMD-(750 pN)₂. H-bond E3(O)S26(H), the one closest to the N-terminus, breaks at 100 ps without reforming. This H-bond has been observed to be the first to break in every constant force or constant velocity SMD simulation. The H-bonds K6(O)E24(H) and E24(O)K6(H) are found to fluctuate between formed and broken states for several times before finally breaking at around 200 ps, i.e., at the time when I27 reaches the plateau region.

View larger version (35K): [in this window] [in a new window]   FIGURE 9   Backbone H-bond energies versus time for H-bonds between -strands A and B. Results from simulations SMD-(750 pN)1 and SMD-(750 pN)2 are plotted in the first and second row, respectively. H-bond E3(O)S26(H) breaks at 100 ps; the other two H-bonds break at 200 ps.

Fig. 10 presents the key events in the unfolding process in SMD-(750 pN)₁ and SMD-(750 pN)₂, the breaking of H-bonds between A' and G, through energy versus time plots for each bond. Before the final concurrent breaking of these bonds at 1000 ps for SMD-(750 pN)₁ and at 1100 ps for SMD-(750 pN)₂, all H-bonds are found to fluctuate between the formed and broken states. The four H-bonds in the middle of the -sheet, N83(O)V11(H), V11(O)K85(H), K85(O)V13(H), and V13(O)K87(H), are more stable in this regard than the two H-bonds at the ends of the sheet, Y9(O)N83(H) and K87(O)V15(H), i.e., during phase II (lifetime of the mechanical unfolding intermediate) the center bonds fluctuate much less frequently than the terminal ones. After the final concurrent breaking of the six H-bonds, the domain is pulled apart and unravels quickly.

View larger version (66K): [in this window] [in a new window]   FIGURE 10   Backbone H-bond energies versus time for individual H-bonds between -strands A' and G. Data from constant force stretching simulation SMD-(750 pN)1 are plotted in the first and the third row; data from simulation SMD-(750 pN)2 are plotted in the second and the fourth row. All H-bonds break at the end of the life span of the mechanical unfolding intermediate, 200 to 1000 ps for SMD-(750)1 and 200 to 1100 ps for SMD-(750)2, except for H-bond Y9(O)N83(H) in SMD-(750)2, which breaks at 600 ps and does not reform.

An exception among the six A'-G interstrand H-bonds is exhibited by H-bond Y9(O)N83(H), the one closest to the N-terminus of the domain. This bond is found to break in the middle of phase II in SMD-(750 pN)₂ without fully reforming. In SMD-(1000 pN), this H-bond also breaks in the middle of phase II and remains in the broken state for an extended period of time, reforming near the end of phase II before it finally breaks together with the other five H-bonds (data not shown). The instability of this H-bond is also recognized in the equilibration run (Fig. 2). We note that this H-bond is next to a -bulge turn, which connects -strands A and A', and which may cause an instability.

An important characteristic describing I27 domain unfolding is the overall extension of the domain and, thus, monitoring the relationship between H-bond energy and domain extension is of interest. Figs. 11 and 12 present the energy versus extension plots for individual H-bonds of -sheet A-B and -sheet A'-G, respectively. For the H-bonds in -sheet A-B, the bond near the N-terminus, E3(O)S26(H), is broken at 9-10 Å; two other bonds, K6(O)E24(H) and E24(O)K6(H), start to fluctuate between formed and broken states at an extensions of 9 Å and eventually break at extensions of 11 to 12 Å. The six H-bonds between -strands A' and G all break at extensions of 13-15 Å. In both simulation SMD-(750 pN)₁ and simulation SMD-(750 pN)₂ two H-bonds at the ends of the parallel -sheet formed by strands A' and G, Y9(O)N83(H) and K87(O)V15(H), start to fluctuate at extensions of 10-11 Å. Four H-bonds in the middle break at an extension between 14 and 15 Å with little fluctuation between the formed and broken states, except for H-bond V11(O)N83(H) in SMD-(750 pN)₁, which starts to fluctuate at 11 Å before the final burst at an extension of 14 Å.

View larger version (32K): [in this window] [in a new window]   FIGURE 11   Energy versus extension for individual backbone H-bonds between -strands A and B. Data from constant force stretching simulations SMD-(750 pN)1 and SMD-(750 pN)2 are plotted in the first and the second row, respectively. The H-bond energy is calculated with the explicit hydrogen bond term using CHARMM parameters. The extension at which H-bond energies increase to 1 kcal/mol and H-bonds are considered broken, ranges from 10 to 12 Å.

View larger version (51K): [in this window] [in a new window]   FIGURE 12   Energy versus extension for individual backbone H-bonds between -strands A' and G. Data from constant force stretching simulation SMD-(750 pN)1 are plotted in the first and the third row; data from simulation SMD-(750 pN)2 are plotted in the second and the fourth row. The extension at which H-bond energies increase to 1 kcal/mol and H-bonds are considered broken, ranges from 13 to 15 Å.

VdW interactions

Another important energetic contribution, the vdW interaction, has also been monitored during the constant-force stretching simulations. Fig. 13 presents the vdW energy versus extension for SMD-(1000 pN), in which case the simulation had been stopped when the domain extension reached 50 Å. The total vdW energy of the protein is found to increase monotonically with extension, while the total vdW energy between protein and water molecules is found to decrease with extension; the sum of the above two terms remains nearly constant. One can recognize in all three plots that the vdW energy is changing monotonically with extension, which shows that vdW interaction does not provide a significant barrier at the extension range between 13 and 15 Å.

View larger version (39K): [in this window] [in a new window]   FIGURE 13   VdW interaction energies versus extension. Data shown are from simulation SMD-(1000 pN), which had been stopped when the domain extension reached 50 Å. The three plots all show a monotonic increase or decrease along the extension. No jump at the burst extension (12-15 Å) is discernible.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The elastic properties of proteins play an important biological role for single cells and tissues and it is desirable to understand these properties at the atomic level. Elastic properties are not only essential in the case of muscle where titin provides passive force and structural order (Maruyama, 1997; Wang et al., 1993), but in many other cellular functions. For example, the extracellular matrix proteins fibronectin and tenascin play roles in controlling cell adhesion, cell movement, and signaling (Hynes, 1990; Oberhauser et al., 1998). Cell adhesion proteins such as cadherin join cells together in a flexible manner that can sustain cell movement (Nagar et al., 1996). The shape and separation of chromosomes are likely controlled by titin (Machado et al., 1998).

Elastic properties of proteins, even though essential for their cellular functions, have only recently become the subject of intense studies when AFM and optical tweezers permitted the stretching and monitoring of single proteins. AFM has demonstrated that under stretching forces, multidomain mechanical proteins experience one-by-one domain unfolding. These domains, usually Ig and FnIII domains, resist unfolding, as observed on experimental time scales of milliseconds for stretching forces of <100 pN. These observations also revealed that in case of Ig the protection against external forces manifests itself at short extension, e.g., at <10% of the maximum extension of the Ig domain (Carrion-Vazquez et al., 1999b).

Naturally, one seeks to know the molecular mechanism involved in the protection against stretching and unfolding. It is difficult for AFM alone to achieve the goal, but in combination with molecular dynamics simulations that are consistent with observation, one can reveal the mechanism (Rief et al., 1997; Lu et al., 1998; Lu and Schulten, 1999a, b; Carrion-Vazquez et al., 1999b; Marszalek et al., 1999). The simulations showed that interstrand hydrogen bonds protect the Ig domains because the architecture of the protein requires that all hydrogen bonds between strands A and B as well as A' and G are broken nearly concurrently. Only when these bonds have been broken does the conventional unfolding scenario occur in which the remaining interstrand hydrogen bonds break in a zipperlike fashion, i.e., one-by-one, and the unfolding and its reverse, folding, is dominated most likely by van der Waals interactions involving the core of the Ig domain. However, the forces that arise in this later phase are weaker than those required to overcome the hydrogen bond protection. Hence, hydrogen bonds control the initial and key phase of unfolding, whereas van der Waals interactions control the initial phase of refolding, the A-B and A'-G hydrogen bonding only establishing a final `lock.' Domain refolding is currently out of reach of computer simulations, due to the long time scales (milliseconds to seconds) of completing the refolding process. A promising approach to better understand this process is to combine AFM results with chemical and thermal folding/unfolding data (Clark et al., 1999) as reported in Carrion-Vazquez et al. (1999b).

Steered molecular dynamics simulations have indeed been consistent with AFM observations. The main force peak observed in SMD's force-extension profile did not only reproduce AFM observations of single force peak per unfolded domain (Rief et al., 1997, 1998; Oberhauser et al., 1998; Carrion-Vazquez et al., 1999b), but also correctly identified the extension at which the force peak and burst of the domain arise. Upon application of a force between its terminal ends, every Ig domain exhibits a pre-burst increase in contour length of 7-10 Å, as shown in observation and simulation (Marszalek et al., 1999). The main force peak corresponds to a potential barrier that separates the folded and unfolded domains. A first passage time analysis of multiple SMD simulations reproduced the potential barrier height as observed both by AFM and chemical denaturation (Carrion-Vazquez et al., 1999b). An AFM study of pulling 127, the latter constructed with five-Gly insertions, demonstrated that the force peak arises before a straightening out of the F-G loop, but after the detachment of -strand A from the domain (Carrion-Vazquez et al., 1999a). The SMD result that the force peak occurs at the time when H-bonds between -strands A' and G are broken, are consistent with the above AFM data.

The question arises if the direction of the force applied to the protein matters, for example, is the stretching scenario altered if one pulls the C-terminus of a domain not diametrically away from the N-terminus as in previous SMD simulations, but rather pulls it in, say an orthogonal, direction to it. We will solve this problem first on the basis of first principles and then on the basis of test simulations.

Approaching the stated problem from a mathematical point of view, one may answer how quickly the protein treated as a rigid body responds to a force applied in the orthogonal direction to the N- to C-terminus direction, ignoring presently the water shell. Permitting a most simple description of a protein as a sphere of uniform density of radius R and mass M, the geometry of applied forces leads to the rotation of the sphere around a point on its surface while a torque is applied at the point diametrically opposite. This arrangement leads to a rotation of the sphere around an axis through the surface point that can be identified with the z-axis and can be assumed not to change in time. If (t) describes the angle around this axis, then the Euler equation governing the time-dependence of this angle is

(1)

where 2Rf_o is the torque due to an applied constant force that is continuously applied at the same surface point tangentially to the sphere, i.e., at a 90° angle to the C- to N-terminus direction. This is a simplification of the force direction in an actual simulation, but should still provide an order of magnitude estimate for the protein's rotation time. The moment of inertia I in Eq. 1, in the present case, is 1.4MR². Assuming (0) =  0 and (0)  =  0, the solution of Eq. 1 is

(2)

The time t_o it takes to rotate the sphere by 90° is then

(3)

Using M = 10⁴ m_H, R = 20 Å, f_o =  750 pN yields the estimate t_o  10 ps. A test simulation of fixing the N-terminus of I27 and pulling the C-terminus of it with a 750 pN force has been performed. The direction of pulling was always perpendicular to the direction of N- to C-terminus as described in the above theoretical analysis. It took 20 ps for I27 to rotate by 90°.

The estimate derived suggests that inertia effects of the protein due to its mass are so small that the protein can adjust its orientation to the forces usually applied in an SMD simulation within a few picoseconds. Even with the typical AFM pulling force, 100 pN, the rotation takes <30 ps. The calculation neglects two important characteristics of immunoglobulin domains: the inherent elasticity and the possibility that hydrogen bonds are broken before the protein adjusts its orientation as described. The role of these characteristics needs to be checked in actual simulations that we have carried out for a test.

SMD simulations with two different pulling directions have been performed. In one case, the pulling direction was perpendicular to the direction of N- to C-terminus and parallel to the plane defined by -sheet A-G (Fig. 14, a-d). In the other case the pulling direction was perpendicular to the direction of N- to C-terminus and perpendicular to the plane defined by -sheet A29-G (Fig. 14, e-h). For each of the pulling directions described above we have performed one 750 pN constant-force SMD and one 0.5 Å/ps constant-velocity SMD. In all four simulations, the I27 first rotates ~ 60°-80° within 20 ps. During this rotation the key H-bonds between -strands A' and G remain formed, and the domain extension is <10 Å. At this point the I27 started the concurrent A'-G H-bond breakage and then unfolded easily.

View larger version (27K): [in this window] [in a new window]   FIGURE 14   Snapshots of I27 unfolding in constant velocity pulling SMD simulations. In both simulations (a-d and e-h) the N-terminus of I27 is fixed and the C-terminus is pulled along the direction perpendicular to the N- to C-terminus direction. (a) Reorientation of I27 when the pulling direction is within the plane defined by -sheet A'-G. This plane is parallel to the paper. (b) I27 first rotates by ~45° without extension, then gradually extends by 10  while it rotates another 15° (c). From (c) to (d) H-bonds between -strands A and B and between A' and G break concurrently. (d) I27 immediately after the key H-bond breakage. (e) Reorientation of I27 when the pulling direction is perpendicular to the plane defined by -sheet A'G. This plane is perpendincular to the paper. (f) I27 first rotates by 60° without extension, then gradually extends 10 Å while it rotates by another 15 Å (g). The concurrent H-bond breakage between A'-G happens from (g) to (h). (h) I27 immediately after the key H-bond breakage.

The test simulations show that when I27 is pulled in directions other than the one from N- to C-terminus, the domain first aligns itself with the pulling direction. This rotation coincides with phase I of the forced unfolding of I27, the pre-burst extension, as described in the Introduction. Phase II of the forced unfolding of I27, concurrent breaking of H-bonds between  strands A'-G, occurs only after the direction from N- to C-terminus has aligned itself (within 10 to 30°) along the pulling direction. Thus, the protection barrier provided by the H-bonds between A'-G still exist on the unfolding pathway. Unzipping the H-bonds between two -strands, which would result in a much lower protection energy barrier as shown for the forced unfolding of C2 domain (Lu and Schulten, 1999a), does not happen in I27 unfolding even when the original pulling direction is not along the N- to C-terminus.

The nanosecond simulations reported in this paper have clearly shown that water molecules play a key role in overcoming the hydrogen bond protection of stretched Ig domains. Water molecules continuously form and break H-bonds with backbone oxygen and hydrogen, as shown in Figs. 3 and 8. During the equilibration, the backbone-backbone H-bonds are strongly favored over water-backbone H-bonds; when water-backbone H-bonds form, they break rapidly as shown, for example, in Fig. 3. Under stretching forces, the protein has a tendency to elongate and backbone H-bonds are less favored. This gives water molecules a better chance to approach backbone atoms and form water-backbone H-bonds. A comparison of Figs. 9 and 10 to Fig. 2 shows that backbone H-bonds break more often under the action of a force than in the force-free case. In nanosecond simulations, as reported here, the number of break/reform cycles of H-bonds, during the time period before the domain is pulled apart, ranges between and 1 and 20. On the physiological unfolding time scale, these H-bond fluctuations will take place many more times than in the simulations. Longer time-scale simulations with better sampling of H-bond fluctuations can further clarify the fluctuations of backbone H-bonds and the role of water molecules.

When a backbone H-bond is broken during stretching, the two bonding partners may be pulled apart, such as is the case for the E3(O)S26(H) bond, as shown in Fig. 7. This prevents the reformation of this H-bond. However, due to the protein architecture, the six H-bonds between -strands A' and G need to break concurrently before the domain can elongate. Indeed, if one of the six backbone H-bonds is broken and a water-protein H-bond forms, the interstrand H-bond has a tendency to reform, even under strong stretching forces. This is shown in Fig. 8 for the case SMD-(750 pN)₁: when the V11(O)K85(H) bond is broken, the remaining H-bonds between -strands A' and G are still formed, the extension of the domain does not increase, and the V11(O)K85(H) bond quickly reforms. Only when a sufficient number of H-bonds between -strands A' and G are broken simultaneously through attack by water molecules (Fig. 8, b and d) does the domain extend, pulling -strands A' and G apart and preventing the reformation of the H-bonds between these two strands.

One can now understand why the simulations of Paci and Karplus yielded results different from those in Lu et al. (1998) and Lu and Schulten (1999b), and in the present paper. These authors assumed an implicit solvent model that, among other contradictions of their simulations with AFM experiments, does not account for the role of individual water molecules competing for hydrogen bonding to backbone atoms. Our simulations using explicit water demonstrated that the key event in the force-induced unfolding of Ig domains is the water-mediated breaking of interstrand hydrogen bonds. Hydrophobic interactions, which are governed by the exposed surface area, have been found in earlier simulations (Lu et al., 1998; Lu and Schulten, 1999a) to exhibit no distinct feature at the force peak extension where the protection barrier arises. The simulations in Lu and Schulten (1999a) revealed the same monotonic change of hydrophobic interactions as in the case of the unfolding of -helical domains. This result agrees with AFM observations that showed that the response of Ig domains and of helical proteins to external forces is very different (Rief et al., 1997, 1999; Lu and Schulten, 1999a), i.e., only the former exhibit a strong force peak due to H-bond protection.

The H-bond protection suggested in this paper and in previous studies (Lu et al., 1998; Lu and Schulten, 1999b) can be tested experimentally. One straightforward test involves mutation of residues in I27 that can disrupt the interstrand H-bonds involved in the protection mechanism. One such test, in which K6 had been changed to a proline in order to break the H-bonds between -strands A and B, has already been performed (Marszalek et al., 1999). The resulting AFM force-extension curve revealed that the pre-burst stretching phase of I27 was indeed abolished, which is consistent with the SMD prediction (Marszalek et al., 1999). Further tests involving other mutations affecting -strands A' and G are presently being carried out (J. Fernandez, private communication).