This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Rousseau, R. Articles by Marx, D. Search for Related Content PubMed PubMed Citation Articles by Rousseau, R. Articles by Marx, D.

Temperature-Dependent Conformational Transitions and Hydrogen-Bond Dynamics of the Elastin-Like Octapeptide GVG(VPGVG): A Molecular-Dynamics Study

Roger Rousseau ^* , Eduard Schreiner ^*, Axel Kohlmeyer ^* and Dominik Marx ^*

^* Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, Bochum, Germany and International School for Advanced Studies, Trieste, Italy

Correspondence: Address reprint requests to Eduard Schreiner, Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, Universitätsstraße 150, 44801 Bochum, Germany. Tel.: 49-234-322-2121; E-mail: eduard.schreiner{at}theochem.ruhr-uni-bochum.de.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS, MODELS, AND TECHNICAL... RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
A joint experimental/theoretical investigation of the elastin-like octapeptide GVG(VPGVG) was carried out. In this article a comprehensive molecular-dynamics study of the temperature-dependent folding and unfolding of the octapeptide is presented. The current study, as well as its experimental counterpart (see companion article in this issue) find that this peptide undergoes an inverse temperature transition (ITT), leading to a folding at 40–60°C. In addition, an unfolding transition is identified at unusually high temperatures approaching the normal boiling point of water. Due to the small size of the system, two broad temperature regimes are found: the ITT regime at 10–60°C and the unfolding regime at T > 60°C, where the peptide has a maximum probability of being folded at T 60°C. A detailed molecular picture involving a thermodynamic order parameter, or reaction coordinate, for this process is presented along with a time-correlation function analysis of the hydrogen-bond dynamics within the peptide as well as between the peptide and solvating water molecules. Correlation with experimental evidence and ramifications on the properties of elastin are discussed.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS, MODELS, AND TECHNICAL... RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Vertebrate elastic fibers, as contained in vascular walls, skin, or lungs, allow for reversible deformations upon mechanical stress. These fibers consist of two main types of protein components: a fibrous component, fibrilen; and an amorphous component, the globular elastin protein (Pasquali-Ronchetti et al., 1995). The essential elasticity is provided by the latter protein, elastin, which is known to have very unique viscoelastic properties in the water-swollen state; see Debelle and Alix (1999), Debelle and Tamburro (1999), Martino et al. (2002), Reiersen and Rees (2001), Urry (1997), and Urry et al. (2002) for recent overviews. The insoluble elastin is an extensively cross-linked polymer of a precursor tropoelastin. Tropoelastin is composed of two types of domains. One of them is rich in lysine and provides cross-linking of the individual monomers resulting in lysinonorleucin, desmosine, and isodesmosine links, which are responsible for its typical yellow color. The other is made predominantly of hydrophobic amino acids with the highly repetitive pentameric repeat unit, VPGVG (amino acids valine, V; proline, P; and glycine, G) which has been found to be crucial to elastin's functionality (Reiersen et al., 1998; Urry, 1993, 1997; Urry et al., 2002).

An interesting peculiarity of both elastin and tropoelastin, is that they undergo folding by increasing the temperature beyond typically 25°C. The term inverse temperature transition (ITT) was coined for this apparently paradoxical change from a disordered (extended) to an ordered (folded) conformation upon heating. A striking manifestation of this phenomena is the ability to grow crystalline solid-state structures from a solution of cyclic elastin-like oligopeptides by heating (Cook et al., 1980; Urry et al., 1978). As such, tropoelastin and its synthetic analogs have been a subject of intense investigation in the field of biopolymers and protein engineering. Elastin-like polymers have the promise of providing materials for biomechanical devices (Nath and Chilkoti, 2001; Urry, 1993, 1997) as well as temperature-dependent molecular switches (Reiersen et al., 1998). Despite the intense research efforts, the detailed structure of tropoelastin still defeats elucidation. Hence a molecular-level picture which relates the protein structure to the viscoelastic properties is still a matter of debate. What is beyond doubt, however, is the decisive role of water as a plasticizer; dry elastin is brittle (Partridge, 1962). Still, the potentially crucial aspect of the protein's hydration water remains largely unexplored with the notable exception of recent molecular-dynamics (MD) simulations (Li et al., 2001b).

