Journal of Chemistry

Volume 2013 (2013), Article ID 407862, 13 pages

http://dx.doi.org/10.1155/2013/407862

Research Article

Toward Structure Prediction for Short Peptides Using the Improved SAAP Force Field Parameters

Kenichi Dedachi,¹ Taku Shimosato,¹ Toshiya Minezaki,² and Michio Iwaoka¹

¹Department of Chemistry, School of Science, Tokai University, Kitakaname, Hiratsuka-shi, Kanagawa 259-1292, Japan

²Laboratory of General Education for Science and Technology, School of Science, Tokai University, Kitakaname, Hiratsuka-shi, Kanagawa 259-1292, Japan

Received 9 November 2012; Revised 5 January 2013; Accepted 9 January 2013

Academic Editor: Yifat Miller

Copyright © 2013 Kenichi Dedachi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Based on the observation that Ramachandran-type potential energy surfaces of single amino acid units in water are in good agreement with statistical structures of the corresponding amino acid residues in proteins, we recently developed a new all-atom force field called SAAP, in which the total energy function for a polypeptide is expressed basically as a sum of single amino acid potentials and electrostatic and Lennard-Jones potentials between the amino acid units. In this study, the SAAP force field (SAAPFF) parameters were improved, and classical canonical Monte Carlo (MC) simulation was carried out for short peptide models, that is, Met-enkephalin and chignolin, at 300 K in an implicit water model. Diverse structures were reasonably obtained for Met-enkephalin, while three folded structures, one of which corresponds to a native-like structure with three native hydrogen bonds, were obtained for chignolin. The results suggested that the SAAP-MC method is useful for conformational sampling for the short peptides. A protocol of SAAP-MC simulation followed by structural clustering and examination of the obtained structures by ab initio calculation or simply by the number of the hydrogen bonds (or the hardness) was demonstrated to be an effective strategy toward structure prediction for short peptide molecules.

1. Introduction

In conformational analysis of short peptides, Monte Carlo (MC) and molecular dynamics (MD) simulation techniques have been widely applied [1, 2], in which a set of potential energy functions (a so-called force field), such as ECEPP [3, 4], AMBER [5–8], CHARMM [9], OPLS [10, 11], and GROMOS [12], is employed to define the relation between the structure and the potential energy. However, it is still not practical to search for the possible conformers comprehensively within a limited computing time without using advanced sampling methods, such as multicanonical [13, 14] and replica-exchange [15–17] methods, which have been developed to overcome the sampling problem. In recent years, various types of coarse-grained force fields, such as AMBER-UA [18], UNRES [19–23], and MARTINI [24, 25], have also been developed in order to reduce time to compute the potential energy. For instance, the UNRES force field successfully simulated the aggregation process and stable structure of β-amyloid protein oligomers that are regarded as pathologic factors of Alzheimer’s disease [26].

On the other hand, we recently discovered the interesting feature of protein structures that Ramachandran-type potential energy surfaces of single amino acid units in water obtained by ab initio molecular orbital calculation are almost identical with statistical distributions of the Ramachandran plots obtained from protein databank (PDB) [27, 28]. In other words, the amino acid residues in proteins seem to statistically follow Boltzmann distributions on the single amino acid potential (SAAP) surfaces in water. Although physicochemical implications of this unexpected similarity are not yet clear, the finding strongly suggests prominent importance of SAAP as a determinant of protein structures. This point would also be supported by the previous experimental result by Dobson and coworkers [29, 30] that and dihedral angle distributions of individual amino acid residues of a polypeptide in a random-coil state can be expressed by use of the Ramachandran plots, as well as by the theoretical approaches by Sakae and Okamoto [31] and Kamiya et al. [32] to the optimization of conventional all-atom force field parameters by using the Ramachandran plots.

A discovery of the similarity between the SAAP in water and the statistical structure of the amino acid residues in folded proteins has prompted us to develop a new force field called SAAP for polypeptide molecules [33, 34]. The SAAP force field (SAAPFF) is entirely different from conventional force fields in that a polypeptide is divided to the amino acid units, not to the atomic units, solvent effects are implicitly included in the parameters, and the atomic charges are not constant but variable depending on the conformation of the amino acid units. In the SAAPFF, the total potential energy () for a polypeptide is expressed by (1) as the following:

where is a sum of potential energies (SAAP) for the individual amino acid units, and are electrostatic and Lennard-Jones potentials between the amino acid units, and is the other correlation term, which is ignored in a current version of SAAPFF. Detailed description of these terms and the simulation program was given previously (see also Supplementary Material available online at http://dx.doi.org/10.1155/2013/407862) [33, 34]. Canonical MC simulation using the SAAPFF reasonably reproduced randomly fluctuating conformation of Met-enkephalin. However, for chignolin the SAAP-MC simulation did not afford the native β-hairpin structure when the simulation was started from the extended structure [34]. Thus, the previous SAAPFF parameters were a little too rough to be applied for structure prediction of the short peptide. In this study, the SAAPFF parameters have been improved in several points, and accuracy of the improved SAAP simulation suite to search for possible conformers of short peptides has been evaluated by performing the molecular simulation for Met-enkephalin and chignolin. The results were then compared with those obtained by the conventional AMBER-MD simulation method using the amber99sb force field [1] combined with the generalized Born solvent model [35–37]. A protocol of SAAP-MC simulation followed by structural clustering is suggested to be an effective strategy toward structure prediction for short peptide molecules.

2. Method of Calculations

2.1. Model Peptides

Met-enkephalin and chignolin without N- and C-terminal protecting groups were employed as model short peptides. Met-enkephalin is a short peptide, which consists of five residues (Tyr–Gly–Gly–Phe–Met) as shown in Figure 1. According to previous experimental [38] and theoretical [39–50] studies, the structure of this short peptide randomly fluctuates in water, while the presence of some characteristic conformations with intramolecular hydrogen bonds was indicated in vacuo [51–53]. Met-enkephalin is a peptide widely used for testing performance of developed molecular simulation methodologies.

Figure 1: Structure of Met-enkephalin.

Chignolin has a unique folded β-hairpin structure in water, although it is a short peptide consisting of ten amino acid residues (Gly–Tyr–Asp–Pro–Glu–Thr–Gly–Thr–Trp–Gly) [54]. Molecular simulation for this miniprotein has been performed actively by using AMBER and OPLS force fields [55–64]. As shown in Figure 2, chignolin folds around Pro4–Glu5–Thr6 to form a β-hairpin structure, which is stabilized by three intramolecular hydrogen bonds, Gly7(N)Asp3(O) (H-bond I), Thr8(N)Asp3(O) (H-bond II), and Asp3(N)Thr8(O) (H-bond III), and a hydrophobic interaction between the aromatic rings of Tyr2 and Trp9.

Figure 2: Structure of chignolin. Three native hydrogen bonds (H-bonds I, II, and III) and a hydrophobic interaction between Tyr2 and Trp9 are indicated by arrows [58].

2.2. Improvement of the SAAPFF Parameters

The SAAPFF parameters are comprised of a potential energy, atomic coordinates, atomic charges, and Lennard-Jones potential parameters for all possible conformations of the single amino acid unit (HCO-Xaa-NH₂). The conformations are defined by dihedral angles, , , , , and so forth, each with an interval of 15 degrees. In this study, the SAAPFF parameters previously developed [33, 34] were improved in the following points.

First, the parameters for main-chain, Pro, and Val in water, which were obtained in the previous version [34] by single-point calculation with the IEFPCM [65–67] solvent model using the structures optimized in vacuo, were replaced with those obtained by geometry optimization in water at the HF/6-31+G(d) level using the IEFPCM model in a Gaussian03 program (rev. E.01) [68]. The similar calculation level has been widely applied for small biomolecules [69, 70]. Relaxation of the potential surfaces in the implicit water model should improve accuracy of the SAAPFF.

Second, atomic charges of the SAAPFF parameters, which are variables as a function of the structure of an amino acid unit, were switched for all amino acids from Mulliken charges to electrostatic potential (ESP) charges [71] that were obtained at the recommended higher MP2/6-31G(d,p) level [72] using the IEFPCM solvent model. The ESP charges are defined to reproduce the electrostatic potential surface and are used in conventional force fields as reliable atomic charges. Therefore, the use of ESP charges would reproduce electrostatic interactions, that is, the term, more appropriately than the Mulliken charges.

Third, the parameters of potential energies for the main-chain unit were further improved as follows. The SAAPFF parameters for most amino acids, except for Gly, Ala, Pro, and Val, are divided into the main-chain and side-chain units according to the side-chain separation approximation method [34]. In this approximation, the main-chain unit was capped with a tert-butyl group, which was fixed eclipsed to the main-chain part during geometry optimization to obtain the SAAPFF parameters. However, this method brought a statistical error that the SAAP energies thus obtained are slightly larger than those obtained without the approximation, due probably to distortion of the bond angles derived from the fixed geometry. Indeed, such a statistical error was evident for the case of Val [34]. Therefore, the energies were replaced with those calculated for the relaxed structures, which were obtained by geometry optimization with the tert-butyl group relaxed in water, while for structural parameters (i.e., atomic coordinates and atomic charges) the geometry with the tert-butyl group fixed eclipsed was still employed in order to avoid abnormally short atomic contacts between the main-chain and side-chain units after the connection.

2.3. Simulation Conditions

Conformational properties of Met-enkephalin and chignolin in water were studied by SAAP-MC and AMBER-MD [73] simulations. The SAAP-MC simulation was carried out at 300 K in water by using the improved SAAPFF parameters, in which solvation effects were implemented implicitly. The conventional Metropolis method [74, 75] and the Mersenne Twister random number generator [76] were employed for the MC simulation. At each MC step, one dihedral angle, randomly chosen from all dihedral angles including main-chain and side-chain dihedral angles, was changed randomly with a maximum displacement angle of ±32 degree. The new structure was constructed by the interpolation of the SAAPFF parameters and subsequent connection of the single amino acid units, and the potential energy was calculated by (1). The potential energy thus obtained was then compared with the previous one. Ten trajectories with total MC steps of 100 million were obtained for Met-enkephalin, while twenty trajectories with total MC steps of 2 billion were obtained for chignolin. As an initial structure of Met-enkephalin for the SAAP-MC simulation, the extended structure from PDB (PDB ID: 1PLW) [77] was employed. On the other hand, extended (defined by all main-chain and dihedral angles of −180 degree) and native folded (PDB ID: 1UAO) [54] structures were selected for chignolin as the initial structures. An Intel Xeon W3550 3.06 GHz processor with 12 GB memory was employed as a platform for calculation. The computing time of SAAP-MC simulation for Met-enkephalin was about 6 hours and for chignolin was about 11 days.

For comparison, ten trajectories of AMBER-MD [73] simulation were also calculated for Met-enkephalin and chignolin, respectively, in water at 300 K. The simulation time was 1 s with a time step of 2 fs, and a cut-off distance of nonbonded interactions was 12 Å. Solvation effects were implicitly considered by using the generalized Born model () [35–37] with a dielectric constant of 78.5, following the recent study by Götz et al. [78]. The AMBER10 program [73] with the amber99sb force field [1] was used. The computing time of the AMBER-MD simulation to obtain one trajectory was about 165 hours for Met-enkephalin, and that for chignolin was about 18 days, on the same calculation platform used for the SAAP-MC simulation.

2.4. Data Analysis

The 10,000 structures of Met-enkephalin that were extracted in every 10,000 step from the SAAP-MC simulation trajectory were classified to twenty structural clusters by using a clustering algorithm called the -means method [79] based on the RMSD for the main-chain C, O, and N atoms or all heavy atoms except for the hydrogen atoms. Similarly, ten structural clusters were obtained for chignolin from the 20,000 structures extracted in every 100,000 step from the SAAP-MC simulation trajectories. On the other hand, 10,000 structures were extracted in every 100 ps from the AMBER-MD trajectories for both Met-enkephalin and chignolin. The obtained structures were classified to ten or twenty clusters by application of the -means method under the same conditions to those applied for the SAAP-MC simulation results.

For chignolin, the RMSD values with respect to the native structure were also calculated for the structures extracted from the SAAP-MC and AMBER-MD trajectories by using the amber-ptraj program [7, 8] based on the main-chain atoms or all heavy atoms. The free-energy surfaces of chignolin were further analyzed by using all trajectories obtained from the SAAP-MC simulation according to (2):

where is a probability of the structures with the corresponding structural parameters, is the probability of the structure that had the lowest energy, is a difference in the free energy, is a gas constant, and is an absolute temperature. Ramachandran-type free-energy surfaces were obtained for each amino acid residue by dividing the surface into 10 × 10 degree units, while main-chain RMSD versus hydrogen bond or hydrogen bond versus hydrogen bond free-energy surfaces were obtained by dividing the surface into 0.1 × 0.1 Å units.

2.5. Ab Initio Calculation

Relative energies of the representative structures obtained for chignolin by the SAAP-MC simulation were calculated in water by the ab initio molecular orbital method. The calculation was performed at the HF/IEFPCM/6-31+G(d,p) level by using a Gaussian03 program [68]. The polarizable continuum model using the integral equation formalism variant (IEFPCM), which is equipped as the default self-consistent reaction field method, was employed for the calculation in water.

3. Results

3.1. Met-Enkephalin by SAAP-MC and AMBER-MD

Energetic trajectories of SAAP-MC and AMBER-MD simulations for Met-enkephalin were obtained in an implicit water model at 300 K (see Figure S1). In SAAPFF, the total energy is basically expressed as a sum of , , and . The values of these energetic terms were maintained stable during the SAAP-MC simulation with large fluctuation for and , while the value of was almost zero, as observed previously [34]: the large dielectric constant of water would make the electrostatic interaction between the amino acid residues ignorable. Similarly, the total energy was maintained stable for the AMBER-MD simulation. However, the distribution of the distance between the terminal Cα atoms () was significantly different for the two simulation methods. The value of dispersed in a range from 5 to 12 Å in the SAAP-MC simulation, while those converged at around 6 Å in the AMBER-MD simulation (Figure 3). The results suggested that the structures obtained from the SAAP-MC simulation contain diverse conformations and are more extended than those obtained from the AMBER-MD simulation on the average. Indeed, when the structures were classified into twenty clusters based on the RMSD values for the main-chain atoms, various types of structures were obtained from the SAAP-MC trajectory. The representative examples are shown in Figure 4. The structures are different from each other with a diverse value, and there is no hydrogen bond in the three representative structures. In contrast, the structures of Met-enkephalin obtained from the AMBER-MD trajectory seemed to be trapped in local energy minimums during the simulation (see Figure S2).

(a)

(b)

Figure 3: Histograms of the distance between Cα atoms of Tyr1 and Met5 for 10,000 structures of Met-enkephalin obtained by the SAAP-MC (a) and AMBER-MD (b) simulations at 300 K in water. The results from ten trajectories were averaged. See the text for details of simulation conditions.

(a)

(b)

(c)

(d)

Figure 4: Representative structures of Met-enkephalin obtained by the SAAP-MC simulation at 300 K in water (trajectory 1). The obtained 10,000 structures were analyzed by a structural clustering algorithm using the -means method based on the main-chain RMSD. A distance of the hydrogen bond is given in a unit of Å.

3.2. Chignolin by SAAP-MC and AMBER-MD

Twenty trajectories were obtained for chignolin by SAAP-MC simulation. The structures in each trajectory were classified to ten clusters based on the main-chain RMSD by using the -means method. According to the main-chain folds, that is, combination of the dihedral angles for the amino acid residues, the obtained clusters were further classified to three representative structures, which are called here A, B, and C (Figure 5). Populations of the three structures are summarized in Table 1. Structure A has a native-like main-chain fold stabilized by three hydrogen bonds, which are the same as those observed in the native structure (i.e., H-bonds I–III). However, the hydrophobic interaction between the aromatic side chains of Tyr2 and Trp9 is not present. Structures B and C have nonnative conformation with a common type II β-turn around Thr6 and Gly7. A weak hydrogen bond exists between Glu5(O) and Thr8(N) with a distance about 4.0 Å. However, the conformations of the C-terminals are different from each other. Other minor structures have various local and global conformations, such as a type I β-turn, a helical conformation, and other conformations.

Table 1: Existence ratios (%) of representative structures of chignolin obtained by the SAAP-MC simulation at 300 K in water.

Trajectory number
A ()
B
C
Other

1
58.5 (11.2)
25.9
2.3
13.3

2
54.7 (11.3)
26.6
6.6
12.1

3
50.4 (8.4)
42.4
3.8
3.4

4
49.3 (11.4)
30.7
2.9
17.1

5
45.4 (10.6)
37.8
6.1
10.7

6
44.3 (8.9)
47.9
0.0
7.8

7
37.9 (19.5)
40.2
2.4
19.5

8
37.6 (8.5)
27.9
0.0
34.5

9
36.5 (11.2)
30.7
6.2
26.6

10
34.7 (6.6)
27.3
4.2
33.8

11
32.1 (8.0)
32.9
5.5
29.5

12
31.2 (9.7)
23.5
0.0
45.3

13
19.8 (3.5)
43.7
5.3
31.2

14
14.2 (2.5)
31.8
0.0
54.0

15
13.6 (0.0)
56.8
6.3
23.3

16
12.8 (2.6)
60.0
7.9
19.3

17
10.5 (0.0)
59.3
0.0
30.2

18
1.1 (0.0)
65.1
0.0
33.8

19
0.0 (0.0)
70.4
5.5
24.1

20
0.0 (0.0)
83.2
10.2
6.6

Average
29.2 (6.7)
43.2
3.8
23.8

Figure 5: Representative structures of chignolin obtained by the SAAP-MC simulation at 300 K in water. The obtained 20,000 structures were analyzed by a structural clustering algorithm using the -means method. Structures (A), (B), and (C) were obtained based on the main-chain RMSD. Structure (A′) was obtained based on the all-atom RMSD. Distances of the hydrogen bonds and the hydrophobic interaction are given in a unit of Å.

In the trajectories from 1 to 12, the existence ratio of structure A was high (31.2 to 58.5%), whereas the ratio was much lower than that of the misfolded structure B in the remaining trajectories. The mean ratio for structure A averaged for all trajectories was 29.2%. On the other hand, when the clustering analysis was performed based on all-atom RMSD, the native-like structure (A′) was obtained from 0 to 19.5% (Table 1). The existence ratio of structure A′ averaged over all trajectories was 6.7%. Structure A′ maintains three native H-bonds as well as a hydrophobic interaction between Tyr2 and Trp9 (Figure 5). The all-atom and main-chain RMSD values of structure A′ from the native structure were 2.6 and 2.0 Å, respectively, suggesting that structure A′ keeps the native fold. The trace of the main-chain RMSD with respect to the native structure is shown in Figure 6(a) for the case of trajectory 10, where structures A and B are involved in almost equal amounts. The RMSD values fluctuated between 1 and 6 Å, and repeats of the structural transitions were observed. Similarly, in the AMBER-MD simulation, the main-chain RMSD value fluctuated between 1 and 6 Å (Figure 6(b)). The structural clustering analysis showed that the native structure with three native hydrogen bonds was involved about 30% (Figure S5).

(a)

(b)

Figure 6: Traces of the main-chain RMSD obtained for chignolin by the SAAP-MC (trajectory 10) (a) and AMBER-MD (trajectory 1) (b) simulations at 300 K in water. The RMSD values were calculated with respect to the native structure.

The convergence of the SAAP-MC simulation is not clear from Table 1. However, the similar relative populations were obtained for structures A (~30% including 6% of A′) and B (~60%) when the simulation was started from the native structure. Moreover, in our preliminary results of the SAAP-MC simulation using ten trajectories, almost the same mean populations were obtained for the representative structures (data not shown). Based on these observations, we employed the twenty trajectories shown in Table 2 for further analysis.

Table 2: Relative energies of structures A–C determined by ab initio calculation at the HF/IEFPCM/6-31+G(d,p) level in water.

Structures
SCF energy

(a.u.)
Relative energy (kcal/mol)

A
−3800.42798
5.52

−3800.43678
0.00

B
−3800.39878
23.85

C
−3800.43002
4.42

The structure with the smallest all-atom RMSD (1.3 Å) obtained from the twenty trajectories is superimposed on the native structure in Figure 7. The main-chain and Cα atom RMSD values of this structure with respect to the native structure were 0.9 and 0.5 Å, respectively, confirming that the structure is almost identical to the native structure.

Figure 7: The structure of chignolin with the smallest all-atom RMSD (1.3 Å) obtained by the SAAP-MC simulation at 300 K in water (white) superimposed on the reference native structure (gray).

The Ramachandran-type free-energy surfaces obtained for Tyr2–Trp9 residues from all trajectories of the SAAP-MC simulation are shown in Figure 8. The profiles for Tyr2, Pro4, Glu5, and Thr6 fit well to the native structures, while the free-energy minimums for Asp3, Gly7, and Thr8 are slightly moved from the positions of the native structures. For Trp9, the stable areas corresponding to α and β regions tend to merge. This is roughly consistent with diversity in the C-terminal conformation of the native structure. The dihedral angles for structures A and A′ are located around the energy minimum points on the free-energy surfaces. In contrast, structure B has different geometry from structures A and A′ for Glu5 and Thr6: the angles are ca. + degree. Structure C has the geometry with angles of Glu5, Thr6, Thr8, and Trp9 different from those of structures A and A′. During the SAAP-MC simulation, mismatch of the main-chain dihedral angles for Glu5, Thr6, and Gly7 to the native structure occurred at high ratios.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8: Ramachandran-type free-energy surfaces for Tyr2–Trp9 residues of chignolin obtained from all trajectories of the SAAP-MC simulation at 300 K in water along with the plots of the native structures determined by NMR [58]. Contour lines are drawn in an interval of 1 kcal/mol.

The free-energy surface of chignolin projected on H-bonds I versus III plane and that projected on the main-chain RMSD versus H-bond II plane obtained from all trajectories are shown in Figure 9. In Figure 9(a), the most stable area was found at  Å and  Å. This energy minimum corresponds to the misfolded structure B. Another potential energy minimum was found at  Å and  Å, corresponding to structures A and A′. Similarly, two stable minimums were located in Figure 9(b);  Å and  Å, corresponding to structure B, and  Å and  Å, corresponding to structure A. The free-energy surfaces showed that structure A maintains three native hydrogen bonds, but it has larger main-chain RMSD from the native structure than misfolded and the most stable structure B. This conflict is ascribed to the N- and C-terminals of structure A, which are significantly apart due to the absence of a hydrophobic interaction between Tyr2 and Trp9 side chains. It should be noted that structure A′ is separated from structure A in Figure 9(b).

(a)

(b)

Figure 9: Free-energy surfaces of chignolin projected on the hydrogen bonds I versus III plane (a) and on the hydrogen bond II versus main-chain RMSD plane (b) obtained from all trajectories of the SAAP-MC simulation at 300 K in water. Contour lines are drawn in an interval of 1 kcal/mol.

3.3. Evaluation of the Relative Energies of Structures A–C by Ab Initio Calculation

To evaluate relative stabilities of structures A, A′, B, and C more accurately, single-point ab initio calculation was carried out. Table 2 lists the SCF and relative energies calculated in the IEFPCM solvent model. Although structure B was most populated in the SAAP-MC simulation, it is significantly unstable compared to the other structures, suggesting that structure B would not be populated in practical solutions. On the other hand, structure A′ with a native fold was found to be the most stable. According to the ab initio calculation results, it is likely that the structures having more hydrogen bonds and/or a hydrophobic interaction between the aromatic rings are more stable.

4. Discussion

4.1. Improvement of SAAPFF Parameters

In this study, we have modified the SAAPFF parameters in the following points. The main-chain, Pro, and Val parameters in water were replaced with those obtained by geometry optimization in water. The atomic charges were replaced from Mulliken to ESP charges. The energies of the main-chain unit were replaced with those obtained for the geometry with the tert-butyl cap being relaxed during the optimization. The effects of these modifications were found to be significant as the SAAP-MC simulation using the new parameters reproduced the native structure of chignolin (structure A′), while it was not obtained when the simulation was performed by using the previous parameters [34].

The clustering of the structures obtained for chignolin by SAAP-MC simulation using the improved parameters showed that structures A and B are dominant in the implicit water model as shown in Table 1. Structure A (29.2%) maintains three native hydrogen bonds, indicating the native main-chain fold, while nonnative structure B (43.2%) is stabilized by a type II β-turn. In addition, the structural clustering based on all-atom RMSD demonstrated that native structure A′, which maintains not only the three native hydrogen bonds but also a hydrophobic interaction between the aromatic side chains of Tyr2 and Trp9, is involved in 6.7% of the total structures. The free-energy analysis clearly indicated that the main-chain dihedral angles of the each amino acid residue and the distances of H-bonds I–III for structures A and A′ are consistent with the native structure (see Figures 8 and 9), although their populations were less than that of the misfolded structure B (vide infra).

In the meantime, randomly fluctuating structure was obtained for Met-enkephalin by using both the previous [34] and improved SAAPFF parameters (see Figures 3 and 4). This would be due to the two Gly residues, whose SAAP parameters have been modified only for the atomic charges in this study, being involved in the mid of the amino acid sequence.

There are a number of reports in the literature on the molecular simulation for chignolin by using conventional force fields to obtain the native structure from the unfolded state [55–64]. In these studies, various state-of-the-art techniques of molecular simulation have been employed, including a multicanonical MD method [55, 56], a replica exchange MD method [57–60], a long-term conventional MD method by using a high-speed computer system [62], and originally developed MHOP and MSES methods [63, 64]. The values of all-atom RMSD for the obtained native structure were in a range from 1.3 to 1.8 Å with respect to the experimentally determined structure, which is smaller than that of structure A obtained as a major native-like structure by the SAAP-MC simulation in this study. Moreover, the AMBER-MD simulation carried out in this study showed that the native structure of chignolin is a little more frequently obtained by using the amber99sb force field than SAAPFF. Thus, it is obvious that further improvement of the SAAPFF parameters is required in future works.

4.2. Efficiency of SAAP-MC Simulation for Conformational Sampling

Previous studies applying an explicit water molecule model [43] or a generalized-ensemble method with RISM theory [39–42] have shown that Met-enkephalin has extended conformations in water. This feature was reasonably reproduced in this study by application of SAAPFF as shown in Figures 3 and 4. Clustering of the structures obtained from the SAAP-MC simulation indicated diverse conformations of Met-enkephalin in water. On the other hand, the AMBER-MD simulation for Met-enkephalin carried out in this study produced rather uniform folded conformations with intramolecular hydrogen bonds (see Figure S2), probably due to insufficient computing time. Thus, the SAAP-MC method would be useful for sampling possible conformations for short peptide molecules in an implicit water model, although the accuracy is not satisfactory in terms of the reproducibility of the relative energies.

The efficiency for conformational sampling arises probably from a less number of variables used in the SAAPFF than that in the conventional force fields. In SAAPFF, the structure of peptides is defined only by the dihedral angles of the each amino acid unit, not by the cartesian coordinates of each atom. Acceleration of conformational sampling by reducing the number of structural parameters is a common strategy of coarse-grained force fields, such as united-atom force fields [18–26, 80, 81]. Efficiency of the SAAP-MC simulation for conformational sampling can also be explained by the fact that the solvent effects are implicitly included in the parameters of SAAPFF. This enables the SAAP-MC simulation in an implicit water model to be performed in the same computing time as that in vacuo.

4.3. Prediction of Native Structures

Although the accuracy of SAAPFF is not yet satisfied, the efficiency for conformational sampling would be advantageous for prediction of the stable structures of peptides in water. Therefore, we subsequently explored the application of the SAAP-MC method to predict the native structure of chignolin.

Indeed, it was found by single-point ab initio calculation in water at the HF/IEFPCM/6-31+G(d,p) level that structure A′ is the most stable one in water among the structures shown in Figure 5 (see Table 2). The procedure should involve a large error in the energy because the structure is not optimized. However, the tendency of the relative energies would be reliable to some extent as previously demonstrated for HCO-Ala₂-NH₂  in vacuo [33]. Thus, a protocol of SAAP-MC simulation followed by structure clustering and examination of the obtained structures by ab initio calculation would be a useful strategy toward the prediction of the native structure for short peptides. The relative energies obtained by ab initio calculation also suggested the possibility that structure B can be produced as an artifact of the force field because the relative energy was quite large (>20 kcal/mol).

As for prediction of the native structure, another simpler method would be possible based on the hardness (i.e., the number of interamino acid interactions or the depth of the potential hole on the free-energy surfaces) of the structures. As seen in Figure 5, structure B of chignolin maintains only one hydrogen bond, while structures A and A′ hold three hydrogen bonds. In terms of the number of hydrogen bonds, it can be assumed that structures A and A′ are more native-like than structure B, although the population of structure B was larger than those of structures A and A′. Indeed, the areas around the free-energy minimums corresponding to structures A and A′ in Figures 9(a) and 9(b) are slightly narrower than those corresponding to structures B and C. Such robust structures on the conformational surface would be predicted to be more native-like than the structures with more conformational flexibility. Thus, we suggest that a protocol of SAAP-MC simulation followed by structural clustering and examination of the obtained structures by the number of hydrogen bonds (or the hardness) may be another simple strategy toward structure prediction for short peptide molecules.

5. Conclusions

The parameters of SAAPFF, which was previously developed to analyze the structures and folding of polypeptides, have been improved in several points in this study. The main-chain, Pro, and Val parameters in water were replaced with those obtained by geometry optimization in the IEFPCM solvent model. The atomic charges were replaced from Mulliken to ESP charges. (3) The energies of the main-chain unit were replaced with those obtained for the geometry with the tert-butyl cap being relaxed during the optimization. To investigate the accuracy of the improved SAAPFF parameters, the SAAP-MC simulation was carried out for Met-enkephalin and chignolin as model peptides. Analysis of the obtained structures revealed that the SAAP-MC simulation reasonably reproduced randomly fluctuating conformation for Met-enkephalin, while three folded structures, one of which corresponds to a native-like structure with three native hydrogen bonds, were obtained for chignolin. The latter result supported that the accuracy of the SAAPFF parameters has been remarkably improved from the previous ones, although the native structure of chignolin could not be assigned to the most stable structure in the conformational space defined by SAAPFF.

In the meantime, efficiency of SAAP-MC simulation for conformational sampling was demonstrated for Met-enkephalin as the SAAP-MC simulation afforded diverse structures. The feature, combined with the structural clustering analysis, was subsequently applied to the structure prediction of chignolin. Among the representative structures obtained by the clustering, structure A′ with a native fold was assigned to the most stable structure according to ab initio calculation. Thus, a protocol of SAAP-MC simulation followed by structural clustering and examination of the obtained structures by ab initio calculation would be a useful strategy toward prediction of the native structure for short peptides. Alternatively, the structure with a maximum number of hydrogen bonds or with the hardest conformation on the potential surface may be considered to be the most native-like among the clustered structures. Applications of this method to structure prediction of other peptide molecules as well as further improvement of the SAAPFF parameters are being conducted by our group.

Acknowledgment

This work was supported by Grant-in-Aid for Scientific Research on Innovative Areas (no. 2120005) from the Ministry of Education, Culture, Sports, Science and Technology. The SAAP force field parameters are available for download at http://saap.sc.u-tokai.ac.jp/.

References

V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, and C. Simmerling, “Comparison of multiple amber force fields and development of improved protein backbone parameters,” Proteins: Structure, Function and Genetics, vol. 65, no. 3, pp. 712–725, 2006.

J. W. Ponder and D. A. Case, “Force fields for protein simulations,” Advances in Protein Chemistry, vol. 66, pp. 27–85, 2003.

F. A. Momany, R. F. McGuire, A. W. Burgess, and H. A. Scheraga, “Energy parameters in polypepltides. VII. Geometric parameters, partial atomic charges, nonbonded interactions, hydrogen bond interactions, and intrinsic torsional potentials for the naturally occurring amino acids,” Journal of Physical Chemistry, vol. 79, no. 22, pp. 2361–2381, 1975.

G. Némethy, K. D. Gibson, K. A. Palmer et al., “Energy parameters in polypeptides. 10. Improved geometrical parameters and nonbonded interactions for use in the ECEPP/3 algorithm, with application to proline-containing peptides,” Journal of Physical Chemistry, vol. 96, no. 15, pp. 6472–6484, 1992.

S. J. Weiner, P. A. Kollman, D. A. Case et al., “A new force field for molecular mechanical simulation of nucleic acids and proteins,” Journal of the American Chemical Society, vol. 106, no. 3, pp. 765–784, 1984.

S. J. Weiner, P. A. Kollman, D. T. Nguyen, and D. A. Case, “An all atom force field for simulations of proteins and nucleic acids,” Journal of Computational Chemistry, vol. 7, no. 2, pp. 230–252, 1986.

D. A. Pearlman, D. A. Case, J. W. Caldwell et al., “AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules,” Computer Physics Communications, vol. 91, no. 1–3, pp. 1–41, 1995.

D. A. Case, T. E. Cheatham III, T. Darden et al., “The Amber biomolecular simulation programs,” Journal of Computational Chemistry, vol. 26, no. 16, pp. 1668–1688, 2005.

B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus, “CHARMM: a program for macromolecular energy, minimization, and dynamics calculations,” Journal of Computational Chemistry, vol. 4, no. 2, pp. 187–217, 1983.

W. L. Jorgensen and J. Tirado-Rives, “The OPLS potential functions for proteins. Energy minimizations for crystals of cyclic peptides and crambin,” Journal of the American Chemical Society, vol. 110, no. 6, pp. 1657–1666, 1988.

W. L. Jorgensen and J. Tirado-Rives, “Potential energy functions for atomic-level simulations of water and organic and biomolecular systems,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 19, pp. 6665–6670, 2005.

W. R. P. Scott, P. H. Hünenberger, I. G. Tironi et al., “The GROMOS biomolecular simulation program package,” Journal of Physical Chemistry A, vol. 103, no. 19, pp. 3596–3607, 1999.

U. H. E. Hansmann, Y. Okamoto, and F. Eisenmenger, “Molecular dynamics, langevin and hybrid Monte Carlo simulations in a multicanonical ensemble,” Chemical Physics Letters, vol. 259, no. 3-4, pp. 321–330, 1996.

N. Nakajima, H. Nakamura, and A. Kidera, “Multicanonical ensemble generated by molecular dynamics simulation for enhanced conformational sampling of peptides,” Journal of Physical Chemistry B, vol. 101, no. 5, pp. 817–824, 1997.

K. Hukushima and K. Nemoto, “Exchange monte carlo method and application to spin glass simulations,” Journal of the Physical Society of Japan, vol. 65, no. 6, pp. 1604–1608, 1996.

Y. Sugita and Y. Okamoto, “Replica-exchange molecular dynamics method for protein folding,” Chemical Physics Letters, vol. 314, no. 1-2, pp. 141–151, 1999.

Y. Okamoto, “Generalized-ensemble algorithms: enhanced sampling techniques for Monte Carlo and molecular dynamics simulations,” Journal of Molecular Graphics and Modelling, vol. 22, no. 5, pp. 425–439, 2004.

L. Yang, C. Tan, M.-J. Hsieh et al., “New-generation amber united-atom force field,” The Journal of Physical Chemistry B, vol. 110, no. 26, pp. 13166–13176, 2006.

A. Liwo, S. Ołdziej, M. R. Pincus, R. J. Wawak, S. Rackovsky, and H. A. Scheraga, “A united-residue force field for off-lattice protein-structure simulations. I. Functional forms and parameters of long-range side-chain interaction potentials from protein crystal data,” Journal of Computational Chemistry, vol. 18, no. 7, pp. 849–873, 1997.

A. Liwo, M. R. Pincus, R. J. Wawak, S. Rackovsky, S. Ołdziej, and H. A. Scheraga, “A united-residue force field for off-lattice protein-structure simulations. II. Parameterization of short-range interactions and determination of weights of energy terms by Z-score optimization,” Journal of Computational Chemistry, vol. 18, no. 7, pp. 874–887, 1997.

A. Liwo, R. Kaźmierkiewicz, C. Czaplewski et al., “United-residue force field for off-lattice protein-structure simulations: III. Origin of backbone hydrogen-bonding cooperativity in united-residue potentials,” Journal of Computational Chemistry, vol. 19, no. 3, pp. 259–276, 1998.

A. Liwo, C. Czaplewski, J. Pillardy, and H. A. Scheraga, “Cumulant-based expressions for the multibody terms for the correlation between local and electrostatic interactions in the united-residue force field,” Journal of Chemical Physics, vol. 115, no. 5, pp. 2323–2347, 2001.

A. V. Rojas, A. Liwo, and H. A. Scheraga, “Molecular dynamics with the united-residue force field: Ab initio folding simulations of multichain proteins,” Journal of Physical Chemistry B, vol. 111, no. 1, pp. 293–309, 2007.

S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. De Vries, “The MARTINI force field: coarse grained model for biomolecular simulations,” Journal of Physical Chemistry B, vol. 111, no. 27, pp. 7812–7824, 2007.

L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman, and S. Marrink, “The MARTINI coarse-grained force field: extension to proteins,” Journal of Chemical Theory and Computation, vol. 4, no. 5, pp. 819–834, 2008.

A. Rojas, A. Liwo, D. Browne, and H. A. Scheraga, “Mechanism of fiber assembly: treatment of Aβ peptide aggregation with a coarse-grained united-residue force field,” Journal of Molecular Biology, vol. 404, no. 3, pp. 537–552, 2010.

M. Iwaoka, D. Yosida, and N. Kimura, “Importance of the single amino acid potential in water for secondary and tertiary structures of proteins,” Journal of Physical Chemistry B, vol. 110, no. 29, pp. 14475–14482, 2006.

M. Iwaoka, M. Okada, and S. Tomoda, “Solvent effects on the φ-ψ potential surfaces of glycine and alanine dipeptides studied by PCM and I-PCM methods,” Journal of Molecular Structure: THEOCHEM, vol. 586, pp. 111–124, 2002.

L. J. Smith, K. M. Fiebig, H. Schwalbe, and C. M. Dobson, “The concept of a random coil : residual structure in peptides and denatured proteins,” Folding and Design, vol. 1, no. 5, pp. R95–R106, 1996.

K. M. Fiebig, H. Schwalbe, M. Buck, L. J. Smith, and C. M. Dobson, “Toward a description of the conformations of denatured states of proteins. Comparison of a random coil model with NMR measurements,” Journal of Physical Chemistry, vol. 100, no. 7, pp. 2661–2666, 1996.

Y. Sakae and Y. Okamoto, “Optimization of protein force-field parameters with the Protein Data Bank,” Chemical Physics Letters, vol. 382, no. 5-6, pp. 626–636, 2003.

N. Kamiya, Y. S. Watanabe, S. Ono, and J. Higo, “AMBER-based hybrid force field for conformational sampling of polypeptides,” Chemical Physics Letters, vol. 401, no. 1–3, pp. 312–317, 2005.

M. Iwaoka and S. Tomoda, “The SAAP force field. A simple approach to a new all-atom protein force field by using Single Amino Acid Potential (SAAP) functions in various solvents,” Journal of Computational Chemistry, vol. 24, no. 10, pp. 1192–1200, 2003.

M. Iwaoka, N. Kimura, D. Yosida, and T. Minezaki, “The SAAP force field: development of the single amino acid potentials for 20 proteinogenic amino acids and monte carlo molecular simulation for short peptides,” Journal of Computational Chemistry, vol. 30, no. 13, pp. 2039–2055, 2009.

G. D. Hawkins, C. J. Cramer, and D. G. Truhlar, “Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium,” Journal of Physical Chemistry, vol. 100, no. 51, pp. 19824–19839, 1996.

G. D. Hawkins, C. J. Cramer, and D. G. Truhlar, “Pairwise solute descreening of solute charges from a dielectric medium,” Chemical Physics Letters, vol. 246, no. 1-2, pp. 122–129, 1995.

V. Tsui and D. A. Case, “Theory and applications of the generalized Born solvation model in macromolecular simulations,” Biopolymers, vol. 56, no. 4, pp. 275–291, 2001.

W. H. Graham, E. S. Carter II., and R. P. Hicks, “Conformational analysis of Met-enkephalin in both aqueous solution and in the presence of sodium dodecyl sulfate micelles using multidimensional NMR and molecular modeling,” Biopolymers, vol. 32, no. 12, pp. 1755–1764, 1992.

A. Mitsutake, M. Kinoshita, Y. Okamoto, and F. Hirata, “Combination of the replica-exchange Monte Carlo method and the reference interaction site model theory for simulating a peptide molecule in aqueous solution,” Journal of Physical Chemistry B, vol. 108, no. 49, pp. 19002–19012, 2004.

A. Mitsutake, M. Kinoshita, Y. Okamoto, and F. Hirata, “Multicanonical algorithm combined with the RISM theory for simulating peptides in aqueous solution,” Chemical Physics Letters, vol. 329, no. 3-4, pp. 295–303, 2000.

M. Kinoshita, Y. Okamoto, and F. Hirata, “First-principle determination of peptide conformations in solvents: combination of Monte Carlo simulated annealing and RISM theory,” Journal of the American Chemical Society, vol. 120, no. 8, pp. 1855–1863, 1998.

M. Kinoshita, Y. Okamoto, and F. Hirata, “Solvation structure and stability of peptides in aqueous solutions analyzed by the reference interaction site model theory,” Journal of Chemical Physics, vol. 107, no. 5, pp. 1586–1599, 1997.

L. Ramya and N. Gautham, “Effects of hydration on the conformational energy landscape of the pentapeptide Met-enkephalin,” Journal of Chemical Theory and Computation, vol. 5, no. 8, pp. 2180–2190, 2009.

K. Y. Sanbonmatsu and A. E. García, “Structure of Met-enkephalin in explicit aqueous solution using replica exchange molecular dynamics,” Proteins: Structure, Function and Genetics, vol. 46, no. 2, pp. 225–234, 2002.

D. Van Der Spoel and H. J. C. Berendsen, “Molecular dynamics simulations of leu-enkephalin in water and DMSO,” Biophysical Journal, vol. 72, no. 5, pp. 2032–2041, 1997.

L. Zhan, J. Z. Y. Chen, and W. Liu, “Comparison of predicted native structures of met-enkephalin based on various accessible-surface-area solvent models,” Journal of Computational Chemistry, vol. 30, no. 7, pp. 1051–1058, 2009.

L. Su and R. I. Cukier, “Hamiltonian and distance replica exchange method studies of met-enkephalin,” Journal of Physical Chemistry B, vol. 111, no. 42, pp. 12310–12321, 2007.

D. A. Evans and D. J. Wales, “The free energy landscape and dynamics of met-enkephalin,” Journal of Chemical Physics, vol. 119, no. 18, pp. 9947–9955, 2003.

M. M. Palian, Boguslavsky, D. F. O'Brien, and R. Polt, “Glycopeptide-membrane interactions: glycosyl enkephalin analogues adopt turn conformations by NMR and CD in amphipathic media,” Journal of the American Chemical Society, vol. 125, no. 19, pp. 5823–5831, 2003.

B. A. Berg and H. Hsu, “Metropolis simulations of Met-Enkephalin with solvent-accessible area parametrizations,” Physical Review E, vol. 69, no. 2, Article ID 026703, 9 pages, 2004.

K. Vengadesan and N. Gautham, “Energy landscape of Met-enkephalin and Leu-enkephalin drawn using mutually orthogonal Latin squares sampling,” Journal of Physical Chemistry B, vol. 108, no. 30, pp. 11196–11205, 2004.

S. G. Itoh and Y. Okamoto, “Effective sampling in the configurational space of a small peptide by the multicanonical-multioverlap algorithm,” Physical Review E, vol. 76, no. 2, Article ID 026705, 2007.

A. Mitsutake, U. H. E. Hansmann, and Y. Okamoto, “Temperature dependence of distributions of conformations of a small peptide,” Journal of Molecular Graphics and Modelling, vol. 16, no. 4–6, pp. 226–238, 1998.

S. Honda, K. Yamasaki, Y. Sawada, and H. Morii, “10 Residue folded peptide designed by segment statistics,” Structure, vol. 12, no. 8, pp. 1507–1518, 2004.

D. Satoh, K. Shimizu, S. Nakamura, and T. Terada, “Folding free-energy landscape of a 10-residue mini-protein, chignolin,” FEBS Letters, vol. 580, no. 14, pp. 3422–3426, 2006.

T. Terada, D. Satoh, T. Mikawa, Y. Ito, and K. Shimizu, “Understanding the roles of amino acid residues in tertiary structure formation of chignolin by using molecular dynamics simulation,” Proteins: Structure, Function and Genetics, vol. 73, no. 3, pp. 621–631, 2008.

D. Van Der Spoel and M. M. Seibert, “Protein folding kinetics and thermodynamics from atomistic simulations,” Physical Review Letters, vol. 96, no. 23, Article ID 238102, 2006.

M. M. Seibert, A. Patriksson, B. Hess, and D. Van Der Spoel, “Reproducible polypeptide folding and structure prediction using molecular dynamics simulations,” Journal of Molecular Biology, vol. 354, no. 1, pp. 173–183, 2005.

W. Xu, T. Lai, Y. Yang, and Y. Mu, “Reversible folding simulation by hybrid Hamiltonian replica exchange,” Journal of Chemical Physics, vol. 128, no. 17, Article ID 175105, 2008.

S. Kannan and M. Zacharias, “Enhanced sampling of peptide and protein conformations using replica exchange simulations with a peptide backbone biasing-potential,” Proteins: Structure, Function and Genetics, vol. 66, no. 3, pp. 697–706, 2007.

X. Dou and J. Wang, “Folding free energy landscape of the decapeptide chignolin,” Modern Physics Letters B, vol. 22, no. 31, pp. 3087–3098, 2008.

A. Suenaga, T. Narumi, N. Futatsugi et al., “Folding dynamics of 10-residue beta-hairpin peptide chignolin,” Chemistry, vol. 2, no. 5, pp. 591–598, 2007.

S. Roy, S. Goedecker, M. J. Field, and E. Penev, “A minima hopping study of all-atom protein folding and structure prediction,” Journal of Physical Chemistry B, vol. 113, no. 20, pp. 7315–7321, 2009.

K. Moritsugu, T. Terada, and A. Kidera, “Scalable free energy calculation of proteins via multiscale essential sampling,” Journal of Chemical Physics, vol. 133, no. 22, Article ID 224105, 2010.

B. Mennucci, E. Cancès, and J. Tomasi, “Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a unified integral equation method: theoretical bases, computational implementation, and numerical applications,” Journal of Physical Chemistry B, vol. 101, no. 49, pp. 10506–10517, 1997.

E. Cancès, B. Mennucci, and J. Tomasi, “A new integral equation formalism for the polarizable continuum model: theoretical background and applications to Isotropic and anisotropic dielectrics,” Journal of Chemical Physics, vol. 107, no. 8, pp. 3032–3041, 1997.

G. Scalmani and M. J. Frisch, “Continuous surface charge polarizable continuum models of solvation. I. General formalism,” Journal of Chemical Physics, vol. 132, no. 11, Article ID 114110, 2010.

M. J. Frisch, G. W. Trucks, H. B. Schlegel, et al., Gaussian 03, Revision C.02, Gaussian, Wallingford, UK, 2004.

Y. K. Kang, “Ab initio and DFT conformational study of proline dipeptide,” Journal of Molecular Structure: THEOCHEM, vol. 675, no. 1–3, pp. 37–45, 2004.

M. Makowski, J. Makowska, and L. Chmurzyński, “Ab initio study of the energetics of protonation, deprotonation and homocomplexed cations and anions formation in systems modeling side chains of biomolecules,” Journal of Molecular Structure: THEOCHEM, vol. 674, no. 1–3, pp. 61–67, 2004.

B. H. Besler, K. M. Merz Jr., and P. A. Kollman, “Atomic charges derived from semiempirical methods,” Journal of Computational Chemistry, vol. 11, no. 4, pp. 431–439, 1990.

V. V. Kostjukov, N. M. Khomytova, A. A. H. Santiago, R. L. Ibarra, D. B. Davies, and M. P. Evstigneev, “Calculation of the electrostatic charges and energies for intercalation of aromatic drug molecules with DNA,” International Journal of Quantum Chemistry, vol. 111, no. 3, pp. 711–721, 2011.

D. A. Case, T. A. Darden, T. E. Cheatham III, et al., AMBER 10, University of California, San Francisco, Calif, USA, 2008.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953.

K. Binder, “Applications of Monte Carlo methods to statistical physics,” Reports on Progress in Physics, vol. 60, no. 5, pp. 487–559, 1997.

M. Matsumoto and T. Nishimura, “Mersenne Twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Transactions on Modeling and Computer Simulation, vol. 8, no. 1, pp. 3–30, 1998.

I. Marcotte, F. Separovic, M. Auger, and S. M. Gagné, “A multidimensional ¹H NMR investigation of the conformation of methionine-enkephalin in fast-tumbling bicelles,” Biophysical Journal, vol. 86, no. 3, pp. 1587–1600, 2004.

A. W. Götz, M. J. Williamson, D. Xu, D. Poole, S. La Grand, and R. C. Walker, “Routine microsecond molecular dynamics simulations with AMBER on GPUs. 1. Generalized Born,” Journal of Chemical Theory and Computation, vol. 8, no. 5, pp. 1542–1555, 2012.

J. Shao, S. W. Tanner, N. Thompson, and T. E. Cheatham III., “Clustering molecular dynamics trajectories: I. Characterizing the performance of different clustering algorithms,” Journal of Chemical Theory and Computation, vol. 3, no. 6, pp. 2312–2334, 2007.

H. A. Scheraga, M. Khalili, and A. Liwo, “Protein-folding dynamics: overview of molecular simulation techniques,” Annual Review of Physical Chemistry, vol. 58, pp. 57–83, 2007.

Y. Sakae and Y. Okamoto, “Optimisation of OPLS-UA force-field parameters for protein systems using protein data bank,” Molecular Simulation, vol. 36, no. 14, pp. 1148–1156, 2010.