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 Jang, H. Articles by Zhou, Y. Search for Related Content PubMed PubMed Citation Articles by Jang, H. Articles by Zhou, Y.

Assembly and Kinetic Folding Pathways of a Tetrameric ß-Sheet Complex: Molecular Dynamics Simulations on Simplified Off-Lattice Protein Models

Hyunbum Jang ^*, Carol K. Hall ^* and Yaoqi Zhou 

^* Department of Chemical Engineering, North Carolina State University, Raleigh, North Carolina 27695-7905; and Department of Physiology and Biophysics, State University of New York at Buffalo, Buffalo, New York 14214

Correspondence: Address reprint requests to Carol K. Hall, E-mail: hall{at}turbo.che.ncsu.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION Models and simulation method RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
We have performed discontinuous molecular dynamic simulations of the assembly and folding kinetics of a tetrameric ß-sheet complex that contains four identical four-stranded antiparallel ß-sheet peptides. The potential used in the simulation is a hybrid Go-type potential characterized by the bias gap parameter g, an artificial measure of a model protein's preference for its native state, and the intermolecular contact parameter , which measures the ratio of intermolecular to intramolecular native attractions. The formation of the ß-sheet complex and its equilibrium properties strongly depend on the size of the intermolecular contact parameter . The ordered ß-sheet complex in the folded state and nonaligned ß-sheets or tangled chains in the misfolded state are distinguished by measuring the squared radius of gyration and the fraction of native contacts Q. The folding yield for the folded state is high at intermediate values of , but is low at both small and large values of . The folded state at small is liquid-like, but is solid-like at both intermediate and large . The misfolded state at small contains nonaligned ß-sheets and tangled chains with poor secondary structure at large . Various folding pathways via dimeric and trimeric intermediates are observed, depending on . Comparison with experimental results on protein aggregation indicates that intermediate values are most appropriate for modeling fibril formation and small values are most appropriate for modeling the formation of amorphous aggregates.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION Models and simulation method RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
We have been engaged in a computational study aimed at understanding the basic physical principles underlying protein aggregation, a cause or associated symptom of a number of fatal protein deposition diseases (Clark and Steele, 1992; Eaton and Hofrichter, 1990; Gallo et al., 1996; Massry and Glasscock, 1983; Moore and Melton, 1997). An example is Alzheimer's disease (Selkoe, 1991; Simmons et al., 1994), which is thought to be caused by the abnormal deposition of insoluble fibril (amyloid) plaques in the extracellular space of brain tissue during aging. The major component of these plaques is a 39- to 43-residue ß-peptide, called ß-amyloid peptide (Aß), whose predominant secondary structure in the fibril is a ß-sheet. Although extensive experimental studies (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) have been conducted on the fibrillization of ß-amyloid peptide (Aß) and other amyloid forming peptides, the causes of amyloid fibril formation are largely unknown and the fundamental mechanisms underlying this aggregation are still under investigation. Recent experimental observations that many different proteins are capable of forming fibrils under slightly denatured concentrated conditions have led leaders in the field to suggest that fibril formation is an intrinsic property of polypeptides (Booth et al., 1997; Chiti et al., 1999). The idea here is that fibril formation is not strongly dependent upon sequence but is instead a consequence of the ability of peptides to experience hydrogen bonding and hydrophobic interactions. This suggests that computer simulations of multipeptide systems that take these types of interactions into account could offer insights into the basic principles behind fibril formation in peptides.

