INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Interactions, delay, and feedback are the three key characteristics of complex fluids or soft matter (de Gennes, 1992). Proteins, whose internal motions are decisive on their folding, stability, and function (Lee and Wand, 2001) are no exceptions. They consist of a structural core, namely the interior, which is mostly dominated by hydrophobic interactions (Makhatadze and Privalov, 1995) and the outer surface, which may be molten (Zhou et al., 1999), composed of hydrophilic residues, and exposed to solvent. Relatively less research efforts have been devoted to the "spacer" region that must provide the necessary conformational degree of freedom for both ends, the hydrophobic interior and the hydrophilic exterior (see, for example, Melchionna et al., 1998). Through the spacer, a remote-control architecture can be created to establish a feedback mechanism between the rigid and more flexible parts (Yilmaz and Atilgan, 2000).

The internal motions of proteins are combinations of the equilibrium transitions among conformational substates and the high-frequency small-amplitude fluctuations about the substates (Frauenfelder and McMahon, 1998). These motions can be recorded effectively by monitoring the positional movements of residues, manifested in Debye-Waller or temperature factors, via molecular dynamics (MD) simulations (Caves et al., 1998; Doruker et al., 2000; Baysal and Atilgan, 2001a,b). The motions confined to each substate or the quadratic envelope of the equilibrium landscape can be determined by simple analytical methods. (Bahar et al., 1997; Atilgan et al., 2001). Though these harmonic vibrations may be too fast to be directly involved in most physiology, through a coupling mechanism between the harmonic motions and the equilibrium transitions (Baysal et al., 1996), these small-amplitude motions may contribute to biological activity.

The presence of conformational substate hierarchy, perceived macroscopically as a complicated heterogeneity of packing densities due to the interior, spacer, and core regions, may delimit the protein-glass analogy (Green et al., 1994), because motions in different regions freeze out at different temperatures (Frauenfelder et al., 1991, 1999). However, as the temperature decreases, the average positional fluctuations of each region shall converge to a common value (Melchionna et al., 1998), because all regions behave like hard, solid materials. At a higher temperature, however, the protein is a soft matter and active, thus, the packing density does not necessarily remain the only parameter that governs the flexibility. Together with the packing order, which constructs a map pointing out the residue pairs in contact, they give a broader view of the nature of collective motions of different regions. Some function-related intrinsic mechanisms may be extracted by analyzing these concerted motions, such as the active sites (Bahar et al., 1999; Atilgan et al., 2001; Baysal and Atilgan, 2001b) and the remotely controlling residues during binding (Baysal and Atilgan, 2001a). Therefore, it is of utmost interest to investigate how flexibility of proteins responds to temperature changes under physiological or extreme conditions (Dvorsky et al., 2000; Vitkup et al., 2000).

There has been a lot of research trying to understand the glassy relaxation behavior of proteins (Bizzarri et al., 2000) and RNA chains (Pagnani, et al. 2000). Recent incoherent neutron-scattering measurements of elastic intensities of proteins (Bicout and Zaccai, 2001) identify the conformational flexibility as essential for enzyme catalysis and that the positional fluctuations are important for proteins' functionality because they behave like a "lubricant" (Zaccai, 2000). These experimental studies have measured a protein dynamics force constant to quantify the molecular resilience and paved the way for the investigation of time-dependent mechanical behavior. So, the investigation of the effect of temperature on the relaxations of the time-delayed correlations of the positional fluctuations may take us one step further in understanding proteins as soft matter. We, therefore, consider the bovine pancreatic trypsin inhibitor (BPTI), which is a well-studied inhibitor, for instance, by normal-mode analysis (Levitt et al., 1985), NMR spin-relaxation measurements and MD simulations (Smith et al., 1995). We monitor the positional fluctuations and their time-delayed correlations at different temperatures with the aid of extensive MD simulations. A stretched exponential with temperature-dependent relaxation time is fitted to the relaxation of time-dependent correlations. This fit leads to two parameters related to the thermodynamic and kinetic aspects of the process, respectively: 1) the stretching exponent is an indicator of the degree of concerted rearrangements among different conformational substates (Frauenfelder et al., 1991), and 2) the temperature dependence of the inverse relaxation time is of the same form as that of polymers at temperatures above but close to the glass transition. It is further shown that the effect of temperature on the time-dependent mechanical behavior is equal to a stretching of the real time for temperatures above or below the glass transition. We thus refer to the behavior of BPTI as thermorheologically simple.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Molecular dynamics simulations

