Extended Thermodynamic Integration: Efficient Prediction of Lambda

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Extended Thermodynamic Integration: efficient prediction of lambda derivatives at non-simulated points Anita de Ruiter, and Chris Oostenbrink J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00458 • Publication Date (Web): 05 Aug 2016 Downloaded from http://pubs.acs.org on August 21, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Extended Thermodynamic Integration: efficient prediction of lambda derivatives at non-‐simulated points Anita de Ruiter & Chris Oostenbrink Institute for Molecular Modeling and Simulation University of Natural Resources and Life Sciences (BOKU) Vienna, Austria

Last modified: 04.08.16

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 27

Abstract Thermodynamic Integration (TI) is one of the most commonly used free-‐energy calculation methods. The derivative of the Hamiltonian with respect to lambda, , is determined at multiple λ-‐points. Because a numerical integration step is necessary, high curvature regions require simulations at densely spaced λ-‐points. Here, the principle of extended TI is introduced, where values are predicted at non-‐simulated λ-‐points. Based on three model systems, it is shown that extended TI requires significantly fewer λ-‐points than regular TI, to obtain similar accuracy.

Introduction Over the last few decades, molecular dynamics (MD) simulations have taken an important role in multiple research fields such as physics, material sciences, biotechnology and structural biology. Whereas experimental studies mostly give macroscopic insights, MD simulations have the great advantage that ligands, proteins, DNA and other (bio)materials can be studied at a molecular level. In combination with free-‐energy calculations, MD simulations become even more powerful. This is especially true for the field of drug design where one can e.g. predict binding affinities of different ligands towards a protein or the influence of a mutation in the protein on ligand binding. The challenge in drug design is to achieve the high accuracy that is required to differentiate the binding affinities of closely related molecules and at the same time keep the computational costs as low as possible, since often a large number of ligands needs to be evaluated.1–4 Also, the methods should be easy to use and require as little human input as possible to allow for automated workflows.4 Although end-‐point and pulling methods are also applied for the calculation of (relative) binding free energies, the rigorous alchemical free-‐energy methods Thermodynamic Integration (TI)5 and Bennett’s acceptance ratio (BAR)6,7 are arguably the most appropriate ones. TI determines the free energy difference between a state A and B by making use of multiple intermediate states as defined by a coupling parameter λ. The Hamiltonian is dependent on λ and results in the Hamiltonian for state A when λ=0 and state B when λ=1. Simulations are performed at several λ values and the ensemble averages of the derivative of the Hamiltonian with respect to λ are integrated numerically to obtain the free energy difference.


2

Page 3 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


! !" ∆𝐺!"

= !

𝑑𝐺 λ = 𝑑λ

!

!

𝜕𝐻 λ 𝜕λ

𝑑λ !

eq. 1

The overall variance of a free energy difference obtained by TI is dependent on the curvature of the as function of λ.8 This means that smoother TI curves lead to more accurate free energy differences. The fact that numerical integration is used for TI also makes clear that regions of high curvature require a more dense spacing of simulated λ-‐points. Typically one starts with an equidistant spacing and adds additional simulations if required, as the shape of the curve is not known a priori.9 BAR determines the free energy difference between two values λ and λ+Δλ by !"# ∆𝐺!,!!!" = 𝑘! 𝑇 ln

𝑓 𝐻 λ − 𝐻 λ + 𝛥λ + 𝐶 !!!! + 𝐶 𝑓 𝐻 λ + 𝛥λ − 𝐻 𝜆 − 𝐶 !

eq. 2

where f x =

1 . 1 + e!/!!!

A self-‐consistent solution towards C is iteratively sought, such that the ensemble averages in eq. 2 are equal to each other, which directly gives the free energy difference. Note that the differences in the Hamiltonian are calculated from simulations of both λ-‐points and that a reevaluation of the trajectory to obtain the Hamiltonian for the other λ-‐point is necessary. The use of multiple simulated states simultaneously slightly improves the efficiency and leads to the definition of MBAR.10 Several papers have compared these two methods in terms of efficiency, accuracy and ease of use.9,11,12 BAR seems to be the most efficient one in terms of the number of simulations that are required, but it requires more disk space and post-‐processing time. This is especially true if additional simulations at intermediate λ-‐points are required, because the simulations of the surrounding λ-‐points need to be reevaluated. On the other side, TI requires more simulations in order to get a smooth curve for numerical integration, but only very little post-‐processing which is, moreover, independent of the neighboring λ-‐points. In this paper we describe a way to predict values for surrounding λ points. This allows the calculation of accurate free-‐energy differences with fewer λ-‐points than was previously possible with TI. Similar predictions were already performed by Tidor et al. to improve the


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 27

quadrature of the TI curve for λ-‐points deviating by 0.025 and 0.05 from the simulated values.13 However, the potential energy function that was used in that work, was a linear combination of those of the two end states, leading to the well-‐known end-‐state problems.14 Li and Yang also describe a reweighting method to obtain values for other λ-‐points, named λ-‐WHAM.15 This method was based on a linear combination of the end states, which makes the calculation of the energy differences between two λ-‐points straightforward. In the current study, however, soft-‐ core potentials are used, which makes the prediction of derivatives at non-‐simulated λ-‐points more difficult. Furthermore, we will make predictions over larger ranges of λ. In the following sections, we will show the theory behind this approach, which we refer to as extended TI, its application to three different model systems in water, as well as one more challenging protein-‐ ligand system. Additionally, the accuracy, efficiency and user friendliness will be discussed, as well as further possibilities of the code that are still to be exploited.

Theory As was briefly pointed out in the introduction, the goal of this study is to optimize TI calculations, such that fewer λ-‐points need to be simulated to achieve reliable free-‐energy differences. The idea behind this is shown in Figure 1. From a simulation at e.g. λ=0.2 we will predict what the ensemble average of the derivative of the Hamiltonian with respect to λ would be at other, non-‐ simulated, λ-‐points. We will first focus on the non-‐bonded potential energy terms, as these functions are the most complex due to the soft-‐core potentials.14 In GROMOS, the coupling parameter (µμ) of each interaction type, x, can be set individually based on a fourth order polynomial of λ:16,17 𝜇! = 𝑎! λ! + 𝑏! λ! + 𝑐! λ! + 𝑑! λ + 𝑒!

eq. 3

As an aside, the individual µμx-‐values, may furthermore be specified separately for interactions between different groups of atoms. At the moment there are twelve interaction types that can be treated with individual µμ-‐values, four of which are relevant for the non-‐bonded potential energy.

µμlj and µμcrf scale the Lennard Jones and Coulomb-‐reaction field interactions, respectively, whereas µμslj and µμscrf scale the corresponding softness. Usually, the parameters of eq. 3 are chosen such that µμx = 0 for λ=0 and µμx=1 for λ=1.


4

Page 5 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The λ dependent Lennard Jones potential energy (VLJ) and Coulomb-‐reactionfield interaction (VCRF), can hence be written as 𝑉!" 𝒓, 𝑡, λ =

1 − 𝜇!"

!

! ! ! 𝑉!" 𝑟!" , 𝑡, 𝜇!"# + 𝜇!" 𝑉!" 𝑟!" , 𝑡, 1 − 𝜇!"#

eq. 4

!,!

and 𝑉!"# 𝒓, 𝑡, λ =

1 − 𝜇!"#

!

! ! ! 𝑉!"# 𝑟!" , 𝑡, 𝜇!"#$ + 𝜇!"# 𝑉!"# 𝑟!" , 𝑡, 1 − 𝜇!"#$

eq. 5

!,!

where the sums run over all interacting atom pairs i,j , rij is the distance between the atoms i and

j and n is an integer as set by the user (power dependence of µμ). 𝑉!"! represents the Lennard Jones interaction energy for either end state A or B: ! 𝑉!" 𝑟!" ; 𝑡; 𝜇!"# =

