Metadynamics As a Tool for Mapping the Conformational and Free

5632

J. Phys. Chem. B 2010, 114, 5632–5642

Metadynamics As a Tool for Mapping the Conformational and Free-Energy Space of Peptides s The Alanine Dipeptide Case Study Jirˇ´ı Vymeˇtal†,‡ and Jirˇ´ı Vondra´sˇek*,‡ Faculty of Science, Charles UniVersity in Prague, AlbertoV 2038/6, 128 00 Prague 2, Czech Republic, and Institute of Organic Chemistry and Biochemistry AS CR, V.V.i., FlemingoVo nám. 2, 166 10 Prague 6, Czech Republic ReceiVed: February 1, 2010; ReVised Manuscript ReceiVed: March 12, 2010

There is a need for a fast, accurate, and reliable method of sampling the conformational space of peptides and proteins in order to obtain a balanced free-energy profile which can lead to our understanding of protein structure. We have utilized metadynamics for the conformational study of the solvated alanine dipeptide molecule, and our results show that the method has proven to be competent as a fast, robust, and reliable method for the conformation free-energy calculations of peptides in an explicit solvent, surpassing traditional methods such as umbrella sampling. We have also addressed the issue of the influence of different water models on the resulting free-energy profile in order to consistently decompose the setting of our simulation. All of the explicit water models for the simulation of biomolecules TIP3P, TIP4P, TIP4P/Ew, TIP5P, and SPCE have exhibited similar effects on the conformational preferences of alanine dipeptide with no significant differences. On the other hand, by comparing the potential energy surface in the gas phase and the freeenergy surface in a water environment, we have shown that the interaction with water molecules is one of the most important structure-driving elements, with a great influence on the free-energy surface (FES) of the solvated peptide and the conformational preferences of the peptide backbone. All of the tested force fields (ff03, ff99SB, opls-aa, and charmm27) appreciably differ in the population of the individual conformers and the barriers between them. Significant divergence was found on both the potential energy surface (PES) as well as free-energy surface (FES) calculated by charmm27. We have therefore concluded that the differences originate dominantly from the parametrization of the peptide backbone in the given force field rather than from a noncovalent interaction with water molecules. Introduction The application of the physical principles governing biomolecules in order to adopt their unique 3D structure implies a need for reasonably accurate methods able to map conformational hyperspace at the Gibbs free-energy level. In the case of proteins, the question is rather complex and includes the issue of at which level the local conformational preferences of amino acids in short sequential context are propagated into the final protein structure. The role of amino acid side chains is not quite clear in the structural context because the backbone preferences define the number of possibilities available for the side chains to interact. This could crucially constrain the spatial possibility of finding the near-energy optimum geometry for each interacting side chain. Fortunately, we can study the role of the backbone at a lower structural level, for example, at the level of small peptides. These studies could provide, in a great detail, a relation between the main-chain and side-chain influences on the resulting structure. We need a full structure-energy picture based on thermodynamically correct roots that explicitly involves the solvent in overcoming the limitation usually brought by potential energy surface mapping. It is now generally accepted that the non-native or denaturated state of a protein is not a random ensemble of all of the possible conformations of the protein molecule but rather a limited set of conformations retaining some of the local structural preferences dictated by * To whom correspondence should be addressed. E-mail: jiri.vondrasek@ uochb.cas.cz. Tel: +420220410324. Fax: +420220410321. † Charles University in Prague. ‡ Institute of Organic Chemistry and Biochemistry AS CR.

the sequence and by the interaction within the local environment.1,2 To succeed in the computational modeling of this task, we first need a method capable of properly exploring the free-energy surface of peptides, which is closely connected with the accuracy of the computational method evaluating the potential energy of the system in every step. In the field of biomolecules and their simulations, we need a reliable empirical potential providing a good estimate of the thermodynamic properties of the molecule. As the relatively accurate methods of computational chemistry are applied to large molecules, there is still a question whether they are generally parametrized well or at which level the canonical ab initio methods can be adopted. The relative tendencies of the peptide molecule to adopt a certain 3D arrangement can be quantified by conformational free energy. However, the calculation of this quality as well as prediction of the most populated conformers for flexible molecules like peptides is notoriously difficult. First of all, the number of all of the possible conformational states grows exponentially with the count of rotatable bonds and the length of the polymer. A direct sampling of conformational space by means of a classical simulation of molecular dynamics suffers from an ineffective crossing of the barriers and other rare events, as a result of which the important regions of conformational space may remain in finite computational time. An effective way to investigate the fundamental aspects of conformational space efficiently relies on the course-graining approach. Fortunately, satisfactory insight into multidimensional conformational space of a molecule (with 3N - 6 degrees of freedom, where N is the number of atoms) can often be acquired

10.1021/jp100950w  2010 American Chemical Society Published on Web 04/02/2010

Metadynamics of Alanine Dipeptide by the introduction of several suitable collective variables (order parameters, reaction coordinates). The wisely selected collective variables (CV) allow a distinction to be drawn between the individual conformers or at least the families of characteristic conformers. The free-energy landscape F(s) as a function of general CV(s) is related to the corresponding probability density F(s) up to the arbitrary constant C

1 F(s) ) - ln〈F(s)〉 + C β where the ensemble-averaged probability density is given by

〈F(s)〉 )