We have used BPTI as the model system, a 58-amino-acid inhibitory protein. The initial structure is that from the Protein Data Bank (PDB) (Berman et al., 2000), PDB code 5pti (Wlodawer et al., 1984). All atoms are treated explicitly. In the PDB structure, two separate locations for the side chains of Glu-7 and Met-52 are reported; the first of each of these locations is selected in the initial structure of the simulations. The PDB structure is soaked in water such that it forms a 6-Å-thick layer around the molecule. Together with the structural water molecules reported in the PDB structure, this treatment leads to a total of 703 solvent molecules. A minimum amount of 0.40 g water per gram protein is required to fully activate protein dynamics and functionality (Gregory, 1995; Bizzarri et al., 2000). The constituents in this study correspond to 1.95 g water per gram protein, which is well above this threshold.

The consistent valance forcefield (Dauber-Osguthorpe et al., 1988) implemented within the Molecular Simulations Inc. InsightII 2000 software package is used in the initial structure refinement and the subsequent MD calculations. Group-based cutoffs are used with a 10-Å cutoff distance. A switching function is used with the spline and buffer widths set to 1.0 and 0.5 Å, respectively. The system is first energy minimized to 7.5 × 10⁵ kcal/mol/Å of the derivative by 4973 conjugate gradients iterations. The minimization takes 1292 s on a Silicon Graphics Origin 200 computer with a 350 MHz CPU. These new coordinates of the system are treated as the initial coordinates in all the MD simulations.

All bonds of the protein and water molecules are constrained by the RATTLE algorithm (Andersen, 1983). Initial velocities are generated from a Boltzmann distribution at the designated temperature. Integration is carried out by the velocity Verlet algorithm. The systems are equilibrated for 200 ps with a time step of 2 fs, while maintaining the temperature by direct velocity scaling. Since the fluctuations increase at higher temperatures, we resort to longer simulations to gather more reliable data as the temperature is increased. Thus, the data-collection stages are of length 2.0-2.8 ns, depending on the temperature: 2.0 ns for T < 230 K, 2.4 ns for 230 < T < 290 K, and 2.8 ns for T > 290 K. Also, second independent MD runs of duration 2.0 ns are made for T > 290 K. At this stage, a time step of 1 fs is used and temperature control is achieved by the extended system method of Nosé. (Nosé, 1984) Data are recorded every 2 ps, and each 400-ps portion of the trajectories is treated as a separate sample. We thus have 5-12 data sets at each temperature, and the calculated quantities are averaged over these.

The fluctuation vector

Throughout this study, we study the time- and temperature-dependent properties of the fluctuation vector attached to the C atoms of the protein. We compute the fluctuation vector, R, as follows. For any given 400-ps piece of the trajectory, we first make a best-fit superposition of the recorded structures to the initial structure by minimizing the root mean square deviations of the C atoms. We then compute the average structure, R(T), from the 200 best-fitted structures. Here, the brackets denote the time average. We finally make another best-fit superposition of the recorded structures to this average structure. Each structure of this final trajectory is denoted by R(t, T), and the coordinates of the ith residue is given by R_i(t, T). In this manner, the resulting trajectory has contributions from the motions of the internal coordinates only. The fluctuation vector for a given residue i at a given time t from a given trajectory obtained at temperature T, R_i(t, T), is thus the difference between the position vectors for the ith residue of the best-fitted and the average structures,

(1)

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The "protein glass transition" and thermal fluctuations

The glass transition is a second-order phase transition, and is therefore indicated by a continuous step in heat capacity. We first start out by computing the temperature range over which the glass transition occurs, which has been experimentally measured to occur at temperatures as low as 150 K for bacteriorhodopsin (Réat et al., 1998) and as high as 220 K (Daniel et al., 1998) for a variety of proteins. For a given MD simulation, the heat capacity at the simulation temperature may be computed from (E  E)²/k_BT². Using this procedure to compute the heat capacity of the system as a function of temperature, the glass transition is found to occur in the temperature range of 190-210 K. Similarly, a sudden drop in the heat capacity is observed at 320 K and is attributed to the onset of unfolding (melting) in the protein.

One might raise the question on how the glassy behavior of the protein is coupled to the glassiness of the hydration water. Above a certain hydration level (Gregory 1995; Bizzarri et al., 2000), which is satisfied in this study as pointed out in the MD simulation details, water and protein form an interacting system with unique properties that would not be found either in the dry protein or in bulk water. In fact, MD simulations and elastic incoherent neutron-scattering experiments on dry and partially hydrated protein on the one hand (Arcangeli et al., 1998; Lehnert et al., 1998), and on bulk water on the other (Sciortino et al., 1996, and references cited therein) corroborate this finding. Moreover, Bizzarri et al. have made a detailed analysis of plastocyanin hydration water, and conclude that the interaction of water molecules with the protein atoms is required to activate protein dynamics. Finally, a recent study (Pal et al., 2002) confirms that the structured water molecules around the protein are needed to form a system that maintains the enzymatic function of the protein subtilisin Carlsberg.

