INTRODUCTION

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We have been engaged in a computational study aimed at understanding the basic physical principles underlying protein aggregation, a phenomenon with serious biomedical problems that include a host of ultimately fatal protein deposition diseases (Clark and Steele, 1992; Eaton and Hofrichter, 1990; Gallo et al., 1996; Massry and Glasscock, 1983; Moore and Melton, 1997; Selkoe, 1991; Simmons et al., 1994). The long-term goal of our work is to reveal molecular-level mechanisms underlying protein aggregation and fibril formation. Because protein aggregation often involves the association of -strands, or the conversion of -helices to -strands (Benzinger et al., 1998, 2000; Burkoth et al., 1998; Esler et al., 1996, 2000; Lazo and Cowning, 1998; Lynn and Meredith, 2000; Sunde et al., 1997; Zhang et al., 2000), we have focused part of our efforts on the development of models of -strand proteins. These models must be capable of capturing the essential physical features of real -strand proteins and at the same time be simple enough to allow the simulation of multi-protein systems with current computer capability. In this paper we introduce one such model, a four-strand protein interacting via a hybrid Go-type potential (Go and Taketomi, 1978, 1979; Taketomi et al., 1975; Ueda et al., 1978) and describe the folding properties of an isolated model chain as a function of temperature. This sets the stage for a later paper in which we examine the folding of several of these model -strand proteins into fibrils.

Low-resolution or simplified protein models are best suited to our purpose because they allow us to study the behavior of single or multi-protein systems over relatively long time scales. These models provide numerous insights into the thermodynamic and kinetic properties of protein folding. The low-resolution models used in simulation studies of protein folding can be divided into two classes: on-lattice models (Bratko and Blanch, 2001; Chan and Dill, 1994; Dill and Stigter, 1995; Go and Taketomi, 1978, 1979; Gupta and Hall, 1997, 1998; Gupta et al., 1999; Harrison et al., 1999, 2001; Kolinski et al., 1995; Lau and Dill, 1989; Miller et al., 1992; Skolnick and Kolinski, 1991; Taketomi et al., 1975; Ueda et al., 1978) and off-lattice models (Dokholyan et al., 1998, 2000; Guo and Thirumalai, 1995, 1996; Guo and Brooks, 1997; Guo et al., 1997; Nymeyer et al., 1998; Pande and Rokhsar, 1998; Shea et al., 2000; Takada et al.,. 1999; Zhou and Karplus, 1997a, 1997b, 1999). The first low-resolution lattice models of proteins were introduced by Go et al. (Go and Taketomi, 1978, 1979; Taketomi et al., 1975; Ueda et al., 1978) in the late 1970s. In their models, the protein was treated as a sequence of beads on a two- or three-dimensional lattice with the native tertiary structure being determined by a set of preassigned attractive residue-residue pairs at the nearest-neighbor lattice sites. Another class of lattice model used to describe protein folding is the HP lattice model proposed by Lau and Dill (1989). To mimic the hydrophobic effect, the protein is modeled as a chain of polar (P) and hydrophobic (H) residues with an attractive potential between nonbonded H beads and zero potential between nonbonded H and P or P and P residues. Simulation results for H-P models have provided useful information on protein folding-refolding pathways for isolated chains (Chan and Dill, 1994; Gupta and Hall, 1997; Gupta et al., 1999; Miller et al., 1992). In recent years these models have been extended to the treatment of protein aggregation for multichain systems (Bratko and Blanch, 2001; Dill and Stigter, 1995; Gupta and Hall, 1998; Harrison et al., 1999, 2001).

Recently, off-lattice models have gained attention as they can more accurately represent the real features of proteins. Guo et al. (Guo and Thirumalai, 1995, 1996; Guo and Brooks, 1997; Guo et al., 1997) have investigated the thermodynamics and kinetics of protein folding using Langevin simulations on off-lattice protein models. Their heteropolymer models consist of N connected beads, which correspond to three types of residues: hydrophobic, hydrophilic, and neutral. They found two characteristic temperatures, a collapse transition temperature T and a folding transition temperature T_f in their thermodynamic studies of -barrel and four-helix bundle proteins. Zhou and Karplus (1997a, 1997b, 1999) recently introduced an off-lattice heteropolymer model for a three-helix bundle protein with a hybrid Go-type potential. They used a bias gap parameter, g, to vary the strength of native contacts relative to that of nonnative contacts. Despite the model's simplicity, thermodynamic studies on this protein-like model exhibit complex protein transitions that have been observed experimentally including a collapse transition, a disordered globule to ordered globule transition, an ordered globule to native transition, and a transition to a surface frozen inactive state (Ptitsyn, 1995).