∫ dX δ(s - s'(X))e-βU(X) ∫ dX e-βU(X)

Here, δ is a Dirac’s delta function utilized for the effective integration of all of the variables except for CVs, which is also a function of the Cartesian coordinates X as well as the potential energy U(X). Although an estimation of F(s) may be obtained from an ordinary unbiased MD simulation, the biased protocols offer incomparably better efficiency because they overcome the sampling problem. Particularly, the methods utilizing the adaptively constructed bias potential have become intensively developed with promising outcomes. One of them, the recently elaborated metadynamics, aims at a reconstruction of the free-energy surfaces (FES) and acceleration of rare events. Since being proposed by Laio and Parrinello,3 metadynamics have been successfully applied to various computational fields, including the folding of small proteins and peptides,4-9 the study of protein-protein interactions,10 or the study of conformational transitions of biomolecules.11 The theoretical background and application of metadynamics and its further variants are comprehensively reviewed in ref 12. In our present study, we have compared the computational efficiency and accuracy of metadynamics with the wellestablished umbrella sampling method enhanced by the weighted histogram analysis method (WHAM).13 From our point of view, metadynamics surpasses the umbrella sampling (US) methodology in several ways. First of all, the sampling of conformational space can be performed within a single simulation, in contrast to US, which confines the exploration of the FES to a priori determined areas delimited by the harmonic bias potential. Every part of the FES is thus in US methodology sampled in an independent simulation of a properly equilibrated system. Afterward, the collected histograms need to be reweighed and merged by the WHAM procedure in order to obtain the final FES. The number of simulations (windows) required for the reconstruction of the FES may be practically fairly high, in a range from tens to hundreds in the case of a multidimensional FES. On the other hand, metadynamics constructs the bias potential continuously with a destabilizing effect on the current conformation and facilitates transitions between regions of conformational space. Essentially, the dynamics of the investigated system evolve in a nonequilibrate, non-Markovian way and dramatically extend the sampling of the CV in comparison to the equilibrium Boltzmann distribution. The bias potential is usually constructed as a sum of the Gaussian functions, which

J. Phys. Chem. B, Vol. 114, No. 16, 2010 5633 are continuously appended in the course of the simulation time. Finally, the estimation of the FES as an equilibrium property is based directly on the bias potential without any reweighing in the postprocessing as

F(s) ) -Vbias(s) The molecular system in question, alanine dipeptide, serves as an excellent benchmark for a wide range of computational methods because of its moderate size and complexity.14-21 There are computational studies at various levels of methods from the empirical potentials and classical MD simulations up to ab initio dynamics in an implicit solvent. Last but not least, alanine dipeptide has also been studied experimentally. Assuming rigidity of peptide bonds and neglecting the rotation of the methyl groups, all of the conformers may be fully characterized by two independent dihedrals. These two parameters represent an obvious choice of the collective coordinates for the construction of the two-dimensional FES in both methods or the potential energy surfaces (PES). Methods Model System. The alanine dipeptide was built by the tLEaP program from amber tools22 as an alanine residue (Ala) terminated by acetyl (Ace) and N-methyl (Nme) capping groups. The molecule was solvated in a cubic box with an edge of 30 Å by approximately 880 water molecules. Force Fields. Several common molecular mechanics force fields for biomolecules were utilized in this study, namely, the (i) amber ff03,23 (ii) Simmerling’s modification of ff99 (ff99SB),24 (iii) charmm 27,25 and (iv) opls-aa.26 The Amber force fields were adopted from Amberport27 (Feb 2009 release), a package of parameters for several amber force fields adapted for Gromacs. Analogously, the port of the charmm27 force field (without CMAP corrections) for Gromacs was obtained from Par Bjelkmar from the Lindahl lab (http://www.dbb.su.se/User: Bjelkmar/Ffcharmm, Nov 2008 release). The parameters for the opls-aa force field (release 2001) were obtained from a standard Gromacs distribution (ver 4.0.4). Water Models. We modeled the solvent using five different widespread explicit water models TIP3P, TIP4P,28 TIP4P/Ew, TIP5P,29 and SPCE.30 The parameters for the TIP4P/Ew water model were adapted from the original paper,31 whereas the others were taken from a standard Gromacs distribution (ver. 4.0.4). Simulation Details. All of the simulations were performed by an “in-house-modified” Gromacs package for molecular dynamics (ver. 4.0.4)32 enhanced by our implementation of metadynamics. The steepest-descent gradient method of geometry optimization was applied prior to the simulations in order to correct any incidental sterical clash and atomic overlaps in the system, followed by two-step equilibration procedures. First, the NVT ensemble was simulated, and the final temperature of 300 K was reached utilizing a Berendsen weak coupling thermostat starting from the initial randomly generated velocities. The convenient pressure was adjusted by a weak coupling barostat in the second phase (NpT ensemble) of equilibration. The length of the both preparative simulations was 500 ps, with a total of 1 ns for the entire equilibration of each system. The production phase of the simulations used for the freeenergy calculations employing both methods, umbrella sampling and metadynamics, was done with the same setup. The simulations were performed at a constant temperature

5634

J. Phys. Chem. B, Vol. 114, No. 16, 2010

of 300 K controlled by a Nose´-Hoover thermostat with the time constant τ set to 0.1 ps and at a constant pressure of 1 bar assured by a Parrinello-Rahman barostat with a time constant of τ )1 ps and an isothermal compressibility of 4.5 × 10-5 bar-1. The periodic boundary conditions for the cubic box were utilized to avoid discontinuities in the simulated systems. The time step of 2 fs was used for the integration of the equations of motion. All of the bond lengths were constrained by means of a fourth-order LINCS algorithm33 with two iterations. The calculation of the nonbonded interactions was optimized by a grid-based neighbor list algorithm updated every 10 steps. The van der Waals forces were controlled by a switch function at a distance of 10 Å, diminishing to zero at 14 Å, as implemented in Gromacs. The dispersion corrections to the total energy and pressure were applied in the simulated systems. The electrostatics were treated by a Particle Mesh Ewald (PME) method of the fourth order with a tolerance of 1 × 10-5, a real space cutoff of 10 Å, and a spacing of the Fourier grid of 1 Å. Umbrella Sampling and WHAM. The umbrella sampling methodology34,35 enhanced by the weighted histogram analysis method (WHAM)13 was applied to a series of simulations in order to get a free-energy profile in the φ and ψ coordinates. The whole dihedral two-dimensional conformational space was regularly divided into 324 windows by 20° in both coordinates. In each window, the system was augmented by a harmonic potential (we utilized the dihedral restraint feature of Gromacs32) with a stiffness constant of 100 kJ mol-1 rad-2, which restrains the values of the φ and ψ dihedrals around the center of the respective window. Molecular dynamics in each window was propagated for 2.4 ns, with the first 0.4 ns being regarded as equilibration with an additional umbrella potential and discarded from further analysis. The dihedral values were recorded every 0.1 ps, resulting in 20 000 samples in the histogram for each window. Consequently, the WHAM procedure was performed for the analysis and concatenation of all of the histograms into the free-energy profile by means of the wham-2d (revision 192) program by Alan Grossfield (available at http://membrane. urmc.rochester.edu/Software/WHAM/WHAM.html), which follows algorithms described in ref 36. The histograms were counted with a bin size of 2 × 2°, and the WHAM equations were settled to a convergence criterion of 1 × 10-6 kcal mol-1. The total simulation time needed for the construction of one profile in our study was 777.6 ns. In order to estimate the confidence intervals of the acquired profiles, all of the simulations and the whole procedure were independently repeated 10 times for different initial conditions. Metadynamics Protocol. The type of metadynamics that we used was the direct version in the two collective coordinates represented by the φ and ψ dihedrals of alanine dipeptide. Our protocol consisted of two parts: The bias potential was constructed in the course of the first part until a reasonable convergence of the free-energy profile was reached. The data for statistical analysis were collected in the second part of the simulation. The total length of each simulation was 60 ns. The first 40 ns were regarded as the first “flooding” phase, whereas the remaining 20 ns served for the estimation of the errors and the levels of confidence in the resulting free-energy profiles. The bias potential was regularly updated at 4 ps intervals throughout the simulation; and 15 000 Gaussians were successively deposited for the reconstruction of each freeenergy profile. The parameters of the Gaussians were chosen

Vymeˇtal and Vondra´sˇek on the basis of our previous experience and the pilot calculations as a compromise between a fine resolution and the computational cost. The width of a Gaussian (σ) determining the spatial resolution was set to 0.3 radians (approximately 17°). The height of the Gaussians was gradually decreased every 20 ns, starting at 0.2 kJ mol-1 for the first period and at 0.1 and 0.05 kJ mol-1 for the consecutive periods, in order to reduce the disruption in the profiles and rather promote their better convergence. Reconstruction of the Free-Energy Profiles in Metadynamics. The bias potential was saved in an external file during the entire simulation as a list of parameters for the individual Gaussians. The reconstruction and statistical analysis of the final free-energy profiles were performed by special programs and scripts developed by us. All of the Gaussians were rendered and continuously summed on a grid with a 500 × 500 point resolution. During the reconstruction of the profile, the ∆F values among individual conformers were calculated and served as an indicator of convergence. This measure of convergence was performed after each 50th Gaussian, and the corresponding profiles (snapshots) were stored. The estimation of the levels of confidence was performed by means of the block average method applied to the last 20 ns of the simulation. The average free-energy profile and relative conformational free-energy differences were reported as a result of a single simulation of metadynamics. The onedimensional projection of the 2D free-energy profile was obtained by a numerical integration of the relation F(x1) ) -(1/β) ln ∫e-βF(x1,x2) dx2, where x1 and x2 refer to the φ or ψ coordinate. The statistical analysis was performed in the same manner as in the previous case. Calculation of Relative Conformational Free Energies. The conformational free energies were determined based on the definition of the individual conformers and the specific allocation of the free-energy landscape. In this study, six distinctive conformational states of the amino acid backbone were distinguished, corresponding to the typical minima found on the obtained profiles and regularly covering all of the conformational space. The nomenclature was adopted from an ab initio quantum study of alanine dipeptide.21 Our definition of the conformers in terms of φ and ψ dihedrals is as follows (Chart 1). The conformational free energies were calculated via a numerical integration over the corresponding area Ω on a grid representing the obtained free-energy profiles

