Generalized Verlet Algorithm for Efficient Molecular Dynamics Simulations with Long-range Interactions

For the purpose of molecular dynamics simulations of large biopolymers we have developed a new method to accelerate the calculation of long-range pair interactions (e.g. Coulomb interaction). The algorithm introduces distance classes to schedule updates of non-bonding interactions and to avoid unnecessary computations of interactions between particles which are far apart. To minimize the error caused by the updating schedule, the Verlet integration scheme has been modified. The results of the method are compared to those of other approximation schemes as well as to results obtained by numerical integration without approximation. For simulation of a protein with 12 637 atoms our approximation scheme yields a reduction of computer time by a factor of seven. The approximation suggested can be implemented on sequential as well as on parallel computers. We describe an implementation on a (Transputer-based) MIMD machine with a systolic ring architecture.     

Binding mechanism of spinosine and venenatine molecules with p300 HAT enzyme: Molecular screening,molecular dynamics and free-energy analysis

The chromatin modification is regulated by the histone acetyltransferase (HAT) and histone deacetyltransferase (HDAC) enzymes; abnormal function of these enzymes leads to several malignant diseases. The inhibition of these enzymes using natural ligand molecules is an emerging technique to cure these diseases. The in vitro analysis of natural molecules, venenatine, spinosine, palmatine and taxodione are giving the best inhibition rate against p300 HAT enzyme. However, the detailed understanding of binding and the stability of these molecules with p300 HAT is not yet known. The aim of the present study is focused to determine the binding strength of the molecules from molecular dynamics simulation analysis. The docking analysis confirms that, the venenatine (−6.97 kcal/mol - conformer 8), spinosine (−6.52 kcal/mol conformer −10), palmatine (−5.72 kcal/mol conformer-3) and taxodione (−4.99 kcal/mol conformer-4) molecules form strong hydrogen bonding interactions with the key amino acid residues (Arg1410, Thr1411 and Trp1466) present in the active site of p300. In the molecular dynamics (MD) simulation, the spinosine retain these key interactions with the active site amino acid residues (Arg1410, Thr1411, and Trp1466) than venenatine and are stable throughout the simulation. The RMSD value of spinosine (0.5 to 1.3 Å) and venenatine (0.3 to 1.3 Å) are almost equal during the MD simulation. However, during the MD simulation, the intermolecular interaction between venenatine and the active site amino acid residues (Arg1410, Thr1411, and Trp1466) decreased on comparing with the spinosine-p300 interaction. The binding free energy of the spinosine (−15.30 kcal/mol) is relatively higher than the venenatine (−11.8 kcal/mol); this increment is attributed to the strong hydrogen bonding interactions of spinosine molecule with the active site amino acid residues of p300.     

Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials

The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.     

Investigation of glucose binding sites on insulin

Possible insulin binding sites for D-glucose have been investigated theoretically by docking and molecular dynamics (MD) simulations. Two different docking programs for small molecules were used; Multiple Copy Simultaneous Search (MCSS) and Solvation Energy for Exhaustive Docking (SEED) programs. The configurations resulting from the MCSS search were evaluated with a scoring function developed to estimate the binding free energy. SEED calculations were performed using various values for the dielectric constant of the solute. It is found that scores emphasizing non-polar interactions gave a preferential binding site in agreement with that inferred from recent fluorescence and NMR NOESY experiments. The calculated binding affinity of -1.4 to -3.5 kcal/mol is within the measured range of -2.0 +/- 0.5 kcal/mol. The validity of the binding site is suggested by the dynamical stability of the bound glucose when examined with MD simulations with explicit solvent. Alternative binding sites were found in the simulations and their relative stabilities were estimated. The motions of the bound glucose during molecular dynamics simulations are correlated with the motions of the insulin side chains that are in contact with it and with larger scale insulin motions. These results raise the question of whether glucose binding to insulin could play a role in its activity. The results establish the complementarity of molecular dynamics simulations and normal mode analyses with the search for binding sites proposed with small molecule docking programs.     

The consistency of large concerted motions in proteins in molecular dynamics simulations.