In this paper we investigate the thermodynamic phase behavior of three different four-strand antiparallel -strand peptides: the -sheet, the -clip, and the -twist. These -strand peptides each have 39 connected residues (beads) and have different native state conformations as shown in Fig. 1. Nonbonded beads can interact through a hybrid Go-type potential (Go and Taketomi, 1978, 1979; Taketomi et al., 1975; Ueda et al., 1978) modeled as a square-well or square-shoulder potential depending on the value of the bias gap parameter g. Discontinuous molecular dynamic (DMD) simulations (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996) were performed on model protein systems with bias gaps ranging from 0.3 to 1.5. The bias gap measures the difference in interaction strength between the native and nonnative contacts: the larger it is the more the native state is favored over the nonnative state. Intermediate values of the bias gap are thought to be the most representative of a real protein in their equilibrium and dynamic behavior. By exploring how variations in our model protein's bias gap, an artificial measure of its "preference" for the native state, influence the types of phase behavior observed, we can get a feeling for how a real protein's preference for its native state, as measured for example by the energy difference between the denatured and native state, is manifested in the protein's phase diagram and vice versa.

View larger version (34K): [in this window] [in a new window]   FIGURE 1   Global energy minimal structures for: (a) the -sheet, (b) the -clip, and (c) the -twist model proteins.

	MODELS AND SIMULATION METHOD

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We consider three different off-lattice protein models whose native states are four-strand antiparallel -strand peptides. The global energy minimal structures for the three different -strand models are shown in Fig. 1. The "zig-zag" shaped native protein in Fig. 1 a is called a -sheet; sometimes it is known as betabellin (Lim et al., 2000). The paper clip and helix-shaped native structures in Fig. 1, b and c are called -clip and -twist, respectively (Kolinski et al., 1995). A qualitative difference between the -sheet structure and the other two structures is immediately apparent. In the native state, the -sheet is two-dimensional (planar), whereas the -clip and -twist have well-defined three-dimensional four-barrel structures. The reduced squared radius of gyration per bead, R/²N, the total number of native contacts, N, and the reduced squared end-to-end distance, r/², for the three global energy minimal structures are summarized in Table 1. Here N is the number of beads and  is the diameter of each bead along the chain. Note that chains with a large number of native contacts are generally more compact and consequently have a small value of the squared radius of gyration.

View this table: [in this window] [in a new window]   TABLE 1   Reduced squared radius of gyration per bead, R/2N, the total number of native contacts, N, and the reduced squared distance between beads number 1 and 39, r/2, in the global energy minimal structure for the three different -strand models

Each model protein contains 39 connected beads. Each bead represents an amino acid residue that can be regarded as being localized at the C atom. Nonbonded beads can interact with each other through a square-well or square-shoulder potential (Zhou et al., 1996, 1997),

(1)

in which  is the bead diameter,  is an energy parameter with  > 0,  is the square-well diameter with  = 1.5 throughout, and r is the distance between residues i and j. The quantity B_ij, the square-well depth or square-shoulder height represents the interaction strength between nonbonded residue pair i and j and is defined as

(2)

in which B_N and B_O are measures of the relative strengths of the energies associated with the native and nonnative pair interactions in this Go-type potential (Go and Taketomi, 1978, 1979; Taketomi et al., 1975; Ueda et al., 1978). Nonbonded pairs of beads that are in contact in the global energy minimal structure experience an attractive interaction, i.e., B_N < 0, when their square-wells overlap. On the other hand, nonbonded pairs of beads that are not in contact in the global energy minimal structure experience either an attractive interaction (B_O < 0) or a repulsive interaction (B_O > 0). The parameter B_O can be either positive or negative depending on the size of the bias gap parameter g,

(3)

in which g is the bias gap (Zhou and Karplus, 1997a, 1997b, 1999). The bias gap measures the ratio of the interaction strength between the native contacts and nonnative contacts. Note that for g > 1, B_O > 0; in this case nonnative contacts are repulsive so that the native state structure is strongly favored over any nonnative state structure. For 0 < g < 1, B_O < 0; all nonbonded contact pairs are attractive, but native contacts are always more favorable than nonnative contacts. For g = 0, native and nonnative contacts are equally favorable, B_O = B_N; the model reduces to a homopolymer. The bias gap is an artificial measure of a model protein's preference for its native state; in a real protein this preference for the native state might be measured for example by the energy difference between the native and nonnative state.

The interaction between two bonded beads, i and i + 1, is given by an infinitely deep square-well potential,

(4)

in which  is a flexibility parameter that controls the bond length. The bond length between neighboring beads can be varied over a small distance between the infinitely high potential barriers at (1  ) and (1 + ). This method for creating a flexible bond length was introduced by Bellemans et al. (1980). The flexible bond-length parameter is set to  = 0.1 through the paper.

