Folding Studies of a Linear Pentamer Peptide Adopting a Reverse

J. Phys. Chem. B 2000, 104, 8023-8034

8023

Folding Studies of a Linear Pentamer Peptide Adopting a Reverse Turn Conformation in Aqueous Solution through Molecular Dynamics Simulation Xiongwu Wu and Shaomeng Wang* Georgetown Institute for CognitiVe and Computational Sciences and Lombardi Cancer Center and Departments of Oncology and Neuroscience, Georgetown UniVersity Medical Center, Research Building, WP05, 3970 ReserVoir Road, Washington, DC 20007 ReceiVed: February 10, 2000; In Final Form: April 27, 2000

Type II reverse turn conformation has been shown to play a role in protein folding. In this study, reversible folding of a linear pentamer peptide YPGDV, which was determined by NMR experiments to have a significant population (50%) of a type II turn conformation in aqueous solution, was simulated in water with atomic detailed representation for both the peptide and solvent molecules at 300 K using the self-guided molecular dynamics (SGMD) simulation. During a 2 nanosecond (ns) SGMD simulation starting from a fully extended conformation, the peptide folds into a type II turn-like conformation and then undergoes unfolding and refolding several times. A cluster analysis method was developed based on conformational population and used to analyze the trajectories of the 2 ns SGMD simulation and of other simulations in this study. Five major conformational clusters were obtained from the 2 ns SGMD simulation and the most populated conformational cluster is a type II reverse turn-like conformation. The structure of the most populated conformational cluster identified through the SGMD simulation, or the “folded” states, is entirely consistent with the NMR data and the estimated relative NMR NOE strengths of proton pairs based upon the SGMD trajectory are in good agreement with the experimental data. Structural analysis of these major conformations indicated that the strong charged-charged interactions of the N-terminal positively charged group with the C-terminal and the side chain of Asp4 negatively charged groups and the hydrogen bonding between the amide group of Asp4 and the carboxyl group of Tyr1 are important for the stability of the “folded” states. In comparison, 2 ns conventional MD simulations starting from either a fully extended or a “folded” conformation failed to achieve folding or unfolding of this peptide, suggesting that longer simulations are needed for adequate conformational sampling using the conventional MD method. Based upon our estimated conformational free energy (calculated as the sum of conformational energy and solvation free-energy), the major conformational cluster obtained through the 2 ns SGMD simulation has the lowest estimated conformational energy and has a conformational free energy approximately 5 kcal/mol lower than any conformational cluster obtained through the 2 ns MD simulation started with the fully extended conformation. Folding studies of other short peptides and small proteins with defined secondary and/or tertiary structures using the SGMD simulation method are currently under way and will be reported in due course.

Introduction Protein structures are composed of elementary secondary structures, including turns, R-helices, and β-sheets. The formation of these secondary structures is believed to play a role in protein folding, which still remains to be one of the major unsolved and most challenging problems in biological science in post-genomic era. An in-depth understanding of folding of these basic secondary structures is an important step toward solving the protein-folding problem and for the design of peptides of specific biological functions with defined secondary structures. Short peptides, which adopt significant populations of secondary structures in solution, are excellent model systems with which to study folding of protein secondary structures. In recent years, experimental1-4 and computational5-19 studies have been carried out to study the thermodynamic and kinetic aspects of protein secondary structure folding. These studies have begun to provide important insights into the protein-folding problem * To whom all correspondence should be addressed. Phone: (202)-6872028. Fax: (202)-687-0617. E-mail: [email protected].

and may be useful for the design of functional peptides with defined secondary structures. Among several computational approaches, molecular dynamics (MD) simulation is a promising approach to study secondary structure folding. Theoretically, by simulating folding of such peptides with detailed atomic description of peptides and solvent (water) molecules, one may be able to predict the folded structure or major conformational states, to access the time scale of peptide and protein-folding events, to analyze factors contributing to the stability of the folded structures, and to uncover the folding mechanism.20 A major difficulty with MD simulation approach is that the time scale that can be accessed even for a relatively short peptide in aqueous solution is still quite limited, typically in the range of a few to a few-tens nanoseconds (10-9-10-8 s) with atomic detailed representation for the peptide and the solvent molecules with current computer power. With such short time-scale, only a limited number of peptide-folding events can be investigated. Another difficulty is that current force fields, such as CHARMM21,22 and AMBER23-25 were developed based upon

10.1021/jp000529i CCC: $19.00 © 2000 American Chemical Society Published on Web 07/28/2000

8024 J. Phys. Chem. B, Vol. 104, No. 33, 2000 small model molecules, and their applicability to peptide and protein folding has yet to be extensively tested. Because of these reasons, few studies have been reported to simulate the folding of peptides in solution with atomic detailed representation for both the peptide and the solvent molecules.18 Recently, Daura et al. performed a series of long MD simulations (50 ns for each simulation) at different temperatures (300-360 K) and successfully observed the reversible folding process of a β-heptapeptide, which was experimentally determined to adopt a stable left-handed-31-helix in methanol.18 The peptide is composed of β-amino acids instead of R-amino acids, and most simulations were conducted above room temperature, in methanol solvent rather than in water, which make the simulation study arguably less relevant to a real biological system. Nevertheless, their study is still highly significant in that for the first time the reversible folding of a peptide in solution with atomic detail was accomplished. To achieve reversible folding of the β-peptide, long simulations were needed not only to sample large conformational space but also, and perhaps more importantly, to ensure frequent interconversion between sets of folded and unfolded conformations. It is of note that, at room temperature (300 K), folding and unfolding of this β-peptide only occurred once during a 50 ns MD simulation in methanol starting from a folded conformation. It is expected that the conformational change of a peptide would be even slower in water because of its higher density, and reversible folding of peptides in water at room temperature through MD simulations would be more computationally demanding. Duan and Kollman have recently completed the very first 1 µs MD simulation of a 36-residue helical protein in water with atomic representation for both the protein and water molecules and have obtained important insights into the early protein-folding events.19 For such long simulation, of course, substantial computing resource was required and the simulation took a 3-month CPU time on a 256 CPU Cray T3E supercomputer with a highly efficiently paralleled AMBER program.19 Reverse turns in proteins have been thought to play a role in initiation of protein folding.27,28 Dyson et al. have designed a series of short peptides and investigated their solution structures.1 They found that while many of these peptides adopt random coiled structures, a few peptides have defined conformations in solution. One such peptide, a pentamer peptide with a YPGDV sequence with a trans-Pro was found to form a high population (50%) of a type II reverse turn conformation in solution based upon NMR experiments. The presence of a turn conformation of this peptide in aqueous solution was also confirmed by circular dichroism (CD) measurement, where the ellipticity had a positive maximum near 206 nm and a negative minimum near 190 nm.1 The folding of this peptide has been previously attempted using a detailed atomic model in water with the conventional MD method,5 but the simulation results obtained were somewhat inconsistent with the experimental results. It was, however, not clear if the inconsistency between the experimental and simulation results can be attributed either to the force field used in the simulation or to the insufficient simulation time (2 ns). We have recently developed a new molecular dynamics simulation method, termed self-guided molecular dynamics (SGMD) simulation method.29,30 We have shown through several model systems that SGMD simulation achieves a much greater conformational searching efficiency than the conventional MD simulation without significantly altering the conformational space of the system. Therefore, the SGMD simulation method enhances our conformational searching capability for