The origin of the viscoelastic properties of elastin is controversially discussed, recalling concepts such as classical rubber elasticity (Hoeve and Flory, 1974), various librational entropy mechanism (Debelle and Tamburro, 1999; Urry, 1993; Wasserman and Salemme, 1990), and multiphase models (Debelle and Alix, 1999; Li et al., 2001b). Much of the experimental data has been summarized in excellent reviews (Li and Daggett, 2002; Martino et al., 2002; Urry, 1997) on both elastin and elastin-like polypeptides. Much data exist that suggest that both elastin and its synthetic mimics display an interesting conformational dynamics in solution. Early nuclear magnetic resonance studies suggest that elastin under physiological conditions is composed of highly mobile chains (Fleming et al., 1980; Torchia and Piez, 1973), which correlated well with the observation from bi-refringence (Aaron and Gosline, 1980) measurements that elastin, in the same state, is isotropically distributed with the chains adopting random conformations. Studies of the thermomechanical properties of water-swollen elastin and tropoelastin are consistent with the interpretation of elastin as a classical rubber with heavy emphasis on the role of entropic contributions to elasticity (Hoeve and Flory, 1974). These findings would support a single-phase model where the structure of elastin is random and undergoes large amplitude fluctuations. Recent single-molecule atomic force microscopy and spectroscopic measurements (Urry et al., 2002) studies have severely challenged this interpretation and suggest that, above the ITT temperature, a structurally ordered model of poly(VPGVG) can also account for these elastic properties due to significant protein librational contributions to the overall entropy as suggested two decades ago (Urry et al.,1983). Other spectroscopic experiments such as Fourier transform infrared (FT-IR), nuclear magnetic resonance (NMR), circular dichroism (CD), and Raman measurements (Debelle et al., 1995; Nicolini et al., 2004; Arad, 1990; Reiersen et al., 1998; E. Schreiner, C. Nicolini, B. Ludolph, R. Ravindra, N. Otte, A. Kohlmeyer, R. Rousseau, R. Winter, and D. Marx, unpublished material; Urry et al., 1985a) on elastin-like polypeptides suggest the presence of ß- and -turns as a common structural motif which may be in dynamic equilibrium and very short and/or distorted antiparallel ß-strands and disordered structures with possible dynamic sheet-coil-turn transitions. Recent NMR studies also point to the high mobility of the elastin chains in the water-swollen elastins (Perry et al., 2002); and in solution, multiple temperature states for poly(GVGVP), each with different dynamical behavior and peptide-water interactions, were observed (Kurková et al., 2003). These later data would instead support a multiphase model of elastin (Debelle and Alix, 1999), where the fine details of the protein structure and, in particular, its interactions with water, play a crucial role in understanding its elasticity. Overall, there is a wealth of data on the structure and properties of these systems, and although the models do not agree on their interpretation, a central ingredient is the intriguing and complex dynamical behavior of elastin-like polymers themselves and the crucial interplay with solvating water molecules.

Atomistic simulations on elastin or elastin-analogs are scarce. For early work, see Debelle and Tamburro (1999). Wasserman and Salemme (1990), in a review, studied the elastic properties of these proteins by performing short MD simulations under constraints to mimic an external pulling force. Here the authors found that entropic forces from reduced librational motion of the protein were indeed a main contributing factor in the elastic behavior. However, the authors were understandably limited to simulation times not exceeding 1 ns. It is noted that, in general, short simulation times generally lead to unreliable statistics, due to insufficient sampling of configuration space, and hence may provide only qualitative insights. A crucial factor also to be considered is the explicit atomistic treatment of the protein/water interface and characterization of its role in determining elastin's properties. As a significant advance along these lines, recent studies on solvated 90-amino-acid polypeptides of sequence (VPGVG)₁₈ have been reported, which explicitly include water (Li et al., 2001b). The importance of the protein/water interface was worked out based on thermodynamic and static structural—as opposed to dynamic—considerations. This work was successfully able to reproduce the ITT, clarifying how water may contribute to the entropic contribution to the elastic forces (Li et al., 2001a) as well as probe the role of elastin mimics as potential molecular machines (Li et al., 2002). Overall this work was successfully able to qualitatively relate its findings to much of the experimentally observed phenomena, although some debate has arisen with regard to the ability of the proposed structures to account for the temperature-dependence of elastin's water content above the ITT (Urry et al., 2002). Despite the significant improvement that these simulations have brought to our understanding of the atomic-level details concerning elastins' behavior, these quite-long biopolymers were simulated on the order of a few nanoseconds for a given temperature. In view of the expected long relaxation times of proteins of that size (Urry et al., 1985b), it would be highly desirable to have access to simulation times that allow for quantitative statistical-mechanical analysis for a given number of amino acids. This is especially true for low-frequency oscillations of the protein backbone, which play a central role in many models of elastin's properties. Furthermore, the dynamics and, in particular, the kinetics of hydrogen bonds, remains unexplored up to now. Thus, questions regarding the coupling of the dynamics of the protein, the protein/water interface, and hydrogen bonding remain open issues.

A significant experimental finding was the recent demonstration (Reiersen et al.,1998) that oligopeptides of the kind GVG(VPGVG)_n display the ITT even in the limit of only one pentameric repeat unit. Using the limiting value n = 1 opens up the possibility to perform MD simulations that allow for fairly long simulation times with sufficiently many solvating water molecules. We therefore have launched a joint experimental/simulation study of the smallest possible elastin model, the octamer GVG(VPGVG) (having only 640 Da); see the companion article, Nicolini et al. (2004). This minimal elastin model was chosen, at the expense of experimental difficulties, to allow for extensive MD simulations, at 12 temperatures between 280 and 390 K, with a thorough conformational sampling of 32 ns each (after equilibration), followed by in-depth statistical mechanical analysis with a focus on dynamics. The current approach will allow us to obtain a crucial link between the multitude of experimental results and the atomistic simulations by directly comparing results obtained on essentially the same system under thermodynamically consistent conditions. Moreover, the detailed statistical mechanical analysis which is possible on the current MD simulations will provide a complementary perspective on simulation work already existing in the literature.

	METHODS, MODELS, AND TECHNICAL DETAILS