! 𝐶!" (𝑖, 𝑗) 1 − 𝐶!! (𝑖, 𝑗) ∙ ! ! ! ! ! 𝛼!" 𝜇!"# 𝐶!"# 𝑖, 𝑗 + (𝑟!" ) 𝛼!" 𝜇!"# 𝐶!"# 𝑖, 𝑗 + (𝑟!" )!

eq. 6

! where 𝐶!" 𝑖, 𝑗 and 𝐶!! 𝑖, 𝑗 denote the Lennard Jones parameters of atom i and j in state X, ! ! 𝐶!"# 𝑖, 𝑗 = 𝐶!" 𝑖, 𝑗 /𝐶!! 𝑖, 𝑗 and αlj is the softness parameter of the Lennard Jones interaction. ! Similarly, 𝑉!"# represents the Coulomb reaction field interaction energy in state A or B and is

defined by ! 𝑉!"# 𝑟!" ; 𝑡; 𝜇!"#$ =

𝑞!! 𝑞!! 4𝜋𝜀!

1 ! 𝛼!"# 𝜇!"#$ +

! 𝑟!"!

!

− !

! !𝐶!" 𝑟!"

! ! 𝛼!"# 𝜇!"#$ + 𝑅!"

! !

−

1 − ! !𝐶!" 𝑅!"

eq. 7

where 𝑞!! and 𝑞!! are the partial charges on atoms i and j in state X, ε0 is the dielectric permittivity of vacuum, αcrf is the softness parameter of the Coulomb-‐reaction field interaction and Crf and Rrf are parameters of the reaction field contribution. We will now slightly rewrite eq. 4 into 𝑉!" 𝒓, 𝑡, 𝜆 = 1 − 𝜇!"

!

! ! 𝑉!" 𝑟!" , 𝑡, 𝜇!"# + 𝜇!" !,!

! 𝑉!" 𝑟!" , 𝑡, 1 − 𝜇!"# .

eq. 8

!,!

The derivative of the Lennard Jones potential energy with respect to λ is then


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

𝜕𝑉!" 𝒓, 𝑡, λ = −𝑛 1 − 𝜇!" 𝜕λ

+

!!!

1 − 𝜇!"

!!! ! 𝑉!" 𝑟!" , 𝑡, 𝜇!"# + 𝑛𝜇!" !,!

! !,!

Page 6 of 27

! 𝑉!" 𝑟!" , 𝑡, 1 − 𝜇!"# !,!

! 𝜕𝑉!" 𝑟!" , 𝑡, 𝜇!"# ! − 𝜇!" 𝜕𝜇!"#

! 𝜕𝑉!" 𝑟!" , 𝑡, 1 − 𝜇!"# !,!

𝜕 1 − 𝜇!"#

𝑑𝜇!" 𝑑λ

eq. 9

𝑑𝜇!"# 𝑑λ

We can now define ! 𝐸!" 𝒓, 𝑡, 𝜇!"# =

𝑉!!! 𝑟!" , 𝑡, 𝜇!"#

eq. 10

!,! ! 𝐸!" 𝒓, 𝑡, 𝜇!"# =

! 𝑉!" 𝑟!" , 𝑡, 1 − 𝜇!"# !,!

! 𝑑𝐸!" 𝒓, 𝑡, 𝜇!"# = !,! ! 𝑑𝐸!" 𝒓, 𝑡, 𝜇!"# = !,!

! 𝜕𝑉!" 𝑟!" , 𝑡, 𝜇!"# 𝜕𝜇!"#

! 𝜕𝑉!" 𝑟!" , 𝑡, 1 − 𝜇!"# 𝜕𝜇!"#

such that eq. 8 and eq. 9 become !

𝑉!" 𝒓, 𝑡, 𝜆 = 1 − 𝜇!"

! ! ! 𝐸!" 𝒓, 𝑡, 𝜇!"# + 𝜇!" 𝐸!" 𝒓, 𝑡, 𝜇!"#

eq. 11

and 𝜕𝑉!" 𝒓, 𝑡, λ = −𝑛 1 − 𝜇!" 𝜕λ +

!!!

1 − 𝜇!"

!

!!! ! ! 𝐸!" 𝒓, 𝑡, 𝜇!"# + 𝑛𝜇!" 𝐸!" 𝒓, 𝑡, 𝜇!"#

! ! ! 𝑑𝐸!" 𝒓, 𝑡, 𝜇!"# − 𝜇!! 𝑑𝐸!" 𝒓, 𝑡, 𝜇!"#

𝑑𝜇!" 𝑑𝜆

eq. 12

𝑑𝜇!"# 𝑑𝜆

Similarly, we can reorganize eq. 5, and define similar quantities for the Coulomb-‐reaction field interaction: ! 𝐸!"# 𝒓, 𝑡, 𝜇!"#$ =

! 𝑉!"# 𝑟!" , 𝑡, 𝜇!"#$

eq. 13

!,! ! 𝐸!"# 𝒓, 𝑡, 𝜇!"#$ =

! 𝑉!"# 𝑟!" , 𝑡, 1 − 𝜇!"#$ !,!

! 𝑑𝐸!"# 𝒓, 𝑡, 𝜇!"#$ = !,! ! 𝑑𝐸!"# 𝒓, 𝑡, 𝜇!"#$ = !,!

! 𝜕𝑉!"# 𝑟!" , 𝑡, 𝜇!"#$ 𝜕𝜇!"#$

! 𝜕𝑉!"# 𝑟!" , 𝑡, 1 − 𝜇!"#$ 𝜕𝜇!"#$

resulting in: 𝑉!"# 𝒓, 𝑡, λ = 1 − 𝜇!"#

!

! ! ! 𝐸!"# 𝒓, 𝑡, 𝜇!"#$ + 𝜇!"# 𝐸!"# 𝒓, 𝑡, 𝜇!"#$

eq. 14

and


6

Page 7 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝜕𝑉!"# 𝒓, 𝑡, λ = −𝑛 1 − 𝜇!"# 𝜕λ +

1 − 𝜇!"#

!!!

!

!!! ! ! 𝐸!"# 𝒓, 𝑡, 𝜇!"#$ + 𝑛𝜇!"# 𝐸!"# 𝒓, 𝑡, 𝜇!"#$

! ! ! 𝑑𝐸!"# 𝒓, 𝑡, 𝜇!"#$ − 𝜇!"# 𝑑𝐸!"# 𝒓, 𝑡, 𝜇!"#$

𝑑𝜇!"# 𝑑𝜆

eq. 15

𝑑𝜇!"#$ 𝑑𝜆

Up to this point, the above formulas have only been slightly reorganized. To allow for an efficient prediction of at other λ-‐points (or other parameterizations of individual µμ), the quantities in eq. 10 will not only be calculated for µμslj according to the simulated λ-‐point, but for a ! range of values, e.g. 𝜇!"# = 0, 0.01, 0.02, … , 1 . These additional calculations can be performed

efficiently because at the time of calculation the pairlist has already been generated and the lookup of several parameters (e.g. C12 and C6) does not have to be repeated. The sets of quantities ! ! ! ! ! ! ! ! 𝐸!" 𝒓, 𝑡, 𝜇!"# and 𝐸!" 𝒓, 𝑡, 𝜇!"# as well as 𝑑𝐸!" 𝒓, 𝑡, 𝜇!"# and 𝑑𝐸!" 𝒓, 𝑡, 𝜇!"# will be written to the

energy and free energy trajectory files, respectively. The total VLJ(r,t,λ) and ∂VLJ(r,t,λ)/∂λ for other λ values (or other parameterizations of the individual µμ values) can be easily reconstructed from these energy and free energy trajectories, since the remaining factors in eq. 11 and eq. 12 are not dependent on the configuration of the system. A very similar procedure is followed for the Coulomb-‐reaction field interactions, by storing the properties of eq. 13 for different values of