The long-term goal of our work is to reveal the molecular-level mechanisms underlying protein aggregation and fibril formation with particular focus on the role played by ß-strands. Because protein aggregation often involves ß-strand association, 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 a model fibril made up of ß-sheets. To study protein aggregation by computer simulation, one must devise protein models that capture not only the basic physical features of real proteins, but also allow the simulation of many proteins with current computer capability. Our approach is to use low-resolution or simplified protein models; such models are best suited for simulating multiprotein systems because they allow us to study folding behavior over relatively long timescales. Investigations based on these models have already provided numerous insights into the thermodynamic and kinetic properties of protein folding for isolated chains (Chan and Dill, 1994; Dokholyan et al., 1998, 2000; Guo and Brooks, 1997; Guo and Thirumalai, 1995, 1996; Gupta and Hall, 1997; Gupta et al., 1999; Jang et al., 2002a,b; Kolinski et al., 1995, 1999; Lau and Dill, 1989; Miller et al., 1992; Nymeyer et al., 1998; Pande and Rokhsar, 1998; Shea et al., 2000; Skolnick and Kolinski, 1991; Zhou and Karplus, 1997a,b; 1999) and of protein aggregation for multichain systems (Bratko and Blanch, 2001; Dill and Stigter, 1995; Dima and Thirumalai, 2002; Gupta and Hall, 1998; Harrison et al., 1999, 2001).

We recently introduced three minimalist models of four-strand antiparallel ß-strand peptides: the ß-sheet, the ß-clip, and the ß-twist (Jang et al., 2002a,b). In our study of the folding thermodynamics of these three models (Jang et al., 2002a), discontinuous molecular dynamics (DMD) simulations (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996) were performed to determine how the thermodynamic properties of an isolated peptide vary with temperature. Despite these models' simplicity, they undergo a complex set of protein transitions similar to those observed in experimental studies on real proteins (Ptitsyn, 1995). Starting from high temperature, these transitions include a collapse transition, a disordered-to-ordered globule transition, a folding transition, and a liquid-to-solid transition. In our study of the folding kinetics for the same models (Jang et al., 2002b), the ß-sheet exhibits a fast-track folding pathway without becoming trapped in any intermediate. In contrast, the ß-clip and ß-twist exhibit multiple folding pathways that include trapping in intermediates and direct folding to the native state. The folding speed of the model proteins with different native state topologies strongly depends on the contact order in the native state (Plaxco et al., 1998). The thermodynamics and kinetics results presented in our previous papers set the stage for this work in which we examine the folding and assembly of several of the model ß-sheet proteins into fibrils.

In this article, we investigate the equilibrium properties and the folding kinetics of a tetrameric ß-sheet complex. The ß-sheet complex consists of four identical four-stranded antiparallel ß-sheet peptides, with each peptide containing 39 connected residues (beads). Nonbonded beads interact both intramolecularly and intermolecularly 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 values of the bias gap parameter g and the intermolecular contact parameter . The bias gap g measures the difference in interaction strength between the intramolecular 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 real proteins in their equilibrium and dynamic behavior. The intermolecular contact parameter measures the ratio of intermolecular to intramolecular native attractions.

Our model of the ß-sheet complex was designed in part to mimic the small Aß oligomers that are observed (albeit indirectly) in the early stages of fibril formation (Harper et al., 1997). These oligomers are now widely believed to serve as the nuclei that seed the growth of the fibrils that characterize Alzheimers's and other amyloid diseases (Pallitto and Murphy, 2001). In fact recent evidence suggests that it is these oligomers, rather than the fully formed fibrils, that are the toxic species in Alzheimers's disease (Kirkitadze et al., 2002). Thus it is of interest to understand the mechanisms associated with oligomer complex nucleation. Our model was designed in part to mimic real Aß oligomer complex formation in several respects. First, our peptides are 39 residues long and Aß is between 39 and 42 residues long. Second, Aß has several extended stretches of hydrophobic residues (Aß(17-21), Aß(32-42)); our peptide has essentially four stretches of residues that act with an external hydrophobic interaction. Third, the Aß peptides in Aß fibrils (and indeed the peptides in all fibrils) experience intramolecular hydrogen bonding within the ß-sheets and intermolecular hydrophobic interactions between the ß-sheets. The model peptides in our ß-sheet complex experience intramolecular interactions within the sheets that are reminiscent of hydrogen bonding; these interactions are characterized by the Go-model bias gap parameter, g. They also experience intermolecular interactions between the sheets that are reminiscent of hydrophobic interactions; these interactions are characterized by the intermolecular contact parameter, . The Go-model bias gap parameter is set to an intermediate value (g = 0.9) because Go models are considered to be the most realistic at these parameter values. The intermolecular contact parameter, , which essentially measures the ratio of the intermolecular (hydrophobic) interaction and the intramolecular (hydrogen bonding) interaction, is selected in the range between 0.2 1. The reasons for considering a range of intermolecular contact parameter values are: 1), the best values for the relative strengths of the hydrophobic and hydrogen bonding interactions in simplified models are a matter of debate (Pace et al., 1998), and 2), intermolecular hydrophobic interaction strengths depend upon the intrinsic conditions of the protein (sequence and number of hydrophobic residues) and on external conditions (pH and temperature) (Fraser et al., 1991a,b; Kowalewski and Holtzman, 1999; Snyder et al., 1994). Furthermore, by exploring how the aggregation mechanism varies with the intermolecular contact parameter we provide insights to other research workers interested in modeling aggregation phenomena.