Simulations were performed using the DMD algorithm (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996). All simulations were started from randomly generated configurations, which were obtained from self-avoiding random walks. The initial velocities were chosen at random from a Maxwell-Boltzmann distribution at the simulation temperature. This DMD algorithm is very efficient for calculating the properties of systems containing hard spheres or square-well spheres (Smith et al., 1996; Zhou et al., 1997) and is significantly faster than standard molecular dynamics on Lennard-Jones chains. The algorithm proceeds by searching for the next collision time and collision pair, advancing all beads in time to the next collision event, and then calculating the velocity change of the colliding pair. The DMD simulation was conducted in the canonical ensemble with a constant number of particles, volume, and temperature. To maintain constant temperature, the Anderson stochastic collision method (Anderson, 1980) was used. In this method the system's particles collide with imaginary or ghost particles, which serve as an effective heat bath. In the collision with the imaginary particles, the velocity of each bead is reassigned from a Maxwell-Boltzmann distribution at the simulation temperature. A reduced ghost particle density of n_g*  n_g³ = 0.1, in which n_g is the number density of ghost particles, was used to ensure that 1% ~ 10% of the collisions were ghost particle collisions (Zhou and Karplus, 1999; Zhou et al., 1997).

To determine the location of the collapse transition, the mean-squared radius of gyration, R_g², was determined in which

(5)

with N equal to the number of beads, x_i, y_i, z_i equal to the coordinates of bead i, x_c, y_c, z_c equal to the center of mass coordinates of the chain, and _cf denoting a configurational average. To save computing time, the summation is taken over configurations sampled every 1000 collisions after equilibration. Other transitions were determined by calculating the reduced specific heat,

(6)

and the reduced internal energy,

(7)

in which k_B is Boltzmann's constant and   denotes an average over all collisions after equilibration. In calculating the specific heat and energy according to Eqs. 6 and 7, the weighted histogram method (Ferrenberg and Swendsen, 1989; Zhou et al., 1997) was used. In this method, temperature-independent degeneracy factors are extracted from the energy probability distribution in simulation runs at different temperatures. Once the degeneracy factors are known, a partition function is obtained for any temperature, and consequently all thermodynamic quantities can be calculated. Details of the method were reported elsewhere (Zhou et al., 1997). By plotting C_v* versus T*, in which T* is the reduced temperature, T*  k_BT/, one can determine the equilibrium transitions from the locations of the peaks and plateaus.

The degree of bead mobility can be characterized by the root-mean-squared pair separation fluctuation, _B, defined by

(8)

in which r_ij denotes the separation between beads i and j, and  _cf denotes the configurational average. To characterize the liquid-to-solid transition of the system, the root-mean-squared fluctuation for bead i, _Li, defined to be

(9)

yields the Lindemann disorder parameter,

(10)

in which bead positions, r_i (i = 1 to N), are obtained by shifting the center of mass coordinates to the origin. The quantity _Li in Eq. 9 defines the Lindemann disorder parameter for an individual bead, and _L in Eq. 10 defines the Lindemann disorder parameter. The Lindemann criterion (Lindemann, 1910) is applied to the system to determine the solid or liquid phase. A system is regarded as a solid if its Lindemann disorder parameter, _L, is in the range of 0.1 to 0.15 (Bilgram, 1987; Löwen, 1994; Stillinger, 1995; Zhou et al., 1999b), whereas a substance with _L > 0.15 is considered liquid.

The progression of a conformation toward the native state can be monitored by introducing the fraction of native contacts formed (Lazaridis and Karplus, 1997; ali et al., 1994), Q, defined by

(11)

in which N_native represents the number of native contacts in a given conformation and N is the total number of native contacts in the global energy minimal structure. The fraction of native contacts Q has a value between 0 and 1. When Q = 1 the model system can be regarded as being in the native structure, whereas when Q  0 the chain becomes a random coil that indicates the denatured state.

Simulations were performed for up to 10⁹ collisions to ensure equilibration. Each simulation was started from a different initial configuration at the temperature of interest. Equilibrium averages were typically taken by discarding the first one-half of the simulation data. All equilibrium results were averaged over at least four independent runs. For simulations at low temperature, a random initial configuration was generated at a relatively high temperature of T* = 1.0 and then gradually annealed toward the target temperature after which the equilibrium simulation was performed. This procedure can prevent the system from becoming trapped in a metastable state, i.e., a local free-energy minimum. Simulations were performed for bias gaps in the range of g = 0.3 to 1.5 at reduced temperatures in the range T* = 0.07 to 5.0.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Thermodynamic properties of -sheet