µμscrf. Note that for the application described in the current paper not all energy and free energy terms of eq. 10 and eq. 13 (and the following terms for the kinetic and bonded interactions) need to be written to disk separately. However for future applications this flexibility will be very useful, see Discussion for more details. Changes of atom masses as well as covalent interactions are easier to implement, since these are only influenced by a single µμx parameter. This also means that the corresponding EX and dEX quantities do not have to be written out separately for state A and B. The kinetic energy contributions are determined by the quantity 𝐸! 𝒑, 𝒓, 𝑡, 𝜇! = !

! 𝑚 ! !

𝜇! (𝐯! (𝑡))!

eq. 16

with 𝑚! 𝜇! = 1 − 𝜇! 𝑚!! + 𝜇! 𝑚!! where 𝑚!! is the mass of perturbed atom i in state X, vi is the velocity of atom i and µμK is the individual µμ for the kinetic energy contributions. However, when calculating the kinetic energy


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 27

contributions for other µμK values as the simulation is performed with, one has to keep in mind ! that the momenta, rather than the velocities, are conserved: 𝒑 = 𝑚 𝜆! 𝐯 ! = 𝑚 𝜆!! 𝐯 ! . Thus, the

predicted kinetic energy contributions for 𝜇!! from a simulation at 𝜇!! is calculated by

!

!

𝑚! 𝜇!! 𝑚! 𝜇!!

! !

𝐸! 𝒑, 𝒓, 𝑡, 𝜇!! =

eq. 17

(𝐯! (𝑡))!

where the velocities (and coordinates) correspond to 𝜇!! . For the derivative with respect to µμK, we have to consider the kinetic energy as a function of the momenta, rather than of the velocities, which results in ! !

𝑑𝐸! 𝒑, 𝒓, 𝑡, 𝜇!! = − !

𝑚!! − 𝑚!!

𝐯! (𝑡) !

eq. 18

Again, for predictions of dEK at µμK values other than the simulated one, the velocities need be scaled, because only the momenta are preserved: 𝑑𝐸! 𝒑, 𝒓, 𝑡, 𝜇! = − !

! !

𝑚! 𝜇!! 𝑚! 𝜇!!

𝑚!! − 𝑚!!

!

eq. 19

𝐯! (𝑡) !

A quartic potential energy term is used to describe the covalent bond stretching interaction 𝐸! 𝒓, 𝑡, 𝜇! = ! ! !

where 𝑘!

! !

! !

1 − 𝜇! 𝑘!

! !

+ 𝜇! 𝑘!

𝑏!! (𝑡) −

1 − 𝜇! 𝑏!!! + 𝜇! 𝑏!!!

! !

eq. 20

is the force constant for perturbed bond n in state X, bn is the current bond length

and 𝑏!!! is the ideal bond length of bond n for state X. The corresponding derivative with respect to µμb is dE! 𝒓, 𝑡, 𝜇! = !

! !

! !

−4 𝑘!

! !

+ 𝜇! 𝑘!

∙ 𝑏!!! + 𝜇! 𝑏!!! − 𝑏!!! ! !

+ 𝑘!

! !

− 𝑘!

! !

− 𝑘!

𝑏!!! − 𝑏!!!

eq. 21

𝑏!! (𝑡) − 𝑏!!! + 𝜇! 𝑏!!! − 𝑏!!!

𝑏!! (𝑡) − 𝑏!!! + 𝜇! 𝑏!!! − 𝑏!!!

! !

!

.

The covalent bond angle bending interactions are described by 𝐸! 𝒓, 𝑡, 𝜇! = !

! !

! !

1 − 𝜇! 𝑘!

∙ cos 𝜃! 𝑡 −

! !

+ 𝜇! 𝑘!

1 − 𝜇! cos 𝜃!!! + 𝜇! cos 𝜃!!!


eq. 22 !

8

Page 9 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


! !

where 𝑘!

is the bond angle force constant for bond n in state X, θn is the current bond angle

and 𝜃!!! is the optimal bond angle for bond n in state X. Its derivative with respect to µμθ is 𝑑𝐸! 𝒓, 𝑡, 𝜇! = !

! !

! !

−2 𝑘!

! !

+ 𝜇! 𝑘!

! !

cos 𝜃!!! − cos 𝜃!!!

− 𝑘!

eq. 23

∙ cos 𝜃! (𝑡) − cos 𝜃!!! − 𝜇! cos 𝜃!!! − cos 𝜃!!! ! !

+ 𝑘!

! !

cos 𝜃! (𝑡) − cos 𝜃!!! − 𝜇! cos 𝜃!!! − cos 𝜃!!!

− 𝑘!

!

.

The dihedral-‐angle torsion interactions are described by ! !

𝐸! 𝒓, 𝑡, 𝜇! =

1 − 𝜇! 𝑘!

! !

1 + cos 𝑚!

𝜑! (𝑡) − 𝜑!!!

eq. 24

! ! !

+ 𝜇! 𝑘! ! !

where 𝑘!

! !

1 + cos 𝑚!

𝜑! (𝑡) − 𝜑!!! (!)!

is the force constant of dihedral n in state X, 𝑚!

is the multiplicity of dihedral n in

state X, φn is the current dihedral angle and 𝜑!!! is the phase shift. We can define the quantities 𝐸!! 𝒓, 𝑡 =

! !

1 + cos 𝑚!

! !

1 + cos 𝑚!

𝑘!

! !

𝜑! (𝑡) − 𝜑!!!

! !

𝜑! (𝑡) − 𝜑!!!

!

𝐸!! 𝒓, 𝑡 =

𝑘!

eq. 25

!

and the derivative with respect to µμφ for dihedral-‐angle torsion interactions is ! !

dE! 𝒓, 𝑡 =

𝑘!

! !

1 + cos 𝑚!

𝜑! (𝑡) − 𝜑!!!

eq. 26

! ! !

− 𝑘!

! !

1 + cos 𝑚!

𝜑! (𝑡) − 𝜑!!!

.

This quantity is independent of µμφ and thus only needs to be calculated once. Similar to the quantities in eq. 10 and eq. 13 the Ex and dEx quantities of the kinetic and covalent interactions are written to the energy and free energy trajectory file, respectively, for 𝜇!! = 0, 0.01, 0.02, … , 1 . We then have for the kinetic and covalent interactions 𝐾 𝒑, 𝒓, 𝑡, λ = 𝐸! 𝒑, 𝒓, 𝑡 𝑉! 𝒓, 𝑡, λ = 𝐸! 𝒓, 𝑡, 𝜇! 𝑉! 𝒓, 𝑡, λ = 𝐸! 𝒓, 𝑡, 𝜇! 𝑉! 𝒓, 𝑡, λ = 1 − 𝜇! 𝐸!! 𝒓, 𝑡 + 𝜇! 𝐸!! 𝒓, 𝑡

eq. 27

and


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 27

𝜕𝐾 𝒑, 𝒓, 𝑡, λ 𝜕𝜇! = 𝑑𝐸! 𝒑, 𝒓, 𝑡, 𝜇! 𝜕λ 𝜕λ 𝜕𝑉! 𝒓, 𝑡, λ 𝜕𝜇! = 𝑑𝐸! 𝒓, 𝑡, 𝜇! 𝜕λ 𝜕λ 𝜕𝑉! 𝒓, 𝑡, λ 𝜕𝜇! = 𝑑𝐸! 𝒓, 𝑡, 𝜇! 𝜕λ 𝜕λ 𝜕𝑉! 𝒓, 𝑡, λ 𝜕𝜇! = 𝑑𝐸! 𝒓, 𝑡, 𝜇! 𝜕λ 𝜕λ

