Altruistic Metadynamics: Multisystem Biased Simulation - The Journal

Article pubs.acs.org/JPCB

Altruistic Metadynamics: Multisystem Biased Simulation Petr Hošek,† Daniela Toulcová,† Andrea Bortolato,‡ and Vojtěch Spiwok*,† †

Department of Biochemistry and Microbiology, University of Chemistry and Technology, Prague, Technická 3, Prague 6, 166 28, Czech Republic ‡ Heptares Therapeutics Ltd., BioPark, Broadwater Road, Welwyn Garden City AL7 3AX, United Kingdom S Supporting Information *

ABSTRACT: We present a new simple extension of multiple walker metadynamics which makes it possible to simulate simultaneously multiple different molecular systems and to predict their free energy surfaces, named Altruistic metadynamics. Similarly to basic metadynamics, it uses a bias potential in the form of hills summed over the simulation. Each system adds a big hill to its “own” bias potential and smaller hills to bias potentials of other systems, hence, each system enhances sampling of other systems. This makes it possible to achieve either faster reaching of the stationary point or higher accuracy of the calculated free energy surfaces. This should be efficient in modeling of series of chemically similar systems, for example, in computational drug screening by metadynamics. The method was tested on model energy surfaces, alanine dipeptide modeled in different force fields and monosaccharides of D-hexopyranose series.

■

INTRODUCTION Enhanced sampling techniques have been applied in many fields of chemistry as an efficient alternative to conventional molecular simulations. Application of these techniques makes it possible to simulate slow or rarely occurring processes and to predict their free energy surfaces with efficient use of computational resources. One of these techniques is metadynamics.1 Metadynamics enhances sampling by a historydependent bias potential, which can be used a posteriori to predict the free energy surface. Metadynamics enhances sampling by adding a low-dimensional bias potential Vbias to the potential energy Vpot. The bias potential accumulates during the simulation and gradually flattens (or “floods” as a “computational sand”) the free energy surface of the system. In the original version of metadynamics,1,2 the bias potential Vbias is defined as a sum of Gaussian hills gt′ each multiplied by height h deposited at different times t′ during the simulation. Each Gaussian hill is a function of a CV s and is defined by width σ (gt′(s) = exp(−(s − s(t′))2/2σ2); the whole concept can be easily generalized for multiple CVs). The free energy surface F is estimated as a negative image of bias potential: F = −Vbias = − ∑ hgt ′

potential of metadynamics and related techniques in drug discovery or protein engineering. However, drug discovery and protein engineering are highthroughput disciplines, that is, they deal with experimental or computational testing of a large series of compounds, protein mutants, and so on. Efficient application of enhanced sampling techniques in these disciplines therefore requires efficient parallelization. We can speculate, for example, that parallel simulation of binding of individual compounds to the target protein would be more efficient than serial simulations provided that information from a simulation of each compound can be exploited to enhance sampling in simulations of other compounds. Metadynamics can be efficiently parallelized as its multiple walker version.8 The system is simulated by this technique in multiple replicas (walkers). These walkers develop independently, except that they share their history-dependent bias potentials. The bias potential is calculated as a sum over N walkers: N

F = −Vbias = −∑ ∑ hgi , t ′ i = 1 t ′≤ t

where i is the index of a walker. Indexes i in multiple walker metadynamics represent different replicas of the same

(1)

t ′≤ t

This technique has been successfully applied to model formation of a protein−ligand complex,3−5 to fold a protein6 and to solve many other problems.7 This illustrates the great © XXXX American Chemical Society

(2)

Received: January 5, 2016 Revised: February 5, 2016

A

DOI: 10.1021/acs.jpcb.6b00087 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

a control. Ramachandran torsions φ and ψ were used as CVs. Height of hill (h in eq 3) was set to 0.1 kJ/mol and width was set to 0.3 rad. Hills were added every 1 ps (50000 hills). Simulations of D-hexopyranoses used Cremer-Pople ring puckering parameters24 as CVs:

molecular system (with different initial conditions) driven by the same collective variables. Multiple walker metadynamics can be efficiently parallelized thanks to the fact that distribution of the bias potential represents a nonintensive communication between computer nodes. However, this method was designed for parallel simulation of multiple replicas of the same system and does not provide convergent free energy surfaces when individual replicas represent chemically different systems. Here we propose an extension of multiple walker metadynamics to simulate a series of similar yet chemically different systems in parallel using the same set of collective variables (CVs). Because each walker “helps” other walkers to improve their sampling, we refer to this technique as Altruistic metadynamics. It was tested on three systems: model energy profiles, alanine dipeptide modeled in different force fields and D-hexopyranoses.

qx =

qy =

qz =

■

6

⎤ ⎡ 2π cos⎢ (j − 1)⎥ ⎦ ⎣ 3

(4)

⎤ ⎡ 2π sin⎢ (j − 1)⎥ ⎦ ⎣ 3

(5)

1 3

∑ zj

1 3

∑ zj

1 6

∑ zj(−1) j − 1

j=1

6 j=1

6 j=1

(6)

(3)

where j is an index of an atom (the order was O5, C1, C2, C3, C4, C5) and zj is a perpendicular distance of an atom from the plane that fits all six atoms. These CVs were implemented to Plumed by MathEval utility. The system contained one monosaccharide molecule (β-D-allose, β-D-altrose, β-D-glucose, β-D-mannose, β-D-gulose, β-D-idose, β-D-galactose, and β-Dtalose) and 861−863 water molecules (TIP3P) model.16 Force field parameters of were generated by Glycam server (http:// glycam.org/, Glycam version 0625). Force field parameters of βD-idose were taken from other monosaccharides and its atomic charges were calculated by RESP method26 from a wave function calculated at HF/6-31G*//HF/6-31G* level of theory in Gaussian 03.27 Height of hill (h in eq 3) was set to 0.1 kJ/ mol and width was set to 0.01 nm. Hills were added every 1 ps (30000 hills).

The parameter α (Altruistic factor) controls the degree of “altruism” and ranges from 0 to 1. Value of α set to 0 corresponds to a series of independent (non-Altruistic) metadynamics simulations with doubled hill heights. For the special case of two replicas, the value α set to 1 corresponds to a multiple walker metadynamics. Simulation Details. Altruistic metadynamics was first tested on model energy profiles. Calculations were carried out using an ad hoc program written in C. The source code can be found as Supporting Information. Simulations of alanine dipeptide (Ace-Ala-Nme) in different force fields and D-hexopyranoses were done in Gromacs 5.09 with a modified version of Plumed 2.1.0.10 Each system was minimized by L-BFGS method.11 This was followed by a production 50 ns (alanine dipeptide) or 30 ns (D-hexopyranoses) metadynamics simulations. Temperature was kept constant (300 K) by Parrinello-Bussi thermostat.12 Pressure was kept constant (100 kPa) by Berendsen barostat.13 Particlemesh Ewald method14 was used to treat electrostatic interactions. Time step was set to 2 fs. LINCS constraints15 were used to all bonds. Classical metadynamics was performed as Altruistic metadynamics with α set to 0. In alanine dipeptide simulation the system contained one molecule of alanine dipeptide and 864 water molecules (TIP3P model16). Seven Amber family force fields were used (Amber94,17 Amber96,18 Amber99,19 Amber99SB,20 Amber99SB-ILDN,21 AmberGS,22 and Amber0323). During the simulations, we revealed that Amber99SB and Amber99SBILDN are identical for alanine dipeptide, because ILDN corrections are applied to amino acids with longer side chains. Nevertheless, we decided to keep both force fields in the set as

RESULTS The new flooding scheme of Altruistic metadynamics was first tested on model 1D energy surfaces (Figure 1). Testing of novel metadynamics schemes on simple potential surfaces is a practical and widely used approach.1,7,28 It is possible to test metadynamics on flooding of a potential energy surface (rather than free energy surface) provided that the CVs are the only relevant degree of freedom and the (fictive) temperature is small. The shape of energy surfaces in Figure 1 mimics a free energy surface of protein−ligand interaction. The distance x may be viewed as an analogy of a distance between the core of the protein and the ligand. Figure 1A,B represents different systems, A is a “bad” and B is a “good” ligand binding to the same protein. The evolution of the coordinate x was simulated by a simple Monte Carlo walker (see Supporting Information for the source code). Altruistic metadynamics was tested on two walkers (N = 2) and α = 0.5. Such a simulation can be described as a biwalker metadynamics with different hill heights. The first walker adds a higher hill (1.5 h) to the first bias potential and lower hill (0.5 h) to the second bias potential (eq 3). It is possible to explain the flooding scheme in the way that the bias potential formed by hills of height h floods the average free energy surface of both systems, whereas the difference between hill heights (1.5 vs 0.5 h) is responsible for compensation of the differences in free energy surfaces of both systems. Value of α = 0 corresponds to serial non-Altruistic metadynamics simulations. This option is efficient when free energy surfaces of studied systems are significantly different. It is possible to improve sampling by setting α to a value from 0 to 1 when the free

METHODS Altruistic Metadynamics. Altruistic metadynamics proposed here is a variant of a multiple walker metadynamics.8 In altruistic metadynamics, the indexes i correspond to different molecular systems driven by the same collective variables (e.g., different ligands binding to the same target protein). For each jth system, the bias potential is defined and the free energy surface can be calculated as F = −Vbias ⎡ α = −⎢(2 − α) ∑ hgj , t ′ + ⎢⎣ N−1 t ′≤ t

N

∑ i = 1, i ≠ j

⎤ ∑ hgi ,t′⎥⎥ t ′≤ t ⎦

■

B


Article


Figure 2. Schematic view of flooding schemes used in this study.

“computational sand”. Figures of corresponding energy surfaces flooded by the bias potential can be obtained in Supporting Information (Figures S1−S3). Figure 3A shows the evolution of errors (i.e., root-mean-square deviations of predicted vs exact energy difference). High errors at the beginning of simulation correspond to flooding of individual minima. Low errors in the second half of simulation correspond to the state when both minima are flooded and the bias potential distributes uniformly. The classical metadynamics reached the stationary point faster (in 600 vs 1000 steps). However, this difference disappeared when all simulations started from the same energy minimum (see Supporting Information, Figure S4). Interestingly, accuracy was increasing with the value of Altruistic factors as shown by boxplots of errors (Figure 3B). Errors clearly show a decreasing trend (P-value < 10−13 for the slope). Analogously we carried out a series of simulation with same values of (2 − α)h in eq 3 (scheme 2 in Figure 2). This means that height of a hill acting of its “own” systems was independent of Altruistic factor. Results are depicted in Figure 4. These simulations reached the stationary point faster for high Altruistic factor (Figure 4A), but their accuracies were comparable and did not show any significant trend (Figure 4B, P-value = 0.466 for the slope). These two experiments show that one can either achieve better accuracy or faster reaching of the stationary point by means of Altruistic metadynamics. Next, the method was tested on conformational modeling of solvated alanine dipeptide (Ace-Ala-NMe) modeled in different force fields. Alanine dipeptide is a popular test case for enhanced sampling techniques due to its small size and energy surface similar to Ramachandran diagram of proteins. Dependency of its free energy surface in the space of Ramachandran torsions on the force field has been documented.29 Instead of modeling multiple different molecular systems, we modeled the same system in seven different force fields, so the potential and free energy surfaces were different yet similar. Simultaneous simulation of alanine dipeptide in six Amber force fields (in fact, the number of force fields was seven, but Amber99SB and Amber99SB-ILDN are identical for alanine dipeptide) by 50 ns Altruistic metadynamics with α set to 0.5 gave an excellent agreement with the conventional metadynamics (Figure 5). Simulations with other values of Altruistic factor can be found in Supporting Information (Figure S5). Similarly to the test on model energy profiles, we observed that the accuracy was increasing with the value of Altruistic factor α (P-value = 0.0008 for the slope, Figure 6A) when the same amount of “computational sand” (scheme 1 in Figure 2) was used.

Figure 1. Pair of model 1D energy surfaces (A, B) and their flooding by altruistic metadynamics with time evolution of coordinates x (C). Phases i−iv are indicated.

energy surfaces are similar. However, in the case of value α = 1 it reaches the stationary point only for identical systems. The random walk started from the minimum at x ∼ 3.5 for both systems. This minimum, which corresponds to the bound state, was flooded in the phase i by both walkers. Then the bad ligand escaped this minimum and migrated to the minimum corresponding to the unbound state. In the phase ii the good ligand was bound and it flooded the bound state minimum (by 1.5 h for the good ligand’s surface and 0.5 for the bad ligand’s surface). Similarly, the bad ligand was unbound and flooded the unbound state minimum (by 1.5 h for the bad ligand’s surface and 0.5 for the good ligand’s surface). The good ligand escaped the minimum in phase iii and starting from phase iv the energy profile was flattened and both ligand can freely migrate by diffusion. This clearly shows that altruistic metadynamics makes it possible to calculate both free energy surfaces simultaneously. In order to further explore advantages of Altruistic metadynamics we performed a similar series of simulations of 25 systems with random depths of energy minima. These simulations were performed with two schemes differing in hill heights (Figure 2). First, it was performed with same heights of hills h in eq 3 (scheme 1 in Figure 2). This means that total sum of bias potential (amount of “computational sand”) was the same independently on Altruistic factor. Second, it was performed so that the values of (2 − α)h were the same (scheme 2 in Figure 2). This means that the maximum height of a hill acting on systems was the same independently on Altruistic factor. The studied model energy profiles were same in both series of simulations. Accuracy of Altruistic metadynamics was evaluated by comparison of predicted and known energy difference between the two minima (energy(6.5) − energy(3.5)). Figure 3 shows evolution of errors (A) and overall accuracy (B) of predicted energy differences for the scheme 1, that is, the same amount of C


Article


Figure 3. Accuracy of Altruistic metadynamics on 25 model energy profiles. Heights of hills h were same for all simulations. Errors were evaluated as root-mean-square deviations of predictions from known energy differences between the two minima (energy(6.5) − energy(3.5)). Evolution of errors of classical (α = 0.0) and Altruistic (α = 0.5 and 0.9) metadynamics (A). Errors were calculated for each system in 200 bins (each having 10 steps). Overall accuracy of classical (α = 0.0) and Altruistic (α = 0.1 to 0.9) metadynamics (B). Errors were calculated for each system in steps 1000−2000. They are depicted as Tukey boxplots for each value of α.

Figure 4. Accuracy of Altruistic metadynamics on 25 model energy profiles. Heights (2 − α)h were the same for all simulations. Errors were evaluated analogously to Figure 3. Evolution of errors of classical (α = 0.0) and Altruistic (α = 0.5 and 0.9) metadynamics (A). Errors were calculated for each system in 200 bins (each having 10 steps). Overall accuracy of classical (α = 0.0) and Altruistic (α = 0.1 to 0.9) metadynamics (B). Errors were calculated for each system in steps 1500−2000. They are depicted as Tukey boxplots for each value of α.

Figure 5. Free energy surface of alanine dipeptide in different force fields simulated by conventional and Altruistic metadynamics (α = 0.0 and 0.5, N = 7).

D


Article


Figure 6. Errors of Altruistic metadynamics of alanine dipeptide (A) and β-D-hexopyranoses (B) with different values of Altruistic factor α. Inaccuracies were calculated for alanine dipeptide as standard deviations of the free energy difference between αR and C5 minima (F(−1.1,−0.3) − F(−2.6,+2.8), in radians) calculated over the last 25 ns. Inaccuracies for β-D-hexopyranose were calculated as standard deviations of the free energy differences between 4C1 and 1C4 minima (F(0,0,0.05) − F(0,0,−0.05), in nm) calculated over last 10 ns. Errors were depicted as Tukey boxplots.

Figure 7. Free energy surface of monosaccharides of β-D-hexopyranose series simulated by conventional and Altruistic metadynamics (α = 0.0 and 0.5, respectively, N = 8). Free energy surfaces are depicted as isosurfaces at +40 kJ/mol relative to a global minimum.

extension of a multiple walker metadynamics and requires relatively small changes in a code for metadynamics simulation and data analysis. In order to show that the method provides accurate free energy estimates it is necessary to show that i. it is capable to reach a stationary point and ii that the bias potential can be used to predict the free energy surface. First, we demonstrated that Altruistic metadynamics evolves toward the stationary point where the sum of energy surface and bias potential is constant over CVs (except for a noise) and the bias potential added after this point distributes uniformly over CVs. Figure 1 shows that Altruistic metadynamics reached the point when the sum of energy surface and bias potential was almost perfectly horizontal. Similar reaching of a stationary point was observed for other simulations of a model and real system (Figures 3−5, 7, S1−S4). Reaching of the stationary point by Altruistic metadynamics was demonstrated in all cases (Figures 1, 3, 4, 7, S1−S4) by computational experiments. A mathematical proof of Altruistic metadynamics convergence toward a stationary point is highly challenging and we have not gathered enough evidence to provide it in this article. Nevertheless, it can be demonstrated that the bias potential distributes uniformly over CVs once the stationary point is reached. At this point the free energy surface is flattened and in turn the CV space is sampled uniformly. Uniform sampling of the CV space (uniform addition of Gaussian hills) leads to uniform addition of the bias potential as apparent from eq 3 (the newly added bias potential is a weighted sum of hills and weighted sum of constants is also

Another molecular test case was a series of saccharides of the β-D-hexopyranose series. Equilibrium between chair, boat and skew-boat forms of these monosaccharides is important for the function of many carbohydrate processing enzymes or structure of glycosaminoglycans. In previous studies these equilibria were modeled by metadynamics with two or three Cremer-Pople coordinates as CVs.24,30 Energetically favorable conformations of a six-membered ring lie on a sphere. North and south poles of the sphere correspond to 4C1 and 1C4 chairs, respectively. Equator corresponds to boat and skew-boat conformations. Here, we used altruistic metadynamics to calculate free energy surfaces of all eight monosaccharides in one parallel calculation. Figure 7 shows a comparison of the results of altruistic (α = 0.5) and conventional metadynamics (α = 0). 3D free energy surfaces are depicted as isosurfaces at +40 kJ/mol relative to the global minimum, that is, the equilibrium probability of points located on the surface of the 3D blobs is exp(−40/kT) (in kJ/ mol) relative to the global minimum. This comparison demonstrates an excellent agreement with serial metadynamics. Again we observed accuracy improvement with increasing Altruistic factor α (P-value = 0.005 for the slope, Figure 5B) when the same amount of “computational sand” was used (scheme 1 in Figure 2).

■

DISCUSSION This paper presents a new enhanced sampling technique which makes it possible to efficiently simulate multiple systems in parallel and to predict their free energy surfaces. It is a simple E


Article

The Journal of Physical Chemistry B constant). Further flooding therefore does not introduce any bumps to the flooded profile. Second, it is desirable to show that the bias potential can be used to predict the free energy surface. We have demonstrated by a computational experiment that free energy surfaces of alanine dipeptide and D-hexopyranoses were predicted by Altruistic metadynamics in a great agreement with the conventional metadynamics, which is a standard method of free energy prediction. Moreover, it can be shown theoretically that the negative value of Altruistic bias potential approximates the free energy surface. Dama and co-workers31,32 formulated a general scheme for evaluation of convergence of metadynamics and the relationship between the bias potential and free energy. Stationary points of metadynamics are states satisfying:

∫ g(s , s′)H(s′)e−[F(s′)+V (s′)]/kT ds′ = C1

decrease of inaccuracy approximately from 1.1 to 0.7 kJ/mol, both without additional computing costs. We see a great potential of Altruistic metadynamics in drug design and protein engineering where a series of systems (drug candidates or protein mutants) can be modeled by metadynamics with the same set of collective variables. Other type of simulation that can benefit from Altruistic flooding are those where a series of systems are simulated with a Hamiltonian gradually altered along some parameter λ (in the spirit of free energy perturbation).33 We also believe that Altruistic metadynamics may inspire development of other flooding schemes working on multiple molecular systems in parallel. The capability of the method to improve sampling or to improve accuracy may be further strengthen by a combination with well tempered metadynamics.28 Development of a well tempered Altruistic biasing scheme is a nontrivial task and is not covered by this manuscript.

(7)

■

,where g(s,s′) is a Gaussian hill centered in s′ (g(s,s′) = h exp(−(s − s′)2/2σ2)), H(s′) is a hill scaling factor (equal to 1 in the classical metadynamics and to exp(−V(s′)/ΔT) in the welltempered metadynamics28), F is free energy, V is bias potential, and C1 is a constant. This condition is satisfied in classical metadynamics when F(s′) = − V(s′) + const, because ∫ g(s,s′)ds′ is a constant function over s. In altruistic metadynamics with two walkers, constant hill heights and altruistic factor 0.5, it is possible to write this condition as

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b00087. Source code and the results of additional Altruistic metadynamics simulations (PDF).

■

1

2

+ 0.5e−[F2(s ′) + 0.5Σg1(s ′) + 1.5Σg2(s ′)]/ kT ⎤⎦ds′ = C2

*E-mail: [email protected]. Notes

(8)

The authors declare no competing financial interest.

■

One possible solution of this equation is that both exponential functions are equal to one, that is, F1(s′) = −1.5 ∑g1(s′) − 0.5 ∑g2(s′) + const and F2(s′) = −0.5 ∑g1(s′) − 1.5 ∑g2(s′) + const. It is also possible to consider the solution that exponential functions are not equal to one and after multiplication by corresponding factors (1.5 and 0.5) they become mutual mirror images:

ACKNOWLEDGMENTS This project was supported by Czech Ministry of Education, Youth and Sport via COST actions GLISTEN (CM1207, LD14133) and Czech Science Foundation (15-17269S). Computational resources were provided by MetaCentrum (LM2010005), CERIT-SC (CZ.1.05/3.2.00/08.0144) and IT4Innovations Centre of Excellence project (CZ.1.05/ 1.1.00/02.0070, LM2011033).

1.5e−[F1(s ′) + 1.5Σg1(s ′) + 0.5Σg2(s ′)]/ kT + 0.5e−[F2(s ′) + 0.5Σg1(s ′) + 1.5Σg2(s ′)]/ kT = C3

(9)

■

However, an analogous equation can be derived also for the second walker:

REFERENCES

(1) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (2) Gervasio, F. L.; Laio, A.; Parrinello, M. Flexible Docking in Solution Using Metadynamics. J. Am. Chem. Soc. 2005, 127, 2600− 2607. (3) Limongelli, V.; Bonomi, M.; Parrinello, M. Funnel Metadynamics as Accurate Binding Free-Energy Method. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 6358−6363. (4) Herbert, C.; Schieborr, U.; Saxena, K.; Juraszek, J.; De Smet, F.; Alcouffe, C.; Bianciotto, M.; Saladino, G.; Sibrac, D.; Kudlinzki, D.; et al. Molecular Mechanism of SSR128129E, An Extracellularly Acting, Small-Molecule, Allosteric Inhibitor of FGF Receptor Signaling. Cancer Cell 2013, 23, 489−501. (5) Mason, J. S.; Bortolato, A.; Weiss, D. R.; Deflorian, F.; Tehan, B.; Marshall, F. H. High End GPCR Design: Crafted Ligand Design and Druggability Analysis Using Protein Structure, Lipophilic Hotspots and Explicit Water Networks. In Silico Pharmacol. 2013, 1, 23. (6) 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, e1000452. (7) Spiwok, V.; Hošek, P.; Šućur, Z. Enhanced Sampling Techniques in Biomolecular Simulations. Biotechnol. Adv. 2015, 33, 1130−1140.

0.5e−[F1(s ′) + 1.5Σg1(s ′) + 0.5Σg2(s ′)]/ kT + 1.5e−[F2(s ′) + 0.5Σg1(s ′) + 1.5Σg2(s ′)]/ kT = C3

(10)

eqs 9 and 10 give: e−[F1(s ′) + 1.5Σg1(s ′) + 0.5Σg2(s ′)]/ kT = e−[F2(s ′) + 0.5Σg1(s ′) + 1.5Σg2(s ′)]/ kT

AUTHOR INFORMATION

Corresponding Author

∫ g(s , s′)⎡⎣1.5e−[F (s′)+1.5Σg (s′)+0.5Σg (s′)]/kT 1

ASSOCIATED CONTENT

S Supporting Information *

(11)

This means that F1(s′) = −1.5∑g1(s′) − 0.5∑g2(s′) + const, and F2(s′) = −0.5∑g1(s′) − 1.5∑g2(s′) + const is the only solution of eq 7. Thus, negative value of the bias potential in each system of Altruistic metadynamics approximates its free energy surface. On the model energy profiles we have demonstrated that Altruistic metadynamics can either improve convergence of calculated free energy differences (i.e., the state when all minima are flooded is reached earlier) or improve their accuracy (i.e., fluctuations of free energy differences are smaller). We have confirmed this for solvated alanine dipeptide, where setting of Altruistic factor to 0.9 lead to decrease of inaccuracy approximately from 1.5 to 0.9 kJ/mol, or for Dhexopyranoses, where setting of Altruistic factor to 0.9 lead to F


Article


Kudin, K. N.; Burant, J. C. et al. Gaussian 03, Revision C.02;Gaussian, Inc.: Wallingford, CT, 2004. (28) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. (29) Vymětal, J.; Vondrásě k, J. Metadynamics as a Tool for Mapping the Conformational and Free-Energy Space of Peptides − the Alanine Dipeptide Case Study. J. Phys. Chem. B 2010, 114, 5632−5642. (30) Davies, G. J.; Planas, A.; Rovira, C. Conformational Analyses of the Reaction Coordinate of Glycosidases. Acc. Chem. Res. 2012, 45, 308−316. (31) Dama, J. F.; Parrinello, M.; Voth, G. A. Well-Tempered Metadynamics Converges Asymptotically. Phys. Rev. Lett. 2014, 112, 240602. (32) White, A. D.; Dama, J. F.; Voth, G. A. Designing Free Energy Surfaces that Match Experimental Data with Metadynamics. J. Chem. Theory Comput. 2015, 11, 2451−2460. (33) Gil-Ley, A.; Bussi, G. Enhanced Conformational Sampling Using Replica Exchange with Collective-Variable Tempering. J. Chem. Theory Comput. 2015, 11, 1077−1085.

(8) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Efficient Reconstruction of Complex Free Energy Landscapes by Multiple Walkers Metadynamics. J. Phys. Chem. B 2006, 110, 3533− 3539. (9) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (10) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185, 604−613. (11) Byrd, R. H.; Lu, P.; Nocedal, J. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM J. Sci. Comput. 1995, 16, 1190−1208. (12) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity-Rescaling. J. Chem. Phys. 2007, 126, 014101. (13) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (14) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N.log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089. (15) 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, 1463−1472. (16) 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, 926−935. (17) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (18) Kollman, P. A. Advances and Continuing Challenges in Achieving Realistic and Predictive Simulations of the Properties of Organic and Biological Molecules. Acc. Chem. Res. 1996, 29, 461−469. (19) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (20) 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., Genet. 2006, 65, 712−725. (21) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (22) García, A. E.; Sanbonmatsu, K. Y. Alpha-Helical Stabilization by Side Chain Shielding of Backbone Hydrogen Bonds. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 2782−2787. (23) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (24) Cremer, D.; Pople, J. A. General Definition of Ring Puckering Coordinates. J. Am. Chem. Soc. 1975, 97, 1354−1358. (25) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2008, 29, 622−655. (26) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269−10280. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; G


Altruistic Metadynamics: Multisystem Biased Simulation - The Journal

Recommend Documents