DMD simulations were performed for the -sheet protein to investigate phase behavior as a function of the size of the bias gap. At bias gaps of g = 0.3, 0.7, 0.9, and 1.3, thermodynamic averages including the squared radius of gyration R_g², the root-mean-squared pair separation fluctuation _B, the specific heat C_v, and the internal energy E were calculated as a function of temperature. These results for g = 0.3, 0.7, 0.9, and 1.3 are shown in Fig. 2 (a-d), respectively. All quantities are described in terms of reduced units in the figures. Statistical averages were obtained over at least four independent simulations. The error bars are the standard deviation in the measured values and are only shown for values larger than the size of the symbol. The results for the specific heat and the energy were obtained using the same weighted histogram method (Ferrenberg and Swendsen, 1989) that was introduced in the study of homopolymers (Zhou et al., 1997). The equilibrium transition temperatures of the model can be identified from the maximal value of the temperature derivative of the squared radius of gyration, dR_g²/dT*, and/or from peaks or plateaus in the specific heat.

View larger version (19K): [in this window] [in a new window]   FIGURE 2   Average values of the reduced squared radius of gyration per bead, Rg2/2N, the root-mean-squared pair separation fluctuation, B, the reduced specific heat per bead, Cv*/N, and the reduced internal energy per bead, E*/N, as a function of the reduced temperature, T*, for the -sheet at four selected values of the bias gap, (a) g = 0.3, (b) g = 0.7, (c) g = 0.9, and (d) g = 1.3.

To determine the collapse transition, the temperature derivative of the squared radius of gyration, dR_g²/dT*, is calculated along a spline fit of the data in R_g² versus T* graphs as shown in Fig. 2. For a small bias gap, g = 0.3, in Fig. 2 a the maximal value of dR_g²/dT* occurs at T* = 1.47, indicating the collapse transition temperature. For g = 0.7, 0.9, and 1.3 in Fig. 2, b, c, and d, the collapse transitions occur at T* = 1.0, 0.95, and 0.88, respectively. The chain has a conformational change from a random coil to the disordered globule state at the collapse transition. The disordered globule has more native contacts than the random coil, but due to the disordered nature of the state, none of the structures appear to be well ordered or native-like. The properties of disordered globule are similar to those observed experimentally for a premolten globule (Uversky and Ptitsyn, 1996a). Additional evidence for the existence of the collapse transition at g = 0.3 is that the _B versus T* curve possesses a broad, but clear, maximum located at a temperature just above T* = 1.47. In addition, a plateau in the specific heat is observed near this temperature. In contrast to the broad peak for g = 0.3 in the _B curve, for g = 0.7, 0.9, and 1.3, there is a well-defined peak in the _B curve that evidently supports the collapse transition. The collapse transition peak in the _B versus T* curve can be explained in the following way. The root-mean-squared pair separation fluctuation _B increases, as expected, as T* increases but then decreases above the collapse transition temperature because the stretched chain in the random coil state is constrained by chain connectivity, resulting in less fluctuations in the bead-bead separations.

The transition from the disordered globule state to the ordered globule state is associated with a distinct peak in the specific heat. This transition is strong with large energy change indicating a large conformational change. For g = 0.3, the peak is found at a temperature of T* = 0.58, far below the collapse transition. This indicates that the disordered globule is a stable state over a wide temperature range, a trend characteristic of smaller values of g. However, for g = 1.3 the peak is found at a temperature of T* = 0.80, only slightly below the collapse transition temperature. The collapse transition will eventually coincide with the disordered-to-ordered transition for g > 1.5. The ordered globule is similar to the molten globule, a liquid-like phase that exists at temperatures between the native and the unfolded coil state. This molten globule is a thermodynamic (not kinetic) intermediate that has a native-like secondary structure but has deficient side-chain interactions (Ptitsyn, 1995; Ptitsyn and Uversky, 1994; Schulman and Kim, 1996; Uversky and Ptitsyn, 1996b).

The folding transition from the ordered globule state to the native state can be identified from the plateaus in the specific heat. At this transition, there are no major changes in structure, but the chain has a subtle conformational change toward the native state. For the smaller gap models, g = 0.3, no evidence for the folding transition was observed. This is because for g < 0.7, a well-ordered native state does not appear to exist; instead kinked or rolled structures were observed. Consequently, no plateaus associated with the folding transition were observed for the smaller gap models. For the g = 0.7 model, although a peak in the specific heat is found at a temperature of T* = 0.35, the peak is not related to the folding transition. It is associated with the transition from the ordered globule state to a more-collapsed ordered globule state that we will call the partially ordered globule state. For the larger gap models, g = 0.9 and 1.3, the plateaus at T* = 0.40 in the specific heat are associated with the folding transition. This folding transition temperature is much lower than the disordered-to-ordered globule transition. This indicates that the ordered globule is a stable state over a wide temperature range for the larger gap models.