Wu and Wang molecular systems with rugged energy surface. In this paper, we investigated the folding of this pentamer peptide (YPGDV) in aqueous solution with atomic detailed representation for both the peptide and the water molecules using the SGMD method and the conventional MD method. Through this study, we hope to answer several questions highly relevant to peptide folding studies. (1) Is current force field accurate enough to study the folding of peptides in aqueous solution? (2) Are we able to achieve reversible folding of this peptide in solution with relatively short simulations with the SGMD method? (3) What is the “folded” state of this peptide in solution? (4) What interactions are important for the stability of the “folded” state? Method and Simulation Setup 1. SGMD Simulation Method. The fundamental idea of the SGMD simulation method is the derivation of a guiding force from its own simulation trajectory (a self-guiding force) and introduction of this self-guiding force into the equation of motion to enhance the systematic motion of the system.29,30 In conventional MD simulation, the equation of motion has the following form:

mir1 i ) fi

(1)

In the SGMD simulation, the equation of the motion is modified, as shown in eq 2.

mir1 i ) fi + λgi

(2)

Where mi and r1 i are the mass and the acceleration of atom i, fi is the force on atom i derived from the force field of the system, λ is the guiding factor, and gi is the guiding force, which can be estimated from the forces the atom experienced previously. Because of the additional guiding force in eq 2, this equation may be called as the self-guided equation of motion. The definition and the estimation for the guiding force gi are detailed below. The introduction of the guiding force, gi, is to enhance the systematic motions in a molecular system, which are defined as the average motions over a long time period and are primarily controlled by nonbonded interactions. But for molecule systems in which atoms are covalent bonded to each other, atoms have high frequency motions due to strong covalent bonding interactions. If these strong covalent interactions were included in the calculations of guiding forces, the guiding forces would be dominated by the high frequency noise of these strong covalent interactions. Two problems would occur by applying such guiding forces to a molecule system of interest. (1) The guiding forces would distort the basic bonded molecular geometry and (2) the guiding forces would not be effective in enhancing the systematic motions of the system. To overcome these problems, we introduced the concept of bonded substructures. Furthermore, we applied the nonbonded interactions evenly to the bonded substructures when calculating the guiding forces, as detailed below. A bonded substructure, Si, is defined for each atom i as the part of the molecule, which includes atom i itself, and all other atoms connecting to atom i through either a bond, a bond angle, or a dihedral angle. In other words, all atoms in a molecule connecting to atom i with three or fewer covalent bonds belong to the bonded substructure, Si. Atom i is called the central atom of Si. Each atom has a bonded substructure and two bonded substructures in a molecule can overlap partially or even totally with each other.

Folding Studies of a Linear Pentamer Peptide

J. Phys. Chem. B, Vol. 104, No. 33, 2000 8025

For each atom i, the nonbonded interactions used for the guiding force calculation are the interactions between the atom and atoms outside its bonded substructures.

fij + λgi ∑ j∉S

G ) Si

GSi

Here, is the contribution of the force acting on atom i to the guiding forces. fij is the interaction between atoms i and j, and the summation ∑j∉Si runs over all atoms outside the bonded substructure, Si. GSi comprises all forces acting on atom i including the actual guiding force, λgi, except bond length, bond angle, dihedral angle, improper dihedral angle, and 1-4 nonbonded interactions as defined in either AMBER23-25 or CHARMM21,22,40 force field. GSi is then distributed to all atoms in the bonded substructure, Si and for each atom, the distributed force is proportional to its mass, mj, so that the same acceleration is resulted for all atoms in the bonded substructure.

mj

∑ k∈S

GSi

(4)

mk

δt δt δt ) ri t + (fi(t) + λgi(t)) 2 2 mi

) (

∑ j∈S

Gi )

i

GSi j )

∑ j∈S i

mi

∑mk

GSj )

∑ j∈S i

k∈Sj

mi

∑mk

(

∑ fjk + λgj)

(5)

k∉Sj

gi(t) )

∫t-tt Gi(τ)dτ

1 tl

(6)

l

gi(t) =

( ) 1-

tl

gi(t - δt) +

δt

χE )

(

Ek t -

δt

∑ t j∈S l

i

δt

∑mk

(

∑ fjk(t) + λgj(t - δt)) k∉S

+

2

δt

fir3 i t -

2

i

+

2

(fi +

λg(S) i )r3 i

(7)

j

k∈Sj

where δt is the time step used in a MD simulation. Equation 7 shows that the guiding force at the current time step, gi(t), can be estimated from the guiding force at the previous time step, gi(t - δt), and the nonbonded interactions at current step, fjk(t), k ∉ Sj. In CHARMM and AMBER programs, since nonbonded forces are calcaluted separately, the guiding force can be easily calculated. The guiding forces obtained from eq 7 are combined with forces derived from the force field to solve the equation of the self-guided motion as shown in eq 2. We used the leapfrog algorithm31,32 to integrate the equation of the self-guided motion:

δt

t-

i

δt 2

δt

)

1/2

(10)

For constant temperature simulation, Berendsen et al.33 proposed a velocity scaling method that uses a velocity rescaling factor

(

χB ) 1 +

))

δt T -1 tT T′

(

1/2

(11)

to maintain a constant temperature. Where T′ is the current kinetic temperature and tT is a preset time constant. In this work, we combined the energy scaling factor, χE, and Berendsen’s velocity rescaling factor χB to obtain a scaling factor, χT, when performing a constant temperature self-guided MD simulation.

(12)

Therefore, the new velocities at t + (δt/2)t are scaled according to the following equation:

(

r3 ′i t +

δt δt ) χTr3 i t + 2 2

)

(

)

(13)

In summary, the detailed implementation of the SGMD method can be outlined as follows: (1) At time step t, calculate the forces between atom i and atoms within and outside its bonded substructure, ∑j∈Sifij(t) and ∑j∉Sifij(t), respectively. The total force acting on atom i is

fi(t) )

mi

(9)

( ) ∑ ( ) ( ) ∑ ( ) Ek t -

χT ) χBχE

As can be seen from eq 6, the guiding force gi(t) is a timedependent parameter and its characteristics also depends on the sampling time, tl, used in the calculation. For computing efficiency, this averaging is approximated by a simple memory function as shown in eq 7.

δt

)

Where r3 i(t) and fi(t) are the velocity and the force of atom i at time t, respectively. Because of the extra guiding forces, the system would not move along an iso-energy surface without additional treatment. To conserve the energy, we scaled the new velocities at t + (δt/2)t with an energy scaling factor, χE, to compensate the perturbation caused by the guiding forces:

k∈Sj

The guiding force for atom i is calculated as a time average of Gi over period of a local sampling time, tl:

(8)

δt δt 2

(

i

GSj i is the distribution of GSi on atom j, which belongs to Si. Because an atom can be included in several bonded substructures, it receives a distributed force from all these bonded substructures and the total distributed force on atom i is

)

ri(t + δt) ) ri(t) + r3 i t +

(3)

i

GSj i )