eq. 28

After determining the contributions of all the interaction types from a simulation at λS, we can reconstruct the total energy 𝐻 𝒓, 𝑡, λ!

!! ,

as well as its derivative 𝜕𝐻 𝒓, 𝑡, λ!

𝜕λ

!! for

any

other value of λ, indicated by λP: 𝐻 𝒓, 𝑡, λ!

!!

= 𝑉!" 𝒓, 𝑡, λ!

!!

+ 𝑉!"# 𝒓, 𝑡, λ!

+ 𝑉! (𝒓, 𝑡, λ! )

!!

+ 𝐾(𝒑, 𝒓, 𝑡, λ! )

!!

+ 𝑉! (𝒓, 𝑡, λ! )

!!

!!

+ 𝑉! (𝒓, 𝑡, λ! )

!!

eq. 29

and 𝜕𝐻 𝒓, 𝑡, λ! 𝜕λ

= !!

𝜕𝑉!" 𝒓, 𝑡, λ! 𝜕λ +

+ !!

𝑑𝑉! (𝒓, 𝑡, λ! ) 𝜕λ

where we use the notation

!! to

!!

𝜕𝑉!"# 𝒓, 𝑡, λ! 𝜕λ

+

𝑑𝑉! (𝒓, 𝑡, λ! ) 𝜕λ

+ !!

!!

+

𝑑𝐾(𝒑, 𝒓, 𝑡, λ! ) 𝜕λ

𝑑𝑉! (𝒓, 𝑡, λ! ) 𝜕λ

eq. 30 !!

!!

indicate that the preceding quantity was obtained from a

simulation at λS. Hence, the estimate 𝜕𝐻 𝒓, 𝑡, λ!

𝜕λ

!!

has been calculated from a

configurational ensemble obtained at λS and should therefore be reweighted to get the ensemble average appropriate for λP:18

𝜕𝐻 𝒓, 𝑡, λ! 𝜕λ

𝜕𝐻 𝒓, 𝑡, λ! 𝜕λ = !! !!

𝑒! !

!!

𝑒! !

!,!,!! !! !,!,!!

eq. 31

!! ! !!

!,!,!! !! !,!,!!

!! !

!!

The focus in this paper will be on the prediction of at λ values not equal to λS, without different parameterizations of the individual µμ values. From the pre-‐computed quantities in eqs. 10, 13, 17 and 19–26 for 𝜇!! = 0, 0.01, 0.02, … , 1 we can thus predict for 100 λP values ranging from 0 to 1 (see Figure 1). The explicit calculation of both the Hamiltonians and its λ-‐derivative, including a softcore potential, is the main difference compared to previously described reweighting schemes.13,15 Both these approaches do not include softcore potentials and


10

Page 11 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Li and Yang estimate the energy difference in the exponential of eq. 31 through the derivative at λS, which is not very accurate for non-‐linear λ-‐dependencies.., The reweighting in eq. 31 is only reliable if the difference between λS and λP is not too large, such that there is sufficient overlap in the phase space.

Combining predictions As pointed out in the previous section, it is in principle possible to predict for all λ values between 0 and 1 from a single simulation. However, the reliability of this prediction is low for large differences between λS and λP because the phase space overlap is likely to be insufficient. The obvious solution for this is to simulate at multiple λ values, such that the phase space overlap between neighboring states is improved. Each simulation can still predict for the full range of λ-‐points or for a subset just up to the neighboring λ-‐points. The predictions in the region λS1 < λP < λS2 will be calculated based on the predictions from the neighboring simulations λS1 and λS2 through a linear weighting scheme: 𝑤!! =

λ! − λ!! λ!! − λ!!

𝑤!! =

λ! − λ!! λ!! − λ!!

eq. 32

where wS1 and wS2 are the weighting factors for the predictions from the simulations at λS1 and λS2, respectively, 𝜕𝐻 𝒓, 𝑡, λ! 𝜕λ

= 𝑤!! !!

𝜕𝐻 𝒓, 𝑡, λ! 𝜕λ

+ 𝑤!! !! (!!! )

𝜕𝐻 𝒓, 𝑡, λ! 𝜕λ

eq. 33

!! (!!! )

Choice of λS values

One of the disadvantages of TI, as already pointed out in the introduction, is that the shape of the curve is not known a priori. This makes the initial choice of λ values to simulate rather arbitrary and can lead to unnecessary simulations in regions of low curvature. On the other hand, in regions where the curvature is high, multiple simulations have to be performed in a very narrow region of λ-‐points in order to capture this shape with the numerical integration. Typically one starts TI with 11 or 21 equidistant λ-‐points and adds points where needed. With extended TI, one can either use a similar scheme as for regular TI (i.e. equidistant points), or determine which λ-‐point would be the best to simulate next based on the predictions from the previously simulated λ-‐points. In both approaches, fewer λ-‐points will be required to reach


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 27

accurate results, because the predictions at intermediate λ-‐points smooth the curve for numerical integration. In an attempt to optimize the choice of λ-‐points for the most accurate free energy differences, we suggest the following scheme: 1.

Simulate at λ=0 and λ=1

2.

From each simulation, predict for all λ-‐points

3.

For every value of λP calculate the absolute difference between the predictions from each of the neighboring λS

4.

Weight the differences by 𝑤 = λ!! − λ!! − λ!! − λ! − λ! − λ!!

5.

Find the absolute maximum of the weighted difference. The corresponding λP-‐point is the next λ-‐point to simulate

6.

Repeat from step 2

The idea behind this approach is to introduce an extra λS where the prediction is the least reliable, i.e. where the neighboring predictions disagree most. The weighting step (4) lowers the significance of differences that are very close to one of the simulated states, because the prediction from the closest simulated state is very reliable and we thus do not need an extra simulation at this point. Other selection schemes could be based on the variance of , where λ-‐points would be added ‘nearby’ the highest variance λ-‐point.19 This requires, however, multiple initial simulations and the actual choice of a ‘nearby’ λ-‐point. The advantage of our scheme is that the curvature at nonsimulated λ-‐points is captured very quickly and additional λ-‐ points can be added accordingly. Whenever non-‐equidistant λ-‐points are mentioned below, these were determined by the above-‐ mentioned scheme, unless explicitly stated differently.

Methods Three different alchemical free energy transformations are performed to test and demonstrate the use of extended TI; (1) Methanol to dummy atoms, (2) phenol to dummy atoms and (3) p-‐ nitrobenzamidine to p-‐methoxybenzamidine (see Figure 2). In addition to these transformations in water, the benzamidine transformation is also performed with the molecule bound to trypsin. All molecular dynamics simulations in this study have been performed with a modified version of GROMOS1120 in combination with the 54a721,22 (methanol and phenol) or the 53a623


12