Having identified the locations of the glass transition and melting temperatures, T_g and T_m, respectively, we now turn our attention to the properties of the fluctuation vectors of the C atoms. The fluctuations of the atoms due to harmonic, quasiharmonic, anharmonic, valley hopping, and other types of motions occurring in the molecule are reflected in, for example, the intensities of the diffraction spots in x-ray diffraction experiments. Thus, the Debye-Waller factors yield the mean square displacements of all nonhydrogen atoms in these experiments, and are given by B = 8²/3u, where u_i is the motion of the ith atom relative to the reference, i.e., it is the atomic displacement. Similarly, incoherent quasielastic neutron scattering experiments provide information on the mean-square displacements of the hydrogen atoms. Such data are very important inasmuch as they yield information on the types of motions operating on the whole molecule or parts of it, and also lead to deriving general concepts of protein dynamics (Frauenfelder and McMahon, 1998).

We present the thermal fluctuations averaged over the C atoms of the protein in Fig. 1 (filled circles). The temperature-dependent behavior of atomic fluctuations in BPTI displays a similar behavior to that observed in other experimental (Tsai et al., 2000; Bicout and Zaccai, 2001) and simulation (Smith et al., 1990; Melchionna et al., 1998; Tarek et al., 2000) studies on various proteins: as the system is heated, the onset of large thermal fluctuations is observed around T_g. Melchionna et al. have shown that the solvent-exposed parts of the protein Cu, Zn superoxide dismutase display larger fluctuations than the residues that have medium or no solvent exposure; the latter two show similar fluctuations (Melchionna et al., 1998). Their interpretation of this observation is that the glass transition is driven by the exterior residues. To corroborate these results, we next analyze the fluctuations in the protein interior and exterior separately (open squares and open circles in Fig. 1, respectively). This separation may be done in one of several ways. Using the solvent-accessible surface areas, for instance, will not let us distinguish between residues just below the protein surface and those in the core of the protein. We thus utilize the depth program (Chakravarty and Varadarajan, 1999), which differentiates between such residues by calculating the depth of a residue from the protein surface. At a residue depth of ~4 Å, the size of the fluctuations for the surface and interior residues converge to the same values below T_g. Above T_g, much larger fluctuations are observed at the surface residues, as in the work of Melchionna et al. (1998), although the fluctuations grow with temperature for the protein interior as well.

View larger version (19K): [in this window] [in a new window]   FIGURE 1   Thermal fluctuations averaged over all residues (filled circles), interior residues (open squares), and exterior residues (open circles). The distinction of exterior and interior residues is described in the text. Line fits to the data in the temperature range Tg < T < Tm gives an estimate of the Kauzmann temperature TK, which is found to be the same (TK  171 K) for the protein interior and exterior.

To answer the question, is the protein interior still glassy above the protein-glass transition, we turn to the theory of glasses. The Kauzmann temperature, T_K, is the temperature that would have been attained at zero entropy had the glass transition phenomena not been operative (Debenedetti and Stillinger, 2001; Stillinger et al., 2001). It may be obtained by extrapolating the fluctuation versus temperature data for the amorphous molecule to the point where the fluctuations cease. At zero fluctuation, a single conformation would exist, and the entropy would be zero. Because T_K > 0 K, violating the third law of thermodynamics, this phenomenon is referred to as the Kauzmann paradox (Stillinger et al., 2001). We make the extrapolation on the fluctuation data for the whole protein and the protein interior and exterior separately at the temperature interval 220-310 K (above T_g, but below T_m), and obtain the same Kauzmann temperature of T_K  171 K in all regions (Fig. 1). Thus, the protein interior goes through the glass-transition phenomena at the same temperature window as the rest of the protein. This conclusion is valid for BPTI only, and should be proven for other proteins before generalization, because there is experimental evidence from elastic incoherent Neutron scattering experiments (Réat et al., 1998) that, in bacteriorhodopsin, T_g is ~150 K when the global motions are explored, but it is slightly above 200 K for the residues that characterize the active site only.

Relaxation phenomena

We now characterize the motion of the fluctuation vector by a relaxation function of time and temperature, C(t, T),

(2)