(

r3 i t +

fij(t) ) ∑fij(t) ∑ j∈S j∉S i

i

(2) Calculate the guiding force, gi(t), according to eq 7 using forces at current time step, ∑j∉Sifij(t), and the guiding forces at previous time step, gi(t - δt). At the beginning, gi(0) is simply set to 0. (3) Forward the velocities to the next half time step, t + (δt/ 2), according to eq 8. For constant energy simulations, the velocities are scaled with the energy-scaling factor calculated according to eq 10. For constant temperature simulations, the velocities are scaled according to eq 13. (4) Forward the positions to the next time step, t + δt, according to eq 9, but with the scaled velocities, r3 ′i(t + (δt/2)). If bond lengths or other internal coordinates are to be fixed, the SHAKE algorithm developed by Ryckaert et al.41 is applied to obtain the constrained positions.

8026 J. Phys. Chem. B, Vol. 104, No. 33, 2000

Wu and Wang

Figure 1. Chemical structure of the 5-mer peptide in this study.

(5) Repeat steps (1) to (4) until the end of the simulation. It should be pointed out that a SGMD simulation trajectory is influenced by the choice of both the local averaging time tl and the guiding factor λ. 2. Conformational Cluster Analysis Method. Conformational clustering is an effective approach to analyze a large number of conformations obtained during long simulations and conformational search,34-39 which reduces the complexity of data while minimizing the loss of information. In this work, we developed a new clustering method, which determines conformational clusters objectively according to their population, as detailed below. Many conformations obtained from a simulation may be redundant or very similar to each other and when performing conformational cluster analysis, it is critical to assess the conformational similarity between two conformations using a similarity function. The backbone φ and ψ dihedral angles of each amino acid residue can well describe the secondary structure of a peptide. Since our major interest in this work is the secondary structure of the peptide a similarity function is defined on the basis of the values of φ and ψ dihedral angles of each residue. For each backbone dihedral angle θi (φ or ψ), we define the similarity function between two conformations, m and n, as follows:

si(m,n) )

C C + [θi(n) - θi(m)]2

(14)

Where C is a constant, whose value determines the dependence of si(m,n) on the dihedral angle difference [θi(n) - θi(m)]. The brackets around θi(n) - θi(m) indicate that dihedral angle periodicity is taken into consideration. For a peptide of M backbone dihedral angles, the similarity function for the entire peptide between two conformations, m and n, is defined as M

S(m,n) )

si(m,n) ∏ i)1

(15)

Thus, S(m,n) has a value of 1 when m and n are identical, and its value approaches zero when m and n are very different. Using the similarity function, we can approximately estimate the

conformational population through a similarity density F(k), which is computed as follows:

F(k) )

1

N

∑ S(k,n) Nn)1

(16)

Where N is the total number of conformations sampled during a simulation. The larger the value of F(k), the more conformations similar to conformation k. F is a continuous function throughout the conformational space, and the value of F(k) is associated with the confomrational population in the vicinity of conformation k. We define a cluster as conformations around a central conformation, c, whose similarity density, F(c), is the highest among them. In other words, on the similarity density surface, F(k), each peak corresponds a cluster. To avoid statistical noise in determining peaks on the F surface, we defined a cluster resolution S0. No two peaks with a similarity value greater than S0 will be distinguished. If there are two peaks with a similarity higher than S0, the two peaks will be combined into a single cluster, and the one with higher F will be identified as the cluster center. The central conformations of any two clusters should have a similarity value smaller than S0. The conformational clustering is performed in two steps: (1) identifying central conformations of clusters, c, which should satisfy the following conditions: F(c) g F(i) for any conformation i when S(i,c) > S0; (2) clustering conformations in such a way that conformation i will be assigned to cluster c if similarity value S(i,c) is higher than similarity values to any other cluster central conformation. In this study, we chose C equal to 300 degree2 in eq 2 and S0 ) ((3/4))M, where M is the total number of backbone dihedral angles of the peptide and equal to 8 for the pentamer peptide studied in this work. This means that on average the difference for each dihedral angle between the central conformations of two clusters is greater than 10°. This clustering method is relatively objective, since it depends only on the definition of the similarity function eq 2 and the cluster resolution S0. By clustering the conformations, a large number of conformations recorded during a simulation are reduced to a limited number of clusters. 3. Simulation Setup. The CHARMM program40 (version 26) was used to perform all conventional MD simulations. The SGMD simulation method29,30 has been implemented in the

Folding Studies of a Linear Pentamer Peptide

J. Phys. Chem. B, Vol. 104, No. 33, 2000 8027

CHARMM program (version 24) and was used for all the SGMD simulations. The chemical structure of YPGDV is shown in Figure 1. The peptide was solvated using 805 TIP3P water molecules25 in a cubic box, and the box size was set to be 29 Å. Cubic periodic boundary conditions were applied in all simulations. CHARMM22 force field was used.21 Nonbonded interactions are smoothly truncated with the method described in CHARMM.40 The cutoff distance for nonbonded interactions was set at 13 Å, and a switch function was applied for both the electrostatic and LJ interactions between 8 and 12 Å to smooth the change across the cutoff. A neighbor list was built and updated every 10 MD steps for nonbonded interaction calculation. The simulation temperature was set to be 300 K. The SHAKE algorithm41 was applied to constrain all bonds so that 2 fs (2 × 10-15 s) time steps can be used in our simulations. The trajectories were recorded every 2 ps for analysis. For the SGMD simulations, the local averaging time tl is set to be 2 ps and the guiding factor was set to be 0.5. Guiding force was calculated and applied to all the atoms in the system. 4. Conformational Free Energy Calculation. The conformational free energy (∆G) for each peptide conformation was estimated by summation of the conformational (Ep) and solvation (∆Gs) energies. This implies that we assume that the conformational entropy difference between two peptide conformations is relatively small and can be ignored. Conformational energy (Ep) for each conformation was obtained during the CHARMM simulation. Solvation energy (Gs) for each conformation was calculated using the AMSOL program (version 6.5.3), running on an SGI workstation. Peptide conformations obtained from the simulations were used directly for solvation energy calculations without further minimization. SM1-AM1 Hamiltonian was used and the peptide was thus assigned a net charge of -1 in all the solvation energy calculations, according to its protonation state at pH equal to 7 in aqueous solution. Results and Discussion 1. Conformational Changes of the Peptide during a 2 ns SGMD Simulation. A 2 ns SGMD simulation was carried out starting from a fully extended conformation, in which two backbone dihedral angles φ and ψ for each amino acid residue were set to be 180°, except for Pro2, in which φ and ψ were set to be -73° and 180°, respectively. During the simulation, the peptide experienced significant conformational change. Within 10 ps, Pro2 quickly reached a conformation with (φ, ψ) around (-40, 90) and remained in this conformational region for the remaining period of the simulation. At 34 ps, the φ and ψ angles of Gly3 changed to (-70, 90), where a hydrogen bond was formed between the NH group of Asp4 and the carbonyl group of Pro2. However, this hydrogen bond existed for only 2 ps and the φ and ψ angles of Gly3 quickly changed to (60, 60), where a new hydrogen bond was formed between the NH group of Asp4 and the carbonyl group of Tyr1 and a type II reverse turn formed. At 38 ps, a strong ionic interaction occurred between the Nterminal, positively charged NH3 group and the negatively charged carboxylate group in the side chain of Asp4. This conformation lasted for about 80 ps. At 120 ps, the backbone conformation underwent a dramatic conformational change, where the ionic interaction between the N-terminal NH3 group and the carboxylate group of Asp4 side chain decreased substantially, but the hydrogen bond between the NH of Asp4 and the carbonyl group of Tyr1 remained during this conformational change. From 120 to 300 ps, major conformational