TOP ABSTRACT INTRODUCTION METHODS, MODELS, AND TECHNICAL... RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
System
The system consists of an 8-amino-acid oligopeptide of sequence GVG(VPGVG) capped by methylamine (-NH-CH₃) and acetyl (-CO-CH₃) groups at its C- and N-termini, respectively. In accordance with the experimental finding (Reiersen et al., 1998) that VPGVG is the minimal repeat unit necessary for observing the ITT irrespective of the end groups, the caps are not expected to change the basic scenario as found in our complementary experimental study on the zwitterionic species (Nicolini et al., 2004; E. Schreiner, C. Nicolini, B. Ludolph, R. Ravindra, N. Otte, A. Kohlmeyer, R. Rousseau, R. Winter, and D. Marx, unpublished material; further checks are discussed below). The octapeptide is solvated with 2127 waters and its center of mass is fixed to the center of a 50 Å-diameter spherical droplet surrounded by stochastic boundary conditions (Brooks and Karplus, 1983). This sphere diameter allows us to maintain at least 3–4 molecules between the peptide and the stochastic boundary even for the configurations where elastin is at its maximum elongation of 27 Å. We employ the all-atom CHARMM force field (MacKerell et al., 1998) for the peptide and the TIP3P (Jorgensen et al., 1983) water potential.

Simulations
Within the EGO molecular-dynamics package (Eichinger et al., 1997, 2000) a weak Berendsen thermostating scheme (coupling time constant of 100 fs) was used to control the temperature using a 1 fs MD time step and the SHAKE algorithm for fixing the bonds between hydrogen and heavy atoms. This conservative choice of parameters for a careful integration of the equations of motion has allowed us to perform stable simulations which included equilibration times of 2–5 ns (depending on temperature) before obtaining 32 ns trajectories at 12 temperatures between 280 and 390 K. An additional run was performed at 400 K (not shown) to stretch the temperature range as much as possible beyond the normal boiling point of water. For many properties, such as hydrogen-bond breaking and reformation-rate constants, the same trend as that obtained from extrapolation of the data from 340 to 390 K was observed. However, average structural parameters in particular showed deviations from this trend, which we interpret tentatively as failures of the simulation approach. Note that 400 K exceeds by 100 K the temperature range for which the employed force fields for both water and the peptide are optimized, i.e., ambient conditions. As a convergence check, the run obtained at 280 K was extended to >40 ns to test the reliability of both average and dynamical properties of the protein; the results at 280 K are quoted for the 32 ns MD run for statistical consistency. Internal pressures inside the droplet were found to be <0.1–0.2 kbar depending on the temperature, which is well below the pressure of the order-of-many-kbars required to suppress protein folding entirely (Nicolini et al., 2004; Tamura et al., 2000).

To cross-check the possible influence of the capping groups we have computed two 25 ns trajectories at 280 and 320 K in larger droplets (60 Å diameter with 3700 water molecules) for both the capped and zwitterionic peptides. We obtained qualitatively identical results to those presented here, indicating that the ITT is neither significantly affected by the presence of the caps (which is in line with experimental findings; see Reiersen et al., 1998), nor by the spherical boundary conditions for the chosen size of the solvating droplet.

Hydrogen-bond analysis
Concerning the definition of hydrogen bonds (HBs) in water, and thus the population variable, h(t), we employ a standard structural criterion based on both the distance (R_OH < 2.6 Å as obtained from including the entire first nearest-neighbor peak in the intermolecular OH radial distribution function) and the angle, _O–HO > 130° for h(t) = 1 (Luzar, 2000; Luzar and Chandler, 1996a,b; Starr et al., 1999; Stillinger, 1975; Xu and Berne, 2001). As a check we also performed the analysis with only the distance criterion, and yet again with the distance criterion, but using a tighter angular criterion of _O–HO > 150°. Consistent with previous findings (Starr et al., 1999; Xu and Berne, 2001) the absolute value of rate constants, k_w and , is sensitive to the actual choice of the parameters of the definition and was found to vary as much as a factor of 2–3 between different definitions. However, we always obtained the same qualitative trends that both HB numbers and rate constants (k_w and ) feature a linear behavior with respect to temperature. Similarly, for peptide/peptide and peptide/water HBs, the same criterion based on distances (R_AH < 2.6 Å where the proton acceptor, A, is either an oxygen atom of a water molecule or of a peptide carbonyl group) and angles (_D–HA > 130° with the proton donor, D, being either a water O-atom or peptide backbone N-atom) was used. As a check we varied the distance to R_AH < 2.4 Å, keeping the angle cutoff unaltered, and performed the analysis with only the distance criterion. Again we found the same general behavior of two temperature regimes based on the dynamics of peptide/water HBs and the temperature-dependence of n_pp (peptide-peptide) and n_bw (bridging water, i.e., peptide-water-peptide) to be preserved. In conclusion, although the HB definitions do change the reported absolute values of the HB numbers and rate constants as expected (Starr et al., 1999; Xu and Berne, 2001), the qualitative behavior remains the same, and thus the trends in the temperature-dependence for these quantities are not affected.