1 F(Ω) ) - ln β

A(φ,ψ)∈Ω e-βF(ψ,φ) dφ dφ

Calculation of Conformational Populations. The corresponding conformational populations were calculated as the derived properties from individual FESs via the numerical integration

P(Ω) )

A(φ,ψ)∈Ω e-βF(ψ,φ) dφ dφ Aφ∈ -180,180)φ∈ -180,180) e-βF(ψ,φ) dφ dφ 〈

〈

Umbrella Sampling versus Metadynamics Comparison of Free-Energy Profiles. The conformational free energy of alanine dipeptide was compared in order to assess the deviations originating from different methods of free-energy calculation or a distinct interaction potential. The comparisons were

Metadynamics of Alanine Dipeptide

J. Phys. Chem. B, Vol. 114, No. 16, 2010 5635

CHART 1

performed by means of two different approaches. First, the relative free-energy differences among the defined conformers were evaluated. The other method introduced by Alonso37 compared energy surfaces by means of a physically adapted form of the correlation coefficient 2 d1,2 ) √(σ21 + σ22)(1 - r1,2 )

Results (1)

where indexes 1 and 2 refer to the energy profiles (F and F(2)) and σ2 refers to the corresponding statistical variances of the values on a given profile N

σ21 )

∑

1 (F(1) - µ1)2 N i)1 i

N

µ1 )

∑

1 F(1) N i)1 i

where N stands for the number of calculated points on an energy surface. r1,2 is the Pearson correlation coefficient between two energy surfaces

r1,2

coordinates were sampled in regular steps of 0.72°, and the final plots were obtained on a grid of 500 × 500 points. Visualization. The two-dimensional free-energy profiles and potential energy profiles were visualized utilizing the MayaVi2 tool kit. The other plots were elaborated in the XmGrace plotting tool (http://plasma-gate.weizmann.ac.il/Grace/).

N (Fi(1) - µ1)(Fi(2) - µ2) 1 ) N i)1 σ1σ2

∑

We exploited d1,2 as an estimator of the average deviation in ∆F for an arbitrary conformational transition in two different potentials. Therefore, d1,2 has the physical meaning of a unit of energy and can be considered as a measure of similarity. Moreover, the significance of the diversions revealed by d1,2 can be related to the magnitude of the thermal noise (RT or kBT). If d1,2 values of two potentials exceed the level of thermal fluctuations, the resulting dynamics and properties of the systems driven by them may differ substantially. Calculation of Potential Energy Profiles. The profiles of the potential energy in the φ and ψ coordinates were constructed by means of a restrained minimization of the pregenerated structures. Molecular geometries with the desired combination of φ and ψ dihedrals were prepared from an extended conformer by rotation around the corresponding bonds and consequently optimized by the conjugant gradient method as implemented in Gromacs with additional harmonic restrains (with a force constant of 100 000 kJ mol-1 rad-1) placed on the φ and ψ dihedrals. Both of the

Convergence of Metadynamics. We examined the properties of the selected metadynamics simulation protocol in order to ensure its reliability and utilization. We then chose a solvated alanine dipeptide molecule with an explicit solvent model (TIP3P) as a model system for this study. The convergence plot of the conformational ∆F is shown in Figure 1. The vertical axes divide the simulation into three parts according to the different heights of the Gaussians used in this time interval. The gradual decrease of their amplitudes obviously facilitated the convergence, as was demonstrated by the almost steady course of ∆F in the last part of the simulation. Comparison of Two Different Methods for the Calculation Profiles of Free Energy: Umbrella Sampling (WHAM) and Metadynamics on the Alanine Dipeptide Molecule. We performed a set of simulations in order to test the accuracy and