Figure 2. Representative conformations during the first 300 ps of the SGMD simulation. From top to bottom are conformations at 0, 34, 38, 120 and 300 ps, respectively.

changes occurred, primarily involving Tyr1, Asp4, and Val5. During this period, a hydrogen bond between the NH group of Val5 and the carbonyl group of Tyr1 switched on and off and hydrogen bonding/ionic interactions occurred between the C-terminal carboxylate group and the NH group of Gly3. Interactions between the carboxylate group of Asp4 side chain and the N-terminal NH3 group were evident in some conformations but transient. At 304 ps, strong ionic interactions occurred between the N-terminal positively charged NH3 group and the C-terminal negatively charged carboxylate group. In this conformation, the side-chains between Tyr1 and Pro2, and between Asp4 and Val5, were in close proximity, suggesting hydrophobic interactions between these groups. This conformation apparently represents the most populated conformational space sampled during the 2 ns SGMD simulation based upon our conformational analysis. A number of representative conformations are shown in Figure 2 to illustrate the conformational changes during the first 300 ps simulation. The peptide underwent continuous and frequent conformational changes during the entire 2 ns simulation. 2. Conformational Cluster Analysis. Using the conformational clustering method described in the Method section, the 1000 conformations recorded during the 2 ns SGMD simulation were analyzed and a total of 5 major clusters were identified. The central conformations of these 5 clusters are the conformations at 1.240 ns for cluster I, at 0.808 ns for cluster II, at 1.058

8028 J. Phys. Chem. B, Vol. 104, No. 33, 2000

Wu and Wang

Figure 3. Five major clusters identified from the trajectory obtained from the 2 ns SGMD simulation (upper row: clusters I, II, and III; lower row: clusters IV and V).

ns for cluster III, at 0.082 ns for cluster IV and at 0.032 ns for cluster V, respectively. The 3D structures of these 5 central conformations are shown in Figure 3. Our analysis showed that conformational cluster I represents 49.5% of total conformations recorded during the simulation, while conformational clusters II and III represent 20.5% and 22.0%, respectively. Conformational clusters IV and V account for 5.5% and 1.2%, respectively. Altogether, these 5 clusters account for 98.7% of the total conformations. As can be seen from Figure 3, conformational cluster I is characterized by two consecutive turns, one from residues 1 to 4 and the other from 2 to 5. The turn from residues 1 to 4 has a hydrogen bond between 1O and 4HN and resembles a type II reverse turn. The turn from residues 2 to 5 has a hydrogen bond between 3HN and 4O but doesn’t resemble any known turn structures in proteins. In addition to these two hydrogen bonds, there are strong charge-charge interactions of the N-terminal positively charged group with the C-terminal and the side chain of Asp4 negatively charged carboxylate groups. Both clusters II and III have a turn conformation that doesn’t resemble any known turn structure in proteins. The major interactions in clusters II and III occur between the N-terminal positively charged group and the C-terminal negatively charged carboxylate group. Cluster IV resembles a type II conformation but is different from that of cluster I. The major interactions in cluster IV are the charged interactions between the N-terminal positively charged group and the negatively charged group of Asp4 side chain and the hydrogen bonding between 1O and 4HN. The crucial difference between cluster I and cluster IV is that cluster IV lacks the charged-charged interactions between the N-terminal positively charged group and the negatively charged carboxylate group of the C-terminal. Cluster V is a typical random coil conformation and has neither intramolecular hydrogen bonding nor charged interactions. The NMR experimental data suggested that the β-turn conformation formed by this peptide may be a type II turn conformation. To compare with the standard reverse type II turn, the backbone dihedral angles of these 5 central conformations,

TABLE 1: NOE Data Obtained from the NMR Experiment and Estimated from the SGMD Simulationa NOE intensity

atom pair types

residue pairs

NMR experiment

SGMD simulation

dRN(i,I+1)

2-3 3-4 4-5 2-4 3-4

very strong moderately strong strong weak medium

8.54 (very strong) 2.24 (strong) 3.38 (strong) 0.20 (medium) 0.88 (medium)

dRN(i,I+2) dNN(i,I+1)

ad R RN(i,j), intramolecular distance between C H and NH on residues i and j; dNN(i,i+1), intramolecular distance between NH and NH on residues i and j.

as well as that of the ideal type II reverse turn, are shown in Table 2. An ideal type II reverse turn in proteins adopts a conformation with φ2 ) -60°, ψ2 ) 120°, φ3 ) +90°, and φ3 ) 0°. As can be seen from Table 2, although none of the major clusters are the ideal type II reverse turn conformation found in proteins, clusters I and IV closely resemble the type II reverse turn conformation. 3. Comparison with NMR Experimental Results. The solution conformation of this pentamer peptide YPGDV has been extensively studied by Dayson et al. using NMR techniques.1 The Asp4 amide proton temperature coefficient for this peptide is only 3.3 ppb/K, considerably lower than that for other pentamer peptides investigated in their studies. This suggested that a high population of β-turn in solution might exist. Indeed, the NOE connectivities observed for this peptide confirmed its β-turn conformation in solution. For example, a strong NOE between the Gly3 and Asp4 amide protons was observed. The sequential cross-peak between the NH of Asp4 and the CRH of Gly3 is reduced in intensity as compared to the Asp4 CRHVal5 NH and Pro2 CRH-Gly3 NH NOE cross-peaks. Of great importance, a cross-peak between the CRH resonance of Pro2 and the NH resonance of Asp4 was observed. Taken together, the very low Asp4 NH temperature coefficient and the dNN(3,4) and dRN(2,4) NOE connectivities provide unequivocal evidence for a relatively high population of a 4-1 hydrogen-

Folding Studies of a Linear Pentamer Peptide

J. Phys. Chem. B, Vol. 104, No. 33, 2000 8029

TABLE 2: Dihedral Angles of the Central Conformation for Each Cluster and the Ideal Type II Turn simulation SGMD

MD (ext.)

MD (folded) ideal type II turn

clusters

ψ1

φ2

ψ2

φ3

ψ3

φ4

ψ4

φ5

1 2 3 4 5 1 2 3 4 1 2

127.41 136.91 171.12 153.33 -177.02 152.82 146.65 141.11 156.09 160.03 178.46 -

-43.00 -38.79 -51.16 -57.52 -51.34 -63.98 -49.69 -52.57 -62.86 -49.32 -37.05 -60

123.47 123.00 125.27 123.03 131.58 116.85 125.77 145.53 128.61 118.18 126.68 120

130.49 90.71 80.60 60.36 -69.28 -76.03 -81.77 -164.01 -71.19 133.61 91.07 90