Page 13 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(benzamidine) force field parameter set. We first describe the simulation settings used for the model systems in water. Cubic simulation boxes were used, together with 800-‐1500 SPC water molecules for the simulations in solution.24 A chloride ion is included in the benzamidine simulations to neutralize the system. The systems are weakly coupled25 to an external bath with a temperature of 300 K (298 K for benzamidine) and a coupling time of 0.1 ps, with separate baths for the solute and solvent degrees of freedom. The pressure is kept constant at 1 atm by isotropic weak coupling25 with a relaxation time of 0.5 ps and a compressibility of 7.624 x 10-‐4 (kJ mol-‐1 nm-‐3)-‐1 for methanol and phenol, and of 4.575 x 10-‐4 (kJ mol-‐1 nm-‐3)-‐1 for benzamidine. The center of mass movement of the solute is removed every 1000 steps. A triple range cutoff scheme is used to calculate the longe range interactions, for which a pairlist is generated every 5th time step. Interactions within 0.8 nm are calculated at every time step, according to the pairlist. The interactions between 0.8 and 1.4 nm are calculated only together with the pairlist update and are kept constant in between. A reaction field contribution26 is added to the electrostatic interactions and forces to account for a homogeneous medium beyond 1.4 nm, using a relative permittivity of 61 (78.5 for phenol), see eq. 5 and eq. 7. SHAKE27 is used to constrain all bond lengths with a geometric accuracy of 1 x 10-‐4. The only exceptions are the two bonds were the bond types are being modified in the benzamidine simulations. These bond-‐stretching interactions are evaluated trough the quartic expression of eq. 20. The softness parameters for atoms that are being perturbed are set to αLJ = 0.5 and αCRF = 0.5 nm2 and the power dependence of the λ coupling (n in eq. 4) is equal to 2 for the methanol simulations and 1 for the other systems. All simulations are performed for 1 ns per λ-‐point. As the first test case, the methanol simulations are performed at 101 equidistantly spaced λ-‐points, to allow direct comparison of the predicted values from extended TI. The other two systems are simulated at 21 equidistantly spaced λ-‐points. For the methanol system, the regular TI data could directly be used as reference data, whereas for the latter two systems, the regular TI curve is smoothened by extended TI based on these 21 λS-‐ points. The phenol simulations are repeated five times with different initial velocities to obtain statistically independent runs and determine the precision of the method. The simulations of the trypsin-‐benzamidine complex have been set up as described in our previous paper.28 Initial simulations showed that benzamidine was moving out of the binding site at certain λ-‐points. Therefore, an additional distance restraint was introduced between the CG


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 27

atom of ASP171 and the CE atom of the benzamidine with a force constant of 500 kJ mol-‐1 nm-‐2. Similar to the benzamidine in water, the trypsin-‐benzamidine complex is simulated at 21 equidistantly spaced λ-‐points. However, these simulations are performed for 3 ns each, to allow more extensive sampling. The quantities described in eqs. 10, 13, 17, 19 and 20-‐23 are calculated for 𝜇!! = {0, 0.01, 0.02, …, 1} every 20 time steps and at the same time written to disk. In contrast, since the quantities described in eqs. 25 and 26 are independent of µμX, only a single value needs to be calculated and written to disk at every 20th time step. After post-‐processing the data, including the reweighting steps as described in eqs. 31-‐33 the trapezoidal rule is used as numerical integration algorithm. For a fair comparison of the BAR and (extended) TI, all methods were performed based on all data points (including correlated data points) and bootstrapping based on 100 bootstrap samples was used as uniform error estimate.

Results Panel A in Figure 3 shows the results of extended TI for the methanol test case. The reference TI data is indicated by the black curve and integration results in a free energy difference of 23.8 kJ/mol. The two regions of high curvature make it necessary to simulate at many λ-‐points when using regular TI. The colored curves in Figure 3A show the extended TI curves, based on predictions from 2 (red), 3 (green), 5 (blue), 6 (orange) or 11 (cyan) equidistant values of λS. When using extended TI, the information of only the two endpoints is already enough to give a rough estimate of the shape of the curve as well as the location of the maximum and minimum (red curve). Furthermore, with as few as 6 equidistant points the predicted curve (orange) can barely be distinguished from the reference TI curve. This also becomes clear from the differences in the ΔG values with respect to the reference data, as shown in Table 1. For comparison these differences are also shown for regular TI and BAR with n λ-‐points. Extended TI with 6 equidistant points only differs by 0.8 kJ/mol from the reference, whereas this is still -‐5.9 kJ/mol for regular TI. The deviations of the BAR results from the reference data are very small even from just 2 λ-‐ points. However, for BAR to obtain reliable and precise free energy differences, sufficient overlap between the forward and backward energy differences from neighboring λ-‐points should be present. It has therefore been suggested that the overlap integral (OI) between each of the


14

Page 15 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


neighboring λ-‐points should at least be 0.01.6,11 Analysis of the OIs shows that BAR results are reliable only with 5 or more λ-‐points. Close agreement between free-‐energy differences may still be a fortuitous result of cancelation of errors. Therefore, Table 1 also shows the mean absolute error (MAE), which is used to compare each point of the predicted curve with the reference curve and thereby prevents cancelation of errors that can be present in the ΔG-‐ΔGREF measure. For the calculation of the MAE of a regular TI, which has fewer data points than the reference data, a linear interpolation is used for the intermediate λ-‐points, as appropriate for integration with the trapezoidal rule. The MAE for extended TI with 6 equidistant points is 1.5 kJ/mol and for TI with the same λ-‐points 16.8 kJ/mol. For 11 points, the deviation of the regular TI curve is still 7.0 kJ/mol. Using non-‐equidistant points as discussed in the theory section improves the performance of extended TI even more. Here, 5 points lead to ΔG-‐ΔGREF of 1.4 kJ/mol and a MAE of 1.5 kJ/mol, whereas for 5 equidistant λS-‐points these values were still 4.0 and 4.1 kJ/mol, respectively. The phenol test system is expected to be more difficult because 13 atoms are being perturbed into dummy atoms, instead of only 3 atoms for methanol. As Figure 3B and Table 2 show, this is especially an issue when only a few λ-‐points are being used. With just two or three λ-‐points, none of the methods result in an accurate free energy difference. BAR seems to have a small deviation from the reference free energy for 5 equidistant λ-‐points (1.6 kJ/mol), but in fact the overlap integral criterium is only met when 6 λ-‐points are used, which then results in a deviation of 1.5 kJ/mol. For extended TI with equidistant points, the predicted curve looks reasonably similar to the reference one between λ = 0.75 and 1, but is far off for all other λ-‐points, especially around λ=0.3 (see Figure 3B). This is predominantly caused by the prediction from the simulation at λ=1. This prediction is extremely difficult, because we try to calculate the (free) energies for the situation where phenol has real atoms instead of dummy atoms, even though that space is still taken up by water molecules. The solvent and solute atoms would thus be clashing, causing very large energy estimates. With 3 equidistant λS-‐points, extended TI is already doing a lot better, but still has a local minimum around λ=0.6 which is not present in the reference TI curve. Starting from 5 equidistant points this additional local minimum is no longer predicted (MAE = 4.9 kJ/mol) and with 6 points the whole curve is represented very well with the MAE equal to 1.7 kJ/mol. Table S1 in the supporting information compares the root-‐mean-‐square deviation from


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 27

the reference data for different λ-‐intervals, using BAR, extended TI and TI. As expected, regular TI cannot handle λ-‐intervals larger than 0.1 reliably, while both BAR and extended TI show better agreement up to Δλ = 0.3 and 0.2, respectively. For some larger λ-‐intervals, no overlap was observed between the states, leading to large deviation for both BAR and extended TI. While BAR can handle larger λ-‐intervals, including the curvature at intermediate λ-‐values as in extended TI significantly enhances the efficiency of TI. Whereas for methanol non-‐equidistant λ-‐points led to more accurate results than equidistant points, the opposite is true in the phenol case. The problem is that the choice of the λ-‐points is based on the difference between the predictions from neighboring simulation points. Since the predicted values from λS = 0 and 1 are so different at low λP, additional λS-‐points are placed relatively close together at low values. For 5 non-‐equidistant points, the proposed simulated λS were 0, 0.25, 0.4, 0.55 and 1. The region around the minimum in the reference curve can only be predicted accurately once a larger λS (0.8) was included as 6th point. The benzamidine case is more complicated since not only non-‐bonded but also covalent changes are being made. Furthermore, this example represents an alchemical change that does not involve the removal of many atoms, but rather a change of character (Figure 2). Similar to the phenol test case, using just 2 λ-‐points is not enough to obtain accurate results with any of the methods for benzamidine in water (see Figure 3C and Table 3). However, the extended TI result is not so far off as in the phenol case. The general shape of the curve is now similar to the reference data, although the values around the minimum around λ=0.3 and the maximum around λ=0.75 are off by 20 and 10 kJ mol-‐1, respectively. Using 3 equidistant λ-‐points results in an accurate free energy difference (error of 0.4 kJ/mol) although there is some cancellation of errors as shown by the MAE of 1.7 kJ/mol. From 5 points onwards also the MAE is low, both using equidistant as well as non-‐equidistant λ-‐points. The overlap integral criterium for BAR is reached with 5 λ-‐points which results in an error of 0.6 kJ/mol with respect to the reference data. An additional challenge is to perform the same transformation of benzamidine when it is bound to the protein trypsin. Similar to the results for benzamidine in water, simulations at two λ-‐ points are not sufficient to obtain accurate results, but the general shape of the curve is captured quite well (Figure 3D). Using data from three simulated λ-‐points greatly improve the results for the region up to λ=0.5. In contrast to the benzamidine in water, the predicted curve does not