The thermodynamic results presented in Fig. 2 differ from those found in real proteins in that the specific heat in the unfolded state (random coil) is lower than in the folded state (ordered globule) even for the larger gap models, g > 0.9; the opposite is true for real proteins (Makhatadze and Privalov, 1995). The low value for the random-coil specific heat is due to the absence of protein-solvent interactions (Privalov, 1992) in this model, which diminishes the size of fluctuations in the energy in the unfolded state and results in a more stable energy state (Zhou et al., 1999a).

The results for the -sheet that we have just described are summarized in Fig. 3, which shows the phases that occur in the space spanned by the reduced temperature T* and the bias gap g. The diagram shows a complex set of protein transitions that is qualitatively similar to those observed for a model three-helix bundle protein (Zhou and Karplus, 1997a, 1997b, 1999) and for the protein-like lattice models (Dinner et al., 1994; Doniach et al., 1996; Pande and Rokhsar, 1998). However, some quantitative differences between the -sheet results and the three-helix bundle results are immediately apparent. This is indicated in the figure by drawing solid lines for transitions that are qualitatively the same transitions as those observed for the three-helix bundle protein and dotted lines for transitions which are not observed in the three-helix bundle protein. The upper line in the figure indicates the collapse transition; it separates the random coil phase from the disordered globule phase. Moving down in the figure, the second line from the top indicates the transition from the disordered globule to the ordered globule; it is obtained from the first peak in the specific heat. The third line from the top is a folding transition that separates the ordered globule from the native state for g > 0.7, and a collapse transition to a partially ordered globule (different from the native state) for 0.3 < g < 0.7. As described later in this paper the partially ordered globule state contains kinked or rolled structures that are distinctly different from their native state conformation. Interestingly, for 0.7 < g < 0.9, we observe a transition from the native state to the partially ordered globule state that is not observed in the three-helix bundle protein. This transition is associated with the specific heat peak for g = 0.7 at T* = 0.35 and the small specific heat peak for g = 0.9 at T* = 0.115. This transition and the partially ordered globule phase are not observed for real proteins. Because they occur at small values of the bias gap, they can be regarded as remnants of the model's homopolymer behavior (as g  0, the model reduces to a homopolymer) in which a transition to a high-density configuration occurs at low temperatures (Zhou et al., 1996, 1997). For g < 0.9 a transition occurs at low temperature to a solid phase; this is obtained through application of the Lindemann criterion as will be discussed later. No solid phase was observed for g > 0.9 at any temperature investigated. Two triple points are found: one at g = 0.3 and T* = 0.58 where the disordered, ordered, and partially-ordered globule phases meet, and one at g = 0.7 and T* = 0.35 where the native, ordered, and partially-ordered globule phases meet.

View larger version (21K): [in this window] [in a new window]   FIGURE 3   Phase diagram for the -sheet model as a function of the reduced temperature, T*, and the bias gap, g.

Folding of the -sheet strongly depends on the bias gap and temperature. A set of sample structures for the -sheet characteristic of the final states observed during the simulations is presented in Fig. 4 a for g = 0.7 at T* = 0.3, Fig. 4 b for g = 0.7 at T* = 0.6, Fig. 4 c for g = 1.3 at T* = 0.5, and Fig. 4 d for g = 1.3 at T* = 0.2. The structure in Fig. 4 a corresponds to the partially ordered globule, the kinked or rolled structure that appears at low temperature for g < 0.7. (Recall that no native state structures appear at any temperatures for g < 0.7.) The partially ordered globule state occurs for smaller gap models because the attraction between nonnative contacts disturbs the tendency to order into the native state at low temperatures. The structures in Fig. 4, b and c, are examples of the ordered globule structure. The ordered globule has more native contacts than the disordered globule. As g increases, the structure corresponding to the ordered globule becomes more native-like (compare Fig. 4, b and c). The structure in Fig. 4 d is the native state structure of the -sheet. Note that for the off-lattice -sheet models with g  0.7 none of the structures correspond to the global energy minimal state, but for 0.7 < g < 1.0 the geometry in the native state is very similar to the global energy minimal structure shown in Fig. 1 a so that the number of native contacts is very close to that of the global energy minimal structure. For highly optimized models, g  1.0, the native state is the global energy minimal state (Zhou and Karplus, 1999).