A detailed investigation is presented into the effect of limited sampling time and small changes in the force field on molecular dynamics simulations of a protein. Thirteen independent simulations of the B1 IgG-binding domain of streptococcal protein G were performed, with small changes in the simulation parameters in each simulation. Parameters studied included temperature, bond constraints, cut-off radius for electrostatic interactions, and initial placement of hydrogen atoms. The essential dynamics technique was used to reveal dynamic differences between the simulations. Similar essential dynamics properties were found for all simulations, indicating that the large concerted motions found in the simulations are not particularly sensitive to small changes in the force field. A thorough investigation into the stability of the essential dynamics properties as derived from a molecular dynamics simulation of a few hundred picoseconds is provided. Although the definition of the essential modes of motion has not fully converged in these short simulations, the subspace in which these modes are confined is found to be reproducible.     

Thermodynamic and kinetic determinants of Thermotoga maritima cold shock protein stability: a structural and dynamic analysis

The cold shock protein (CSP) from hyperthermophile Thermotoga maritima (TmCSP) is only marginally stable (DeltaG(T(opt)) = 0.3 kcal/mol) at 353 K, the optimum environmental temperature (T(opt)) for T. maritima. In comparison, homologous CSPs from E. coli (DeltaG(T(opt)) = 2.2 kcal/mol) and B. subtilis (DeltaG(T(opt)) = 1.5 kcal/mol) are at least five times more stable at 310 K, the T(opt) for the mesophiles. Yet at the room temperature, TmCSP is more stable (DeltaG(T(R)) = 4.7 kcal/mol) than its homologues (DeltaG(T(R)) = 3.0 kcal/mol for E. coli CSP and DeltaG(T(R)) = 2.1 kcal/mol for B. subtilis CSP). This unique observation suggests that kinetic, rather than thermodynamic, barriers toward unfolding might help TmCSP native structure at high temperatures. Consistently, the unfolding rate of TmCSP is considerably slower than its homologues. High temperature (600 K) complete unfolding molecular dynamics (MD) simulations of TmCSP support our hypothesis and reveal an unfolding scheme unique to TmCSP. For all the studied homologues of TmCSP, the unfolding process first starts at the C-terminal region and N-terminal region unfolds in the end. But for TmCSP, both the terminals resist unfolding for consistently longer simulation times and, in the end, unfold simultaneously. In TmCSP, the C-terminal region is better fortified and has better interactions with the N-terminal region due to the charged residues, R2, E47, E49, H61, K63, and E66, being in spatial vicinity. The electrostatic interactions among these residues are unique to TmCSP. Consistently, the room temperature MD simulations show that TmCSP is more rigid at its N- and C-termini as compared to its homologues from E. coli, B. subtilis, and B. caldolyticus.     

Calculation of the relative binding free energy of 2'GMP and 2'AMP to ribonuclease T1 using molecular dynamics/free energy perturbation approaches

We present a calculation of the relative changes in binding free energy between the complex of ribonuclease T1 (RNase Tr) with its inhibitor 2'-guanosine monophosphate (2'GMP) and that of RNase T1-2'-adenosine monophosphate (2'AMP) by means of a thermodynamic perturbation method implemented with molecular dynamics. Using the available crystal structure of the RNase T1-2'GMP complex, the structure of the RNase T1-2'AMP complex was obtained as a final structure of the perturbation calculation. The calculated difference in the free energy of binding (delta delta Gbind) was 2.76 kcal/mol. This compares well with the experimental value of 3.07 kcal/mol. The encouraging agreement in delta delta Gbind suggests that the interactions of inhibitors with the enzyme are reasonably represented. Energy component analyses of the two complexes reveal that the active site of RNase T1 electrostatically stabilizes the binding of 2'GMP more than that of 2'AMP by 44 kcal/mol, while the van der Waals' interactions are similar in the two complexes. The analyses suggest that the mutation from Glu46 to Gln may lead to a preference of RNase T1 for adenine in contrast to the guanine preference of the wild-type enzyme. Although the molecular dynamics equilibration moves the atoms of the RNase T1-2'GMP system about 0.9 A from their X-ray positions and the mutation of the G to A in the active site increases the deviation from the X-ray structure, the mutation of the A back to G reduces the deviation. This and the agreement found for delta delta Gbind suggest that the molecular dynamics/free energy perturbation method will be useful for both energetic and structural analysis of protein-ligand interactions.     