-63.24 -86.93 -82.80 40.42 103.93 -54.26 -109.18 -53.89 -55.19 -65.85 5.17 0

-75.04 -77.72 -59.98 -66.55 -77.23 -73.38 -61.29 -82.08 -171.82 -84.37 -144.15 -

-123.31 85.02 157.91 -34.64 90.54 127.43 -50.48 104.08 -55.13 -128.78 -144.82 -

-126.59 38.51 51.24 -125.95 69.52 -138.70 -122.94 -166.26 -97.25 -111.46 -71.38 -

bonded β-turn in this peptide. Formation of β-turn by this peptide was also confirmed independently by CD measurements.1 The presence of a significant population in water solutions of this peptide was indicated by the positive ellipticity with a maximum near 206 nm and by the minimum near 190 nm. We first examined if the central conformations of these major clusters are in agreement with the NMR experimental results. As can be seen from Figure 3, clusters I and IV have a strong 1-4 hydrogen bond between the carbonyl group of Tyr1 and the amide group of Asp4, in good agreement with the NMR experimental results.1 These two clusters account for a significant population of the total conformations (55%). On the basis of the temperature coefficient of the amide group of Asp4, Dayson et al. estimated1 that the turn population of this peptide in aqueous solution was approximately 50% and our simulation results appear to be in good agreement with the experimental results. In NMR experiments, NOE effects represent average structural information in solution, rather than a specific conformation. To provide a more vigorous validation, we estimated the average NOE effects from our simulation conformations. Because NOE depends on the sixth power of distance, the ensemble average effect of NOE may be estimated by the following equation:

NOEij ) Cij 〈r-6 ij 〉

(17)

where NOEij represents the strength of the NOE between atoms i and j; Cij is a constant whose value depends on the atomic type of both i and j; rij is the distance between atom i and j in a conformation. The symbol 〈 〉 indicates an ensemble average. Since Cij depends on experimental conditions and atomic types, it is thus difficult to calculate precisely a quantitative value of NOEij based upon simulation trajectories. To facilitate the calculations, we assumed that Cij adopts the same value regardless of the atomic types of i and j. It was also assumed that a strong NOEij signal is being observed when rij is equal to or less than 3 Å and when rij is equal to 3 Å, NOEij has a relative intensity of 1 and Cij thus assumes a value of 36 (729). The relative NOE intensities for all the atom pairs within a molecule may be calculated from the simulation trajectory according to the following equation:

NOEij ) 729〈r-6 ij 〉

(18)

With the above equation, the relative NOE intensity between two protons can be quantitatively estimated. It is generally assumed if the distance between two protons is within 3 Å, their NOE signal is strong (NOEij g 1), if between 3 and 4 Å, their NOE signal has a medium intensity (1 > NOEij g 0.18), and if between 4 and 5 Å, their NOE signal is week (0.18 > NOEij g 0.04).

Using all the 1000 conformations recorded during the 2 ns SGMD simulation, the ensemble average NOE effects between each relevant proton pair are calculated and the results are listed in Table 1. As can be seen, the agreement between our calculated and the experimental NOE effects is excellent. 4. Conventional MD Simulations. It would be instructive to carry out MD simulations under same conditions (force field, temperature, solvent, and other simulation conditions) using the conventional MD simulation method to validate the results obtained using the SGMD method and to investigate the improvement of the conformational searching efficiency of the SGMD method. For these reasons, we have carried out two simulations using the conventional MD simulation method. One started with the same fully extended conformation as that in the SGMD simulation, and the other started from the central conformation of cluster I (the “folded” conformation) obtained from the SGMD simulation. Starting from the fully extended conformation, a bend appeared in the peptide during the first 240 ps, involving residues 2-4. The distance between the NH group of Asp4 and the carbonyl group of Pro2 oscillated between 1.7 and 4.6 Å, indicating the formation of a weak hydrogen bond. After 240 ps, another bend occurred, involving residues 3-5. This bend propagated to become a type I turn, involving residues 2-5. In some conformations, a hydrogen bond was found between the NH group of Gly3 and the C-terminal carboxylate group. At ca. 700 ps, hydrophobic contacts were made between the side chains of Val5 and Pro2, as well as between Tyr1 and Pro2. This type I turn appeared to be quite stable because it persisted from 300 to 1000 ps. At 1000 ps, a new bend occurred, involving residues 1 and 3 and the peptide adopted an S-shape conformation. This conformation remained until the end of the 2 ns simulation. It is of note that none of these conformations sampled during the 2 ns simulation are in agreement with the experimental data for this peptide because of the lack of the 1-4 hydrogen bond between the Tyr1 CO and Asp4 NH groups, as clearly indicated by the NMR results. Furthermore, the NMR results indicated1 that an electrostatic interaction involving the negatively charged carboxylate group of the Asp4 side chain may play a role to stabilize the turn but this electrostatic interaction was not found in any conformation obtained from the 2 ns MD simulation. In Figure 4a,b, the conformational space sampled during the 2 ns conventional MD simulation and during the 2 ns SGMD simulation are plotted using the 8 backbone φ and ψ angles of the peptide. As can be seen, a major conformational space sampled during the 2 ns SGMD simulation was not accessed during the 2 ns conventional MD simulation. The 2 ns SGMD simulation sampled a much larger conformational space than the 2 ns conventional MD simulation. The most populated conformational cluster I, obtained from the 2 ns SGMD

8030 J. Phys. Chem. B, Vol. 104, No. 33, 2000

Wu and Wang

Figure 4. (a) Backbone dihedral angle space sampled during the 2 ns MD simulation starting from a fully extended conformation; (b) backbone dihedral angle space sampled during the 2 ns SGMD simulation starting from a fully extended conformation; (c) backbone dihedral angle space sampled during the 2 ns MD simulation starting from cluster I obtained in the SGMD simulation.

simulation and agreed with all the NMR experimental data, was not accessed during the 2 ns conventional MD simulation. To investigate whether the major conformational cluster obtained from the 2 ns SGMD simulation is a stable conformation, a conventional MD simulation starting from the central conformation of cluster I was carried out for 2 ns. The conformational space sampled during the 2 ns conventional MD simulation is plotted in Figure 4c. As can be seen, a highly localized conformational space for the peptide was sampled during the 2 ns simulation, indicating that the peptide essentially stayed in a local conformational space around cluster I. This highly localized sampling of the conformational space also suggested that the conformational cluster I is a very stable conformation, but this conformation was not accessed during the 2 ns conventional MD simulation starting from an extended conformation. Previously, two simulations were performed with this peptide for 2.2 ns under similar conditions using the conventional MD method.5 In one simulation, which started with the peptide in a fully extended conformation, the peptide folded into a transient type II turn before unfolding. The peptide remained unfolded for another 0.9 ns before folding into a type I turn involving residues 2-5. Our simulation results using the conventional MD simulation method, started with the peptide with a fully extended conformation, were basically in good agreement with the simulation results.5 However, it was found that the type I turn reported previously was somewhat different from a true type I turn, because the hydrogen bond between Pro 2 and Val 4 is lacking. In addition, the type II turn, which was found to have a transient existence, was also somewhat different from an ideal type II turn because of the lack of hydrogen bonding. Conformational analysis showed that the results obtained from the 2 ns conventional MD simulations started from an extended conformation did not agree with the experimental NMR results.1 For example, NMR data indicated that an intramolecular hydrogen bond for the amide group of Asp 4 should be formed. However, this crucial hydrogen bond was not found in any conformation obtained from current and previous 2 ns simulation started with a fully extended conformation. Thus, in order to properly sample the solution conformational space of this peptide, simulations longer than 2 ns are required using the conventional MD method. It should be noted that performance of conventional MD simulations much longer than 2 ns even on this relatively small peptide system in aqueous solution was prohibitively expensive when previous MD simulations were performed due to insufficient computing power.