View larger version (36K): [in this window] [in a new window]   FIGURE 4   Final chain conformations of the -sheet (a) for g = 0.7 at T* = 0.3 (the partially ordered globule state), (b) for g = 0.7 at T* = 0.6 (the ordered globule state), (c) for g = 1.3 at T* = 0.5 (the ordered globule state), and (d) for g = 1.3 at T* = 0.2 (the native state).

It is interesting to consider the variation in structural properties as folding progress. Fig. 5 shows the average value of the fraction of global energy minimal contacts formed (Lazaridis and Karplus, 1997; ali et al., 1994), Q, as a function of temperature T * and bias gap g. The collapse transition occurs at Q ~ 0.2 for g = 0.3 and at Q ~ 0.3 for g = 1.3. The transition from the disordered globule to the ordered globule occurs at Q values in the range of Q = 0.4 ~ 0.5 for all values of g. In the native state, we find that Q > 0.9 for g > 0.7, whereas for the partially ordered globule, Q has values in the range from 0.6 to 0.9 depending on the size of g. In the three-dimensional graph, the flat surface at the top of the rear edge corresponds to the native state with Q ~ 1 for g > 0.7. 

View larger version (58K): [in this window] [in a new window]   FIGURE 5   Average value of the fraction of global energy minimal contacts, Q, for the -sheet as a function of the reduced temperature, T *, and the bias gap, g.

Further information on the nature of the transition at low temperature can be obtained by calculating the Lindemann disorder parameter (Lindemann, 1910). This parameter is often used to analyze the liquid-to-solid transition (Bilgram, 1987; Löwen, 1994; Stillinger, 1995; Zhou et al., 1999b). The solid phase corresponds to the inactive glassy state that can be observed in many proteins (Ferrand et al., 1993; Frauenfelder and McMahon, 1998; Green et al., 1994; Reat et al., 1998; Tilton et al., 1992). Fig. 6 a shows the average value of the Lindemann disorder parameter, _L, as a function of temperature for two selected values of the bias gap, g = 0.3 and 1.3. Each result for _L at given temperature is a sum of the average of the root-mean-squared fluctuations over all beads. A form of the Lindemann criterion for melting is adopted here in which systems with _L < 0.15 are considered solid-like, whereas systems with _L > 0.15 are considered liquid. The dotted line in Fig. 6 a indicates the criterion for the liquid-to-solid transition. It can be seen that for g = 1.3, as was seen in the phase diagram, no solid phase is found at any investigated temperatures because _L > 0.15, whereas, for g = 0.3 the liquid-to-solid transition occurs at T* = 0.19. This agrees well with the temperature at which the inflection in the root-mean-squared pair separation, _B, curve is found in Fig. 2 a.

View larger version (13K): [in this window] [in a new window]   FIGURE 6   (a) Lindemann disorder parameter, L, for the -sheet as a function of the reduced temperature, T*, for two selected values of the bias gap, g = 0.3 and 1.3. (b) Lindemann disorder parameter for an individual bead, Li, for the -sheet as a function of the reduced distance of bead i from the origin, |ri|/, for g = 0.3 at T* = 0.1, 0.3, and 0.6. (c) Same as (b) but for g = 1.3 at T* = 0.2, 0.6, and 0.8.

Fig. 6, b and c, shows the Lindemann disorder parameter for the individual beads i, _Li, as a function of the reduced distance of bead i from the origin, |r_i|/, for Fig. 6 b, g = 0.3 at T* = 0.1, 0.3, and 0.6, and Fig. 6 c, for g = 1.3 at T* = 0.2, 0.6, and 0.8. Note that the bead positions, r_i, are defined with respect to the origin at the chain's center of mass. In Fig. 6 b, which corresponds to g = 0.3, the three selected temperatures correspond to the three different phases of the chain: the solid-like phase, partially ordered globule, and disordered globule, respectively. At high temperature, T* = 0.6, a large degree of freedom in the bead motion is observed. At low temperature, T* = 0.1, the small bead fluctuations and small bead-bead separations give rise to a solid-like phase for the chain. In Fig. 6 c, which corresponds to g = 1.3, the three selected temperatures correspond to the three different phases of the chain: the native state, ordered globule, and disordered globule, respectively. At high temperature, T* = 0.8, the large bead fluctuations in the _Li graph indicate that the chain undergoes a transition from the disordered globule to the ordered globule. At low temperature, T* = 0.2, the results of the individual bead fluctuations correspond to the native state. The native protein is regarded as a surface-molten solid, which means that the interior of native protein is solid-like, but the surface is liquid-like phase (Zhou et al., 1999b). However we find that the native -sheet at T* = 0.2 is not a surface-molten solid, because all the values of _Li are greater than 0.15, which means that the system is considered liquid. Therefore, the native -sheet is more like a molten globule because all of the beads exhibit liquid-like motion. The reason that the native -sheet is not a surface-molten solid is that it is planar, and consequently core beads have considerable freedom in their motion; they can move perpendicular to the molecular plane even in the native state. This can be regarded as a case in which the topology determines the phase behavior.