where the bar and the brackets denote the average over all residues, and the time average, respectively. If the relaxation function may be characterized by a single phenomenon, it can be modeled with a simple exponential with a single rate constant. However, usually C(t, T) will have contributions from many different homogeneous processes with different relaxation times, and their collective effect on relaxation will be observed as heterogeneous dynamics (Deschenes and Vanden Bout, 2001), represented by

(3)

The above equation may be approximated by the stretched exponential, or the Kohlrausch-Williams-Watts function (Kohlrausch, 1874; Williams and Watts, 1970)

(4)

where k(T) is an inverse relaxation time representing the average contribution of all the processes affecting the relaxation of the fluctuation vector. (T) reflects the complexity of the processes involved (0    1). For one simple process characterized by a simple exponential decay,  = 1.  tends to decrease from 1, not only if a larger number of contributing processes (n > 1) is at play, but also if these processes have time or length scales spanning different orders of magnitudes. It has recently been shown by incoherent quasistatic scattering experiments on partially hydrated lysozyme and ribonuclease A powder samples, that  obtained from relaxation data on proton fluctuations decreases with increasing temperature, indicating the emergence of multiple relaxation modes (Tsai et al., 2001). In the same study,  is reported to be constant at 1 for dehydrated samples, indicating completely harmonic motions on the timescale of the measurements.

We have examined the stretched exponential fits to the decays of C(t, T) data at all the temperatures studied. We find that the data is well approximated by the fits at all time scales slower than ~15 ps, especially at temperatures above the transition. The inset in Fig. 2 shows the relaxation of the positional fluctuations at 150, 180, and 300 K, and the curves fitted to these data using Eq. 4. At the time scales faster than 15 ps, there is a sharp decay that is not represented by the fit. This may very well be because a power-law behavior governs the dynamics in this regime. In fact, there is evidence for such behavior in the sub-picosecond regime for the intermediate scattering of water oxygen atoms around plastocyanin (Bizzarri et al., 2000). However, in the present study, we investigate the relaxation of the positional fluctuations, which probes the collective dynamics of the protein. As such, the time scales of interest are above 20 ps.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   Temperature dependence of the exponent  (Eq. 4). The curve fit is provided to guide the eye. In the inset, the relaxation functions, C(t, T), at the temperatures of 150 K (), 180 K () and 300 K () are displayed. The curves are best fits to the data using Eq. 4.

The temperature dependence of  is displayed in Fig. 2. The general trend is that   0.2 for the protein glass, shows an increase during the transition period, and   0.4 above T_g. However, the general expectancy would be to have a decreasing trend in  with increasing temperature as in the above experimental example. In contrast, the relaxation of the C fluctuation vector used in this study is on a much larger scale than the relaxation defined by the fluctuations of protons. The relaxation function is averaged over all C atoms. Also, the trajectories of the fluctuation vectors are defined for the protein conformations sampled around the average structure obtained from a portion of the trajectory (see Methods). Thus, C(t, T) has contributions from all kinds of processes, such as harmonic motions between bonded atoms, solvent effects, transitions between conformational substates of the backbone and side chains, and larger cooperative fluctuations occurring in, for example, the loop regions. A small  at low temperatures, therefore, does not mean that there is a greater variety of processes at low temperatures, or that each of these processes is more complex (Ediger and Skinner, 2001). Rather, because the conformations are trapped in various substates in the glass, there is no communication between these modes that spans a large range of frequencies. At higher temperatures, some paths are created between previously disconnected modes, leading to cooperativity and ultimately a simpler overall dynamics reflected in the higher  value. This viewpoint does not preclude the emergence of new processes at higher T that did not occur in the glassy state.

To check on these assumptions, we have carried out the same analysis on the relaxation of the NH vector of the Arg-1 residue. Because this vector belongs to the outermost part of the protein, its local environment is the least crowded. Therefore, its motion is the one that is the most decoupled from the rest of the protein, having heavy contributions from simpler processes such as those due to solvent bombardment. Although our MD simulations may be assumed to have been carried out on a fully hydrated protein as opposed to partially hydrated samples, and the experiments have the advantage of filtering the time scales such that the time window examined corresponds mainly to the time scales invoked by H atoms moving with the amino-acid side chains of the protein, the relaxation behavior of this NH vector from simulations should nevertheless be more similar to the relaxation of proton fluctuations observed in experiments. Indeed, for this NH vector, we find that  decreases monotonically from 0.63 at 150 K to 0.22 at 320 K.