Figure 1. The convergence of our protocol for metadynamics. The relative free-energy differences (∆F) among all of the considered conformers plotted as a function of the number of deposited Gaussians in the course of one representative simulation (ff03, TIP3P water). The heights of the Gaussians are 0.2, 0.1, and 0.05 kJ mol-1 in the drawn fields, respectively.

5636


Vymeˇtal and Vondra´sˇek

Figure 2. A comparison of the free-energy profiles provided by the umbrella sampling method and metadynamics. The two-dimensional plot of the conformational free energy of alanine dipeptide (ff03 with TIP3P water) constructed in the φ and ψ dihedrals as the parameters obtained by metadynamics (A) and umbrella sampling (B). The contours are drawn at each 2 kJ mol-1. The one-dimensional projections of the former plots on the individual dihedrals are also presented (C, D) for both methods.

efficiency of both of the methods commonly used for free-energy calculations, that is, metadynamics and umbrella sampling. The free-energy profiles were obtained in the φ and ψ dihedral coordinates, which describe the conformational state of the molecule quite well. The two-dimensional profiles are shown in Figure 2a and c, respectively, together with their onedimensional projection on the individual dihedral coordinates for an explicit comparison in Figure 2b and d. Both profiles were constructed as an average of 10 freestanding, mutually independent profiles, representing 10 independent simulations with the same setup. Such an approach allows the estimation of the confidence level of the obtained values and the reproducibility of the individual free-energy profiles obtained from single simulations. Both of the two-dimensional free-energy landscapes as well as their one-dimensional projections calculated by metadynamics (Figure 2a-d) and the umbrella sampling method showed a significant level of similarity. The overall shapes are the same, and the positions of the topologically important points (minima, barriers, and transition pathways) seem to be perfectly located at the same places. The minima are evident on both 2D plots. The conformational free-energy differences among these minima and their calculated populations are listed in Table 1a-c. They support our utilization of metadynamics and serve as a proof of validity for our implementation of metadynamics. The three most favorable regions, β, RR, and C5 (also known as polyproline helix 2-like, R-helical and extended, respectively), were found to be almost isoenergetic (umbrella sampling versus metadynamics), and moreover, both methodologies provided their identical order. A similar conclusion can be made also for

the less favorable minima (R′, RL, and RD). However, the global minimum on the free-energy surface in the ff03 force field for alanine dipeptide is not possible to determine unambiguously from a single simulation. Both metadynamics and the umbrella sampling protocols yield the β and RR regions as the most stable conformers. The standard deviations for the conformational ∆F reported in Table 1 show metadynamics as a method of more uniform and reproducible results. To support such a conclusion, we tested it by means of the physically rationalized correlation coefficient d1,2. The mutual similarity of the individual freeenergy profiles measured by dA,B is quoted in Table 2a and b. The value can be interpreted as the averaged unsigned differences of the conformational Gibbs free energy ∆F provided that all of the possible conformational transitions take places on two arbitrary energy landscapes. The similarity of the independent metadynamics runs seems to be higher than that obtained by a comparison of independent umbrella sampling runs according to the lesser and more uniformly scattered values of d1,2. Finally, the d1,2 correlation coefficient between the averaged free-energy profile provided by metadynamics and that from the umbrella sampling method is 1.77 kJ/mol-1. This corresponds to the lower similarity than that found for the independent simulations in both methods. Such a difference may also indicate the different ways of reconstructing the free-energy profile. The landscapes in metadynamics are rendered from smooth Gaussian functions, but umbrella sampling and the successive WHAM analysis work with sharply binned histograms. However, the value of 1.77 kJ mol-1 is still lower than the magnitude of thermal noise (RT300K ) 2.49 kJ mol-1), and one can hence assume the physical meaning of such



TABLE 1: Free-Energy Differences (in kJ mol-1) between Individual Conformers As Obtained by Metadymics (a) and Umbrella Sampling with WHAM (b) and the Corresponding Conformational Populations (c) (%)a (a) conformers

RR

R′

β

C5

RL

RR R′ b C5 RL RD