Principal component analysis
Principal component analysis (PCA) is a standard projection method in statistical data analysis, feature extraction, and data compression (Hyvärinen et al., 2001). It also turned out to be a very useful tool for analyzing the complex conformational behavior of biomolecules as obtained from molecular-dynamics trajectories (see Amadei et al., 1993, Balsera et al., 1996, and García, 1992, for concepts, details, and limitations). For the PCA carried out here, all heavy atoms of the peptide backbone were included. In addition, a united-atom representation was introduced for each side chain by assigning a mass equivalent to that of all its constituting atoms, which was located at its center of mass, reducing the number of sites to N = 38. Explicit tests indicate that PCA, including all atoms, provides an essentially identical picture for the low-order modes that are of interest, which is in line with previous findings (Aalten et al., 1997). The dependence of the resulting motions upon trajectory length was also carefully evaluated and results were found to be well converged for simulations on the order of 20 ns in duration. For example, the dot product between the slowest motion, i.e., {m₁}, computed over two distinct 20 ns trajectories, and that from a 30- to 40 ns trajectory at the lowest temperature, 280 K, is 0.98—therefore indicating that the modes are essentially converged. Moreover, the dot products between the vectors {m₁} obtained from PCA at all temperatures always exceed 0.98, indicating that this motion, as such, is not temperature-dependent. Finally, the modes discussed in this investigation are found to exist also when using the zwitterionic octapeptide, and in the simulations with a larger solvation sphere around the peptide, as discussed above.

As an alternate approach to using Cartesian coordinates, we have also applied PCA in internal coordinates (see e.g., Aalten et al., 1997). Natural internal coordinates for protein backbone conformational changes are the usual dihedral angles and , defining rotations around the C-C_carbonyl and C-N_amide bonds of each amino acid, respectively. Concerning the octapeptide, we confirm qualitatively that opening/closing and librational modes in the PCA analysis also exist in / space. However, at variance with PCA in Cartesian coordinates (see below for a detailed discussion), there is no single eigenvector that can be used as a reaction coordinate, the reason being that the motions represented by the eigenvectors in dihedral space turn out to be very local in the sense of involving mainly individual angles and thus decoupled rotations. Therefore, a linear combination of many angular eigenvectors must be used to describe sufficiently collective correlated motion. This is also confirmed by the fact that many more such eigenvectors are needed to represent the variance of the trajectory as compared to performing PCA in Cartesian coordinates, where the first mode already accounts for 40% of the total variance.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS, MODELS, AND TECHNICAL... RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Time evolution and temperature-dependence of structural parameters
To get a first impression of the dynamics of the solvated octapeptide GVG(VPGVG), we consider the time evolution of a few selected quantities describing the overall structure. The first parameter of interest is the distance from the C- to the N-terminus, r_C–N, which is a measure of the extension of this short polymer chain. A typical time-evolution profile of this quantity is presented in Fig. 1 a at a temperature of 330 K. Judging from this parameter, there are two distinctly different types of peptide structure: the majority of the time r_C–N oscillates approximately a distance of 17 Å (corresponding to open structures), whereas a second, less probable configuration is found when r_C–N 5 Å (stemming from closed structures). Representative configurations, as obtained at 330 K for each type of conformation, are depicted in Fig. 2. The presence of two types of peptide conformations is also mirrored by the radius of gyration of the biopolymer, r_gyr, which is a standard measure of a polymer's overall size. This parameter is strongly correlated with the end-to-end distance, r_C–N (see Fig. 1 b), such that for the extended structures, r_gyr 8 Å, whereas for the closed ones, r_gyr 4 Å. As a simple measure of the proximity of hydrophobic side chains of the peptide we also consider the radius of gyration of the valine side chains only, (see Fig. 1 c). Again, this parameter is found to be strongly, although not perfectly, correlated with the end-to-end distance and, as expected, the radius of gyration of the entire peptide.

View larger version (43K): [in this window] [in a new window]   FIGURE 1  Time evolution of selected structural quantities of GVG(VPGVG) in water at 330 K. (a) Distance from C- to N-terminus, rC–N; (b) radius of gyration, rgyr; (c) radius of gyration of the valine side chains, ; (d) number of peptide-peptide HBs, npp; (e) number of water molecules bridging two amino acids of peptide by HBs, nbw; (f) number of HBs between peptide and solvating water molecules, nps; and (g) projection of the trajectory onto the first eigenvector of the covariance matrix, ; see text for definitions. For presentation purposes only, the functions were denoised using the Savitzky-Golay procedure (Savitzky and Golay, 1964), where a time window of 128 ps and a polynomial of sixth degree was used; the analysis was, however, carried out based on the original data sets.

View larger version (20K): [in this window] [in a new window]   FIGURE 2  Graphical representation of typical protein conformations at 330 K including only bridging waters; (a) open state and (b) closed state.

As another measure of the solvation environment of the peptide we consider the number of water molecules within the first solvation shell. This quantity may be broken down into two subclasses: 1), the number of waters solvating the hydrophilic backbone atoms, n_bb (backbone); and 2), those that solvate the hydrophobic side chain carbons, n_sc (side chain). The water molecules of the first solvation shell are defined to be within 4.5 Å of the peptide heavy atoms as obtained from the first nearest-neighbor peak of the radial distribution function. For the hydrophobic side chains the first nearest-neighbor peak of the radial distribution function of water oxygen around peptide heavy atoms is within 4.5 Å. For the backbone atoms, the radial distribution function has two maxima in this region, one stemming from the waters associated with the backbone via HB within 3.0 Å, and one belonging to waters which are in the vicinity of the backbone but without an HB at 4.5 Å. Therefore in both cases the reasonable cutoff of 4.5 Å was chosen. For the current small peptide n_bb and n_sc are roughly equal to 50–70 water molecules, depending on the temperature. We plot in Fig. 3 the time evolution of these two quantities along with the gyration radius, r_gyr, at 300 K. At this temperature we again find that r_gyr shows the signature of both open and closed structures, although with pronounced lower occurrence of the latter. As expected, n_bb is strongly correlated with the opening and closing, dropping in value by 20 water molecules upon entering the closed state. This is in stark contrast with the number of waters around the side chains: n_sc is largely insensitive to the opening and closing motion. A similar happening is found to occur at all temperatures, which strongly suggests that rearrangement of the water structure around the hydrophobic side chains does not play a dominant role in the energetics of forming open and closed structures.