16

Page 17 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


match the reference data so well for λ > 0.5. This also becomes clear when comparing the MAEs for the simulations in water and when bound to trypsin, as shown in Table 3 and Table 4, respectively. In water, the MAE was 1.7 kJ/mol for 3 λ-‐points whereas this is 4.1 kJ/mol for the benzamidine bound to trypsin. For BAR, the overlap integral criterium is met with at least 5 λ-‐ points, leading to an error of -‐1.3 kJ/mol. As expected, the transformations of a molecule bound to a protein are more difficult to predict from fewer λ-‐points than the same transformations in solution. However, with as few as 6 equidistant λ-‐points the MAE is down to 1.1 kJ/mol. The scheme to automatically determine the next λ-‐point to simulate actually resulted in equidistant points up to 5 λ-‐points. Using 6 non-‐equidistant λ-‐points gives a slightly higher MAE (2.0 kJ/mol) than the equidistant scheme (1.1 kJ/mol). Overall, extended TI is giving very good results using only 6 equidistant λ-‐points (for benzamdine in water even with 5 λ-‐points), whereas regular TI usually requires around 21 λ-‐ points to obtain the same level of accuracy. Hereby, extended TI comes very close to the possibilities of BAR which also required 5 or 6 λ-‐points to get reliable results for our test systems. While we do not claim that extended TI is more accurate than BAR, one of the advantages of extended TI over BAR is that no additional post-‐analysis on coordinate trajectories is required when simulations at intermediate λ-‐points are needed. A further advantage of the current implementation of extended TI in the GROMOS code is moreover, the possibility to write out additional energy and free energy terms which not only allows the performance of extended TI but also much faster post-‐analysis for application of BAR. So far we have only discussed the accuracy of extended TI, but the precision is important just as well. To gain more insight in the latter, we have repeated the phenol simulations 5 times with different initial velocities. Additionally, for each of these 5 simulations we test the robustness of extended TI by reducing the amount of data included in the analysis. We use the full 1 ns, just the first 500 ps, or data collected at 80 fs (every second stored energy) and 200 fs (every fifth stored energy) intervals. The MAEs averaged over the 5 repeats as well as the standard deviations over the repeats are shown in Table 5 and Table 6. As was shown above for the single simulation of phenol using just the endpoints leads to very large errors, even when using all data (1 ns and


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 27

using every frame). The average MAE over the 5 repeats is around 3 x 102 kJ/mol with a standard deviation of 2 x 102 kJ/mol. Including λ=0.5 greatly improves the average MAE and also the standard deviation is much lower with a value of 2.9 kJ/mol (see Table 5). Using only the first half of the simulation causes an increase of the average MAE of 1 x 102 kJ/mol for 2 λ-‐points. For 3 equidistant points this is difference is smaller but still around 4 kJ/mol. In both cases the standard deviation was also slightly larger for the predictions based on 500 ps of simulations. From 5 equidistant points onwards the differences between the predictions based on half or the full data are not so pronounced any more. Using every second frame from the simulation data results in the same number of frames as using the first 500 ps, but the results are rather different. In contrast to the results from 500 ps, the average MAEs for every second frame are hardly different from the results of the full data. Using every fifth frame, which corresponds to using data from every 100th time-‐step or every 0.2 ps, again causes higher average MAEs when only few λ-‐points are used, although the standard deviation is lower. With increasing number of λ-‐ points, the results become similar to the full data results. Overall, if 5 or more equidistant λ-‐ points are used the standard deviation is always ≤ 1 kJ/mol, independent of which frames of the simulations were used. Whether results will be influenced by reducing the amount of data included in the analysis will depend on the correlation times of the quantities in question. In the case of TI, practically all correlation times of ∂H/∂λ are larger than 0.2 ps (which equals every 5th frame). This indicates that the TI results will not change much upon including only every 5th frame in the analysis (Table 5). The correlation time of extended TI is more complicated to determine because of the reweighting in eq. 31 and because from every simulation, 100 predictions are made which all have their own correlation times. For simplicity, we only evaluate the correlation time of the numerator in eq. 31, which should be lower than that of the denominator because of the multiplication with ∂H/∂λ. The vast majority of the evaluated correlation times are smaller than 0.2 ps and it can be observed that larger differences between λP and λS tend to give smaller correlation times. This is also reflected in the results of Table 5, where the differences in MAE upon reducing the data to every 5th frame are larger when fewer λ-‐points are used. Interestingly, the accuracy of TI improves more from including more λ-‐points than from longer simulations, whereas this is not necessarily the case for extended TI. Considering a constant total


18

Page 19 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulation time of around 11 ns for each of the five repeats of the phenol simulations, either 11 points of each 1 ns or 21 λ-‐points of each 0.5 ns can be used. TI improves from an average MAE of 9.7 to 2.5 kJ/mol when including more λ-‐points, whereas the results for extended TI hardly change with average MAEs of 1.0 and 0.8 kJ/mol, respectively (not all data shown). As discussed above, the scheme to determine the non-‐equidistant λ-‐points to simulate (see theory section) can give non-‐optimal results when the perturbation is large. This is also shown in Table 6 for the 5 repeats of the phenol system. Using three non-‐equidistant λ-‐points gives much larger average MAEs than the equidistant points and in particular the standard deviations are also very large. Small differences in the simulated end points can determine a different λ-‐point to simulate next and this is what causes the large differences in the MAEs. If the same (albeit non-‐ equidistant) points would be used among the repeats, the standard deviation would be much smaller (data not shown). Similar to the equidistant points, there is not much difference between using the full data or every second frame. The first half of the simulation gives higher average MAEs than the full simulation, which is mostly caused by a single repeat with a very large MAE. Using every fifth frame gives even higher average MAEs, mostly because there are multiple repeats with large errors. This is predominantly caused by the choice of λS-‐points. For example, using 6 non-‐equidistant points the average MAE equals 13 ± 9 kJ/mol, whereas this would be only 5 ± 4 kJ/mol if the same λS-‐points were used as determined from all data (data not shown).

Discussion The results in Table 6 showed that the scheme to determine the next λ point to simulate is not very reliable for large perturbations. However, the scheme can be modified such that at step 1, we simulate e.g. 3 equidistant points instead of only the end points and then continue at step 2 as indicated. This greatly improves the results for the phenol simulations as shown in Table S2. In principle one can start with any choice of λ-‐points, equidistant or not, and then use the scheme to find the λ points that would be the best to add to further refine the curve. This is probably one of the most efficient strategies, because to determine the next λ-‐point to simulate with the above-‐ mentioned scheme, the simulation at the previous λ-‐point must have finished. On the other hand, when equidistant points are used (or any other set of predefined λ-‐points), the simulations can be performed in parallel.


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 27