Dynamics of fluoroquinolones induced resistance in DNA gyrase of Mycobacterium tuberculosis

DNA gyrase is a validated target of fluoroquinolones which are key components of multidrug resistance tuberculosis (TB) treatment. Most frequent occurring mutations associated with high level of resistance to fluoroquinolone in clinical isolates of TB patients are A90V, D94G, and A90V–D94G (double mutant [DM]), present in the larger subunit of DNA Gyrase. In order to explicate the molecular mechanism of drug resistance corresponding to these mutations, molecular dynamics (MD) and mechanics approach was applied. Structure-based molecular docking of complex comprised of DNA bound with Gyrase A (large subunit) and Gyrase C (small subunit) with moxifloxacin (MFX) revealed high binding affinity to wild type with considerably high Glide XP docking score of ?7.88 kcal/mol. MFX affinity decreases toward single mutants and was minimum toward the DM with a docking score of ?3.82 kcal/mol. Docking studies were also performed against 8-Methyl-moxifloxacin which exhibited higher binding affinity against wild and mutants DNA gyrase when compared to MFX. Molecular Mechanics/Generalized Born Surface Area method predicted the binding free energy of the wild, A90V, D94G, and DM complexes to be ?55.81, ?25.87, ?20.45, and ?12.29 kcal/mol, respectively. These complexes were further subjected to 30 ns long MD simulations to examine significant interactions and conformational flexibilities in terms of root mean square deviation, root mean square fluctuation, and strength of hydrogen bond formed. This comparative drug interaction analysis provides systematic insights into the mechanism behind drug resistance and also paves way toward identifying potent lead compounds that could combat drug resistance of DNA gyrase due to mutations.     

High performance computing in biology: multimillion atom simulations of nanoscale systems

Computational methods have been used in biology for sequence analysis (bioinformatics), all-atom simulation (molecular dynamics and quantum calculations), and more recently for modeling biological networks (systems biology). Of these three techniques, all-atom simulation is currently the most computationally demanding, in terms of compute load, communication speed, and memory load. Breakthroughs in electrostatic force calculation and dynamic load balancing have enabled molecular dynamics simulations of large biomolecular complexes. Here, we report simulation results for the ribosome, using approximately 2.64 million atoms, the largest all-atom biomolecular simulation published to date. Several other nano-scale systems with different numbers of atoms were studied to measure the performance of the NAMD molecular dynamics simulation program on the Los Alamos National Laboratory Q Machine. We demonstrate that multimillion atom systems represent a 'sweet spot' for the NAMD code on large supercomputers. NAMD displays an unprecedented 85% parallel scaling efficiency for the ribosome system on 1024 CPUs. We also review recent targeted molecular dynamics simulations of the ribosome that prove useful for studying conformational changes of this large biomolecular complex in atomic detail.     

Effects of deformability and thermal motion of lipid membrane on electroporation: by molecular dynamics simulations

Effects of mechanical properties and thermal motion of POPE lipid membrane on electroporation were studied by molecular dynamics simulations. Among simulations in which specific atoms of lipids were artificially constrained at their equilibrium positions using a spring with force constant of 2.0 kcal/(mol Ų) in the external electric field of 1.4 kcal/(mol Å e), only constraint on lateral motions of lipid tails prohibited electroporation while non-tail parts had little effects. When force constant decreased to 0.2 kcal/(mol Ų) in the position constraints on lipid tails in the external electric field of 2.0 kcal/(mol Å e), water molecules began to enter the membrane. Position constraints of lipid tails allow water to penetrate from both sides of membrane. Thermal motion of lipids can induce initial defects in the hydrophobic core of membrane, which are favorable nucleation sites for electroporation. Simulations at different temperatures revealed that as the temperature increases, the time taken to the initial pore formation will decrease.     