View larger version (33K): [in this window] [in a new window]   FIGURE 3  Time evolution of selected structural quantities of GVG(VPGVG) in water at 300 K. (Top) Radius of gyration, rgyr; (middle) number of water molecules in first solvation shell of hydrophilic backbone atoms, nbb; and (bottom) number of water molecules in first solvation shell of hydrophobic side chain atoms, nsc. See Fig. 1 for smoothing procedure.

View larger version (16K): [in this window] [in a new window]   FIGURE 6  Dynamics of the m1 eigenmode. (a) Autocorrelation function, cm(t), at 280 K, ; 330 K, ; 360 K, ; and 390 K, . (b) Relaxation time m as a function of temperature; here the symbol size reflects approximate error bars, . (c) Peptide backbone quasiharmonic entropy, S, , as a function of temperature; see text for definitions.

View larger version (26K): [in this window] [in a new window]   FIGURE 4  Average quantities as a function of temperature. (a) Distance from C- to N-terminus, rC–N, •; (b) radius of gyration, rgyr, ; (c) radius of gyration of valine side chains, , ; (d) number of HBs between peptide and solvating water molecules, nps, ; (e) number of internal HBs, nIHB, ; and (f) number of water molecules solvating a backbone heavy atom, nbb, ; and number of water molecules solvating a hydrophobic side-chain atom, nsc, ; see text for definitions. Dashed lines are linear connections of the data to guide the eye. Solid lines are linear least-squares fit.

Overall, these averaged structural quantities provide a suggestive, though by no means comprehensive, picture of the atomic-level details of the peptide structural transitions. The presence of two types of structures, open and closed, correlates extremely well with the observation of an isodichroitic point in the CD spectra (Nicolini et al., 2004; Reiersen et al., 1998) and supports the description of this small peptide as behaving like a two-state system. The increase of the contribution of closed structures, at temperatures up to 330 K and its subsequent decrease at higher temperatures, is consistent with the temperature-dependence of more compact structures containing - and ß-turns, as observed in the FT-IR spectra of the complementary experimental investigation (Nicolini et al., 2004). The differential scanning calorimetry (DSC) and pressure perturbation calorimetry (PPC) thermodynamic measurements presented in the joint experimental study may also be interpreted as the peptide ultimately reverting to a structurally and thermodynamically similar state at temperatures near the normal boiling point of water. In addition, the above scenario is qualitatively in accord with that presented in Li et al. (2001b), which describes the ITT using a similar but much longer model peptide. Specifically this study also finds an increasing number of peptide/peptide HBs, an overall average contraction of the peptide at the ITT, a decreasing number of interfacial water molecules, and an increasing number of bridging waters which may be thought of as the internal waters in our small elastin-like analog. The authors of Li et al. (2001b) attribute this shrinkage to a hydrophobic collapse, which is in agreement with our observed temperature-dependence of n_ps and .

However, the current picture suggests the subtle, but crucial, difference that the observed structural transitions in capped GVG(VPGVG) are not caused by a classical textbook hydrophobic collapse. Rather, there is a net exchange of peptide/peptide and bridging-water intramolecular HBs for peptide/water interfacial HBs; the energetic implications of the various HB types are addressed in the section on HB dynamics and kinetics. For the octapeptide, our finding and its interpretation has a striking resemblance to the findings concerning self-association of methanol in water in Dixit et al. (2002), in addition to being in qualitative accord with some concepts of hydrophobic solvation as reviewed in Pratt and Pohorille (2002). We must stress, however, that the peptide model in this previous study has a solvation structure where 80% of the solvating water molecules are located around hydrophobic side chains in contrast to the 50% observed in our much smaller system. Thus, it remains an open issue to address these questions in the case of longer GVG (VPGVG)_n peptides. In the following two subsections these suggestive findings, which are based on simple descriptors, will be further scrutinized by a detailed statistical-mechanical analysis.

Principal component analysis and order parameter
To probe the above findings in a systematic but simple way one would ideally wish to describe these phenomena by a single parameter, i.e., by an order parameter or reaction coordinate. We have employed the PCA to explore this possibility (Amadei et al., 1993; Balsera et al., 1996; García, 1992; Hyvärinen et al., 2001) for presentations of this technique. In this approach one computes, by time averaging, ..., over the entire MD trajectory the covariance matrix , where x(t) are the Cartesian coordinates of the peptide atoms at time t in a frame of reference where the overall translations and rotations of the polymer have been subtracted. The 3N–6 normalized eigenvectors, {m_i}, with nonzero eigenvalues, {_i}, provide a basis in which the complex polypeptide motion may be decomposed. This then provides a theoretical framework which can be thought of as being analogous to harmonic normal modes that are of ubiquitous use to explain spectroscopic observations of small molecules. In practice, this allows us to disentangle relevant peptide backbone folding motions from small-amplitude vibrations via projecting the deviations of the MD trajectories from the average structure onto the eigenvectors, i.e., . Ideally, for uncorrelated small-amplitude vibrations the corresponding projections will simply fluctuate around an average value yielding an approximately Gaussian distribution function, ; the width might be temperature-dependent akin to harmonic molecular vibrations. Strong deviations from such a unimodal behavior would single out interesting modes that can be associated to collective changes involving many atoms. Upon changing the temperature the various modes are expected to mix, particularly the higher-order modes. Thus, a common basis, {m_i}, has to be chosen to ensure a consistent analysis. In the following, projections will be performed onto the basis of eigenmodes of the low-temperature trajectory at 280 K. In addition, the mass weighted covariance matrix C' = M^1/2CM^1/2 provides a crude and very approximate method (Andriciolaei and Karplus, 2001; Schlitter, 1993) for estimating the (vibrational/librational) entropy contribution of each mode by assuming a superposition of independent one-dimensional quantum harmonic oscillators, as