5.0 ( 0.2 -0.2 ( 0.3 2.0 ( 0.2 14.5 ( 0.3 17.5 ( 0.3

-5.2 ( 0.2 -3.0 ( 0.2 9.4 ( 0.2 12.5 ( 0.3

2.3 ( 0.2 14.7 ( 0.2 17.7 ( 0.3

12.4 ( 0.3 15.5 ( 0.3

3.1 ( 0.3

RD

(b) conformers

RR

R′

β

C5

RL

RR R′ B C5 RL RD

4.4 ( 0.2 -0.1 ( 0.2 1.5 ( 0.5 15.3 ( 0.5 17.9 ( 0.2

-4.5 ( 0.3 -2.9 ( 0.5 10.9 ( 0.6 13.5 ( 0.3

1.6 ( 0.3 15.4 ( 0.5 18.0 ( 0.3

13.8 ( 0.7 16.4 ( 0.5

2.6 ( 0.5

RD

(c) conformers/method

RR

R′

β

C5

RL

RD

metadynamics umbrella sampling

37.5 ( 2.6 36.3 ( 2.5

5.0 ( 0.3 6.2 ( 0.6

40.8 ( 2.4 37.3 ( 0.9

16.6 ( 0.9 20.2 ( 2.7

0.11 ( 0.01 0.08 ( 0.02

0.03 ( 0.0 0.03 ( 0.0

The values are reported as x ( σ, where the mean (x) and standard deviation (σ) were calculated from 10 independent simulations for both methods. a

TABLE 2: Similarity of the 2D Free-Energy Landscapes Measured by the d1,2 Correlation Coefficient (in kJ mol-1)a (a) simulation

average

#1

#2

#3

#4

#5

#6

#7

#8

#9

average #1 #2 #3 #4 #5 #6 #7 #8 #9 #10

0.40 0.42 0.42 0.39 0.40 0.40 0.37 0.40 0.41 0.40

0.58 0.61 0.56 0.62 0.55 0.57 0.61 0.66 0.61

0.62 0.58 0.64 0.60 0.63 0.63 0.60 0.62

0.60 0.64 0.63 0.61 0.55 0.64 0.64

0.57 0.60 0.59 0.61 0.59 0.61

0.59 0.56 0.58 0.59 0.58

0.52 0.62 0.65 0.63

0.57 0.60 0.56

0.63 0.57

0.53

simulation

average

#1

#2

#3

#4

#5

#6

#7

#8

#9

average #1 #2 #3 #4 #5 #6 #7 #8 #9 #10

0.70 1.30 0.76 0.77 0.57 0.67 0.70 0.70 0.72 0.66

1.79 0.72 1.17 1.00 1.11 1.06 0.90 0.70 1.10

1.92 1.60 1.33 1.22 1.52 1.77 1.81 1.23

1.10 1.06 1.13 1.09 0.91 0.68 1.19

1.02 1.10 1.02 1.00 1.22 1.06

0.93 0.87 1.00 1.00 0.82

1.06 1.15 1.12 0.84

0.98 1.03 1.14

0.89 1.06

1.13

#10

(b) #10

a

Every free-energy landscape that has resulted from one simulation is compared to that of the others, including the average profile over all of the simulations.

a difference, just like its impact on the modeling of experimentally accessible quantities. We have to stress that our comparison of both methods for free-energy calculations is strictly connected to the simulation protocols and settings used. One can harmonize both the accuracy and efficiency of the two methods by the selection of optimal parameters. We prefer metadynamics also due to the

better computational efficiency (length of simulations) for the parameters considered as “optimal”. Effect of Different Explicit Water Models on the FreeEnergy Landscape of Alanine Dipeptide. To evaluate the effects of different explicit water models on a free-energy landscape, the metadynamics of alanine dipeptide with different water models was performed. The conformational free-energy

5638



Figure 3. The plots of the conformational free energy for alanine dipeptide (ff03) solvated with distinct explicit water models (SPCE, TIP3P, TIP4P, TIP4P/Ew, and TIP5P). One-dimensional projections were constructed for the φ (a) and ψ (b) dihedrals.

TABLE 3: Selected Conformational Transitions and Their ∆F (in kJ mol-1) for Different Water Models (a) and Percentage of the Population of Individual Conformers (b)a (a) SPCE

TIP3P

TIP4P

RR f R′ 5.5 ( 0.1 4.9 ( 0.1 5.3 ( 0.1 RR f β -0.1 ( 0.1 -0.1 ( 0.1 -0.5 ( 0.1 2.4 ( 0.1 1.9 ( 0.1 1.6 ( 0.1 RR f C5 RR f RL 15.1 ( 0.2 14.8 ( 0.1 14.8 ( 0.3 RR f RD 17.1 ( 0.2 17.6 ( 0.2 17.3 ( 0.3 β f C5 2.5 ( 0.1 2.0 ( 0.1 2.1 ( 0.1

TIP4P/Ew

TIP5P

5.4 ( 0.2 0.1 ( 0.2 2.3 ( 0.1 14.8 ( 0.3 17.2 ( 0.1 2.2 ( 0.1

5.1 ( 0.2 -0.6 ( 0.2 2.2 ( 0.2 15.1 ( 0.2 18.4 ( 0.4 2.8 ( 0.2

TABLE 4: Similarity of the Free-Energy Profiles Obtained by Metadynamics for Alanine Dipeptide (ff03) Accompanied by Different Water Modelsa water model

SPCE

TIP3P

TIP4P/Ew

TIP4P

SPCE TIP3P TIP4P/EW TIP4P TIP5P

0.65 0.61 0.81 1.22

0.70 0.77 1.11

0.70 1.08

1.14

TIP5P

a The reported values correspond to d1,2 correlation coefficients in kJ mol-1.

(b) water model/ conformer

SPCE

TIP3P

TIP4P

TIP4P/Ew

TIP5P

RR R′ B C5 RL RD

39.4 ( 1.2 4.4 ( 0.2 41.1 ( 1.2 15.0 ( 0.5 0.09 ( 0.01 0.04 ( 0.00

37.5 ( 2.6 5.0 ( 0.3 40.8 ( 2.4 16.6 ( 0.9 0.11 ( 0.01 0.03 ( 0.00

34.9 ( 1.0 4.1 ( 0.1 42.6 ( 0.8 18.2 ( 0.5 0.09 ( 0.01 0.03 ( 0.00

40.5 ( 1.7 4.6 ( 0.2 38.7 ( 1.6 16.2 ( 0.4 0.11 ( 0.01 0.04 ( 0.00

35.4 ( 1.9 4.6 ( 0.2 45.3 ( 1.7 14.6 ( 0.9 0.08 ( 0.00 0.02 ( 0.00

a The reported uncertainties were estimated as standard deviations by the block average method from one simulation.

landscapes were reconstructed in two dimensions, with the dihedral angles φ and ψ taken as collective variables. We studied the effects of the five different explicit water model parameters, namely, SPCE, TIP3P, TIP4P, TIP4P/Ew, and TIP5P, on the free-energy landscapes. The resulting free-energy plots did not provide any apparent topological differences between the examined water models; the presence and positions of all of the minima, barriers, and saddle points remained well preserved (not shown). Both the qualitative and quantitative agreement is demonstrated in the one-dimensional plots in Figure 3. Both plots are aligned to zero for the global minimum on the curves; therefore, their optimum fit has not been attained (see Figure 3b). The free-energy differences among the selected conformers are listed in Table 3. The relative preference of the individual conformers remained unaltered, but the order of the two mostfavored conformers was not possible to distinguish at the relevant statistical level. Moreover, the standard deviations of the resulting values estimated by the block average method seem to be underestimated in comparison with those calculated from independent simulations (see the table in the previous section). The overall alignment of the resulting free-energy landscapes by means of the d1,2 correlation coefficient is presented in Table 4. According to the obtained values, TIP5P water provided the most distinct free-energy landscape. The differences produced by other models were ascribed to the intrinsic uncertainty of

our protocol as estimated in the previous section (0.6 kJ mol-1). However, the largest found deviation (1.22 kJ mol-1) exceeded this limit by two times, which is still two times smaller than the magnitude of thermal noise (2.49 kJ mol-1). Comparison of Force Fields. Four independent metadynamics runs were performed for a comparison of the four allatomic force fields available in the Gromacs package, namely, ff03, ff99SB, opls-aa, and charmm27. We tested their performance on the alanine dipeptide benchmark molecule in an explicit solvent, namely, TIP3P water. The free-energy landscapes in the φ and ψ coordinates served as indicators of similarity or diversity between individual force field parametrizations. The resulting landscapes are shown in Figure 4. The topological characteristics generally agreed reasonably well for ff03, ff99SB, and opls-aa. The remarkable exception is the charmm27 landscape, which differs in the position and shape of the minima regions in the upper left-hand corner of the calculated profile. The extended conformer (C5) does not seem to be separated by a barrier from the β conformer (polyproline helix-like) in the charmm27 landscape shown. The striking differences between the force fields were found for their amplitudes and the shape of the energy barriers and also for the saddle points, as arises from the two-dimensional plots (A,B,C,D) in Figure 4. The same trend is apparent on the effective free-energy profile along the individual dihedrals φ and ψ on plots (E) and (F). As a result, the minimum-energy paths between the found minima can be distinguished for every investigated force field. For example, the transition barrier between RR and β is lower along the ψ coordinate around -120° for charmm27, contrary to the others, where the other barrier around 90° is more feasible to overcome (Figure 4F). Further differences involve the minimal-energy paths to other higherenergy regions RR and RD passing through distinctive saddle points in all of the force fields. The representative conformational free-energy differences are listed in Table 5. Conformers RR and β were found to be almost



Figure 4. The free-energy landscapes of alanine dipeptide as obtained from metadynamics with different force fields (all with TIP3P water). (A-D) Two-dimensional plots of the conformational free energy depending on the φ and ψ dihedrals for ff03, ff99SB, opls-aa, and charmm27, respectively. The contours are drawn at each 2 kJ mol-1. The one-dimensional profiles along individual dihedrals are provided in (E, F).

energetically indistinguishable in all of the force fields, but corresponding probabilities on the free-energy profiles show differences more sensitively (see the population in Table 5). Conformer β was slightly more preferred to RR in ff99SB and opls-aa. The relative conformational energy of the other minima regions differs substantially in all of the force fields, especially for RR and RD, as demonstrated in Table 5. The mutual overall comparison of the obtained landscapes by means of the coefficient d1,2 is reported in Table 6. Interestingly, the highest similarity among free-energy surfaces in the different force fields was found for opls-aa and ff03 (5.58 kJ mol-1). The values of d1,2 obtained for the other free-energy landscapes

compared were approximately 2 times larger. However, all of the values of similarity coefficients exceeded the level of thermal noise by approximately 2-4 times and also surpassed the effects coming from the usage of different water models. Because of the significant differences in the free-energy profiles obtained from the distinctive force fields, investigation of inherent profiles of the potential energy in the absence of solvent was performed. The relaxed PESs of alanine dipeptide as a function of the φ and ψ dihedrals for all four of the studied force fields are presented in Figure 5. Apparently, the obtained PESs correspond to a gas-phase environment with well-resolved C5, C7ax, and C7eq minima in all of the cases. The numerical

5640



TABLE 5: Selected Conformational Transitions and Their ∆F (in kJ mol-1) for Different Force Fields (all with TIP3P water) (a) and Percentage of the Corresponding Population of Conformers (b)a (a) RR f R′ RR f β RR f C5 R R f RL R R f RD β f C5

ff03

ff99SB

opls-aa

charmm27

4.9 ( 0.1 -0.1 ( 0.1 1.9 ( 0.1 14.8 ( 0.1 17.6 ( 0.2 2.0 ( 0.1

3.2 ( 0.1 -1.7 ( 0.2 -0.1 ( 0.1 5.6 ( 0.1 11.8 ( 0.1 1.6 ( 0.1

2.4 ( 0.1 -1.3 ( 0.1 1.5 ( 0.2 9.5 ( 0.1 14.7 ( 0.2 2.8 ( 0.1

5.9 ( 0.2 0.5 ( 0.2 5.0 ( 0.2 21.9 ( 0.2 17.6 ( 0.2 4.5 ( 0.1

(b) force field/ conformer

ff03

ff99SB

opls-aa

charmm27

RR R′ B C5 RL RD

37.5 ( 2.6 5.0 ( 0.3 40.8 ( 2.4 16.6 ( 0.9 0.11 ( 0.01 0.03 ( 0.00

22.7 ( 0.8 6.4 ( 0.3 44.6 ( 1.3 23.8 ( 0.8 2.39 ( 0.1 0.20 ( 0.01

27.6 ( 1.2 10.5 ( 0.2 46.1 ( 0.7 15.1 ( 0.9 0.6 ( 0.03 0.08 ( 0.01

48.7 ( 2.2 4.5 ( 0.2 40.2 ( 1.9 6.6 ( 0.3 0.01 ( 0.00 0.04 ( 0.00

a The reported uncertainties were estimated as standard deviations by the block average method from one simulation.

TABLE 6: Similarity of the Free-Energy Profiles Obtained by Metadynamics for Alanine Dipeptide (in TIP3P water) Described by Different Force Fieldsa force field

ff03

ff99SB

opls-aa

ff03 ff99SB opls-aa charmm27

10.59 5.58 9.22

11.38 12.98

11.77

charmm27

a

The reported values correspond to d1,2 correlation coefficients in kJ mol-1.

comparison of the potential energy profiles by d1,2 is provided in Table 7. The presented values illustrate the distinctness of the individual force fields, in good agreement with the analogous mutual d1,2 coefficients for solvated free-energy profiles in Table 6. These results clearly indicate that the observed dissimilarity in the FESs of a solvated molecule originates from the different treatment of the intramolecular interactions rather than that of the intermolecular ones with molecules of solvent. In contrast, the latter strongly determine the important character of the freeenergy profile of a solvated molecule. This is obvious from a comparison of the presented PESs and FESs as has already been noticed and reported.38 Surprisingly, the different explicit model of water molecules has a rather negligible impact on the resulting conformational preferences of the studied dipeptide, as was pointed out in the previous section. Discussion The performance of the different water models continues to be one of the most important questions in the field of biomolecular simulations. We performed a direct comparison of the TIP3P, TIP4P, TIP4P/Ew, TIP5P, and SPCE explicit water models and evaluated their impact on the conformational free-energy surface of alanine dipeptide. We have found, surprisingly, that the effect of the model is rather negligible regarding the resulting FES. The results clearly suggest that the interactions of water molecules with alanine dipeptide are described consistently in all of the tested models. This is quite a positive finding because the nonbonded interactions between

the peptide and water are the major agents responsible for the FES profile of the solvated peptide, as is evident from a comparison of the presented PES of the nonsolvated alanine dipeptide. It should be noted that other thermodynamic (e.g., ∆G of solvation) and kinetic properties of the tested water models may differ significantly, as has been reported elsewhere.39 One of the important tasks that we wanted to assess was the accuracy and reliability of our metadynamics simulations in descriptions of the free-energy surface of alanine dipeptide. The commonly used umbrella sampling protocol, which served as a benchmark for our calculations, did not present more consistent and reproducible results than metadynamics. Moreover, one run of metadynamics was able to build up a complete FES in a considerably smaller amount of total computer time than the approximately 100 short umbrella sampling simulations needed for WHAM. The outputs of both methods can be easily compared, and the agreement is reasonably good when the estimated uncertainty of metadynamics and WHAM is taken into account. The quality of the force field essentially underlies the applicability of molecular dynamics and modeling; principally, the reliable and well-balanced parameters are required for the study of peptides and proteins. Alanine dipeptide as a model molecule allows the direct targeting of the intrinsic conformational preferences of the peptide backbone. In our study, we have included the two recently most used amber force fields (ff03 and ff99SB) and two other force fields with a different scheme of parametrization, opls-aa and charmm27. The freeenergy surfaces obtained by the ff03, ff99SB, and opls-aa force fields generally bear a resemblance to each other in the positions of energy minima and barriers. The most significant deviation from these similar surfaces was obtained by the use of charmm27. The similar behavior of charmm27 and the corresponding differences was also found on the potential energy surface. This suggests that it is the distinct force field parametrization of the peptide backbone in charmm27 rather than nonbonded interaction of the peptide with water molecules that is responsible for the manifested dissimilarities. Both amber force fields and opls-aa produced an abundant population in the three major regions of the Ramachandran plot, the extended, helical, and polyproline helix-like (ppII), respectively. The experimental studies16 of short alanine peptides identified a major population in the ppII region and hardly any population in the R-helical region. The currently published computational study on short peptides found almost all of the force fields to be too helical40 while suggesting corresponding corrections. A similar insufficiency was also found by the utilization of many semiempirical methods.19 Surprisingly, despite the relatively poor performance on short peptides, recent force fields have proven to be able to fold small proteins correctly such as Trpcage, villin headpiece, sh3 domain, and many others.41,42 The observed population of the conformers from our metadynamics simulations of alanine dipeptide agrees well with the published results from other laboratories, but a direct and quantitative comparison is still questionable. The conformational preferences derived from a simulation of longer peptides inevitably add the contribution from intermolecular hydrogen bonding to the intrinsic backbone preferences. Such a gentle interplay between individual components should be correctly balanced in a physically correct force field. We believe that modern robust and fast methods for freeenergy calculation such as metadynamics may facilitate the



Figure 5. The potential energy surface (PES) of alanine dipeptide as a function of the φ and ψ dihedrals for ff03 (a), ff99SB (b), opls-aa (c), and charmm27 (d). The color scale is the same as that on the previous plots for the sake of better comparison; the values greater than 40 kJ mol-1 have been left uncolored. The contours are drawn at regular intervals of 2 kJ mol-1.

TABLE 7: Similarity of the Potential Energy Profiles Obtained by the Relaxed Scan of Alanine Dipeptide in Different Force Fieldsa force field

ff03

ff99SB

opls-aa

ff03 ff99SB opls-aa charmm27

16.87 7.35 10.54

16.13 16.01

13.49

charmm27

a

The reported values correspond to d1,2 correlation coefficients in kJ mol-1.

further development of force fields toward a more precise description of peptides and proteins as well as the related phenomena. Conclusions We utilized metadynamics for a conformational study of the solvated alanine dipeptide molecule. Metadynamics proves to be competent as a fast, robust, and reliable method for the conformation free-energy calculations of the peptide in an explicit solvent, surpassing the traditional methods such as umbrella sampling. These results suggests that metadynamics

can be used for detailed mapping of FES for small peptides. Nevertheless, significant progress should be made in choosing the appropriate collective variables for the systems bigger than dipeptides where the two dihedral angles represent the only important degrees of freedom. The most popular explicit water models for the simulation of biomolecules, TIP3P, TIP4P, TIP4P/Ew, TIP5P, and SPCE, showed similar effects on the conformational preferences of alanine dipeptide without any significant differences. However, interaction with water molecules is one of the most important structure-generating elements having a great influence on the FES of the solvated peptide and the conformational preferences of the peptide backbone. All of the tested force fields (ff03, ff99SB, opls-aa, and charmm27) appreciably differ in the populations of individual conformers and the barriers between them. The most discrepancies were found for charmm27, whereas the other mentioned force fields provided a basically similar FES. A significant divergence was also found on the PES calculated by charmm27. Therefore, we can conclude that the differences originate dominantly from the parametrization of the peptide backbone in the given force field rather than from noncovalent interaction with water molecules.

5642


Acknowledgment. This work was supported by Grant No. LC512 from the Ministry of Education, Youth and Sports (MSMT) of the Czech Republic. It was also a part of the research projects No. Z40550506 and MSM6198959216. References and Notes (1) Ozkan, S. B.; Wu, G. A.; Chodera, J. D.; Dill, K. A. Protein folding by zipping and assembly. Proc. Natl. Acad. Sci. U.S.A. 2007, 104 (29), 11987–11992. (2) Camilloni, C.; Sutto, L.; Provasi, D.; Tiana, G.; Broglia, R. A. Early Events in Protein Folding: Is There Something More than Hydrophobic Burst. Protein Sci. 2008, 17 (8), 1424–1433. (3) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (20), 12562–12566. (4) Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free-Energy Landscape for Beta Hairpin Folding from Combined Parallel Tempering and Metadynamics. J. Am. Chem. Soc. 2006, 128 (41), 13435–13441. (5) Piana, S.; Laio, A. A Bias-Exchange Approach to Protein Folding. J. Phys. Chem. B 2007, 111 (17), 4553–4559. (6) Babin, V.; Roland, C.; Darden, T. A.; Sagui, C. the Free Energy Landscape of Small Peptides As Obtained from Metadynamics with Umbrella Sampling Corrections. J. Chem. Phys. 2006, 125 (20), 1-9. (7) Barducci, A.; Chelli, R.; Procacci, P.; Schettino, V.; Gervasio, F. L.; Parrinello, M. Metadynamics Simulation of Prion Protein: Beta-Structure Stability and the Early Stages of Misfolding. J. Am. Chem. Soc. 2006, 128 (8), 2705–2710. (8) Marinelli, F.; Pietrucci, F.; Laio, A.; Piana, S. A Kinetic Model of Trp-Cage Folding from Multiple Biased Molecular Dynamics Simulations. PLOS Comput. Biol. 2009, 5 (8), 1-18. (9) Todorova, N.; Marinelli, F.; Piana, S.; Yarovsky, I. Exploring the Folding Free Energy Landscape of Insulin Using Bias Exchange Metadynamics. J. Phys. Chem. B 2009, 113 (11), 3556–3564. (10) Fiorin, G.; Pastore, A.; Carloni, P.; Parrinello, M. Using Metadynamics to Understand the Mechanism of Calmodulin/Target Recognition at Atomic Detail. Biophys. J. 2006, 91 (8), 2768–2777. (11) Spiwok, V.; Lipovova, P.; Kralova, B. Metadynamics in Essential Coordinates: Free Energy Simulation of Conformational Changes. J. Phys. Chem. B 2007, 111 (12), 3073–3076. (12) Laio, A.; Gervasio, F. L. Metadynamics: A Method to Simulate Rare Events and Reconstruct the Free Energy in Biophysics, Chemistry and Material Science. Rep. Prog. Phys. 2008, 71 (12), 1-22. (13) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. Multidimensional Free-Energy Calculations Using the Weighted Histogram Analysis Method. J. Comput. Chem. 1995, 16 (11), 1339–1350. (14) Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the Equilibrium Boltzmann Distribution from Well-Tempered Metadynamics. J. Comput. Chem. 2009, 30 (11), 1615–1621. (15) Daggett, V.; Kollman, P. A.; Kuntz, I. D. A Molecular-Dynamics Simulation of Polyalanine s An Analysis of Equilibrium Motions and Helix Coil Transitions. Biopolymers 1991, 31 (9), 1115–1134. (16) Graf, J.; et al. Structure and Dynamics of the Homologous Series of Alanine Peptides: A Joint Molecular Dynamics/NMR Study. J. Am. Chem. Soc. 2007, 129 (5), 1179–1189. (17) Kaminsky, J.; Sebek, J.; Bour, P. Molecular Dynamics with Restrictions Derived from Optical Spectra. J. Comput. Chem. 2009, 30 (6), 983–991. (18) Kaminsky, J.; Mata, R. A.; Werner, H. J.; Jensen, F. The Accuracy of Local MP2 Methods for Conformational Energies. Mol. Phys. 2008, 106 (15), 1899–1906. (19) Seabra, G. D.; Walker, R. C.; Roitberg, A. E. Are Current Semiempirical Methods Better Than Force Fields? A Study from the Thermodynamics Perspective. J. Phys. Chem. A 2009, 113 (43), 11938– 11948. (20) Strodel, B.; Wales, D. J. Free Energy Surfaces from an Extended Harmonic Superposition Approach and Kinetics for Alanine Dipeptide. Chem. Phys. Lett. 2008, 466 (4-6), 105–115. (21) Wang, Z.-X.; Duan, Y. Solvation Effects on Alanine Dipeptide: A MP2/cc-pVTZ//MP2/6-31G** Study of (Phi, Psi) Energy Maps and Conformers in the Gas Phase, Ether, And Water. J. Comput. Chem. 2004, 25 (14), 1699–1716.

Vymeˇtal and Vondra´sˇek (22) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Kenneth, M. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26 (16), 1668–1688. (23) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24 (16), 1999–2012. (24) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Bioinform. 2006, 65 (3), 712–725. (25) Mackerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, 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.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102 (18), 3586–3616. (26) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105 (28), 6474–6487. (27) Sorin, E. J.; Pande, V. S. Exploring the Helix-Coil Transition via All-Atom Equilibrium Ensemble Simulations. Biophys. J. 2005, 88 (4), 2472–2493. (28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926–935. (29) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112 (20), 8910–8922. (30) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91 (24), 6269– 6271. (31) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Gordon, T. H. Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120 (20), 9665–9678. (32) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435–447. (33) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18 (12), 1463–1472. (34) Mezei, M. Adaptive Umbrella Sampling s Self-Consistent Determination of the Non-Boltzmann Bias. J. Comput. Phys. 1987, 68 (1), 237– 248. (35) Torrie, G. M.; Valleau, J. P. Non-Physical Sampling Distributions in Monte-Carlo Free-Energy Estimation s Umbrella Sampling. J. Comput. Phys. 1977, 23 (2), 187–199. (36) Roux, B. The Calculation of the Potential of Mean Force Using Computer-Simulations. Comput. Phys. Commun. 1995, 91 (1-3), 275–282. (37) Alonso, J. L.; Echenique, P. A Physically Meaningful Method for the Comparison of Potential Energy Functions. J. Comput. Chem. 2006, 27 (2), 238–252. (38) Avbelj, F.; Grdadolnik, S. G.; Grdadolnik, J.; Baldwin, R. L. Intrinsic backbone preferences are fully present in blocked amino acids. Proc. Natl. Acad. Sci. U.S.A. 2006, 103 (5), 1272–1277. (39) Hess, B.; van der Vegt, N. F. A. Hydration Thermodynamic Properties of Amino Acid Analogues: A Systematic Comparison of Biomolecular Force Fields and Water Models. J. Phys. Chem. B 2006, 110 (35), 17616–17626. (40) Best, R. B.; Buchete, N. V.; Hummer, G. Are Current Molecular Dynamics Force Fields Too Helical. Biophys. J. 2008, 95 (1), L07–L09. (41) Daggett, V. Protein Folding-Simulation. Chem. ReV. 2006, 106 (5), 1898–1916. (42) Scheraga, H. A.; Khalili, M.; Liwo, A. Protein-Folding Dynamics: Overview of Molecular Simulation Techniques. Annu. ReV. Phys. Chem. 2007, 58, 57–83.

JP100950W

Metadynamics As a Tool for Mapping the Conformational and Free

Recommend Documents