Previous conventional MD simulation started with the peptide in an ideal type II reverse turn involving the first four residues showed that the turn unfolded after about 1.4 ns.5 This result appears to be different from our present study since we found that the type II reverse turn-like conformation remained “folded” during the entire 2 ns conventional MD simulation. This is thus intriguing why two conventional MD simulations obtained different results. Close examination showed that the starting conformations of these two simulations are similar but with some crucial difference. While both conformations had the crucial 4-1 hydrogen bond, the ideal type II reverse turn conformation used in previous simulation did not have the charged-charged interactions between the N-terminal positively charged group and the C-terminal or the side chain of Asp4 negatively charged carboxylate groups. On the basis of our 2 ns SGMD simulation results, it was found that these chargedcharged interactions play an important role to the stability of the type II reverse turn-like conformation (clusters I), which was supported by the NMR experimental results. Because of the lack of the charged-charged interactions in the ideal type II turn conformation, this conformation seemed to be less stable than the type II turn-like conformation (cluster I) obtained from the 2 ns SGMD simulation. Therefore, it is not surprising that the ideal type II turn conformation unfolded after 1.4 ns simulation while the type II turn-like conformation remained in such conformation during the entire period of 2 ns simulation in our study. We have shown that the type II turn-like conformations (cluster I) are entirely consistent with the experimental NMR data, suggesting that the formation of the 4-1 hydrogen bond in this peptide doesn’t necessarily require the formation of an ideal type II turn conformation. In fact, the ideal type II turn conformation is only partially in agreement with the experimental NMR data, while the type II turn-like conformations (cluster I) are entirely consistent with the experimental NMR data. Taken together, these results suggest that the type II turn-like conformations but not the ideal type II turn conformations are most likely the folded structure of this peptide in water based upon our simulation results and the experimental data. 5. Estimation of Conformational Free Energy. The “folded” state identified through the SGMD simulation is consistent with the NMR experimental data and the 2 ns conventional MD simulation started with the peptide in its “folded” state also showed that this “folded” conformation is very stable for the peptide remained in that conformation during the entire 2 ns simulation. However, there is a possibility that the agreement

Folding Studies of a Linear Pentamer Peptide

J. Phys. Chem. B, Vol. 104, No. 33, 2000 8031 TABLE 3: Average Conformational Energies, Solvation Free Energies, and Estimated Conformational Free Energies of the 5 Major Conformational Clusters

Figure 5. Conformational free energy profile of the 5-mer peptide against the simulation time.

between the SGMD simulation results and the NMR experimental results is merely incidental. To further examine whether cluster I is the “folded” state, we have estimated the conformational free energy for each conformation obtained from both the SGMD and the MD simulations. The relative conformational free energy may be divided into three components: (1) the conformational energy of the peptide, (2) the solvation free energy of the peptide, and (3) the conformational entropy of the peptide (solute entropy). While the conformational energy and solvation free energy can be directly calculated, an accurate estimation of the solute entropy contribution is quite difficult.44 If we assume that the solute entropy difference between two conformational states is small for this relatively short peptide, the conformational free energy may be roughly estimated by the summation of the conformational energy and solvation free energy of the peptide. Care should be taken in interpreting the results, but such estimates should provide us with a semi-qualitative measure on the relative stability between two different conformational clusters. The estimated conformational free energy of the peptide for each conformation obtained from the 2 ns SGMD simulation and two conventional MD simulations is plotted against the simulation time in Figure 5. The conformational free energy of the peptide in the MD simulation starting from the “folded” conformation is around -330 kcal/mol. However, the lowest conformational free energy reached during the 2 ns MD simulation started from a fully extended conformation is -325 kcal/mol. During the 2 ns SGMD simulation, the conformational states with a conformational free energy around -330 kcal/ mol were repeatedly searched. It is of note that these conformations with a conformational free energy around -330 kcal/mol belong to conformational cluster I and are consistent with the NMR experimental data. The results suggested again that the 2 ns SGMD simulation was able to search a much larger conformational space and identify the “folded” state with a free energy approximately 5 kcal/mol more stable than any conformation searched during the 2 ns MD simulation started from the fully extended conformation. In Table 3, the average conformational energy, solvation free energy, and estimated conformational free energy for the 5 major conformational clusters obtained from the 2 ns SGMD simulation are provided. As can be seen, cluster I, whose structure was found to be consistent with all NMR data, has the lowest estimated conformational free energy among all the 5 clusters, supporting that this conformational cluster may be indeed the “folded” state

clusters

Ep kcal/mol

∆Gs kcal/mol

∆G kcal/mol

δ(Ep + ∆Gs) kcal/mol

I II III IV V