Discontinuous molecular dynamics simulations (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996) were performed on the model systems at an intermediate value of the bias gap (g = 0.9) for different intermolecular contact parameters (0.2 1). At each value of , at least 50 independent trajectories starting from different initial configurations were monitored and used in the averaging. Four separated random coils generated at high temperature were used as the initial configuration in the simulations. The equilibrium properties of folded and misfolded ß-sheet complexes were investigated. Thermodynamic averages were taken by accumulating the last quarter of the simulation data because the ß-sheets assembled into the complex structure early in the simulation and were stable once they were assembled. The kinetics of fibril formation was investigated by monitoring the progress parameters measuring the complex formation, including the squared radius of gyration of the entire system, the squared radius of gyration for the nth chain, the fraction of total native contacts, Q_total, the fraction of intramolecular native contacts, Q_intra, the fraction of intermolecular native contacts, Q_inter, and the fraction of intermolecular native contacts for the nth chain,

Highlights of our results are the following. We find that the equilibrium properties and the kinetic formation of the ß-sheet complex strongly depend on the size of the intermolecular contact parameter . For small , the folded state is an ordered ß-sheet complex in a liquid-like phase. Most folding trajectories follow the path: four monomers dimer and two monomers trimer and monomer or two dimers tetramer. The folding yield is low. Secondary structure begins to form before the ß-sheet complex assembles. Many trajectories fold into a misfolded state that contains two well-aligned eight-stranded ß-sheets, albeit a highly ordered structure in terms of fibril growth. For intermediate , the folded state is a highly ordered ß-sheet complex whose interior ß-sheets are solid-like, and whose outer ß-sheets have both solid-like and liquid-like characteristics, with the entire structure being essentially solid-like. Most folding trajectories follow the path: four monomers dimer and two monomers trimer and monomer tetramer. The folding yield is very high. For large , the folded state is a highly ordered ß-sheet complex in a solid-like phase that is physically almost the same as the global energy minimum structure. Most folding trajectories follow the path: four monomers dimer and two monomers two dimers tetramer. The folding yield is very low. The misfolded state contains tangled chains with poor secondary structure and is observed in a trajectory with early trimer formation. Comparison with experimental results on protein aggregation indicates that intermediate values (where the majority of folding trajectories lead to the highly ordered ß-sheet complex) are most appropriate for modeling fibril formation, whereas small values (where the majority of folding trajectories lead to the misfolded state) are most appropriate for modeling the formation of amorphous aggregates.