(1)

where the (quasiharmonic) frequencies are obtained from the eigenvalues '_i of C'; k_B is the Boltzmann constant.

Upon projecting the MD trajectories onto the eigenvector with the largest eigenvalue, ₁, it is found that this mode, m₁, actually describes the opening and closing motion of the octapeptide, and is present at all temperatures (see Fig. 5 a for a graphical illustration). In Fig. 1 g, we depict this projection, , as a function of time at 330 K where it can be inferred to be perfectly correlated with the other parameters already discussed in the previous section. The octapeptide thus exists in a closed state, , with end-to-end distances of 5 Å; and an extended state, , with end-to-end distances of 17 Å. In addition, the projection process yields for a bimodal distribution function at low temperatures, whereas it adopts a unimodal, but very broad and skewed shape upon heating (see Fig. 5 c).

View larger version (47K): [in this window] [in a new window]   FIGURE 5  Results of the principal component analysis. (a) Visualization of large amplitude opening and closing mode, m1. (b) Visualization of large amplitude librational mode, m2. (c) Distribution function of the projection at 280, 330, 360, and 390 K. (d) Distribution function of the projection at 280, 330, 360, and 390 K. (e) Relative free energy along the projection . (f) Relative free energy along the projection . In a and b, an artificial trajectory using the full first m1 (a), and second m2 (b), eigenvector was synthesized for a suitable graphical presentation. In e–f: 280 K, ; 330 K, ; 380 K, ; and 390 K, . For presentation purposes the free energy profiles for 330, 380, and 390 K were shifted in energy by 5, 10, and 15 kJ mol-1, respectively.

Having obtained a proper order parameter distribution function, , allows one to readily define an effective relative free energy profile (or potential of mean force) according to as a function of temperature. The normalizing reference value is chosen such that the arbitrary value at the minimum of the free energy profile is set to zero for convenience (see Fig. 5 e). At 280 K, the free energy shows two minima, separated by a small barrier, related to both open (broad global minimum at ) and closed (local minimum at ) structures with the latter 3 kJ/mol higher in energy then the former. By 330 K the minimum associated with the closed structure, as well as the barrier, become much lower in relative free energy; thus the entire profile of F broadens out into a flat potential energy landscape. Above 330 K the free energy of the closed minimum gradually increases again and eventually disappears by 370–390 K. Performing an identical analysis for the remaining modes provides approximately harmonic potentials of mean force (see Fig. 5 f for the potential derived from , in agreement with our use of only as an order parameter). The trend in reflects the presence of the temperature regimes observed from the average parameters, i.e., increased folding from 280 to 320 K, the observation of maximum folded structures 330 K, and the unfolding at T 340 K. The presence of two minima in in conjunction with the temperature-induced changes naturally explains the existence of the isosbestic and isodichroitic points observed in FT-IR and CD experiments (Nicolini et al., 2004), respectively, and gives us confidence in the reliability of our simulations.

At 330 K, fluctuations and thus statistical errors are larger than at other temperatures due to the broader configuration space which must be sampled possibly in conjunction with slower relaxation. This effect is similar to the slowing-down problem occurring for statistical convergence at second-order phase transitions where all quantities, and in particular the order parameter itself, undergo large fluctuations. However, the finite system size does not allow for a true phase transition with a well-defined transition temperature but rather, one observes, in agreement with spectroscopic and thermodynamic measurements (Nicolini et al., 2004; Reiersen et al., 1998; E. Schreiner, C. Nicolini, B. Ludolph, R. Ravindra, N. Otte, A. Kohlmeyer, R. Rousseau, R. Winter, and D. Marx, unpublished material) the two temperature regimes which are bracketed around the temperature of maximal folding, T 330 K.

The dynamical motion of the peptide itself can be analyzed by examining the time-dependence of the reaction coordinate as described by the correlation function (see Fig. 6 a). This time-autocorrelation function may be well approximated in the short time limit, t < 100 ps, by an exponential function, c_m(t) exp [-t/_m], with an associated relaxation time, _m, of the peptide's opening and closing motion. The temperature-dependence of this important dynamical parameter is depicted in Fig. 6 b. For temperatures <330 K, _m 300–350 ps is roughly constant within the accuracy of our statistics. Above 330 K, there is a progressive speedup (_m 130–150 ps) which reflects the steepening of the free energy surface in the region of closed structures, i.e., for . From 370–390 K, where the free energy surface contains essentially a single minimum, , with approximately constant curvature, the relaxation time levels off to a plateau value of _m 150 ps. Temperature variations of _m thus reflect three dynamical scenarios corresponding to folding, unfolding, and unfolded, which is slightly different from the observations made on the average structural quantities. Note, however, that the last of the three temperature regimes represents a limiting case of the unfolding transition, where a minimum for closed structures no longer effectively exists.