Computational aided mechanistic understanding of Camellia sinensis bioactive compounds against co-chaperone p23 as potential anticancer agent

Co-chaperon p23 has been well established as molecular chaperon for the heat shock protein 90 (Hsp90) that further leads to immorality in cancer cells by providing defense against Hsp90 inhibitors, and as stimulating agent for generating overexpressed antiapoptotic proteins, that is, Hsp70 and Hsp27. The natural compounds such as catechins from Camellia sinensis (green tea) are also well known for inhibition activity against various cancer. However, molecular interaction profile and potential lead bioactive compounds against co-chaperon p23 from green tea are not yet reported. To this context, we study the various secondary metabolites of green tea against co-chaperon p23 using structure-based virtual screening from Traditional Chinese Medicine (TCM) database. Following 26 compounds were obtained from TCM database and further studied for extra precision molecular docking that showed binding score between −10.221 and −2.276 kcal/mol with co-chaperon p23. However, relative docking score to known inhibitors, that is, ailanthone (−4.54 kcal/mol) and gedunin ( 3.60 kcal/mol) along with ADME profile analysis concluded epicatechin (−7.013 kcal/mol) and cis-theaspirone (−4.495 kcal/mol) as potential lead inhibitors from green tea against co-chaperone p23. Furthermore, molecular dynamics simulation and molecular mechanics generalized born surface area calculations validated that epicatechin and cis-theaspirone have significantly occupied the active region of co-chaperone p23 by hydrogen and hydrophobic interactions with various residues including most substantial amino acids, that is, Thr90, Ala94, and Lys95. Hence, these results supported the fact that green tea contained potential compounds with an ability to inhibit the cancer by disrupting the co-chaperon p23 activity.     

Molecular dynamics simulations of charged and neutral lipid bilayers: treatment of electrostatic interactions

Molecular dynamics (MD) simulations complement experimental methods in studies of the structure and dynamics of lipid bilayers. The choice of algorithms employed in this computational method represents a trade-off between the accuracy and real calculation time. The largest portion of the simulation time is devoted to calculation of long-range electrostatic interactions. To speed-up evaluation of these interactions, various approximations have been used. The most common ones are the truncation of long-range interactions with the use of cut-offs, and the particle-mesh Ewald (PME) method. In this study, several multi-nanosecond cut-off and PME simulations were performed to establish the influence of the simulation protocol on the bilayer properties. Two bilayers were used. One consisted of neutral phosphatidylcholine molecules. The other was a mixed lipid bilayer consisting of neutral phosphatidylethanolamine and negatively charged phosphatidylglycerol molecules. The study shows that the cut-off simulation of a bilayer containing charge molecules generates artefacts; in particular the mobility and order of the charged molecules are vastly different from those determined experimentally. In the PME simulation, the bilayer properties are in general agreement with experimental data. The cut-off simulation of bilayers containing only uncharged molecules does not generate artefacts, nevertheless, the PME simulation gives generally better agreement with experimental data.     

Molecular Dynamics/Free Energy Perturbation Studies of the Thermostable V74I Mutant of Ribonuclease HI

Abstract

Recent site-directed mutagenesis and thermodynamic studies have shown that the V74I mutant of Escherichia coli ribonuclease HI (RNase HI) is more stable than the wild type protein [Ishikawa et al., Biochemistry 32, 6171 (1993)]. In order to clarify the stabilization mechanism of this mutant, we calculated the free energy change due to the mutation Val 74→Ile in both the native and denatured states by free energy perturbations based on molecular dynamics (MD) simulations. We carried out inclusive MD simulations for the protein in water; i.e., fully solvated, no artificial constraints applied, and all long-range Coulomb interactions included. We found that the free energy of the mutant increased slightly relative to the wild type, in the native state by 1.60 kcal/mol, and in the denatured state by 2.25 kcal/mol. The unfolding free energy increment of the mutant (0.66 ± 0.19 kcal/mol) was in good agreement with the experimental value (0.6 kcal/mol). The hysteresis error in the free energy calculations, i.e., forward and reverse perturbations, was only ±0.19 kcal/mol. These results show that the V74I mutant is stabilized relative to the wild type by the increased free energy of the denatured state and not by a decrease in the free energy of the native state as had been proposed earlier based on the mutant X-ray structure. It was found that the stabilization was caused by a loss of solvation energy in the mutant denatured state and not by improved packing interactions inside the native protein.     