In the following section, a full description of the model and the simulation method are presented. Next is the Results section, which is divided into seven subsections; the first presents the results for the equilibrium properties of the ß-sheet complex, the second presents the results for the folding kinetics, the third through fifth present the results for the kinetics of the formation of the ß-sheet complex at = 0.2, 0.5, and 0.8, respectively, the sixth describes the misfolded states, and the seventh discusses the folding pathways for = 0.2 and 0.5. The article concludes with a discussion of the key findings.

	Models and simulation method

TOP ABSTRACT INTRODUCTION Models and simulation method RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
We consider an off-lattice protein model of a ß-sheet complex that consists of four identical four-stranded antiparallel ß-sheet peptides. The global energy minimum structure for the model ß-sheet complex is shown in Fig. 1. Each color represents a different ß-sheet chain: chain 1 (blue), chain 2 (red), chain 3 (green), and chain 4 (yellow). A topology diagram representing a top view of the complex is shown at the bottom of the figure. In the topology diagrams, thick lines indicate chain connectivity at the top of the strands, and thin lines indicate chain connectivity at the bottom of the strands. It can be seen from the topology diagram that the ß-sheet is not perfectly planar because it has an inflection in the middle. The ß-sheet complex consisting of the kinked ß-sheets has more intermolecular native contacts than that consisting of planar ß-sheets. In the global energy minimum state, the total number of intramolecular native contacts, the total number of intermolecular native contacts, and the reduced squared radius of gyration for the entire system, for the ß-sheet complex are summarized in Table 1. For a monomer ß-sheet in the global energy minimum state, the table also shows the number of intramolecular native contacts for the nth chain, where n is the chain number (n = 1-4), the number of intermolecular native contacts for the nth chain, and the reduced squared radius of gyration for the nth chain, Note that but because the intermolecular native contacts for some of the chains are shared, i.e., the interior ß-sheets (red and green) have twice the number of intermolecular native contacts as the outer ß-sheets (blue and yellow) in the global energy minimum state.

View larger version (28K): [in this window] [in a new window]   FIGURE 1  Bead (a) and topology (b) diagrams of the global energy minimum structures for the model ß-sheet complex.

View this table: [in this window] [in a new window]   TABLE 1  Total number of intramolecular, and intermolecular, native contacts and the reduced squared radius of gyration, for a ß-sheet complex in the global energy minimum state

(1)

where 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)

where B_O, and are measures of the relative strengths of the energies associated with the intramolecular native, nonnative, and intermolecular native pair interactions in this Go-type potential (Go and Taketomi, 1978, 1979; Taketomi et al., 1975; Ueda et al., 1978). Within a chain, nonbonded pairs of beads that are in contact in the global energy minimum structure experience an attractive interaction, i.e., when their square-wells overlap. On the other hand, within a chain and in different chains, nonbonded pairs of beads that are not in contact in the global energy minimum 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)

where g is the bias gap (Jang et al., 2002a,b; Zhou and Karplus, 1997a,b; 1999). The bias gap measures the ratio of the interaction strength between the intramolecular 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 the intramolecular native contacts are always more favorable than the nonnative contacts. For g = 0, the intramolecular native and nonnative contacts are equally favorable, ; 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. In different chains, bead pairs that are in contact in the global energy minimum structure experience an attractive interaction, i.e., when their square-wells overlap. The parameter measures the interaction strength between beads in different chains and depends on the size of the intermolecular contact parameter ,

(4)

where is the intermolecular contact parameter. The quantity measures the ratio of the interaction strength between the intermolecular and intramolecular native contacts. For = 1, the intramolecular and intermolecular native contacts are equally favorable, i.e., For < 1, the intramolecular native contacts are more favorable than the intermolecular native contacts, i.e., However, nonnative contacts are always less favorable than the intramolecular and intermolecular native contacts, so that for g > 0. In this model, the intramolecular native interaction, intermolecular native interaction, and nonnative interaction might be regarded as mimicking hydrogen bonding, hydrophobic interactions, and van der Waals interactions, respectively, because in real fibrils, the molecular interactions along the sheets (generally hydrogen bonds) are different from the molecular interactions between the sheets (generally hydrophobic interactions).

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