These findings correlate well with larger peptide backbone fluctuations observed in NMR spectroscopy (Perry et al., 2002) for water-swollen elastins. Unfortunately, the most directly comparable NMR measurements of ¹³C relaxation times for poly(GVGVP) in solution are only obtained for the temperatures below the ITT (Kurková et al., 2003). Nonetheless, the authors report a mean effective correlation time for polymer motion on the order of 500 ps at 298 K in astonishing agreement with the 300-ps timescale reported here at the same temperature. These facts suggest that temperature-dependent measurements of NMR ¹³C relaxation times for this small system may serve as an independent verification of our findings.

Previous MD studies on much larger poly(VPGVG) chains find qualitatively that the peptide becomes "slightly more dynamic," or less rigid at higher temperatures (Li et al., 2001b). This also matches our quantified speedup of the dynamics of the peptide's backbone opening and closing motion. Thus, despite the increase in number of peptide/peptide HBs, there must be an overall increase of freedom in the motion of the peptide backbone. The question then arises: does this extra motion lead to favorable entropic contribution to the ITT from the peptide? To address this, we consider the entropy change of the peptide backbone motion as estimated from a quasiharmonic approximation based on the principal component analysis (Andriciolaei and Karplus, 2001; Schlitter, 1993). Note that this analysis is based on a quasiharmonic approximation (Andriciolaei and Karplus, 2001; Schlitter, 1993), whereas the lowest-order mode, m₁, which contributes 40% to the overall variance of the peptide backbone, is strongly anharmonic and even bimodal at low temperatures. However, its contribution to the total peptide backbone entropy is fairly small, with much larger contributions arising from the more harmonic higher-order modes. Fig. 6 c shows the entropy of the peptide in a united-atom representation as a function of temperature. We find that despite the volume contraction of the peptide <330 K the entropy is increasing with rising temperature, i.e., the entropy of the peptide itself increases upon folding. The associated entropy change is also consistent in magnitude with the experimentally measured entropy increase of 0.15 kJ/mol^-1 K^-1 and 0.11 kJ/mol^-1 K^-1 (at 298 K and 1 bar) obtained from a van't Hoff analysis of equilibrium constants of both CD and FT-IR data (Nicolini et al., 2004; Reiersen et al., 1998; E. Schreiner, C. Nicolini, B. Ludolph, R. Ravindra, N. Otte, A. Kohlmeyer, R. Rousseau, R. Winter, and D. Marx, unpublished material). These findings are similar in spirit to various librational entropy models employed to explain viscoelastic properties of elastin and its synthetic mimetica (Debelle and Tamburro, 1999; Urry, 1993; Urry et al., 2002). Our analysis does not eliminate the possibility that stabilizing entropic terms from the water are also important contributors to the folding behavior; however, our simulations clearly demonstrate that our elastin mimic is a very dynamical entity and that the thermodynamic consequences of this motion, in particular the increasing entropy, are key components in understanding its structural transitions.

Hydrogen-bond dynamics and kinetics
To address the role of water, and in particular, the influence of hydrogen bonding on the peptide's conformational transitions, we analyze the dynamics of the HB network. To this end, we employ the intermittent HB autocorrelation function (Chandra, 2000; Luzar and Chandler, 1996a,b; Starr et al., 1999; Stillinger, 1975; Xu and Berne, 2001), c(t) = h(0)h(t)/h. As usual the HB population variable h(t) is defined to be unity, h(t) = 1, if a particular HB exists at time t and zero h(t) = 0 otherwise (the criteria to define an HB are compiled and discussed in Methods, Models, and Technical Details); and the average ... is taken over all HBs that were present at time t = 0. This function describes the probability that an HB, which was intact at t = 0, is intact at time t independent of possible breakings (and reformation) in the interim time. Similar to the average number of HBs introduced above, the following classes were defined: HBs between water molecules in the bulk c_w(t) (Fig. 7 a); between water molecules which are both located within the first solvation shell of the hydrophobic side chains c_sc(t); between first shell protein solvation water and bulk water c_sw(t) (Fig. 7 b); between the peptide and first shell solvation water c_ps(t) (Fig. 7 c); and direct HB contacts of the peptide with itself c_pp(t) (Fig. 7 d). These functions provide us with a wealth of insights into the complex dynamics and the associated timescales of the solvated peptide in relation to, e.g., relaxation dynamics in bulk water. Comparison of c_w(t) and c_ps(t) nicely shows the well-known trend (Tarek and Tobias, 1999, 2002a,b) of slower dynamics for interfacial peptide water HBs than those within bulk water: depending on the temperature, c_w(t) and c_ps(t) decay to zero between 10–20 ps and 20–50 ps, respectively, at temperatures <330 K. At the highest temperatures both functions decay to zero within 10 ps with the most pronounced speedup of this process occurring at the peptide/water interface. The functions c_w(t) and c_sc(t) are virtually identical, and hence the latter is not shown. Conversely, the function c_sw(t) decays to zero in only a few picoseconds at all temperatures, indicating that only fast HB breaking dynamics occurs for this class of bonds. A curiosity is the function c_pp(t), which instead of decaying to zero still retains a finite value for upwards of 100 ps or more. This behavior arises from the fact that peptide/peptide HBs, once broken, cannot diffuse away from each other inasmuch as they are formed involving the peptide backbone; i.e., they have a finite probability to reform even after long elapsed times due to topology.