The results in Table 5 and 6 show that (at least for this system) it might not be necessary to write out the energies and free energies every 20 time steps, but this can be reduced to at least every 40 time steps. This would reduce both the run time of the simulations as well as the required disk space. Additional reduction of resources can be achieved when the energy contributions are not calculated for the whole range of λ-‐points. Especially when using equidistant λS, this can be easily achieved because the predictions only need to be made up to the neighboring λS-‐points. For example, if 6 equidistant λ-‐points are simulated, then for λS = 0.2 predictions are only required for λP ranging between 0.0 and 0.4. The additional simulation time due to the calculations for extended TI is determined by a number of factors. These include how often the additional calculations are performed (every X timesteps), the number of λP for which the interactions are determined, the number of perturbed atoms and the total size of the system. We used the benzamidine-‐in-‐water system as a test case, where the additional calculations are performed every 20 time steps for 101 λP-‐points and four atoms, two bond lengths, one angle and one dihedral angle are perturbed. The total number of atoms in the system was 2467 and on a single cpu, extended TI takes around 10% more simulation time than a regular TI simulation. For larger systems this percentage would be significantly lower, because the percentage of interactions involving perturbed atoms is much smaller than in a small system. Extended TI thus results in a 69% gain of efficiency since only around 6 λ-‐points need to be simulated instead of the most commonly used 21 λ-‐points with an additional cost of only 10% or less per λ-‐point. The above-‐described application of predicting the TI curve from fewer λ-‐points is not the only possible application of the current code. The energy terms of eqs. 10, 13, 17 and 19-‐26 are written out separately because this leaves much more flexibility in terms of post-‐analysis. Up to now we have only predicted for other λP-‐points than the λS-‐points at which simulations were performed with, but in principle predictions can be made for changes in any of the parameters that influence the free energy pathway, as long as the intrinsic outcome of the calculation does not change. Systematically searching this parameter space can lead to insight in the most efficient pathway and thereby further improve the efficiency of TI. Naden et al. have also explored the parameter space to find the optimal path, by making use of linear combinations of λ-‐independent


20

Page 21 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


basis functions describing the non-‐bonded interactions.19,29 However, their approach circumvents the optimization of softcore parameters. In our current approach, parameter changes can e.g. involve separate changes of parameterization of the individual µμx values (including softness parameters) according to eq. 3 or of the power dependency of µμ, n. As these parameters can have a strong influence on the shape of the TI curve,9,29 they directly influence the efficiency of the TI because high curvature regions require closely spaced λ-‐points in order to integrate accurately. When using BAR one could explore the influence of simulation parameters on efficiency by using non-‐Boltzmann Bennett, but this would require the expensive re-‐ evaluation of the coordinate trajectories.30 The current implementation of extended TI allows the reconstruction of the total energies for each of the parameter sets from the pre-‐computed values from the simulations of a single parameter set and this should therefore be extremely fast. These possibilities will be explored in future work.

Supporting information The Supporting Information is available free of charge on the ACS Publications website. -‐ Root mean square deviations over all possible λ-‐intervals comparing BAR, extended TI and TI for the phenol system -‐ Average MAE values for different initial λ-‐points for the non-‐equidistant λ-‐points scheme.

Acknowledgements The authors thank Drazen Petrov for the initial setup of the phenol system. Financial support by the Vienna Science and Technology Fund (WWTF; grant number LS08-‐QM03) is gratefully acknowledged.


21


Figures

λ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 27

1

0 λ

Figure 1. Schematic overview of extended TI. The red dots represent λ-‐values that are explicitly simulated and for which is obtained. The resulting TI curve is indicated by the red line. By predicting 𝜕𝐻/𝜕λ !! at intermediate, non-‐simulated λ-‐values, the black curve may be obtained.

Figure 2. The three model systems used in this paper; methanol to dummy atoms (top), phenol to dummy atoms (middle) and p-‐nitrobenzamidine to p-‐methoxybenzamidine (bottom). D indicates a non-‐interacting dummy atom.


22

Page 23 of 27

A 200

100

0

B 200 100

[kJ/mol]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 -100 C

-100

-120

-140

-160 D

-100

-120

-140

-160 0

0.2

0.4

0.6

0.8

1

λ

Figure 3. Results for extended TI using equidistant points for methanol (A), phenol (B), benzamidine in water (C) and benzamidine bound to trypsin (D). Different colors show the number of λ points used for extended TI: red: 2 points, green: 3 points, blue: 5 points, orange: 6 points and cyan: 11 points. The thick black curve with error bars shows the data from regular TI.


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 27

Table 1. Deviations from the reference data using different numbers of simulations for methanol (reference: ΔGREF = 23.8 ± 0.0 kJ/mol). The error estimates given on the free energy differences with respect to the reference data are based on bootstrapping analysis and do not include the error estimate of the reference data. MAE extTI TI extTI extTI (non.eq) (eq) (non.eq) 2 59.0 ± 0.1 -‐1.3 ± 2.8 9.1 ± 0.3 9.1 ± 0.3 70.5 10.2 10.2 3 13.8 ± 0.1 -‐0.2 ± 0.3 8.6 ± 0.2 3.9 ± 0.2 28.5 9.2 5.1 5 -‐7.9 ± 0.1 -‐0.6 ± 0.2 4.0 ± 0.3 1.4 ± 0.1 15.9 4.1 1.5 6 -‐5.9 ± 0.1 -‐0.1 ± 0.1 0.8 ± 0.1 1.0 ± 0.1 16.8 1.5 1.1 11 -‐3.2 ± 0.1 0.1 ± 0.1 0.1 ± 0.0 0.6 ± 0.0 7.0 0.8 0.7 n refers to the number of simulated λ-‐points. MAE is the mean absolute error with respect to the reference TI curve. extTI (eq) and extTI (non.eq) refer to the results of extended TI with equidistant and non-‐ equidistant points, respectively. n

TI

ΔG-ΔGREF BAR extTI (eq)

Table 2. Deviations from the reference data using different numbers of simulations for phenol (reference: ΔGREF = 68.0 ± 0.1 kJ/mol). The error estimates given on the free energy differences with respect to the reference data are based on bootstrapping analysis and do not include the error estimate of the reference data. MAE extTI extTI (eq) (non.eq) 2 2 22.5 ± 0.1 (-‐5 ± 1) x 10 -‐243 ± 8 -‐243 ± 8 53.8 251.4 251.4 3 11.3 ± 0.1 -‐24.3 ± 3.8 3.1 ± 0.4 -‐72.1 ± 2.0 47.2 23.2 86.3 5 -‐2.5 ± 0.1 1.6 ± 0.2 3.9 ± 0.2 5.5 ± 0.3 34.7 4.9 15.9 6 -‐12.0 ± 0.1 1.5 ± 0.1 1.6 ± 0.2 1.7 ± 0.1 25.1 1.7 2.7 11 -‐2.7 ± 0.1 0.5 ± 0.1 0.4 ± 0.0 0.5 ± 0.0 9.8 0.5 0.5 n refers to the number of simulated λ-‐points. MAE is the mean absolute error with respect to the reference TI curve. extTI (eq) and extTI (non.eq) refer to the results of extended TI with equidistant and non-‐ equidistant points, respectively. n

TI

ΔG-ΔGREF BAR

extTI (eq)

extTI (non.eq)

TI