(5)

where 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.

The assembly and folding kinetics of the ß-sheet complex are simulated using the discontinuous molecular dynamics (DMD) algorithm (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996). This algorithm is very efficient for calculating the properties of systems containing hard or square-well spheres or chains (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. During any collision process, the total energy of the system is conserved, so that the change in kinetic energy is equivalent to the negative of the change of potential energy. The DMD simulations are conducted in the canonical ensemble with constant number of particles, volume, and temperature; periodic boundary conditions are applied in all directions. The length of the cubic box is twice the length of a fully extended chain, preventing that chain from interacting with more than one image of another chain at one time. To maintain constant temperature, the Andersen stochastic collision method (Anderson, 1980) is used. In this method the system's particles collide with imaginary or ghost particles that 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 where n_g is the number density of ghost particles, was used. As a result, 2% of the collisions were ghost particle collisions. Because the value of the ghost particle density can alter the kinetic results, we use the same value of the ghost particle density for all kinetic simulations.

The system is quenched to the temperature of interest from the denatured or random coil state equilibrated at a high temperature. The initial configurations for the kinetic simulations were generated from equilibrium simulations of the random coil state at T^* = 2.0. The temperature used in this paper is a dimensionless reduced temperature, T^* k_BT/, where k_B is Bolzmann's constant, T is the absolute temperature, and is the energy parameter. The coordinates and velocities of the beads in the random coil state were obtained every 10⁵ collisions during the equilibrium simulation to create independent initial conformations for the quenching studies. This large collision interval ensures that the sampled configurations are independent. After quenching, the four random coils undergo a conformational change toward an energy minimum, i.e., the highly ordered ß-sheet complex (native state) or a local energy minimum, i.e., a misfolded state. All chains are identical, so that any chain can fold into any location within the ß-sheet complex. We investigated the time series of folding behavior from the random coil state that was generated at T^* = 2.0 to the native state, keeping the temperature constant at T^* = 0.2. Our previous thermodynamic study on the isolated ß-sheet models (Jang et al., 2002a) showed that at T^* = 2.0 and T^* = 0.2, the ß-sheet peptide is in a random coil state and the native state, respectively.

Simulations were performed for intermolecular contact parameters in the range of = 0.2 to 1.0 at a bias gap g = 0.9. The intermediate value of the bias gap g = 0.9 is thought to be the most representative of a real protein in its equilibrium and dynamic behavior (Jang et al., 2002a,b). At this bias gap, all nonbonded pairs of beads are attracted with square-well depths B_O = -0.1, and for the intramolecular native contact, nonnative contact, and intermolecular native contact pairs, respectively. For each value of , a number of quenching simulations were performed over at least 50 independent trajectories that start from different initial configurations. The total time elapsed during a simulation of a given number of collisions strongly depended on the size of the intermolecular contact parameter . The time used in the simulations is expressed in terms of reduced time (Zhou and Karplus, 1999), where t is the real time, is the energy parameter, m is the mass of the bead, and is the bead diameter. For all values of , the simulations were performed for 4 x 10⁸ collisions, which is approximately equivalent to t^* = 3.4 x 10⁵, 1.9 x 10⁵, and 1.6 x 10⁵ for = 0.2, 0.5, and 0.8, respectively. This indicates that the simulation speed for large is faster than for small .

The equilibrium properties of the ß-sheet complex were also investigated once oligomerization was achieved. Oligomers, whether folded or misfolded, were typically seen after 1 x 10⁸ collisions for all values of , indicating that the ß-sheets assembled into the oligomer early in the simulation and were stable once they were assembled. Equilibrium averages were taken by accumulating the simulation data over the last quarter of the collisions (i.e., after 3 x 10⁸ collisions) to ensure system equilibration. Sampling was taken over the folding or misfolding trajectories separately. For the folding trajectories, the sampled equilibrium ensembles are well within a Gaussian distribution with a narrow width, because the folded oligomers are structurally similar to each other. For the misfolding trajectories, the sampled equilibrium ensembles are also well within a Gaussian distribution with a wide width, because the misfolded oligomers are structurally slightly different from each other but the majority of the misfolding trajectories end up in some common misfolded state.

The progression of a conformation toward the native state was monitored by introducing the time-dependent fraction of total native contacts formed (Lazaridis and Karplus, 1997; ali et al., 1994), Q_total(t), defined by

(6)

where N_total(t) is the total number of native contacts at a given time, with N_total(t) N_intra(t) + N_inter(t), and is the total number of native contacts in the global energy minimum state, with The fraction of total native contacts is a weighted sum of the fraction of intramolecular native contacts, Q_intra(t), and the fraction of intermolecular native contacts, Q_inter(t), by

(7)

where and is a constant, The assembly of the complex was monitored by introducing the time-dependent fraction of intermolecular native contacts for the nth chain, defined by

(8)

where is the number of intermolecular native contacts for the nth chain at a given time. The value used for the denominator in Eq. 8 is that for the interior ß-sheets at n = 2 or 3, because the ß-sheet has a maximum number of intermolecular native contacts when it is located in the interior of the ß-sheet complex (see Table 1).

To determine the separation between the chains before they assemble, the size of the entire system after oligomerization, and the compactness of the individual chains during the simulations, the time-dependent squared radius of gyration for the entire system, and for the nth chain, were monitored where

(9)

and

(10)

with M_n equal to the number of beads for each ß-sheet monomer. The time-dependent bead positions, r_i(t) (i = 1 to M), were obtained by shifting the center of mass coordinates of the whole system to the origin. Likewise, the time-dependent bead positions for the nth chain, r_nk(t) (k = 1 to M_n), were obtained by shifting the center of mass coordinates of each chain to the origin. To determine the degree of secondary structure formation in the ß-sheet complex, the time-dependent squared radius of gyration for the nth chain, is averaged over the N chains, where N is the total number of the chains, i.e., N = 4, to give the mean squared radius of gyration for the individual chains, defined by

(11)

We have studied the equilibrium properties of the ß-sheet complex after oligomerization. The equilibrium averages were determined by accumulating the last one-quarter of the quenching simulation data, because a ß-sheet complex occurs early in the simulation and is stable once it is assembled. The reduced internal energy can be calculated from

(12)

where _cf denotes the configurational average over the last one-quarter of the total collisions. The degree of bead mobility can be characterized by the root-mean-squared (rms) pair separation fluctuation, _B, defined by

(13)

where r_ij denotes the separation between beads i and j in the ß-sheet complex. To characterize the phase of the ß-sheet complex, the Lindemann disorder parameter for bead i, _Li, is defined to be the root-mean-squared fluctuation for bead i,

(14)

and the Lindemann disorder parameter, _L, is defined to be

(15)

where r_i is the position of bead i (i = 1 to M) and is the mean-squared fluctuation. Because the entire structure of the assembled ß-sheet complex continuously translates and rotates over space even in equilibrium, this causes a large value of the root-mean-squared deviation. In the calculation of the mean-squared fluctuation, the translations are removed by shifting the system's center of mass coordinates to the origin. The rotations are removed by rotating the system so as to minimize the root-mean-squared deviation of the new coordinates from the coordinates set at the beginning of the calculation (Kabsch, 1976). The Lindemann disorder parameter is often used to analyze the liquid-to-solid transition. A form of the Lindemann criterion (Lindemann, 1910) for melting is applied to the systems here to determine if they are in the solid phase or the liquid phase. A system is regarded as a solid if its Lindemann disorder parameter, _L, is in the range of 0.1–0.15, whereas a system with _L > 0.15 is considered to be liquid (Bilgram, 1987; Löwen, 1994; Stillinger, 1995; Zhou et al., 1999).

	RESULTS

TOP ABSTRACT INTRODUCTION Models and simulation method RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
We performed DMD simulations to investigate the equilibrium properties and the folding kinetics of the tetrameric ß-sheet complex. Four separate random coils were quenched from a temperature T^* = 2.0 to T^* = 0.2, the temperature at which the peptides are in their ß-sheet native state. After quenching, the four monomer chains undergo a conformational change toward the highly ordered ß-sheet complex (folded state) or to a partially folded or disordered (misfolded) state. The equilibrium properties of the assembled systems were investigated as a function of the intermolecular contact parameter , in the range 0.2 10. At each value of , at least 50 independent trajectories starting from different initial configurations were monitored and used in the averaging. All simulations were performed at a bias gap g = 0.9.

The dependence of the folding yield on the intermolecular contact parameter, , is shown in Fig. 2. The folding yield is defined to be the ratio of the number of trajectories leading to the folded state and the total number of trajectories at each value of . The folded and misfolded states were determined by visualizing the assembled tetramer structures with graphic tools and by measuring the values in the tetramer state. In particular, the values of for the two inner chains in the folded state should be twice the size of for the two outer chains as was seen in the global energy minimum state. It can be clearly seen that the folding yield is high at intermediate values of ( = 0.4, 0.5, 0.6, and 0.7), with a maximum at = 0.5. At these conditions, a large number of the trajectories fold into the highly ordered ß-sheet complex. However, at small ( = 0.2 and 0.3) and large ( = 0.8, 0.9, and 1.0) values of , the folding yield is low, indicating that only a few trajectories successfully fold into the highly ordered ß-sheet complex. The minimum value of the folding yield occurs at = 1.0.

View larger version (15K): [in this window] [in a new window]   FIGURE 2  Folding yield as a function of the intermolecular contact parameter, .

View larger version (23K): [in this window] [in a new window]   FIGURE 3  (a) The fraction of total native contacts, Qtotal, (b) the fractions of intramolecular native contacts, Qintra, and intermolecular native contacts, Qinter, and (c) the reduced internal energy per bead, E*/M, for the ß-sheet complex as a function of for the folded (solid symbols) and misfolded (open symbols) states.

For the ß-sheet complex, the size of the system after oligomerization and the compactness of the individual chains are determined by measuring the radius of gyration. The reduced squared radius of gyration per bead for the entire system, the reduced mean squared radius of gyration per bead for the individual chains, and the root-mean-squared (rms) pair separation fluctuation, _B, are shown in Fig. 4 as a function of for the folded (solid symbols) and the misfolded (open symbols) states. Here, denotes the average over the folding and misfolding trajectories, separately, and _n denotes the average over the n chains. For the folded state, the value in Fig. 4 a, which characterizes the system size, is small at = 0.2 and 0.3, abruptly increases at = 0.4, and then remains nearly the same for intermediate and large values of . For small , the value is small because the value of the ß-sheet monomer is small. This can be seen from Fig. 4 b, which shows that the mean value for the individual ß-sheet chains, is smallest at = 0.2 and 0.3 and then rapidly increases at = 0.4. The low values of occur at small , because, at these conditions the ß-sheet strands are not fully elongated even in the highly ordered state. Fig. 4 c shows the root-mean-squared (rms) pair separation fluctuation, _B, which measures the degree of bead fluctuations. The _B value is very high at the smaller , gradually decreases as increases, and is stable at the larger . This suggests that the ß-sheet strands are not fully stiff at the smaller due to the strong bead fluctuations. As increases, the intermolecular forces get stronger, preventing bead fluctuations and squeezing the ß-sheets together. The ß-sheet strands become fully elongated and each ß-sheet monomer in the complex structure recovers its intrinsically long molecular shape as seen in the global energy minimum structure in Fig. 1. This causes the values of the whole system and of the individual chains to increase.

View larger version (22K): [in this window] [in a new window]   FIGURE 4  (a) The reduced squared radius of gyration per bead, (b) the reduced mean squared radius of gyration per bead for the individual chains, and (c) the root-mean-squared (rms) pair separation fluctuation, B, for the ß-sheet complex as a function of for the folded (solid symbols) and misfolded (open symbols) states.

Folding kinetics of the ß-sheet complex
The folding kinetics of the ß-sheet complex was investigated. The time-dependent fraction of intermolecular native contacts, Q_inter(t), is shown in Fig. 5 as a function of the reduced time, t^*, at = 0.2, 0.5, and 0.8 for trajectories that end in: a), the folded state, and b), the misfolded state. Here, denotes the average calculated by sampling every 1000 collisions before 1 x 10⁷ collisions, and every 100,000 collisions after 1 x 10⁷ collisions over the folding and misfolding trajectories separately. To highlight the assembly events that occur early in the simulation, a logarithmic timescale is used in the figures to display early assembly behavior of the chains. The simulation results are presented for 4 x 10⁸ collisions, which is approximately equivalent to t^* = 3.4 x 10⁵, 1.9 x 10⁵, and 1.6 x 10⁵ for = 0.2, 0.5, and 0.8, respectively. To document the assembly process of the chains, the figure presents only Q_inter values for clarity. It is immediately apparent that the Q_inter values for both the folded and misfolded states begin to increase at t^* > 100. In contrast, the Q_intra values (not shown in figure) are saturated at the native value by t^* > 100. This is consistent with our previous kinetic study for isolated chain systems (Jang et al., 2002b), which showed that at g = 0.9 the ß-sheet is a fast folding protein and its native state is observed by t^* 100. At low , the assembly process for four chains goes through two distinct stages. In the first stage, the initial random coils collapse to the native ß-sheet at t^* 100, and in the second stage, the four ß-sheets assemble at t^* > 100. However, for large , the individual chains are more likely to collapse and assemble concurrently. For the folded state in Fig. 4 a, the Q_inter values increase in a series of steps after t^* = 100, whereas for the misfolded state in Fig. 4 b, the Q_inter values increase monotonically after t^* = 100 for all values of . The steps in the increase in Q_inter for the folded state correspond to assembly events in which the ß-sheets become aligned in parallel facing each other. This causes an increase in the number of intermolecular native contacts and hence a stepwise increase in Q_inter. In contrast, the continuous increase in Q_inter for the misfolded state indicates that the ß-sheets are assembled but not well aligned.