Recognition Mechanism of siRNA by Viral p19 Suppressor of RNA Silencing: A Molecular Dynamics Study

The p19 protein (p19) encoded from Tombusvirus is involved in various activities such as pathogenicity and virus transport. Recent studies have found that p19 is a plant suppressor of RNA silencing, which binds to short interfering RNAs (siRNAs) with high affinity. We use molecular dynamics (MD) simulations of the wild-type and mutant p19 protein (W39 and W42G) binding with a 21-nt siRNA duplex to study the p19-siRNA recognition mechanism and mutation effects. Our simulations with standard MD and steered molecular dynamics have shown that the double mutant structure is indeed much less stable than the wild-type, consistent with the recent experimental findings. Comprehensive structural analysis also shows that the W39/42G mutations first induce the loss of stacking interactions between p19 and siRNA, Trp⁴²-Cyt1 (Cyt1 from the 5′ to 3′ strand) and Trp³⁹-Gua′19 (Gua19 from the 3′ to 5′ strand), and then breaks the hydrophobic core formed by W39-W42 with nucleotide basepairs in the wild-type. The steered molecular dynamics simulations also show that the mutant p19 complex is “decompounded” very fast under a constant separation force, whereas the wild-type remains largely intact under the same steering force. Moreover, we have used the free energy perturbation to predict a binding affinity loss of 6.98 ± 0.95 kcal/mol for the single mutation W39G, and 12.8 ± 1.0 kcal/mol loss for the double mutation W39/42G, with the van der Waals interactions dominating the contribution (∼90%). These results indicate that the W39/42G mutations essentially destroy the important p19-siRNA recognition by breaking the strong stacking interaction between Cyt1 and Gua′19 with end-capping tryptophans. These large scale simulations might provide new insights to the interactions and co-evolution relationship between RNA virus proteins and their hosts.     

Conformational equilibria of valine studied by dynamics simulation.

The conformational probability distribution of a valine residue in the valine dipeptide and of the valine side chain in an alpha-helix, as well as the change in helix stability for replacing alanine with valine, has been calculated by molecular dynamics simulations of explicitly hydrated systems: dipeptide, tetrapeptide and 10-, 14- and 18-residue oligoalanine helices. All computed free-energy differences are means from at least eight separate slow-growth simulations, four in each direction and are reported with their root-mean-square deviations. Different values for the change in free energy of folding (delta delta G degrees) have been calculated with the use of forcefields having an all-atom and a central-atom representation of methyl groups, etc. The value obtained with the all-atom forcefield agrees well with new experimental values (3 kJ/mol = 0.7 kcal/mol). Furthermore, the most stable valine side-chain rotamer in the helix is different for these two representations. The most stable rotamer for the all atom conformation is the same one that predominates for valines in alpha-helices in proteins of known conformation. The lower conformational freedom of the valine side chain in the helix contributes 1 kJ/mol to the difference in stability computed with the all-atom potential; unfavorable interactions of the side chain with helix, even in the most stable conformation, further increase delta delta G degrees.     

Elastic properties, Young's modulus determination and structural stability of the tropocollagen molecule: a computational study by steered molecular dynamics

The aim of this report is to investigate at microscopic level the elastic properties of a tropocollagen-like molecule submitted to linear traction along its longitudinal axis. For this purpose, we performed steered molecular dynamics (SMD) simulations for a wide range of spring constants in order to test the molecular response based on a two-spring model connected in series. An elastic behavior was observed in an elongation range of 2.5-4% of the molecular length, estimating an "effective molecular elastic constant" of 1.02+/-0.20 kcal/mol A2 in this region. Accordingly, a Young's modulus for the tropocollagen molecule of Y=4.8+/-1.0 GPa was calculated. The complex hydrogen bond network was traced along molecular dynamics (MD) and SMD simulations revealing a rearrangement of these interactions preserving the integrity of the molecular structure when submitted to traction. No evidence of the significant role attributed to water bridges for structural stability was detected, on the contrary facts pointed out that the hydrogen bond network might be the responsible.     

