J. Phys. Chem. B 2002, 106, 1081-1089
1081
Protein Denaturation in Vacuo. Behavior of the Four-r-Helix Bundle of Apocytochrome c′ under Centrifugal Unfolding Conditions Gustavo A. Arteca* De´ partement de Chimie et Biochimie, Laurentian UniVersity, Ramsey Lake Road, Sudbury, Ontario P3E 2C6, Canada
O. Tapia Department of Physical Chemistry, Uppsala UniVersity, Box 532, Uppsala S-751 21, Sweden ReceiVed: July 12, 2001; In Final Form: NoVember 7, 2001
Experiments and simulations on dry proteins in a low-pressure gas suggest that unfolding and refolding transitions are possible even in the absence of water. Here, we use computer simulations of molecular dynamics to gain further insight into the microscopic mechanisms underlying protein denaturation in vacuo. Recent studies of neutral lysozyme provide an initial guide into the principles at play during the in vacuo unfolding of a multidomain protein. In this work, we contrast its behavior with that of cytochrome c′; this protein has the same chain length as lysozyme, but it is organized as a four-R-helix bundle. Using an ensemble of molecular dynamics trajectories for cytochrome c′, we discuss the qualitative response that follows the imposition of unfolding conditions. The unfolding bias is caused by small, yet systematic, “centrifugal” forces associated with unrestrained spinning in vacuo. The emerging mechanisms appear to depend on the native fold. Whereas lysozyme exhibits a fast denaturing transition, cytochrome c′ exhibits restricted deformations within the same time scale. Two characteristic behaviors are found. In the more common response, the protein undergoes a limited reduction in secondary structure, and relaxation leads back to quasi-native states. In the less common response, cytochrome c′ unfolds partially by losing compactness and tertiary contacts in a correlated manner. Our approach illustrates the usefulness of centrifugal unfolding in vacuo as a tool to detect large-scale configurational motions accessible to the native state. The method provides insight into (a) possible pathways associated with the initial steps of unfolding and (b) a qualitative survey of regions of configurational space connecting the native state attractor with the denatured “state”. In this regard, we show for the first time how the same denatured “state” can be reached not only by direct unfolding of the native state but also by reunfolding of partly folded intermediates.
Introduction Dehydrated proteins provide an important reference system to uncover fundamental aspects of the folding behavior intrinsic to a polypeptide chain. As a result of advances in techniques for transferring, and subsequently characterizing, proteins in a near vacuum,1-6 it is becoming clear that anhydrous proteins exhibit an array of phenomena resembling, to some extent, those seen in solution. Not only does it emerge that the native state can be metastable in vacuo, but also it emerges that unfolded structures can relax into compact nativelike conformers with no need of water. In this context, molecular dynamics (MD) simulations are an important tool to develop a microscopic model for the processes at play. Here, we employ simulation conditions that imitate those prevalent when proteins tumble in a gas migration tube. We study the effect of these conditions over an anhydrous R-helix bundle. Our goal is to map the pattern of configurational and molecular shape transitions accessible to an all-helical native state under the same conditions that lead to rapid unfolding in a multidomain R/β-protein native state, such as that of lysozyme. Our protocol uses “in vacuo boundary conditions”, which allow for unrestricted rotation and translation and expose a * To whom correspondence should be addressed. Telephone: +1 (705) 675-1151, ext. 2117. Fax: +1 (705) 675-4844. Electronic mail: gustavo@ nickel.laurentian.ca.
globally neutral protein to mild denaturing conditions. In contrast to other methods, we avoid the use of biased forces,7-9 residue charging,10 or high temperature.11-13 Instead, the method mimics the effect of thermalization within a low-pressure gas by changing the global coupling to a simulated bath. The procedure, known as “centrifugal unfolding”,14 has recently been used to study the unfolding and relaxation of lysozyme,15-19 a small, 129-residue protein with β-sheet and R-helix content. The emerging pattern of configurational transitions suggests the presence of an unfolding mode in which the loss of compactness is correlated with a reduction in chain entanglement.16,17,20 It could be conjectured that the details of this behavior should depend on the structure of the initial native state. Here, we explore this idea by exposing the R-helix bundle of cytochrome c′ to centrifugal unfolding. We discuss the pattern of configurational transitions that accompany the formation and relaxation of partly unfolded transients. Our results provide insight into the mechanism leading to the denatured state of an R-helix bundle, as well as into essential differences among unfolding and refolding paths.21,22 Computer Simulation of Unfolding and Molecular Shape Characterization It is known that centrifugal unfolding conditions cause rapid denaturation in disulfide-reduced and disulfide-oxidized lysozyme
10.1021/jp012692z CCC: $22.00 © 2002 American Chemical Society Published on Web 01/10/2002
1082 J. Phys. Chem. B, Vol. 106, No. 5, 2002 at room temperature with no directed external bias.14 From the present viewpoint, it is sufficient to regard the method as a tool that reveals structural responses intrinsic to a native fold that is exposed to weak disruptive forces. Nevertherless, it is possible that aspects of the simulated dynamics are akin to a realistic denaturing process in vacuo. For example, lysozyme simulations produce unfolded conformers the mean lengths and cross sections of which resemble those observed experimentally for the protein in flight.23 In addition, in vacuo denaturing takes place in a much shorter scale than processes in solution yet still in time scales comparable to the fastest unfolding processes known for small proteins in solution (ca. 10 ns).24 Detailed simulations of unfolding extracted over a large ensemble of molecular dynamics (MD) trajectories are available so far only for lysozyme.20 Here, we provide a needed comparison by applying the same simulation protocol to the analysis of cytochrome c′ dynamics. The procedure used for the MD trajectories and the subsequent analysis of conformational snapshots by large-scale molecular shape descriptors is discussed briefly below. A. Molecular Dynamics Simulations. The simulation protocol produces an ensemble of constant-temperature MD trajectories for the protein cytochrome c′ in anhydrous form. The trajectories correspond to perturbations on the protein, that is, they are nonequilibrium runs in which we analyze the ability of a macromolecule to change states under denaturing conditions. The set of trajectories is not equivalent to a canonical distribution of accessible configurations at a given temperature, T0. All simulations start from the crystal structure of one monomer of cytochrome c′. This protein is a type-2 cytochrome c involved in reactions with small oxidants and reactants in the purple phototropic bacterium Rhodobacter capsulatus (PDB access code 1cpr).25 In the cell, this protein functions as a dimer of two identical 129-residue chains; each of the monomers is a left-twisted four-R-helix bundle and contains a heme group. Helices are packed in two pairs. Within each pair, the helices are nearly parallel; in turn, each pair is rotated approximately 45° with respect to each other. There are no disulfide or salt bridges. In our case, we consider the apo form of a single monomer (i.e., the heme cofactor is removed) as the starting structure. The heme is not essential to achieve the correct native fold.25 A representative snapshot of the native state appears on the left-hand side of Figure 1 (denoted as “native 1cpr”). This fold can be recognized as distinct in terms of shape descriptors with respect to the native states for other proteins with the same chain length (n ) 129).16,26 The protein is modeled with the Gromos 87 force field27,28 using the 37D4 parameter set. Within this force field, amino acids are described by effective charges placed on polar atoms, but the sum of these charges adds up to zero over each residue.28 Within this potential energy, hydrogen bonding emerges as a result of electrostatic and van der Waals interactions. Because the Gromos 87 parameters for the protein do not depend on the environment, we use this potential energy function for in vacuo simulations. We analyze the behavior emerging from an ensemble of 2-nslong MD trajectories. The length of these trajectories is appropriate for the specific goals of this work, namely, the early response of the protein to the presence of unfolding conditions. Regarding other simulation details, we use no periodic boundary conditions and a 13 Å cutoff in nonbonded (van der Waals) interactions. The “in vacuo boundary conditions” are achieved by not removing global rotations and translations. Two en-
Arteca and Tapia
Figure 1. Selected conformational snapshots of cytochrome c′. The left-hand side structure corresponds to the native state in the crystal. In the central snapshot, we indicate the pattern of global rotation observed during the early part of an MD trajectory in vacuo. The righthand side snapshots show two typical mode-I deformations, in which the native state is essentially conserved in a more stretched R-helix bundle (see text).
sembles of trajectories were generated, one using an integration step of ∆t ) 2 fs and another with ∆t ) 1 fs in the Verlet algorithm. Trajectories with ∆t ) 2 fs lead to similar configuration transitions, but temperature regulation is more unstable. For this reason, we base our conclusions on the results obtained for trajectories with the smaller integration step. Independent trajectories are produced by changing the distribution of initial velocities consistent with the desired temperature for the simulated bath (here, T0 ) 293 K). Centrifugal unfolding is modulated, but not directed, by the presence of a strong coupling to the thermal bath (see below). As in previous works, we use the Berendsen thermostat29 to expose the protein to either unfolding or relaxation conditions. The bath coupling is regulated by the single parameter τ, corresponding to the relaxation constant in Newton’s cooling law, written as: dT/dt ) (T0 - T)/τ, where T0 is the bath temperature and T the instantaneous temperature of the protein. The well-known fact that the Berendsen thermostat does not produce a canonical (NVT0) ensemble distribution29,30 is unimportant for the goals of the present work. What is important is that, even when using a strong coupling to the bath (i.e., a small τ value), the thermostat can change the state of the molecule by a smooth dissipation mechanism with little effect on local dynamics.29 This simulation protocol allows the transfer of energy from internal degrees of freedom to collective vibrorotations, through the coupling of torsional and bending modes.20 This possibility is a key feature that makes centrifugal unfolding effective as an approach to perturb the native state. Within this procedure, global rotations are allowed to develop over time, their magnitude and persistence depending on the bath relaxation τ. In a study of lysozyme, we showed that the resulting rotational energy can be as large as one third of the total kinetic energy.14a As a result of sustained nonrigid spinning about an inertial axis, a protein may be subject to “centrifugal forces” that, although weak, can affect systematically a subset of collective motions. Depending on the molecule, some of these modes may involve large-scale unfolding.
Protein Denaturation in Vacuo
J. Phys. Chem. B, Vol. 106, No. 5, 2002 1083
B. Large-Scale Molecular Shape Features. As was done for lysozyme, we monitor the dynamic behavior of cytochrome c′ in terms of large-scale molecular shape descriptors. For better discrimination, we employ descriptors that recognize differences in chain spatial organization for conformers with the same level of compactness or molecular size. The main part of the present analysis is based on the evolution of absolute and relative chain entanglements, as well as a measure of chain anisometry. All of these descriptors are translationally and rotationally invariant and require no optimized superpositions, as is the case for the standard root-mean-square deviations (rmsd). Within the context of geometrical31-34 and topological35,36 characterizations of 3D curves, a chain is said to exhibit entanglements whenever its orthogonal planar projections produce bond-bond crossings (overcrossings). For a given conformation, K, the probability distribution {AN(K)} will denote the probability of observing N g 0 overcrossings, averaged over a large number (m) of random projections.33 The first moment of the distribution, or mean overcrossing number N h , is a useful absolute descriptor of entanglement:37 Nmax
mN
∑N m mf∞ N)0
N h ) lim
(1)
where mN is the number of projections, among the total m, that yield N overcrossings. For a linear chain with n monomers, we have Nmax ) (n - 2)(n - 3)/2. As a relatiVe descriptor of entanglement, we use the root-mean-square deviation of overcrossing probabilities between the instantaneous conformer K and a reference one, K*:26 Nmax
σA(K,K*) ) [
(K*) 2 1/2 ∑ (A(K) N - AN ) ]
(2)
N)0
In our particular case, we can take K* as the native state of cytochrome c′ in the crystal. When a native chain unfolds, N h decreases up to values that depend on the residual secondary structure. For a fully unfolded elongated conformation, we have N h ≈ 0. Similarly, the relative descriptor σ will increase from typical values corresponding to the statistical fluctuations in a native state and its mutations, estimated as σ(K,K*) ≈ 0.04 ( 0.02.26 Because the overcrossing probability distribution of a rod is A0(K) ) 1 and AN(K) ) 0 for N g 1, the root-squaredeviation in overcrossings will increase as a chain unfolds, reaching up to a maximum of σ(Krod,Kcompact) J 1, when comparing a rodlike and a maximally compact conformation for a very long chain. As a complementary tool, we use the chain asphericity,38 given in terms of the principal moments of inertia, {Ii}, as Ω ) {(I1 - I2)2 + (I1 - I3)2 + (I3 - I2)2)}/{2(I1 + I2 + I3)2}. In spheroidal native states, we expect that Ω ≈ 0, and Ω ≈ 1/4 in rodlike conformers. As in previous works, we will use the twodimensional shape map (N h ,Ω) as a convenient space to monitor changes in chain organization and compactness during unfolding-refolding transitions.16-20 Partial Unfolding of 1cpr-Cytochrome c′: Effect of a Strong Coupling to the Bath A strong coupling to the bath is achieved by using a relaxation constant of τ ) 5 fs. In the case of lysozyme, this coupling leads systematically to unfolding into elongated structures by a very sharp transition at t* ) 600 ( 100 ps for every tested trajectory.20 The behavior emerging from an ensemble of
Figure 2. Three representative trajectories represented in terms of the evolution of entanglement complexity. The trajectories correspond to a mode-I (top), mode-II (middle), and mode-III (bottom) deformation. Only the latter shows a substantial amount of unfolding. Overcrossings numbers provide a convenient tool for geometrical and topological analysis of polymers. In the case of knotted DNA,35,36 overcrossings give rise to topological quasi-invariants. Within the context of this work, overcrossings characterize the geometry and connectivity for chains with transient entanglements. It must be noted that the distribution of overcrossings33 characterizes the three-dimensional organization without direct reference to any conventional definition of what constitutes secondary structure. When convenient, we use conventional displays of secondary structure, but only for the sake of illustration. For clarity, the labels “N” and “C” denote the location of the N-terminus and C-terminus, respectively, in selected snapshots.
trajectories of cytochrome c′ is quite different. The results shown below are based on 50 independent MD runs lasting 2 ns, computed with an integration step of 1 fs. The behavior resulting from 35 additional MD runs with an integration step of 2 fs is qualitatively similar. Initially, cytochrome c′ reacts to the application of centrifugal unfolding conditions in a manner that resembles that of lysozyme.20 The effect is depicted in Figure 1. The left-handside snapshot corresponds to the initial crystal structure. The snapshot in the center is representative of those appearing during the first 300 ps of simulation. After a short induction period (t < 100 ps), the main global motion corresponds to a rotation about the axis indicated in the figure. The amino acid side chains located near the outer protein surface “feel” the effect of the resulting centrifugal forces. As a result, the initially uncorrelated distribution of side chains on the surface adopts a more organized distribution in which these chains tend to orient nearly perpendicular to the main helical axis. If this motion persists over a significant amount of time, a pattern of three typical responses can be elicited, as shown by the evolution of chain entanglements in Figure 2. For convenience, we label the “unfolding modes” in terms of the content of secondary and tertiary structure. However, the molecular shape is analyzed only in terms of large-scale descriptors. A. Mode I: A Limited Deformation of the Native Fold. In many trajectories, the motion corresponding to a rotation about the axis of the helical bundle persists for the entire 2 ns run. The typical results obtained at the end of the trajectory are illustrated by the two snapshots on the right-hand side of Figure 1. The deformation pattern involves very limited unfolding with helices becoming nearly parallel and elongated along the main axis of the bundle. This stretching may conserve the original helical content (right-hand side, top) or involve some loss (bottom). The fact that these are the most prevalent responses to the unfolding conditions indicates the intrinsic stability of the R-helix bundle. (Given that the helix is resistant to
1084 J. Phys. Chem. B, Vol. 106, No. 5, 2002
Figure 3. Evolution of the relative deviation in overcrossing probability for the trajectories highlighted in Figure 2. Deviations in overcrossing distribution are measured with respect to the native state in the crystal. The quasi-native structures correspond to mode-I and mode-II deformations. The unfolded structure is obtained by a mode-III trajectory. (For the termini at some of the snapshots, see Figure 2.)
unwinding or to loop-enlargement by rotation, the only dissipative motion apparently left accessible is a stretching along the main helical axis.) The effect on molecular shape is clear in Figure 2, which shows the evolution of the mean overcrossing number for three selected MD runs. In all cases, the N h value increases with respect to that of the native structure in the crystal. This effect arises as a consequence of electrostriction in a vacuum, that is, a compaction of the chain resulting from the lack of screening of electrostatic interactions in absence of solvent.14 This intermediate can be considered as a metastable “in vacuo native state” (IVNS). The “mode-I” deformation (Figure 2, top trajectory) conserves N h at the value corresponding to the IVNS during 2 ns. Figure 3 complements the analysis by showing the relative overcrossing descriptor σ measured with respect to the crystalline native state. The “mode-I” oscillations observed near σ ≈ 0.03 are typical of equilibrium fluctuations about the native fold.26 B. Mode II: Partial Melting of Secondary Structure. The changes involved in mode-II deformations (as well as those for mode-III below) are signaled by the occurrence of a change in the dominant rotation. In these cases, the molecule passes from a rotation along the main helical axis of the bundle to a rotation in a perpendicular direction. The typical outcome of this new deformation mode is illustrated by the middle trajectory in Figure 2. In this case, the geometry of the native state is essentially conserved for ca. 800 ps, with structures that resemble those found in mode I. However, mode-II trajectories undergo further unfolding as a multistep transition in the overcrossing number. As shown by the snapshot at the end of the trajectory, this transition generates a family of structures resembling bundles of four parallel helices but with a substantial loss of helical content at both N- and C-termini. The significance of this change in molecular shape emerges when looking at relative entanglements (Figure 3). The mode-II deformation shows a drift in σ values during the first half of the trajectory with a maximum at t ≈ 500 ps. During this period, the structure is locally distorted with respect to the native fold. After the sudden change in N h (cf. Figure 2), the σ descriptor returns to values typical of a mode-I trajectory. In other words, the overall shape of a four-R-helix bundle (e.g., the tertiary contacts among
Arteca and Tapia
Figure 4. Three representative snapshots illustrating the first transition in overcrossing numbers during the mode-III deformation in Figure 2. The bottom row corresponds to the same snapshots but viewed from the bottom. Note the sharp rearrangement of the four-helix bundle. (For clarity, the labels “N” and “C” denote the location of the termini in selected snapshots.)
Figure 5. Three representative snapshots illustrating the second transition in overcrossing numbers during the mode-III deformation in Figure 2. Note the loss of helical content during the rearrangement. (The label “N” denotes the location of the N-terminus.)
helices) is maintained, despite the reduction in entanglement that accompanies the partial melting of secondary structure. C. Mode III: Partial Unfolding with a Substantial Secondary Structure. This unfolding mode is illustrated by the bottom trajectory in Figure 2, which undergoes a two-step transition toward less-entangled backbones. This partial-unfolding transition leads to structures that retain part of the initial R-helical content yet involve substantial changes in helical packing. Figure 3 shows a sharp contrast between mode-III unfolding and the other two previous modes in terms of relative deviations of overcrossing probability. The two transitions observed in the mode-III trajectory are associated with σ ≈ 0.10 and σ ≈ 0.08; these values resemble those for the partial unfolding of lysozyme.20 (Extensive unfolding corresponds to σ ≈ 0.20.) The snapshots in Figures 4 and 5 illustrate the nature of the partial-unfolding steps. During the first transition, located at t ≈ 900 ps, the main rearrangement corresponds to the breaking of contacts between pairs of parallel helices. As shown at the bottom of Figure 4, the original bundle geometry persists still at t ≈ 800 ps, but it is broken at t ≈ 900 ps at which three helices are nearly coplanar. Figure 5 shows that the second transition (centered t ≈ 1400 ps) involves only a further loss of secondary structure within the disorganized helical bundle. The final structure obtained at 2000 ps is not fully unfolded yet not trivially distorted. The root-mean-square deviation (rmsd) between the native backbone and that of the final partly unfolded transient is 10.3 Å. The above pattern of partial unfolding is representative of the general behavior observed within the ensemble of MD trajectories. Figure 6 collects all data for the change in entanglement complexity. On the right-hand side, we give the distribution of N h values corresponding to last conformational snapshots for each trajectory (at t ) 2 ns). The dominant effect
Protein Denaturation in Vacuo
Figure 6. Results for a set of 50 independent MD trajectories for centrifugal unfolding of cytochrome c′ using τ ) 5 fs. The histogram at the right-hand side represents the population distribution of N h values for the last conformation of every trajectory (arbitrary units). The majority of the trajectories produce very limited unfolding.
Figure 7. Results for a set of 50 independent MD trajectories for centrifugal unfolding of cytochrome c′ using τ ) 5 fs in terms of relative deviations in overcrossing probabilities. The structures indicated as “quasi-native” correspond to mode-I and mode-II deformations.
of imposing centrifugal-unfolding conditions yields mode-I deformations, leading to values of entanglement complexity typical of a slightly compressed form of the native state (N h ≈ 49). Only 10% of the trajectories produce a reduction in mean overcrossing number below that of the native state (N h ≈ 45). Even then, not all of these deformations involve substantial changes with respect to the native structure, as shown by the data for the relative descriptor σ in Figure 7. From Figure 7, we can conclude that most of the trajectories produce conformers with quasi-native shape features (σ e 0.05); only two of the trajectories involve substantial unfolding. The MD runs generated with an integration step of ∆t ) 2 fs show a comparable unfolding yield (data not shown). In conclusion, we estimate a probability of unfolding not larger than 5% for the response of the native state of cytochrome c′ to the presence of centrifugal unfolding conditions. This is a much lower value than the one observed in lysozyme (nearly 100% unfolding yield).20 Unfolding Pattern in (N h ,Ω) Shape Space Further insight into the mechanism of centrifugal unfolding can be obtained from the pattern of molecular shapes spanned by the configurational transitions. To this end, we use the entanglement-and-anisometry space (N h ,Ω), which allows us to describe MD trajectories in terms of changes in compactness and large-scale geometrical organization of the backbone (or,
J. Phys. Chem. B, Vol. 106, No. 5, 2002 1085
Figure 8. Anisometry-and-entanglement map for two typical modeIII (black line) and mode-II (gray line) trajectories. Both trajectories exhibit a degree of correlation between these two properties during the part of the trajectory associated with significant unfolding events. Nevertheless, the two processes can be clearly separated.
loosely speaking, the “protein topology”). These two properties are associated with Ω and N h , respectively. As shown below, this analysis clearly separates the shape pattern for mode-III unfolding from the other behaviors observed under unfolding conditions. Figure 8 illustrates a typical trend in the (N h ,Ω) map. The black line corresponds to the mode-III partial unfolding, whereas the gray line corresponds to a mode-II deformation. The figure also indicates the position of the native state, the final snapshots for the corresponding trajectories, and an intermediate conformer obtained during the mode-III trajectory (snapshot on the righthand side, with N h ≈ 50). From these two trajectories, the following observations can be made: (i) The mode-II relaxation of the native state (gray line) is initiated by an increase in entanglement at nearly constant anisometry. After a persistent transient at N h ≈ 49 and Ω ≈ 0.07, the protein undergoes a partial-unfolding transition characterized by a high correlation between anisometry and entanglement. Along this trajectory, there is a substantial population of intermediates. Despite the difference in the extent of the denaturation, these two features resemble the pattern of molecular shape changes observed during the initial unfolding steps of lysozyme.20 (ii) The mode-III relaxation (black line) is initiated by a similar deformation as the one seen in mode II (leading to transients with N h ≈ 49 and Ω ≈ 0.09). Later, this trajectory undergoes a fast oscillation in asphericity and entanglement. Note that, prior to unfolding, the trajectory generates a transient the molecular shape of which is similar to that of the native state in the crystal. Further along the trajectory, unfolding proceeds with a simultaneous increase in Ω and decrease in N h. The trajectory follows a pattern that resembles the one for mode II, although shifted toward lower values of chain entanglement. There are also two recognizable populations of intermediates along the path in shape space. A continuation of the latter trajectory indicates that the final partly unfolded state can survive up to 10 ns. Nevertheless, as discussed in the next section, it is possible that other mode-III trajectories may show further unfolding. Figure 9 collects the data for all trajectories. Trajectories are indicated in gray, except the mode-III trajectory in Figure 8,
1086 J. Phys. Chem. B, Vol. 106, No. 5, 2002
Arteca and Tapia Relaxation from Partly Unfolded Transients of Cytochrome c′
Figure 9. Anisometry-and-entanglement map for all trajectories under denaturing conditions. The region in gray corresponds to mode-I and mode-II deformations. The black line corresponds to the mode-III run depicted in Figure 8. Structures above the SAW line are noncompact (see text).
which is again highlighted as a black line. The white circle stands for the crystalline native state. The diagram also includes the values for the shape descriptors, N h and Ω, averaged over the accessible configurations of self-avoiding walks (SAWs) with variable excluded-volume interaction and the same backbone length as lysozyme and cytochrome c′.16 These values, indicated by the dotted “SAW line”,16,39 provide a convenient reference to establish whether conformations are “folded” or “unfolded”. Given that the SAWs include only repulsive interactions, conformations below the dotted line are more compact than random walks and thus dominated by attractive interactions. Similarly, conformations above the line are noncompact (likely unfolded) conformers. Figure 9 shows a clear difference between mode-III unfolding and all other processes. Although mode-I and mode-II deformations can produce elongated structures above the SAW line, they span a region of shape space that is distinct from the one covered by the snapshots in Figures 4 and 5. Only the early transients of a mode-III unfolding path are in the shape region corresponding to standard deformations. The unfolding pattern of cytochrome c′ on shape space is akin to that of lysozyme.20 Both proteins show correlated changes in chain compactness and entanglement. To some extent, the response of cytochrome c′ can be compared with that of lysozyme under weak couplings to the thermal bath or near the critical charge for unfolding.39 However, the native fold of cytochrome c′ is more resilient; a strong coupling to the thermal bath does not elicit full unfolding. The pattern for modeIII unfolding of cytochrome c′ and the unfolding of lysozyme also fit with results in other proteins. Simulations on the engrailed homeodomain protein show that unfolding paths in water are dominated by an organized loss of tertiary contacts with a smaller distortion on local secondary structure.24 Despite these differences, the pattern in Figure 9 suggests the occurrence of a well-defined process by which the cytochrome c′ responds to unfolding conditions in vacuo. The ability to undergo these structural rearrangements would therefore appear to be intrinsic to the polymer and would not require the presence of a guiding solvent.
The structures observed in the previous sections could be regarded as typical of those accessible to the native state of cytochrome c′ during the early stages of in vacuo denaturation. The properties of these conformers can be tested by subjecting them to relaxation conditions that allow for refolding or further denaturation. Previous work on lysozyme indicates qualitatively distinct populations of unfolded transients.17 Conformers that are unfolded below a “critical” (yet nontrivial) extent are able to refold reVersibly to natiVelike structures. In contrast, unfolded transients above the critical value relax to long-liVed compact intermediates with nonnatiVe structural features. Here, we discuss whether a similar process is at play for the unfolding of cytochrome c′. Relaxation conditions are handled as done previously for lysozyme.15-18 To this end, an unfolded transient is chosen as the initial structure to restart a simulation. Global rotation and translation are stopped, the initial velocities are rerandomized at temperature T0 ) 293 K, and finally the protein is recoupled to the thermal bath while reimposing in vacuo boundary conditions of unrestrained rotations and translations. We have considered two cases. On one hand, a weak coupling to the bath (τ ) 100 ps) that should favor relaxation to compact structures is considered. On the other hand, a strong coupling (τ ) 5 ps), equal to the one used to produce the initial unfolded transient, is used. Consistent with the quasi-native character of the transients produced by mode-I and mode-II deformations, these structures relax toward the native state under conditions of weak bath coupling. In other words, the cytochrome c′ structures produced by the majority of the MD unfolding runs belong to the catchment region of the native state. For this reason, we focus on the relaxation behavior of the more interesting mode-III transients. Using the last conformational snapshot obtained during the mode-III run with partial unfolding (cf. Figure 8), we have produced an ensemble of trajectories with both weak and strong coupling to the bath. For better comparison, we used the same ensembles of initial velocities employed during the initial unfolding trajectories. The general trends are illustrated by the results in Figure 10. Here, we display two trajectories beginning from the initial partly unfolded transient of the mode-III run (white circle). The following observations can be made: (i) The top trajectory indicates the evolution of the mean overcrossing number in the case of weak coupling. In this case, relaxation takes the initial unfolded transient to values consistent with the metastable form of the in vacuo native state of cytochrome c′. The figure displays the final conformer obtained during the new 2-ns-long trajectory. The resemblance with the native state is not only visually striking; the rmsd between this structure and the crystalline native state is only 4.03 Å. The result is more remarkable when remembering that the initial conformer has a completely different helical packing and an overall rmsd with the native state of 10.3 Å. Although generic potential energy functions have the ability to optimize a given fold for many sequences, near-native refolding is not trivial unless the initial relaxing conformer is in the close vicinity of the native-state attractor. Given the strong differences between initial structure and the native fold in this case, we do not believe that the result is an artifact of the force field. (ii) The bottom trajectory in Figure 10 shows the change in N h when using a strong coupling. In this case, relaxation does not lead to quasi-native conformers. For clarity, we have
Protein Denaturation in Vacuo
Figure 10. Evolution of chain entanglements during the relaxation dynamics of an initial partly unfolded transient for a mode-III trajectory. The partly unfolded initial structure is indicated by the white circle. The black circles denote the position for the selected snapshots. The results show refolding to nativelike structures when using a weak coupling relaxation (top) and reversible re-unfolding to the initial transient when using a strong coupling (bottom).
Figure 11. Relaxation pattern for an ensemble of independent trajectories initiated from the transient in Figure 10. Weak coupling relaxations are in gray, whereas strong coupling relaxations are in black. The symbols to the right indicate (a) relaxation to quasi-native structures, corresponding to mode-I final transients, (b) relaxation to slightly deformed nativelike structures, corresponding to mode-II transients, (c) reversible re-unfolding to the partly unfolded initial structure, corresponding to a mode-III transient, and (d) relaxation leading to full unfolding.
indicated a number of selected snapshots (black circles). As shown, re-unfolding sets in at t > 1000 ps, after the formation of structures with intermediate compactness. The final snapshot displayed, corresponding to t ) 2000 ps, has an rmsd of 10.7 Å with respect to the native structure but an rmsd of 2.42 Å with respect to the initial transient from which the trajectory was started. In summary, the results in Figure 10 show clearly the occurrence of two distinct outcomes: (i) refolding to nativelike structures when using a weak coupling and (ii) reVersible reunfolding toward mode-III partly unfolded transients when using a strong coupling. The pattern over the entire ensemble of trajectories appears in Figure 11. Relaxation runs produced with a weak coupling appear in gray, whereas those with a strong coupling are in black. Relaxation with a weak coupling leads
J. Phys. Chem. B, Vol. 106, No. 5, 2002 1087 either to native structures (denoted by (a)) or quasi-native structures (denoted by (b)). These outcomes match closely the final structures obtained during mode-I and mode-II deformations, respectively. In contrast, the relaxations under strong coupling produce re-unfolding to the initial mode-III transient (denoted by (c)) or further unfolding (denoted by (d)). In the latter case, there is a very sharp transition at t ≈ 800 ps, followed by a sequence of several steps with extensive unfolding yet still retaining some R-helical structure. Figure 12 shows a selection of snapshots to illustrate the nature of the latter process. Up to 700 ps, the snapshots are typical of re-unfolding to the initial mode-III transient. However, a fast transition takes place between 800 and 900 ps, whereby a pair of R-helices moves away perpendicularly to the other pair. This result is triggered by the occurrence of a sustained rotation about a direction perpendicular to the helical axes. Finally, very elongated conformers appear for t g 1000 ps. Even then, we find a residual secondary structure, including the local structure for a pair of helices. Most remarkably, the structure develops some β-sheet content during the simulation that is not present in the native structure (cf. the snapshots for 900, 1000, and 2000 ps in Figure 12). A similar phenomenon of spontaneous β-sheet organization was observed during the relaxation dynamics of a model of lysozyme with reduced attraction.19 The sequence of events in Figure 12 resembles qualitatively, in a reverse order, a hierarchical mechanism of folding. In this case, the organization of local secondary structure and tertiary contacts would occur in a different time scale. Figure 13 completes our discussion with the pattern of molecular shape transitions in the anisometry-and-entanglement space. Starting from the partly unfolded transient (white circle), all trajectories exhibit the same initial response of increasing entanglements and chain compactness. In the case of weak coupling (gray area), the resulting conformers span a region similar to the one observed in Figure 9 for mode-I and mode-II deformations. In the case of strong coupling (black lines), the region is expanded by the occurrence of both partial and full unfolding. The latter fully unfolded structures obtained when relaxing under strong denaturation conditions resemble qualitatively those observed during the full unfolding of lysozyme.14,20 In both cases, the final conformers resemble rods with ca. 120 Å span with residual secondary structure at either one or both ends. Cytochrome c′ did not reach such structures within the initial ensemble of 50 MD unfolding runs with integration step ∆t ) 1 fs. A similar pattern did appear within an additional ensemble of 35 MD runs with integration step of ∆t ) 2 fs. For this reason, we believe that extensive unfolding of cytochrome c′ starting from the native structure cannot be ruled out under strong unfolding conditions. However, it may not be a very likely outcome during the initial step of unfolding as modeled here. Given the diversity of relaxation behaviors accessible to partly unfolded (mode-III) transient, these structures would appear to be located in a “bottleneck” region of the folding landscape, connecting the native state attractor with the unfolded continuum. Conclusions The results in this work provide a first exploration of the pattern of in vacuo unfolding accessible to a “sturdy” native state, the four-R-helix bundle. As discussed before, our approach can be regarded as a tool to introduce an unfolding bias. We
1088 J. Phys. Chem. B, Vol. 106, No. 5, 2002
Arteca and Tapia
Figure 12. Selected snapshot showing the transition toward full unfolding, starting from the mode-III partly unfolded transient. The final structures are highly elongated yet retain a degree of secondary structure. (The label “N” denotes the location of the N-terminus.)
Figure 13. Overall pattern of relaxations in the anisometry-andentanglement space. Trajectories with weak coupling to the thermal bath are in gray; those with a strong coupling are in black. The structures highlighted as “partial re-unfolding” correspond to the return toward mode-III transients. The trajectories leading to full unfolding occupy regions observed during unfolding simulations of lysozyme.17,20
do not know whether the weak “centrifugal” forces that accompany global spinning under a strong coupling to the thermal bath might correspond to realistic motions of a protein diffusing in a low-pressure gas (or protein ion under a weak electric field). However, even if this were not the case, the important notion emerging from our work is that the simulated effects of centrifugal unfolding allow one to “prepare” consistently a molecule toward a state in which it then can develop the intrinsic large-scale motions that can trigger substantial unfolding. In the case of cytochrome c′, direct full unfolding appears to be a rare event. Under the same conditions, the lysozyme fold produces a sharp transition into fully unfolded conformers. However, we have shown here that nontrivial partial unfolding is indeed possible in cytochrome c′, and the resulting structures appear at the threshold of the unfolded continuum. Our results indicate that partly unfolded cytochrome c′ (corresponding to mode-III trajectories) can relax to the native state or re-unfold reversibly or undergo full unfolding depending on the strength of the denaturing conditions. No other intermediate structure found shares these possibilities.
The present analysis suggests therefore that fast spinning does not automatically elicit unfolding. Rotational activation of an unfolding mode may depend on protein tertiary structure. In the case of lysozyme, a protein with a multidomain structure and a slightly ellipsoidal native state, spinning produces fast unfolding. In the case of cytochrome c′, a protein with the same chain length but with strong tertiary contacts between helices, the effect is less dramatic. Nevertheless, the accessible deformations detected in both cases share some common features. The pattern emerging points towards an initial step of unfolding dominated by a correlated change in entanglements (or tertiary contacts) and chain compactness with little effect on local secondary structure. These characteristics would appear to distinguish the unfolding paths from relaxation trajectories that start from the fully unfolded continuum both in solution and in vacuo.21,22,24 In the latter case, the initial step is mostly dominated by a rapid increase of compactness, with little development of nativelike chain entanglements.17,18 Acknowledgment. We thank Prof. Curt Reimann (Lund) for valuable discussions and Koradi et al. for use of their protein visualization program Molmol.40 G.A.A. acknowledges support from NSERC (Canada), and O.T. is grateful to Vetenskaprådet (Council for the Sciences, Sweden) for financial support. References and Notes (1) Shelimov, K.; Jarrold, M. F J. Am. Chem. Soc. 1996, 118, 10313. (2) Gross, D.S.; Schnier, P. D.; Rodrı´guez-Cruz, S. E.; Fagerquist, C. K.; Williams, E. K. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 3143. (3) (a) Shelimov, K. B.; Clemmer, D. E.; Hudgins, R. R.; Jarrold, M. F.; J. Am. Chem. Soc. 1997, 119, 2240. (b) Valentine, J. S.; Anderson, J. G.; Ellington, A. D.; Clemmer, D. E. J. Phys. Chem. B 1997, 101, 3891. (c) McLafferty, F. W.; Guan, Z.; Haupts, U.; Wood, T. D.; Kelleher, N. L. J. Am. Chem. Soc. 1998, 120, 4732. (4) Jarrold, M. F. Acc. Chem. Res. 1999, 32, 360. (5) Hoaglund-Hyzer, C. S.; Counterman, A. E.; Clemmer, D. E. Chem. ReV. 1999, 99, 3027. (6) Jarrold, M. F. Annu. ReV. Phys. Chem. 2000, 51, 170. (7) Marchi, M.; Ballone, P. J. Chem. Phys. 1999, 110, 3697. (8) Ferrara, P.; Apostolakis, J.; Caflisch, A. J. Phys. Chem. B 2000, 104, 4511. (9) Isralewitz, B.; Gao, M.; Schulten, K. Curr. Opin. Struct. Biol. 2001, 11, 224. (10) Reimann, C. T.; Vela´zquez, I.; Tapia, O. J. Phys. Chem. B 1998, 102, 9344. (11) Alonso, D.; Daggett, V. Protein Sci. 1998, 7, 860.
Protein Denaturation in Vacuo (12) Kazmirski, S. L.; Daggett, V. J. Mol. Biol. 1998, 277, 487. (13) Bruscolini, P.; Casetti, L. Phys. ReV. E 2000, 61, R2208. (14) (a) Reimann, C. T.; Vela´zquez, I.; Tapia, O. J. Phys. Chem. B 1998, 102, 2277. (b) Arteca, G. A.; Tapia, O. J. Chem. Phys. 2001, 115, 10557. (15) Vela´zquez, I.; Reimann, C. T.; Tapia, O. J. Am. Chem. Soc. 1999, 121, 11468. (16) Arteca, G. A.; Vela´zquez, I.; Reimann, C. T.; Tapia, O. Phys. ReV. E 1999, 59, 5981. (17) Arteca, G. A.; Vela´zquez, I.; Reimann, C. T.; Tapia, O. J. Chem. Phys. 1999, 111, 4774. (18) Arteca, G. A.; Paulino, M.; Reimann, C. T.; Tapia, O. Phys. Chem. Chem. Phys. 2000, 2, 5259. (19) Arteca, G. A.; Reimann, C. T.; Tapia, O., J. Phys. Chem. B 2000, 104, 11360. (20) Arteca, G. A.; Tapia, O. J. Mol. Graphics Modell. 2001, 19, 102. (21) Finkelstein, A. V. Protein Eng. 1997, 10, 843. (22) Dinner, A. R.; Karplus, M. J. Mol. Biol. 1999, 292, 403. (23) Reimann, C. T.; Sullivan, P. A.; Axelsson, J.; Quist, A. P.; Altman, S.; Roepstorff, P.; Vela´zquez, I.; Tapia, O. J. Am. Chem. Soc. 1998, 120, 7608. (24) Mayor, U.; Johnson, C. M.; Daggett, V.; Fersht, A. R. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 13518. (25) (a) Tahirov, T. H.; Misaki, S.; Meyer, T. E.; Cusanovich, M. A.; Higuchi, Y.; Yasuoka, N. J. Mol. Biol. 1996, 259, 467. (b) Scott, R. A., Mauk, A. G., Eds. Cytochrome c: A Multidisciplinary Approach; University Science Books: Sausalito, CA, 1996.
J. Phys. Chem. B, Vol. 106, No. 5, 2002 1089 (26) Arteca, G. A.; Tapia, O. Int. J. Quantum Chem. 2000, 80, 848. (27) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular Simulation (GROMOS) Library Manual; Biomos: Groningen, Holland, 1987. (28) Åqvist, J.; van Gunsteren, W. F.; Leijonmarck, M.; Tapia, O. J. Mol. Biol. 1985, 183, 461. (29) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (30) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1989. (31) Arteca, G. A.; Mezey, P. G. Biopolymers 1992, 32, 1609. (32) Janse van Rensburg, E. J.; Sumners, D. W.; Wasserman, E.; Whittington, S. G. J. Phys. A 1992, 25, 6557. (33) Arteca, G. A. Biopolymers 1993, 33, 1829. (34) Arteca, G. A.; Nilsson, O.; Tapia, O. J. Mol. Graphics 1993, 11, 193. (35) Katritch, V.; Bednar, J.; Michoud, D.; Scharein, R. G.; Dubochet, J.; Stasiak, A. Nature 1996, 384, 142. (36) Vologodskii, A. V.; Crisona, N. J.; Laurie, B.; Pieranski, P.; Katritch, V.; Dubochet, J.; Stasiak, A. J. Mol. Biol. 1998, 278, 1. (37) Arteca, G. A. J. Chem. Inf. Comput. Sci. 1999, 39, 550. (38) Baumga¨rtner, A. J. Chem. Phys. 1993, 98, 7496. (39) Arteca, G. A.; Vela´zquez, I.; Reimann, C. T.; Tapia, O. Chem. Phys. Lett. 2000, 327, 245. (40) Koradi, R.; Billeter, M.; Wu¨thrich, K. J. Mol. Graphics 1996, 14, 51.