Table 3. Deviations from the reference data using different numbers of simulations for benzamidine in water (reference: ΔGREF = -‐131.8 ± 0.0 kJ/mol). The error estimates given on the free energy differences with respect to the reference data are based on bootstrapping analysis and do not include the error estimate of the reference data. MAE extTI TI extTI extTI (non.eq) (eq) (non.eq) 2 24.3 ± 0.1 2.2 ± 3.6 -‐1.7 ± 0.5 -‐1.7 ± 0.5 25.5 9.9 9.9 3 8.6 ± 0.2 1.1 ± 0.4 0.4 ± 0.2 0.7 ± 0.2 12.7 1.7 2.6 5 1.4 ± 0.1 0.6 ± 0.1 0.0 ± 0.1 -‐0.7 ± 0.1 7.8 0.6 1.3 6 0.8 ± 0.1 0.1 ± 0.1 -‐0.6 ± 0.0 -‐0.4 ± 0.0 7.0 0.7 0.9 11 0.4 ± 0.1 -‐0.1 ± 0.1 0.0 ± 0.0 -‐0.1 ± 0.0 2.1 0.2 0.3 n refers to the number of simulated λ-‐points. MAE is the mean absolute error with respect to the reference TI curve. extTI (eq) and extTI (non.eq) refer to the results of extended TI with equidistant and non-‐ equidistant points, respectively. n

TI

ΔG-ΔGREF BAR extTI (eq)


24

Page 25 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 4. Deviations from the reference data using different numbers of simulations for benzamidine when bound to trypsin (reference: ΔGREF = -‐124.1 ± 0.0 kJ/mol). The error estimates given on the free energy differences with respect to the reference data are based on bootstrapping analysis and do not include the error estimate of the reference data. MAE extTI extTI (eq) (non.eq) 2 14.8 ± 0.1 -‐7.9 ± 1.8 -‐8.5 ± 0.3 -‐8.5 ± 0.3 15.8 9.0 9.0 3 4.3 ± 0.1 -‐3.5 ± 0.4 -‐3.4 ± 0.1 -‐3.4 ± 0.1 12.8 4.1 4.1 5 -‐0.5 ± 0.1 -‐1.3 ± 0.1 -‐1.9 ± 0.1 -‐1.9 ± 0.1 6.0 2.0 2.0 6 1.3 ± 0.0 0.3 ± 0.0 0.2 ± 0.0 -‐2.0 ± 0.0 5.7 1.1 2.0 11 0.5 ± 0.0 0.1 ± 0.0 0.1 ± 0.0 -‐0.4 ± 0.0 2.0 1.0 0.8 n refers to the number of simulated λ-‐points. MAE is the mean absolute error with respect to the reference TI curve. extTI (eq) and extTI (non.eq) refer to the results of extended TI with equidistant and non-‐ equidistant points, respectively. n

TI

ΔG-ΔGREF BAR

extTI (eq)

extTI (non.eq)

TI

Table 5. Average and standard deviation of the MAE over 5 repeats of the phenol simulations, while using all or reduced amount of data. 2, 3, 5,6 and 11 equidistant spaced λ-‐points are used. 1.0 ns 0.5 ns nd every 2 th every 5

2 2 (3 ± 2) x 10 2 (4 ± 2) x 10 2 (3 ± 2) x 10 2 (6 ± 1) x 10

3 22 ± 3 26 ± 4 23 ± 3 24 ± 2

5 4.5 ± 0.7 4.9 ± 1.0 4.6 ± 0.8 5.0 ± 0.9

6 2.1 ± 0.6 2.4 ± 0.6 2.4 ± 0.7 2.6 ± 0.4

11 1.0 ± 0.3 1.2 ± 0.5 1.1 ± 0.3 1.1 ± 0.3

Table 6. Average and standard deviation of the MAE over 5 repeats of the phenol simulations, while using all or reduced amount of data. 2, 3, 5,6 and 11 non-‐equidistant spaced λ-‐points are used. 1.0 ns 0.5 ns nd every 2 th every 5

2 2 (3 ± 2) x 10 2 (4 ± 2) x 10 2 (3 ± 2) x 10 2 (6 ± 1) x 10

3 1 (10 ± 6) x 10 1 (12 ± 7) x 10 1 (10 ± 6) x 10 1 (17 ± 7) x 10

5 17 ± 10 17 ± 2 14 ± 7 27 ± 13

6 5 ± 4 6 ± 4 5 ± 4 13 ± 9


11 0.7 ± 0.3 1.2 ± 0.3 0.8 ± 0.3 0.7 ± 0.4

25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 27

References (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30)

Shirts, M. R.; Mobley, D. L.; Chodera, J. D. In Annual Reports in Computational Chemistry; D.C. Spellmeyer and R. Wheeler, Ed.; Elsevier, 2007; Vol. Volume 3, pp 41–59. Shirts, M. R.; Mobley, D. L.; Brown, S. P. In Drug Design: Structure-‐ and Ligand-‐Based Approaches; Merz, K. M., Ringe, D., Reynolds, C. H., Eds.; Cambridge University Press: New York, NY, 2010; pp 61–86. Mobley, D. L.; Klimovich, P. V. J. Chem. Phys. 2012, 137 (23), 230901. Christ, C. D.; Fox, T. J. Chem. Inf. Model. 2014, 54 (1), 108–120. Kirkwood, J. G. J. Chem. Phys. 1935, 3 (5), 300. Bennett, C. H. J. Comput. Phys. 1976, 22 (2), 245–268. Shirts, M. R.; Pande, V. S. J. Chem. Phys. 2005, 122 (14), 144107–144116. Pham, T. T.; Shirts, M. R. J. Chem. Phys. 2011, 135 (3), 034114. de Ruiter, A.; Boresch, S.; Oostenbrink, C. J. Comput. Chem. 2013, 34 (12), 1024–1034. Shirts, M. R.; Chodera, J. D. J. Chem. Phys. 2008, 129 (12), 124105–124110. Bruckner, S.; Boresch, S. J. Comput. Chem. 2011, 32 (7), 1303–1319. Paliwal, H.; Shirts, M. R. J. Chem. Theory Comput. 2011, 7 (12), 4115–4134. Tidor, B.; Karplus, M. Biochemistry 1991, 30 (13), 3217–3228. Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Chem. Phys. Lett. 1994, 222 (6), 529–539. Li, H.; Yang, W. Chem. Phys. Lett. 2007, 440 (1–3), 155–159. Riniker, S.; Christ, C. D.; Hansen, H. S.; Hünenberger, P. H.; Oostenbrink, C.; Steiner, D.; van Gunsteren, W. F. J. Phys. Chem. B 2011, 115 (46), 13570–13577. Hritz, J.; Oostenbrink, C. J. Chem. Phys. 2008, 128 (14), 144121. Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23 (2), 187–199. Naden, L. N.; Shirts, M. R. J. Chem. Theory Comput. 2015, 11 (6), 2536–2549. Schmid, N.; Christ, C. D.; Christen, M.; Eichenberger, A. P.; van Gunsteren, W. F. Comput. Phys. Commun. 2012, 183 (4), 890–903. Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C. J. Chem. Theory Comput. 2012, 8 (10), 3705–3723. Reif, M. M.; Winger, M.; Oostenbrink, C. J. Chem. Theory Comput. 2013, 9 (2), 1247–1264. Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. J. Comput. Chem. 2004, 25 (13), 1656–1676. Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Reidel: Dordrecht, The Netherlands, 1981; Vol. 11, pp 331–342. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81 (8), 3684–3690. Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. J. Chem. Phys. 1995, 102 (13), 5451–5459. Ryckaert, J.-‐P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23 (3), 327–341. de Ruiter, A.; Oostenbrink, C. J. Chem. Theory Comput. 2012, 3686–3695. Naden, L. N.; Pham, T. T.; Shirts, M. R. J. Chem. Theory Comput. 2014, 10 (3), 1128–1149. König, G.; Boresch, S. J. Comput. Chem. 2011, 32 (6), 1082–1090.


26

Page 27 of 27

For Table of Contents Only

λ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

λ

1


27

Extended Thermodynamic Integration: Efficient Prediction of Lambda

Recommend Documents