It is of interest to determine the extent to which this model displays thermodynamic two-state cooperativity (Ptitsyn, 1995). The most direct way to do this is to examine the shape of the free energy distribution in the transition region to see if it is bimodal. Accordingly, the energy distributions in the vicinity of transitions displayed in Fig. 3 have been investigated. Fig. 7 shows the energy distribution in Fig. 7 a for g = 0.9 at T* = 0.7949 (the disordered-to-ordered globule transition), in Fig. 7 b for g = 0.9 at T* = 0.3998 (the folding transition), and in Fig. 7 c for g = 0.3 at T* = 0.5848 (the disordered-to-partially-ordered globule transition). Fig. 7 a seems to indicate that the g = 0.9 model has a bimodal free energy distribution. As pointed out by Zhou and Karplus (1999), the logarithm of the energy distribution as calculated via the weighted histogram method is directly proportional to the free energy, and the existence of a bimodal shape in the free energy distribution indicates the presence of a two-state transition (Zhou et al., 1996, 1997). Nevertheless, to further verify the existence of a two-state transition, it is useful to also apply the calorimetric criterion for chain models of protein. The criterion says that to have two-state cooperativity in protein folding, H_vH/H_cal  1 around the transition, where H_vH is the van't Hoff enthalpy around the peak of the specific heat after baseline subtraction and H_cal is the calorimetric enthalpy change associated with the entire transition (Jackson and Brandts, 1970; Freire, 1995). We evaluate the criterion for the collapse transition of the large bias gap models using two different equations for the van't Hoff enthalpy, H_vH(T) = 2T (Chan, 2000; Kaya and Chan, 2000) and H(T) = 4k_BT²C_p(T)/H_cal (Zhou et al., 1999a). For the transition in Fig. 7 a, we find that H_vH/H_cal ~ 0.48 and H/H_cal ~ 0.23 at the midpoint temperature. These ratios are significantly smaller than unity, indicating that either the collapse transition is not calorimetrically two-state or the baseline for the specific heat curve is not accurate (Zhou et al., 1999a). The ratios increase with the size of the bias gap, indicating that the collapse transition becomes more cooperative, but none of them becomes equal to or close to unity for the -sheet models we investigated. We suspect that for highly optimized models, i.e., even larger bias gaps, the transition will become more protein-like, displaying the two-state "all-or-none" character associated with temperature-induced protein denaturation.

View larger version (15K): [in this window] [in a new window]   FIGURE 7   Energy distribution in the transition region for the -sheet as a function of energy (a) for g = 0.9 at T* = 0.7949 (the disordered-to-ordered globule transition), (b) for g = 0.9 at T* = 0.3998 (the folding transition), and (c) for g = 0.3 at T* = 0.5848 (the disordered-to-partially-ordered globule transition).

The folding transition in Fig. 7 b for g = 0.9 is not two-state because no bimodal energy distribution occurs. This indicates that the ordered globule and native states of the -sheet are geometrically similar to each other. Hence, the folding transition is weak with small energy change and is continuous between the two states. Note that this transition is associated with the plateaus in the specific heat as seen in Fig. 2, c and d. For the transition from the disordered globule state to the partially ordered globule state in Fig. 7 c at g = 0.3, the transition is continuous even though it is associated with a distinct peak in the specific heat as seen in Fig. 2 a. The lack of a bimodal energy distribution for g = 0.3 suggests that the disordered and partially ordered globule structures are similar to each other.

Thermodynamic properties of -clip and -twist

The thermodynamic averages, R_g²/²N, _B, C_v*/N, and E*/N for the -clip and -twist models, are shown as a function of T* in Fig. 8, a and b (-clip) and Fig. 8, c and d (-twist). For brevity the figures only show results for two selected values of the bias gap, g = 0.3 and 1.3. For the small bias gap, g = 0.3, the results in Fig. 8 a for the -clip and Fig. 8 c for the -twist are qualitatively similar to those in Fig. 2 a for the -sheet model at g = 0.3. However, for the large bias gap, g = 1.3, the thermodynamic averages for both the -clip model and the -twist model (Fig. 8, b and d, respectively) are qualitatively different from those of the -sheet model at g = 1.3 in Fig. 2 d. The differences are that kinks in the R_g²/²N and _B curves at low temperatures and a second peak in the specific heat in Fig. 8, b and d are observed for both the -clip model and the -twist model, but not for the -sheet model in Fig. 2 d. The small peak in the specific heat is related to the folding transition, whereas the kinks in the R_g²/²N and _B curves at low temperatures are connected with the liquid-to-solid transition. The reason for the increases in R_g²/²N at low temperatures for both the -clip model and the -twist model at g = 1.3 is that the R_g² value of the native protein is larger than that of the ordered globule, indicating that the collapsed chain undergoes a subtle conformational change from the ordered globule to the native state by elongating its strand.