We next analyze the temperature dependence of k (Eq. 4) in Fig. 3. Inasmuch as k is an inverse relaxation time, averaged over the many processes affecting the observed relaxation, a very fast relaxation behavior dominates below T_g. This is because large-scale motions are too slow to be observed on the time scale of the simulations due to the frozen-in substates, and the dynamics is dominated by fast modes on the sub-picosecond timescale. As the glass transition sets in, slower relaxation phenomena on the order of picoseconds begin, reflected in the sharp drop in the k values. Above glass transition, the average relaxation time scale is still on the order of picoseconds, but there is an increase trend in the values as the temperature rises due to the cooperativity of different modes. At 320 K, where unfolding has started, the relaxation time drops once more, because the cooperativity due to the folded nature of the protein is lost.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Temperature dependence of the inverse relaxation time, k in Eq. 4, expressed as log k versus 1000/T. The curve fit is provided to guide the eye. Inset: The same data at physiological temperatures (Tg < T < Tm) expressed in the Vogel-Tammann-Fulcher form for all residues (filled circles), and interior (open squares), and exterior residues (open circles) separately.

For temperatures higher than T_g and lower than T_m, which corresponds to the physiological temperatures of proteins, the temperature dependence of k may be characterized by the Vogel-Tammann-Fulcher expression (Angell, 1995; Debenedetti and Stillinger, 2001),

(5)

This relationship suggests a connection between thermal fluctuations, a thermodynamic property, and the relaxation phenomena, a kinetic property. The physical meaning attributed to T₀, the Vogel temperature, is that it is the temperature at which infinite viscosity is reached in a material during the cooling process. In reality, the glass transition sets in before T₀ may be reached; i.e., T_g > T₀. Moreover, it has frequently been observed that T₀ coincides with T_K (Angell, 1995), which is computed to be ~171 K for all regions of the present system (see Fig. 1). Note that the description of T_K is based on a thermodynamic view of the glass transition, whereas that of T₀ finds its roots in a dynamic picture, and their correspondence suggests a connection between kinetics and thermodynamics. Adam and Gibbs have formulated this connection, with the underlying assumption that the slowing down in the dynamics close to T_g originates from the decrease in the number of configurations the system can sample (Debenedetti, 1996). A fit of log k versus 1/(T  T₀) is given in the inset to Fig. 3 in the temperature range 220-310 K. The preexponential factor, A, which itself is a weak function of temperature, gives an estimate of the average collision frequency for the relaxation phenomena described by k(T). Thus, A is on the order of ~400 ns¹ for BPTI in the temperature range of interest. The interior residues are active at a higher frequency than the exterior residues, with values of A on the order of 500 and 300 ns¹, respectively. However, the fact that we can fit the data for the interior and exterior residues equally well to the Vogel-Tamman-Fulcher form (R²  0.7, R²  0.9 excluding the data point corresponding to T = 277 K) suggests that the protein interior, though rigid, is not glassy in this temperature interval.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We have previously shown that the librational motions around the torsional degrees of freedom of even a simple polymer chain carries information on the conformational jumps experienced by the torsional angles (Baysal et al., 1996). There is also supporting experimental evidence from neutron-scattering data that the fast (nanosecond-to-picosecond) thermal motions are essential for the slower (millisecond) relaxations underlying conformational changes (Lehnert et al., 1998). Time-delayed correlations of residual fluctuations, C(t, T), may be viewed as the "dynamic flexibility" characterizing the sampled protein conformations. It is also a function of librational motions, carrying information on the large conformational rearrangements that take place within the protein. A close inspection of the form of C(t, T) (Eq. 4) reveals that time and temperature effects are superposable, leading to thermorheologically simple behavior. Below T_g, protein conformations are trapped in local minima (  0.2), whereas, at physiological temperatures, sampling of different minima is possible (  0.4). We thus associate  with configurational entropy. The inverse relaxation time, k(T), characterizing the kinetic phenomena and  characterizing thermodynamic properties, stretch time and length scales, respectively.

At physiological temperatures, the protein interior and exterior are both mobile, suggesting that glass properties do not prevail. The difference lies in the size of the fluctuations and the time scales of motion operating in these regions. Yet, it is this subtle distinction between the flexibilities of the protein interior/exterior that maintain the integrity of the protein while allowing the conformational freedom necessary for biological activity (Zaccai, 2000).

In view of biological and biotechnological applications, it is also of interest to study the transition region in more detail. The fact that librations lead to large conformational motions above the transition temperature, as opposed to the lack of such mechanisms at low temperatures, as well as the broad nature of the transition of  from 0.2 to 0.4 (Fig. 2), suggests that different relaxation phenomena are invoked at different temperatures. The onset of transitions between conformational substates may further be affected by the nature of the solvent. It seems plausible that these mechanisms may be controlled using subtle temperature modulations or different solvents.