The molecular basis for the inhibition of human cytochrome P450 1A2 by oroxylin and wogonin

In our previous kinetics studies the natural products oroxylin and wogonin were shown to have strong biological affinity for, and inhibitory effects against, human cytochrome P450 1A2, with IC₅₀ values of 579 and 248 nM, respectively; this might lead to the occurrence of drug–drug interactions when co-administered clinically. However, their inhibitory mechanisms against 1A2 remain elusive. In this study, molecular docking and molecular dynamics simulations were performed to better understand the molecular basis of their inhibitory mechanisms towards 1A2. Structural analysis revealed that oroxylin has a different binding pattern from wogonin and another very strongly binding inhibitor α-naphthoflavone (ANF, IC₅₀ = 49 nM). The O₇ atom of oroxylin forms hydrogen bonds with the OD1/OD2 atoms of Asp313, which is not observed in the 1A2–wogonin complex. Because of energetically unfavorable repulsions with the methoxy group at the 6 position of the oroxylin ring, significant conformational changes were observed for the sidechain of Thr118 in the MD simulated model. As a result, the larger and much more open binding-site architecture of the 1A2–oroxylin complex may account for its weaker inhibitory effect relative to the 1A2–ANF complex. Energy analysis indicated that oroxylin has a less negative predicted binding free energy of −19.8 kcal/mol than wogonin (−21.1 kcal/mol), which is consistent with our experimental assays. Additionally, our energy results suggest that van der Waals/hydrophobic and hydrogen-bonding interactions are important in the inhibitory mechanisms of oroxylin whereas the former is the underlying force responsible for strong inhibition by ANF and wogonin.     

The Fe-CO bond energy in myoglobin: a QM/MM study of the effect of tertiary structure

The Fe-CO bond dissociation energy (BDE) in myoglobin (Mb) has been calculated with B3LYP quantum mechanics/molecular mechanics methods for 22 different Mb conformations, generated from molecular dynamics simulations. Our average BDE of 8.1 kcal/mol agrees well with experiment and shows that Mb weakens the Fe-CO bond by 5.8 kcal/mol; the calculations provide detailed atomistic insight into the origin of this effect. BDEs for Mb conformations with the R carbonmonoxy tertiary structure are on average 2.6 kcal/mol larger than those with the T deoxy tertiary structure, suggesting two functionally distinct allosteric states. This allostery is partly explained by the reduction in distal cavity steric crowding as Mb moves from its T to R tertiary structure.     

Equilibrium structure and folding of a helix-forming peptide: circular dichroism measurements and replica-exchange molecular dynamics simulations

We have performed experimental measurements and computer simulations of the equilibrium structure and folding of a 21-residue alpha-helical heteropeptide. Far ultraviolet circular dichroism spectroscopy is used to identify the presence of helical structure and to measure the thermal unfolding curve. The observed melting temperature is 296 K, with a folding enthalpy of -11.6 kcal/mol and entropy of -39.6 cal/(mol K). Our simulations involve 45 ns of replica-exchange molecular dynamics of the peptide, using eight replicas at temperatures between 280 and 450 K, and the program CHARMM with a continuum solvent model. In a 30-ns simulation started from a helical structure, conformational equilibrium at all temperatures was reached after 15 ns. This simulation was used to calculate the peptide melting curve, predicting a folding transition with a melting temperature in the 330-350 K range, enthalpy change of -10 kcal/mol, and entropy change of -30 cal/(mol K). The simulation results were also used to analyze the peptide structural fluctuations and the free-energy surface of helix unfolding. In a separate 15-ns replica-exchange molecular dynamics simulation started from the extended structure, the helical conformation was first attained after approximately 2.8 ns, and equilibrium was reached after 10 ns of simulation. These results showed a sequential folding process with a systematic increase in the number of hydrogen bonds until the helical state is reached, and confirmed that the alpha-helical state is the global free-energy minimum for the peptide at low temperatures.