View larger version (19K): [in this window] [in a new window]   FIGURE 8   Average values of the reduced squared radius of gyration per bead, Rg2/2N, the root-mean-squared pair separation fluctuation, B, the reduced specific heat per bead, Cv*/N, and the reduced internal energy per bead, E*/N, as a function of the reduced temperature, T*, for the -clip model at (a) g = 0.3 and (b) g = 1.3, and for the -twist model at (c) g = 0.3 and (d) g = 1.3.

The Lindemann disorder parameters are shown as a function of temperature and the bias gap for the -clip and -twist in Fig. 9, a and b, respectively. It can be seen that for both models the locations of the kinks in the _L curves at low temperatures and large g are consistent with those at low temperatures found in the R_g²/²N and _B curves as seen in Fig. 8, b and d. The solid-like phase with _L < 0.15 is prevalent at low temperature for the -clip and -twist models. Sharp increases at high temperatures in the _L curves for the large gap models are also consistent with those found in the R_g²/²N curves and with the peaks in the _B curves in Fig. 8, b and d, indicating the collapse transition.

View larger version (74K): [in this window] [in a new window]   FIGURE 9   Lindemann disorder parameter, L, as a function of the reduced temperature, T*, and the bias gap, g, for (a) the -clip and (b) the -twist.

Fig. 10, a and b, show the Lindemann disorder parameter for an individual bead, _Li, a for the -clip model with g = 0.9 at T* = 0.1, 0.26, and 0.5, and b for the -twist model with g = 0.7 at T* = 0.1, 0.28, and 0.6. For both the -clip and -twist models, the three selected temperatures correspond to the three different phases of the chain: the solid-like phase, the native state, and the ordered globule, respectively. In the native states near the liquid-to-solid transition (-clip with g = 0.9 at T* = 0.26 in Fig. 10 a and -twist with g = 0.7 at T* = 0.28 in Fig. 10 b), we find that the interiors of the chains are solid-like with _Li < 0.15, but the surfaces are liquid-like with _Li > 0.15. Thus the native proteins for the -clip and -twist are surface-molten solids. At low temperature, T* = 0.1, for both the -clip and -twist models, the bead motions are frozen with all _Li < 0.1, indicating an inactive solid-like phase.

View larger version (15K): [in this window] [in a new window]   FIGURE 10   Lindemann disorder parameter for an individual bead, Li, as a function the reduced distance of bead i from the origin, |ri|/, for (a) the -clip with g = 0.9 at T* = 0.1, 0.26, 0.5, and (b) the -twist with g = 0.7 at T* = 0.1, 0.28, 0.6.

Fig. 11 a for the -clip and Fig. 11 b for the -twist summarize the protein phases that occur in the space spanned by the reduced temperature T* and the bias gap g. The qualitative differences between the phase behavior results for the -clip and -twist models and those for the -sheet model (Fig. 3) are immediately apparent. For example, the native states for the -clip and -twist models exist at low temperatures for g > 0.5, whereas the native state for the -sheet exists at low temperatures for g > 0.7. For the -clip and -twist models at even lower temperatures, the liquid-to-solid transition exists for all values of g, but for the -sheet it exists only for g < 0.9. For the -clip and -twist models at g  0.9, although not shown in the figures, the energy distribution for the collapse transition is bimodal, indicating that the transition is two-state-like. However, the calorimetric criterion shows that at g = 1.3, H_vH/H_cal ~ 0.8 for the -clip and -twist models, a ratio that is still smaller than the value of H_vH/H_cal ~ 0.96 observed in experimental studies on small proteins (Privalov, 1979). This indicates that the collapse transition is only weakly cooperative.

View larger version (18K): [in this window] [in a new window]   FIGURE 11   Phase diagram as a function of the reduced temperature, T*, and the bias gap, g, for (a) the -clip model and (b) the -twist model.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We have studied the folding thermodynamics for three different types of off-lattice four-strand antiparallel -strand protein models: the -sheet, the -clip, and the -twist. DMD simulations of these models have been performed for different sizes of the bias gap g, a measure of the strength of the native contacts relative to that of the nonnative contacts. Although all three -strand model proteins are simple, they undergo a complex set of protein transitions as observed in experimental studies on real proteins (Ptitsyn, 1995