-230.0 ( 0.8 -213.4 ( 0.7 -207.7 ( 1.2 -150.0 ( 0.7 -67.8 ( 0.3

-97.3 ( 0.7 -107.8 ( 0.6 -115.7 ( 1.0 -166.1 ( 0.7 -253.0 ( 0.3

-327.3 ( 0.2 -321.2 ( 0.3 -323.4 ( 0.3 -316.1 ( 0.2 -320.8 ( 0.2

7.3 7.1 6.8 7.1 7.9

of the peptide in aqueous solution. Although cluster IV also satisfies some of the NMR data, such as the hydrogen bonding between 1O and 4HN, this cluster has the highest conformational free energy. Cluster V, which represents a typical random coil conformation, has the second highest conformational free energy. It is of note that cluster V has the lowest (most negative) solvation free energy and highest conformational energy. This is consistent with it structural characteristic for the structure has neither intramolecular hydrogen bonding nor charged interactions and all the backbone and side chain polar groups are exposed to solvent. Clusters II and III have an estimated conformational free energy 4 to 6 kcal/mol higher on average than cluster I. Therefore, despite the caveat in the our conformational free energy calculations, the estimations still provided a semiquantitative measure of the relative stability of peptide conformations in solution and suggested that cluster I is the most stable conformational state, in agreement with MD and SGMD simulation results. 6. Factors Contributing to the Folding and the Stability of the Type II Turn Conformation. The 2 ns SGMD simulation achieved a reversible folding and identified major conformational clusters that are entirely consistent with the experimental results. This provided us the opportunity to analyze factors that are important to the folding of this peptide in water. This pentamer peptide consists of a Pro at position 2. The influence of Pro residue to peptide and protein conformations has been extensively studied and documented.45-48 Formation of a reverse turn conformation restricts a number of rotatable bonds and results in a loss in conformational entropy. Because the side chain of the Pro is a fused ring, the accessible conformational space of φ and ψ angles of Pro is limited. This characteristic of Pro makes its entropy loss to be minimal when the turn is formed. This peptide consists of a Gly at position 3. Gly is a unique residue because of the lack of a side chain and its φ and ψ angles can adopt certain values that are prohibited or unfavorable for other residues. The Ramachandran plot for Gly obtained from the 2 ns SGMD simulation showed that in the most populated conformation, the φ and ψ angles adopted values of (120, -60), which locates in a unfavorable conformational region for other residues. Indeed, when Gly was substituted by other residues experimentally, these substituted pentamer peptides all displayed a significantly decreased β-turn conformation population in solution.1 Hydrogen bond/ionic interactions may play an important role for the formation of protein secondary structures. A number of hydrogen bonds/ionic interactions are presented in the turn conformations obtained from the SGMD simulation. Furthermore, there appeared to be van der Waals contacts between the side chains of Tyr1 and Asp4. To analyze their roles of these interactions for the solution conformations of this peptide, we have plotted 8 atomic distances against simulation time in Figure 6a-h. As can be seen from Figure 6a, the distance between the oxygen of the backbone carbonyl group of Tyr1 and the

8032 J. Phys. Chem. B, Vol. 104, No. 33, 2000

Wu and Wang

Figure 6. Distance profile between two atoms against simulation time. (a) distance between the oxygen of the backbone carbonyl group of Tyr1 and the hydrogen of the amide group of Asp4; (b) distance between the oxygen of the backbone carbonyl group of Tyr1 and the hydrogen of the amide group of Val5; (c) distance between the oxygen of the backbone carbonyl group of Pro2 and the hydrogen of the amide group of Val5; (d) distance between the beta carbon (Cβ) on the side chain of Tyr1 and the beta carbon (Cβ) on the side chain of Asp4; (e) distance between one of the hydrogen atoms of the ammonium group (HNT1, HNT2, or HNT3 in Figure 1) of Tyr1 and the oxygen of the backbone carbonyl group of Asp4; (f) distance between one of the hydrogen atoms (HNT1, HNT2, or HNT3 in Figure 1) of the ammonium group (NH3) of Tyr1 and one of the oxygen atoms of the side chain carboxylic group (OD1 or OD2 in Figure 1) of Asp4; (g) distance between one of the hydrogen atoms of the ammonium group (HNT1, HNT2, or HNT3 in Figure 1) of Tyr1 and one of the oxygen atoms of the carboxylic group of Val5; (h) distance between the hydrogen atom of the backbone amide of Gly3 and one of the oxygen atoms of the carboxylic group (OCT1 or OCT2 in Figure 1) of Val5.

hydrogen of the backbone amide group of Asp4 oscillated around 2 Å during the simulation, indicating the formation of a strong hydrogen bond and the formation of β-turn conformations during most of the simulation time. But the hydrogen bond was also broken several times and re-formed after certain time period. Other intramolecular hydrogen bonds formed during the simulation included the hydrogen bonds between the backbone carboxyl group of Pro 2 (2O) and the backbone amide group of Val5 (5HN), as shown in Figure 6c, between the backbone carbonyl group of Tyr1 and the backbone amide group of Val5, as shown in Figure 6b, and between the amide group of Gly3 and the C-terminal negatively charged carboxylate group, as shown in Figure 6h. However, for these three hydrogen bonds, although in some conformations, the hydrogen bonding distances became sufficiently short for the formation of strong hydrogen bonds, the distances fluctuated frequently within a wide range, indicating on average, these hydrogen bonds were generally weak. Other major interactions observed that may play a role for the solution conformations of this peptide included the ionic interactions/hydrogen bond by the positively charged N-terminal NH3 group with other groups in the peptide. This positively charged group formed a strong ionic interaction with the C-terminal, negatively charged carboxylate group. A close distance (2 Å) between any of the three NH3 hydrogen atoms and any of the two carboxylate oxygen atoms is an indication of this strong interaction. As can be seen from Figure 6f, after the first 300 ps, the distance between one of NH3 hydrogen

atoms and one of the two carboxylate oxygen atoms became very short (2 Å) with little fluctuation, indicating a very strong interaction between these two groups. This strong interaction was found during the entire remaining simulation period except for at 1565 ps, a small burst occurred and the interaction between these two groups became weak. But the interaction was quickly reestablished. The N-terminal positively charged NH3 group also formed a hydrogen bond with the backbone carbonyl group (4O) of Asp4. As can be seen from Figure 6e, during the 2 ns SGMD simulation, a hydrogen bond was formed between the two groups in approximately 40% of conformations. In addition, a strong ionic interaction was also found between the N-terminal NH3 group and the side chain carboxylate group of Asp4 and existed in approximately 30% of the conformations (Figure 6f). As can be seen from Figure 6e,f, these two interactions seemed to have a strong anti-correlation, i.e., the N-terminal NH3 group switched back and forth between the backbone carbonyl and the side chain carboxylate groups of Asp4. It was found that during the SGMD simulation, there were some van der Waals contacts between the side chains of Tyr1 and Asp4. In Figure 6d, the distance between two β carbon atoms of Try1 and Asp4 was plotted against the simulation time. As can be seen, the distance decreased to less than 6.6 Å, shorter than the diameters of a carbon atom plus that of a water molecular. During the 2 ns SGMD simulation, this distance became as short as 4 Å in many conformations. The close proximity precludes insertion of water molecules between these

Folding Studies of a Linear Pentamer Peptide

J. Phys. Chem. B, Vol. 104, No. 33, 2000 8033

Figure 7. Correlation between estimated conformational free energy and the distances specified in Figure 6a-h.

two groups, indicating that these two groups may have van der Waals contacts in some conformations. To gain a better understanding on their roles of these hydrogen bonding, electrostatic/charged interactions and van der Waals contacts to the stability of the folded conformations (the conformational states with low conformational free energy), we have plotted the estimated conformational free energy shown in Figure 5 against the 8 atomic distances shown in Figure 6ah, as shown in Figure 7a-h. If we take conformations with an estimated conformational free energy below -330 kcal/mol as the “folded” states, as evident from Figure 7a, a strong hydrogen bond formed between the backbone carbonyl group of Tyr1 and the amide group of Asp4 always exists in those “folded” conformational states. Similarly, as can be seen clearly from Figure 7g, a strong charged-charged interaction between the

ammonium group of Tyr1 and the carboxylic group of Val5 exists in all “folded” states, suggesting that this charge-charge interaction may play a role for the stability of the “folded” states for this peptide in solution, in agreement with the experimental observation.1 However, a recent computational study showed that this type of strong charge-charge (or salt-bridge) interaction may not be important for the stability of protein structures due to strongly unfavorable desolvation free energy penalty upon formation of a salt bridge.49 One explanation for the discrepancy between our present and previous studies may be due to the fact that in this small 5-mer peptide, formation of the salt-bridge between the positively charged N-terminal and the negatively charged C-terminal still allows partial exposure of these two highly charged groups to solvent molecules. Thus, the desolvation free energy penalty is not as large as that for salt-bridges

8034 J. Phys. Chem. B, Vol. 104, No. 33, 2000 buried in the interior of proteins. Other interactions seem to play less important roles for the stability of the “folded” conformations since they mainly exist in other solution conformations with an estimated conformational free energy larger than -330 kcal/mol. Conclusion In this paper, we demonstrated that the reversible folding of a pentamer peptide in aqueous solution at room temperature with atomic detailed representation for both the peptide and the solvent molecules can be achieved through the use of a new MD method with an efficient conformational searching ability using moderate computer power (a workstation) and a relatively short simulation time (2 ns). This pentamer peptide was found to form a type II reverse turn based upon the experimental NMR results. The “folded” state identified through our SGMD simulation is consistent with all the experimental data. Structural analysis showed that hydrogen bonding and charged-charged interactions play an important role for the stability of the type II turn conformation in aqueous solution. Two 2-ns conventional MD simulations started either from the fully extended conformation or from the “folded” state failed to achieve folding or unfolding. Conformational free energy estimation showed that the “folded” state has the lowest conformational free energy among all the 5 major clusters identified through the SGMD simulation. Our results demonstrated that the CHARMM force field employed in our current study may be of sufficient accuracy to study the folding of protein secondary structures in solution. Indeed, in addition to our current study, to date, we have successfully completed the folding of peptides, which formed helices or β-hairpins in aqueous solution, with atomic detailed representation for both the peptides and the solvent molecules. These studies pave the way to study the folding of other peptides or even small proteins with defined secondary or tertiary structures and may prove to be valuable for the design of peptides with defined structures and specific biological functions. Acknowledgment. This work is supported by a NIH grant GM-59188 and by funding from the Department of Defense (DOD DAMD17-93-V-3018) to the Georgetown Institute for Cognitive and Computational Sciences. We like to thank the Georgetown University Molecular Modeling Center for the use of its computer facility funded from the National Science Foundation (NSF-CHE-9601976). The SGMD code as implemented in the AMBER and CHARMM programs is available upon request from the authors. References and Notes (1) Dyson, H. J.; Rance, M.; Houghten, R. A.; Lerner, R. A.; Wright, P. E. J. Mol. Biol. 1988, 201, 161. (2) Blanco, F. J.; Jime´nez, M. A.; Herranz, J.; Rico, M.; Santoro, J.; Nieto, J. L. J. Am. Chem. Soc. 1993, 115, 5887. (3) Blanco, F. J.; Jime´nez, M. A.; Pineda, A.; Rico, M.; Santoro, J.; Nieto, J. L. Biochemistry 1994, 33, 6004. (4) Munoz, V.; Eaton, A.; Thompson, P. A.; Hofrichter, J.; Eaton, W. A. Nature 1998, 390, 196. (5) Tobias, D. J.; Mertz, J. E.; Brooks, C. L., III. Biochemistry 1991, 30, 6054. (6) Tobias, D. J.; Sneddon, S. F.; Brooks, C. L., III. J. Mol. Biol. 1992, 227, 1244. (7) Tirado-Rives, J.; Maxwell, D. S.; Jorgesen, W. L. J. Am. Chem. Soc. 1993, 115, 11590.

Wu and Wang (8) Daggett, V.; Levitt, M. J. Mol. Biol. 1992, 223, 1121. (9) Smythe, M. L.; Huston, S. E.; Marshall, G. R. J. Am. Chem. Soc. 1993, 115, 11594-11595. (10) DiCapua, F. M.; Swaminathan, S.; Beveridge, D. L. J. Am. Chem. Soc. 1990, 113, 6145-6155. (11) Sung, S. S.; Wu, X. W. Proteins 1996, 25, 202-214. (12) Sung, S. S.; Wu, X. W. Biopolymers 1997, 42, 633-644. (13) Wu, X. W.; Sung, S. S. Proteins 1999, 34, 295-302. (14) Hirst, J. D.; Brooks, C. L., III. Biochemistry 1995, 34, 7614-7621. (15) Sheinerman, F. B.; Brooks, C. L., III. J. Am. Chem. Soc. 1995, 117, 10098-10103. (16) Young, W. S.; Brooks, C. L., III. J. Mol. Biol. 1996, 259, 558. (17) Shirley, W. A.; Brooks, C. L., III. Proteins 1997, 28, 59. (18) Daura, X.; Jaun, B.; Seebach, D.; van Gunstern,W. F.; Mark, A. E. J. Mol. Biol. 1998, 280, 925-932. (19) Duan, Y.; Kollman, P. A. Science 1998, 282, 740-744. (20) Karplus, M.; Petsko, G. A. Nature 1990, 347, 631. (21) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fisher, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (22) MacKerell, A. D., Jr.; Wiorkiewicz-Kuczera, J.; Karplus, M. J. Am. Chem. Soc. 1995, 117, 11946. (23) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. A. J. Am. Chem. Soc. 1984, 106 (6), 765. (24) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230. (25) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (26) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (27) Lewis, P. N.; Momany, F. A.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1971, 68, 2293. (28) Zimmerman, S. S.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1977, 74, 4126. (29) Wu, X. W.; Wang, S. J. Phys. Chem. B 1998, 102, 7238. (30) Wu, X. W.; Wang, S. J. Chem. Phys. 1999, 110, 9401. (31) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136. (32) Potter, D. Computational Physics; Wiley: New York, 1972; Chapter 5. (33) Benrendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Di Nola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (34) Levitt, M. J. Mol. Biol. 1983, 168, 621. (35) McKelvey, D. R.; Brooks, C. L., III.; Mokotoff, M. J. Protein Chem. 1991, 10, 265. (36) Polinsky, A.; Goodman, M.; Willians, K. A.; Deber, C. M. Biopolymers 1992, 32, 399. (37) Rooman, M. J.; Rodriguez, J.; Wodak, S. J. J. Mol. Biol. 1990, 213, 327. (38) Unger, R.; Harel, D.; Wherland, S.; Sussman, J. L. Proteins 1989, 5, 355. (39) Karpen, M. E.; Tobias, D. J.; Brooks, C. L., III. Curr. Opin. Struct. Biol. 1993, 3, 412. (40) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Jaun, B.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (41) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (42) Hawkins, G. D.; Lynch, G. C.; Giesen, D. J.; Rossi, I.; Storer, J. W.; Liotard, D. A.; Cramer, C. J.; Trular, D. G. QCPE Bull. 1996, 16, 11. (43) Cramer, C. J.; Trular, D. G. J. Am. Chem. Soc. 1991, 113, 8305. (44) Srinivasan, J.; Cheatham, T. E., III; Cieplak, P.; Kollman, P. A.; Case, D. A. J. Am. Chem. Soc. 1998, 120, 3401. (45) Kopple, K. D.; Go, A.; Pilipauskas, D. R. J. Am. Chem. Soc. 1975, 97, 6830. (46) MacArthur, M. W.; Thornton, J. M. J. Mol. Biol. 1991, 218, 397. (47) Yun, R. H.; Anderson, A. G.; Hermans, J. Proteins: Struct,, Funct., and Genet. 1991, 10, 219. (48) Polinsky, A.; Goodman, M.; Willians, K. A.; Deber, C. M. Biopolymers 1992, 32, 399. (49) Hendsch, Z. S.; Sindelar, C. V.; Tidor, B. J. Phys. Chem. B 1998, 102, 4404.