View larger version (33K): [in this window] [in a new window]   FIGURE 7  Temperature- and time-dependence of various HB autocorrelation functions, c(t): HBs between (a) water molecules in the bulk, cw(t); (b) first shell solvation water molecules and bulk water molecules, csw(t); (c) the peptide and first shell solvation water molecule, cps(t); and (d) direct HB contacts of the peptide with itself, cpp(t) (see text for definitions). The choice of contour lines is coded in each panel as a, b, and c, where a and c denote the lowest and the highest contours, respectively, and b defines the relative spacing.

View larger version (27K): [in this window] [in a new window]   FIGURE 8  Temperature- and time-dependence of various HB autocorrelation functions, n(t): HBs between (a) water molecules in the bulk, nw(t), and (b) the peptide and first shell solvation water molecules, nps(t). See text for definitions. The choice of contour lines is coded in each panel as a, b, and c, where a and c denote the lowest and the highest contours, respectively, and b defines the relative spacing.

(2)

where k is the rate constant of HB breaking and is that for HB reformation subsequent to breaking; note that = 1/k defines an average HB lifetime. For bulk water (see Fig. 9 a and Table 1), both time constants k_w and show a simple Arrhenius or classical transition-state-theory behavior, i.e., , at all temperatures (with and kJ/mol). Thus, both HB breaking and reformation increase in rate constantly with temperature in the bulk, as expected for a simple thermal process. Note that the bulk water value of k_w at 300 K is close to data from simulations of pure water at room temperature (Luzar, 2000; Luzar and Chandler, 1996a,b; Xu and Berne, 2001), whereas our value of is smaller by a factor of approximately one-half. Although some of this discrepancy may arise from systematic errors as already discussed, it is more likely due to the difference in the dynamical properties of TIP3P water from SPC or SPC/E models (see van der Spoel et al., 1998, and references therein) used in previous studies (Luzar, 2000; Luzar and Chandler, 1996a,b; Xu and Berne, 2001). In particular, TIP3P water has a larger diffusion coefficient (van der Spoel et al., 1998) that allows the water molecules to diffuse away at a different rate after an HB breaking event. Although the bulk water behaves in an expected manner the same cannot be said for HBs that connect the peptide to its surrounding water shell. As with the average structural parameters and the peptide backbone relaxation time, both k_ps and show two distinct Arrhenius temperature regimes bracketed 330 K. Below 320 K the activation energy to break an interfacial HB is higher than that in the bulk, whereas above 330 K it is much lower (4 kJ/mol). Concurrently, HB reforming at the interface (11 kJ/mol) is close to that for breaking, but at 330 K the apparent activation energy ( kJ/mol) vanishes. In terms of rates (see Fig. 9), both processes to break and reform HBs at the peptide/water interface effectively slow down upon heating (with regard to a linear Arrhenius-extrapolation of the low temperature behavior) once a critical threshold, 320–330 K, is surmounted. Hence, the rate to make interfacial HBs essentially reaches a plateau value above 330 K, whereas the one to break these bonds continues to grow with temperature (see Fig. 9).

View larger version (16K): [in this window] [in a new window]   FIGURE 9  Arrhenius plot of rate constants for HB breaking, k, and reformation, , for HBs between (a) first shell protein solvation water molecules and bulk water molecules, ksw,; first shell protein solvation water molecules around the hydrophilic side chains, ksc, ; water molecules in the bulk, kw, , and , ; (b) the peptide and first shell solvation water molecules, kps, ; and , ; (c) the peptide with itself, kpp, . The symbol size covers the error bars and the lines are linear (i.e., Arrhenius) fits; see text for definitions and Table 1 for the resulting fit parameters.

Conversely, the observation of a profound change in the dynamics of the peptide/water interface is in qualitative accord with thermodynamic measurements (Nicolini et al., 2004; E. Schreiner, C. Nicolini, B. Ludolph, R. Ravindra, N. Otte, A. Kohlmeyer, R. Rousseau, R. Winter, and D. Marx, unpublished material; Urry et al., 2002). In particular, the behavior of the thermal expansion coefficient also indicates an overall decrease in the strength of these interactions. This is also consistent with the proposed hydrophobic collapse that occurs for the ITT of elastins (Li et al., 2001b). However, the most intriguing aspect of the present observation is that the peptide/water HB dynamics, which occur on a timescale of 1–10 ps, mirrors the dynamics of the peptide which occur on timescales two orders-of-magnitude larger.

We now investigate the dynamics of the HBs between interfacial water and bulk water as well as those involving the bridging waters. For the latter class a correlation function c_bw(t) is introduced based on an HB population function h_bw(t) which is unity only if a water molecule is simultaneously bound to two or more different residues of the peptide. For the former of these two HB classes the second term in Eq. 2 may be neglected due to the relatively rapid decay of c_sw(t) to zero. We find that k_sw features essentially an Arrhenius behavior at all temperatures, with a similar value of as bulk water (see Fig. 9 a and Table 1). This picture is entirely in accord with the proposition that these surface water molecules have a hindered orientation which strains the HBs that are formed with bulk water. This does not allow for HB reformation after breaking, hence the process of fast HB reformation does not occur for this class of HBs. In passing, we note that it was discovered recently by femtosecond spectroscopy that the addition of simple salts (such as Mg(ClO)₄, Na(ClO)₄, and Na₂SO₄) had no influence on the rotational dynamics of water molec