View larger version (22K): [in this window] [in a new window]   FIGURE 5  The time-dependent fraction of intermolecular native contacts, Qinter(t), as a function of the reduced time, t*, at = 0.2, 0.5, and 0.8 (a) for the folded state and (b) for the misfolded state.

View larger version (21K): [in this window] [in a new window]   FIGURE 6  The time-dependent reduced mean squared radius of gyration per bead for the individual chains, as a function of the reduced time, t*, at = 0.2, 0.5, and 0.8 (a) for the folded state and (b) for the misfolded state.

The = 0.2 model

View larger version (52K): [in this window] [in a new window]   FIGURE 7  Snapshots showing kinetic ß-sheet complex formation for = 0.2 at selected reduced times, t* = 0, 98, 976, 7.13 x 103, 2.12 x 104, and 1.72 x 105.

View larger version (81K): [in this window] [in a new window]   FIGURE 8  For = 0.2 (a) the time-dependent fractions of intramolecular and intermolecular native contacts, Qintra(t) and Qinter(t), (b) the time-dependent fraction of intermolecular native contacts for the nth chain, (c) the time-dependent reduced squared radius of gyration per bead, (t), (d) the time-dependent reduced squared radius of gyration for the nth chain, all as a function of the reduced time, t*, and (e) the Lindemann disorder parameter for an individual bead, Li, in the ß-sheet complex as a function of the reduced distance of bead i from the origin, |ri|/.