Accelerated Enveloping Distribution Sampling: Enabling Sampling of

Subscriber access provided by UNIV OF DURHAM

B: Biophysical Chemistry and Biomolecules

Accelerated Enveloping Distribution Sampling: Enabling Sampling of Multiple End-States While Preserving Local Energy Minima Jan Walther Perthold, and Chris Oostenbrink J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b02725 • Publication Date (Web): 18 Apr 2018 Downloaded from http://pubs.acs.org on April 19, 2018

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 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 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.

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

The Journal of Physical Chemistry

Accelerated Enveloping Distribution Sampling: Enabling Sampling of Multiple End-States While Preserving Local Energy Minima Jan Walther Perthold and Chris Oostenbrink* Institute for Molecular Modeling and Simulation, Department for Material Sciences and Process Engineering, University of Natural Resources and Life Sciences (BOKU), Vienna, Muthgasse 18, 1190 Vienna, Austria Corresponding Author *[email protected]

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 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 Enveloping distribution sampling (EDS) is an efficient approach to calculate multiple free-energy differences from a single molecular dynamics (MD) simulation. However, the construction of an appropriate reference state Hamiltonian that samples all states efficiently is not straightforward. We propose a novel approach for the construction of the EDS reference state Hamiltonian, related to a previously described procedure to smoothen energy landscapes. In contrast to previously suggested EDS approaches, our reference-state Hamiltonian preserves local energy minima of the combined end-states. Moreover, we propose an intuitive, robust and efficient parameter optimization scheme to tune EDS Hamiltonian parameters. We demonstrate the proposed method with established and novel test systems and conclude that our approach allows for the automated calculation of multiple free-energy differences from a single simulation. Accelerated EDS promises to be a robust and user-friendly method to compute free-energy differences based on solid statistical mechanics.

KEYWORDS Molecular Dynamics Simulation, Free-Energy Calculation, Free-Energy Perturbation, Alchemical Change, GROMOS


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


Introduction In computational drug design, molecular dynamics (MD) simulations are increasingly used to calculate precise binding free energy estimates for novel drug candidates and aid the development thereof.1 With alchemical methods it is possible to calculate free-energy differences between different chemical compounds in solution and when bound to a receptor. Using thermodynamic cycles, those free-energy differences can be translated into relative binding free energies.2 A plethora of methods is available to perform alchemical processes in MD simulations.3 Notably, the free energy perturbation (FEP) approach proposed by Zwanzig4 more than six decades ago is still heavily being used today. The free energy difference ∆𝐹 between two states 1 and 2, defined by their Hamiltonians 𝐻1 and 𝐻2 , can be obtained from an ensemble average over the positions 𝑟⃑ for one of the states, here over state 1 (𝑅 is the gas constant and 𝑇 is the absolute temperature): ∆𝐹2,1 = −𝑅𝑇 𝑙𝑛⟨𝑒 −(𝐻2 (𝑟⃑)−𝐻1 (𝑟⃑))⁄𝑅𝑇 ⟩1

(1)

Clearly, FEP can be extended to multiple Hamiltonians that were not simulated (like 𝐻2 ). However, the perturbation formula only gives meaningful results if the sampled phase-space overlaps significantly with the phase-spaces of the non-simulated Hamiltonians, or in the limit of infinite sampling. A possible answer to this problem is to split the perturbation process into multiple steps, where combinations of the Hamiltonians are used to sample intermediate steps. This, however, diminishes the potential for multi-state perturbations and hence reduces efficiency. A different approach is the creation of a reference state 𝑅, which ideally samples the phase-space of all desired end-states and thereby allows for efficient multi-state perturbations. One such approach is the one-step perturbation (OSP) method. Here, a reference Hamiltonian 𝐻𝑅 (𝑟⃑) is constructed prior to the simulation. Different reference Hamiltonians have been proposed and are still actively being discussed in literature.5-15 Depending on the molecules of interest, different


3


Page 4 of 27

reference Hamiltonians might need to be chosen which is not always straightforward. An elegant approach to create a reference Hamiltonian is enveloping distribution sampling (EDS).16 Here, 𝑛 different Boltzmann-weighted end-state Hamiltonians 𝐻𝑖 (𝑟⃑) are combined in the reference Hamiltonian 𝐻𝑅 (𝑟⃑), which partition function is now the sum of the partition functions of the endstates. Additionally, to ensure equal sampling of the end-states, constant free-energy offset parameters ∆𝐹𝑖𝑅 are added:16-19 𝑅

𝐻𝑅 (𝑟⃑) = −𝑅𝑇 𝑙𝑛 (∑𝑛𝑖=1 𝑒 −(𝐻𝑖 (𝑟⃑)−∆𝐹𝑖

)⁄𝑅𝑇

(2)

)

Here, we call the offset parameters ∆𝐹𝑖𝑅 free-energy offset parameters, because they correspond to the relative free-energies of the end-states if all states are sampled with equal probability. However, Equation (2) was found to result in high energy-barriers, preventing efficient sampling of the minima of all end-states. An additional smoothening parameter 𝑠 was proposed to smoothen the energy landscape of the reference Hamiltonian:17, 20-23 𝐻𝑅𝑠 (𝑟⃑) = −

𝑅𝑇 𝑠

𝑅

𝑙𝑛 (∑𝑛𝑖=1 𝑒 −𝑠(𝐻𝑖 (𝑟⃑)−∆𝐹𝑖

)⁄𝑅𝑇

(3)

)

It is noted that in the special case of equation (3), the offset parameters ∆𝐹𝑖𝑅 do not necessarily correspond to the relative free-energies of the end-states if all states are sampled with equal probability. Moreover, the 𝑠 and ∆𝐹𝑖𝑅 parameters strongly influence each other, complicating the search for their optimal values. Different EDS parameter search schemes and more elaborate combination rules have been proposed to optimize the reference Hamiltonian with parameters 𝑠 and ∆𝐹𝑖𝑅 .23-27 While EDS has been successfully applied for the study of conformational transitions in macromolecules28-31 and for constant pH MD simulations,32-34 the search for reference Hamiltonian parameters in multi-state EDS still represents a difficult problem. Very recently, a generalized replica-exchange MD (REMD) method was proposed to reduce the complexity of the parameter-search problem, albeit with the drawback of high computational cost.35,

36

By


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


introducing additional intermediate states in the EDS equation, combined with an automated update scheme for the offset parameters for each of these states, the construction of an energy landscape with reduced barriers could be simplified as well.37 Theory EDS reference Hamiltonians according to equation (3) with 𝑠 < 1, do not preserve the energy minima of the energy landscape as given in equation (2). Hence, sampling of important energy minima is hindered and inappropriate choices of 𝑠 render some states inaccessible in the reference Hamiltonian. These properties of the EDS reference Hamiltonian with a smoothening parameter 𝑠, and the difficulties in determining 𝑠 and ∆𝐹𝑖𝑅 represent significant drawbacks. Similarly, in the approach using intermediate states proposed by Mori, et al.,37 the reference potential is optimized towards a flat-bottom potential within the bounds of the original energy minima, largely increasing the sampling of less relevant phase-space on the integrated energy surface. Here, we propose an alternative construction of an EDS reference Hamiltonian in which all endstates are accessible and, importantly, in which the energy minima of the different end-states are preserved. Similar to the recently proposed Gaussian accelerated MD (GAMD) method,38 which is a variant of the accelerated MD approach,39 we smoothen the reference Hamiltonian 𝐻𝑅 (𝑟⃑) given in equation (2) by modification with a harmonic potential energy function. However, in contrast to the GAMD approach, we do not pull the energy landscape up to a certain chosen maximum energy parameter, but take the opposite approach; i.e. we pull the energy landscape down to a minimum energy parameter 𝐸𝑚𝑖𝑛 . Moreover, we restrict the modification of 𝐻𝑅 (𝑟⃑) to a region between 𝐸𝑚𝑖𝑛 and a maximum energy parameter 𝐸𝑚𝑎𝑥 . This preserves the energy landscape around local energy minima in greater detail and allows for intuitive tuning of the parameters 𝐸𝑚𝑖𝑛 and 𝐸𝑚𝑎𝑥 . We define the continuous accelerated EDS (A-EDS) reference Hamiltonian 𝐻𝑅∗ (𝑟⃑) as:


5


𝐻𝑅 (𝑟⃑) − 𝐻𝑅∗ (𝑟⃑) =

𝐻𝑅 (𝑟⃑) − 2 (𝐸

𝐸𝑚𝑎𝑥 −𝐸𝑚𝑖𝑛

𝑚𝑎𝑥 −𝐸𝑚𝑖𝑛 )

{

𝐻𝑅 (𝑟⃑) ≥ 𝐸𝑚𝑎𝑥

2

1

Page 6 of 27

(𝐻𝑅 (𝑟⃑) − 𝐸𝑚𝑖𝑛 )2

𝐸𝑚𝑖𝑛 < 𝐻𝑅 (𝑟⃑) < 𝐸𝑚𝑎𝑥

𝐻𝑅 (𝑟⃑)

(4)

𝐻𝑅 (𝑟⃑) ≤ 𝐸𝑚𝑖𝑛

It is possible to define the end-state 𝑖 which is currently being sampled in a simulation of the EDS reference state as the state with the lowest energy 𝐻𝑖 (𝑟⃑) − ∆𝐹𝑖𝑅 . We make use of this definition to propose a simple optimization scheme for the A-EDS parameters 𝐸𝑚𝑖𝑛 and 𝐸𝑚𝑎𝑥 . By keeping track of the unmodified average value of 𝐻𝑅 (𝑟⃑) in each state, 𝐸̅𝑖 , and the average of the maximum transition energy between states within a state round-trip (we define a round-trip as † having visited all states at least once), 𝐸̅𝑚𝑎𝑥 , a maximum energy barrier ∆𝐸𝑚𝑎𝑥 between the end-

states in the unmodified EDS reference Hamiltonian 𝐻𝑅 (𝑟⃑) can be defined: 𝐸̅𝑖 = 〈𝐻𝑅 (𝑟⃑)〉𝑅,𝑖 † 𝐸̅𝑚𝑎𝑥 = 〈max (𝐻𝑅† (𝑟⃑))〉𝑅,𝑟𝑜𝑢𝑛𝑑−𝑡𝑟𝑖𝑝𝑠 † ∆𝐸𝑚𝑎𝑥 = 𝐸̅𝑚𝑎𝑥 − min(𝐸̅𝑖 )

(5)

𝐻𝑅† (𝑟⃑) is the value of 𝐻𝑅 (𝑟⃑) when a new state is sampled. The unmodified EDS Hamiltonian 𝐻𝑅 (𝑟⃑) can now be accelerated by setting 𝐸𝑚𝑖𝑛 and 𝐸𝑚𝑎𝑥 such that the maximum energy barrier ∗ between the states, ∆𝐸𝑚𝑎𝑥 , is reduced to a user-defined parameter ∆𝐸𝑚𝑎𝑥 in 𝐻𝑅∗ (𝑟⃑). 𝐸𝑚𝑎𝑥 simply † corresponds to 𝐸̅𝑚𝑎𝑥 and 𝐸𝑚𝑖𝑛 is calculated as: † 𝐸𝑚𝑎𝑥 = 𝐸̅𝑚𝑎𝑥

𝐸𝑚𝑖𝑛 = {

𝐸𝑚𝑖𝑛 = 𝐸𝑚𝑎𝑥 ∗ ) − 𝐸𝑚𝑎𝑥 𝐸𝑚𝑖𝑛 = 2 (min(𝐸̅𝑖 ) + ∆𝐸𝑚𝑎𝑥

∗ ∆𝐸𝑚𝑎𝑥 ≥ ∆𝐸𝑚𝑎𝑥 ∗ ∆𝐸𝑚𝑎𝑥 < ∆𝐸𝑚𝑎𝑥

(6)

If the resulting value of 𝐸𝑚𝑖𝑛 is less than min(𝐸̅𝑖 ), it is adjusted to: 𝐸𝑚𝑖𝑛 =

1 1 2 ∗ 𝐸𝑚𝑎𝑥 ∆𝐸𝑚𝑎𝑥 +𝐸𝑚𝑎𝑥 min(𝐸̅𝑖 ) − 𝐸𝑚𝑎𝑥 2 − min(𝐸̅𝑖 )𝑖 2

2

∗ ∆𝐸𝑚𝑎𝑥

(7)


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


The proposed optimization scheme for the A-EDS parameters can be used in a non-equilibrium parameter search simulation in which the user only specifies the desired maximum barrier height, ∗ ∆𝐸𝑚𝑎𝑥 . No initial values for 𝐸𝑚𝑖𝑛 and 𝐸𝑚𝑎𝑥 have to be set, as they are adjusted while the reference

state explores the energy landscape. To enable sampling of all states in the beginning of the search simulation, the value of 𝐸𝑚𝑎𝑥 is set to the maximum value of the unmodified Hamiltonian 𝐻𝑅 (𝑟⃑) as the system explores the energy landscape until all states have been visited at least once. After all states have been visited at least once, 𝐸𝑚𝑎𝑥 is calculated as described above. In Figure 1, the EDS Hamiltonian of two harmonic potentials is shown. To a good approximation, 𝐻𝑅 (𝑟⃑) according to equation (2) is given by min(𝐻1 , 𝐻2 − ∆𝐹2𝑅 ). In panel A, the reference state Hamiltonian with a parameter 𝑠 according to equation (3) is shown. Note that both energy minima are shifted on the reaction coordinate and regarding their relative energy levels. Moreover, the choice of a lower 𝑠 parameter to reduce the barrier further would lead to the disappearance of the energy minimum on the left side. In panels B and C, the reference state Hamiltonian was constructed according to the proposed A-EDS formulation given in equation (4). This reference state Hamiltonian 𝐻𝑅∗ (𝑟⃑) preserves the energy minima with respect to their position on the reaction coordinate and their relative energy levels and the energy barrier between the minima is reduced to a small, desired value. The choice of ideal free-energy offset parameters ∆𝐹𝑖𝑅 for the A-EDS Hamiltonian is not influenced by the A-EDS parameters, which is not the case for the EDS Hamiltonian formulated with a parameter 𝑠.


7


Page 8 of 27

A

B

C

Figure 1. EDS reference state Hamiltonians (black solid lines) constructed either according to equation (3) with a parameter of 𝑠 = 0.012 (A) or according to the proposed A-EDS formulation given in equation (4) with parameters 𝐸𝑚𝑎𝑥 = 126 and 𝐸𝑚𝑖𝑛 = 10 (B) or 𝐸𝑚𝑎𝑥 = 126 and 𝐸𝑚𝑖𝑛 = −500 (C). The original harmonic potential energy functions are shown with grey solid and brown dotted lines, respectively. The Hamiltonian shown with brown solid lines was shifted using a freeenergy offset parameter of ∆𝐹2𝑅 = 40 for the construction of the EDS reference state Hamiltonians. Energy minima of the original potential energy functions and the constructed EDS Hamiltonians are shown with dashed lines. In panels B and C, A-EDS parameters 𝐸𝑚𝑎𝑥 and 𝐸𝑚𝑖𝑛 , ∗ and the target maximum energy barrier height ∆𝐸𝑚𝑎𝑥 are shown with red solid and dashed lines, respectively.


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


In addition to the A-EDS parameters 𝐸𝑚𝑖𝑛 and 𝐸𝑚𝑎𝑥 , the free energy offset parameters ∆𝐹𝑖𝑅 can be optimized simultaneously during the parameter search simulation. The free energy offset parameter of the first state is arbitrarily set to zero and the other free energy offset parameters are calculated relative to the first state by simply using the perturbation formula: ∆𝐹1𝑅 = 0 𝑅 ∆𝐹𝑖≠1 = ∆𝐹𝑖,1 = ∆𝐹𝑖,𝑅 − ∆𝐹1,𝑅 =

= −𝑅𝑇 𝑙𝑛 ⟨𝑒 −(𝐻𝑖

(𝑟⃑)−𝐻𝑅∗ (𝑟⃑))⁄𝑅𝑇

⟩ + 𝑅𝑇 𝑙𝑛 ⟨𝑒 −(𝐻1 𝑅

(𝑟⃑)−𝐻𝑅∗ (𝑟⃑))⁄𝑅𝑇

(8)

⟩ 𝑅

As sufficient sampling of each state is required to calculate the free energy offset parameters, the offset parameters are only updated after the state has been sampled for a certain amount of simulation time 𝑡𝑚𝑖𝑛 . Before the free energy parameters are optimized, they are all set to zero or to the values of the relative energy differences between the states in the beginning of the simulation. Methods We implemented the A-EDS reference Hamiltonian 𝐻𝑅∗ (𝑟⃑) and the non-equilibrium parameter search in the GROMOS MD simulation engine.40 To illustrate the A-EDS method, different test systems which partially have already been used earlier14, 16, 23, 37 were simulated. The test systems were: 1. Water to methanol conversion (SPC-MeOH): a SPC water molecule41 and a methanol molecule42 were simulated in a cubic box (edge length 3.31 nm) with 1179 SPC water molecules as solvent. In each state, non-bonded interactions of one molecule were turned on and the other molecule was represented by non-interacting dummy atoms.16, 23 2. Full dipole inversion (FDI): The charges of a GROMOS 53A643 ethane dipole with charges +1e and -1e solvated in a cubic box (edge length 3.42 nm) were inverted.16, 23, 37


9


Page 10 of 27

3. Single water solvation (1x SPC): as system SPC-MeOH, but just a single SPC water molecule was either interacting or non-interacting with the solvent (new system). 4. Multiple (dis-)appearing water molecules (5x SPC): as above, but in each state one of five SPC water molecules was interacting with the other 1175 solvent SPC water molecules.23 5. Charge inversion (CI): the charge of a GROMOS 53A6 Chloride ion was changed from +1e to -0.68e in a cubic box (edge length 3.9 nm) solvated with 1926 SPC water molecules.16, 23 6. Full charge inversion (FCI): same as above, but with a charge change from +1e to -1e (new). 7. Sodium ion to Chloride ion conversion (Na-Cl): a GROMOS 53A6 Sodium ion with charge +1e was perturbed to a Chloride ion with charge -1e (new system, preparation as above). 8. Phenylalanine to Alanine conversion (Phe-Ala): a N- and C-capped GROMOS 53A6 Phenylalanine residue was converted to Alanine (new system, preparation as above). 9. Relative binding free energies of five positive allosteric modulators of the glutamate receptor A2 ligand-binding domain (GRA2-M): alchemical conversions of five different modulators of the 3,4-dihydro-2H-1,2,4-benzothiadiazine 1,1-dioxide (BTD) family were performed both in solution and in the protein to calculate the relative binding free energies of the modulators.14 Note that some simulated systems are quite challenging: a charge difference between two states as big as 2e (systems FCI and Na-Cl) requires extensive changes in solvent orientation, while the reversible destruction and formation of a very large cavity requires big reorganization of the solvent (system Phe-Ala). Moreover, the calculation of relative binding free energies of the GRA2 modulators requires sufficient sampling of the possibly large conformational space of the proteinbound state. The largest system of the test cases is undoubtedly the glutamate receptor A2. This dimeric structure of two times 262 amino acids binds two allosteric modulators of the BTD family


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


at the interface. The five modulators differ in the fluorination pattern of the benzene ring of the 4cyclopropyl-3,4-dihydro-2H-1,2,4-benzothiadiazine 1,1-dioxide core structure (non-fluorinated, four different fluorination substitutions at positions 5, 6, 7 or 8).14 Because the simulations are started from the X-ray structure in complex with two very similar modulators (PDB-ID 3TDJ),44 no further (allosteric) conformational changes of the proteins are expected. All MD simulations were performed using the GROMOS MD++ 1.4.0 simulation engine40 to which the A-EDS functionality was added. All systems were simulated in canonical NVT ensembles at 300 K using Nosé-Hoover chains45 with 5 coupled thermostats and relaxation times of 0.1 ps. The solute and solvent degrees of freedom were coupled to different heat-baths. The equations of motion were solved using a leap-frog integration scheme46 with a time-step of 2 fs. Covalent bond lengths were constrained using the SHAKE47 and SETTLE48 algorithms. The center-of-mass motions of the systems were removed every 10000 steps. Neighbour searching was performed every 5 steps using a group-based cutoff scheme within a cutoff sphere of 1.4 nm. For the calculation of non-bonded interactions, a twin-range cutoff scheme was employed. Within a cutoff of 0.8 nm, non-bonded interactions were calculated every step. Between 0.8 nm and 1.4 nm, non-bonded interactions were evaluated every 5 steps, and held fixed between the updates. For the calculation of electrostatic interactions, a reaction-field contribution49 with a relative dielectric permittivity of 6150 beyond the cutoff sphere was added. A-EDS nonequilibrium parameter search simulations and A-EDS equilibrium simulations were both performed for 50 ns, except for systems GRA2-M, for which they were performed for 1 ns and 10 ns, respectively. GRA2-M system preparation and equilibration was done as described earlier.14 Free-energy offset parameters ∆𝐹𝑖𝑅 were set to 0 kJ∙mol-1 at the start of the parameter search simulations, except for systems GRA2-M, for which they were set to the values of the initial


11


Page 12 of 27

energy differences between the state Hamiltonians. Simulations for reference free-energy calculations using BAR were performed by modifying the non-bonded interactions gradually in 21 equally spaced equilibrium simulation points. Every simulation point was equilibrated for 1 ns and then sampled for 10 ns simulation time. For the modification of the nonbonded interactions, soft-core parameters51 of 0.5 for the Lennard-Jones interactions and 0.5 nm2 for the electrostatic interactions were used. Error estimates for all computed free-energy differences were calculated according to Chodera, et al.52 Reference simulations of the GROMOS 53A6 Na+ and Cl- ions were performed for 25 ns. Results and Discussion Two simulations were performed for all test systems. First, a non-equilibrium parameter search was run using the algorithms outlined above, followed by an equilibrium simulation to determine the final free-energy differences. In Table 1, A-EDS optimization parameters (maximum energy ∗ barrier height, ∆𝐸𝑚𝑎𝑥 and sampling time until the free-energy offset parameters are estimated,

𝑡𝑚𝑖𝑛 ) are reported together with the resulting A-EDS parameters. After the parameter search simulation, equilibrium A-EDS simulations were performed to obtain free-energy differences between the states. To demonstrate the robustness of the method, the values for 𝐸𝑚𝑎𝑥 , 𝐸𝑚𝑖𝑛 and ∆𝐹𝑖𝑅 were rounded in these simulations. Moreover, extensive reference free-energy calculations using Bennett’s Acceptance Ratio (BAR)20 were performed to compare the resulting free-energy values, except for system GRA2-M, for which reference data was taken from reference 14. The data for the equilibrium simulations are summarized in Table 2. Except for system FCI, all freeenergy estimates agree to the reference values calculated with BAR within 3.0 kJ∙mol-1, with a root-mean-square difference of 1.3 kJ∙mol-1. The largest remaining difference of 3.0 kJ∙mol-1 is for system Phe-Ala which possibly does not sample quite sufficient switches between the states yet


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


(see Figure S9). The 6th column in Table 2 reports the number of unique visits of the different states. For most systems a significant number of switches between the states take place, except for systems FCI and Phe-Ala. Even during the relatively short simulation of the 2x GRA2-M (bound) system a large number of transitions is observed. Table 1. A-EDS simulation parameters and results for the nonequilibrium parameter search simulations. Energy units are kJ∙mol-1 and time units are ps. A-EDS Parameter Search Simulation Input Parameters

Resulting Data

System

∆𝑬∗𝒎𝒂𝒙

𝒕𝒎𝒊𝒏

𝑬𝒎𝒂𝒙

𝑬𝒎𝒊𝒏

∆𝑭𝑹𝒊≠𝟏

SPC-MeOH

25

10

-17.6

-79.1

10.6

FDI

25

10

-22.1

-207.0

0.1

1x SPC

50

10

-28.5

-37.7

26.4 0.3 0.3

5x SPC

25

10

15.7

-83.8 0.3 0.0

CI

50

100

-81.1

-381.5

48.1

FCI

50

100

-64.8

-440.8

-166.7

Na-Cl

50

100

-77.7

-743.6

16.3

Phe-Ala

40

10

-18.2

-25.1

0.9 -111.0 22.2

1x GRA2-M (unbound)

100

10

42.0

(42.0) 8.6 -54.2 -231.0 55.6

2x GRA2-M (bound)

100

10

151.7

-17.2 18.3 -102.4


13


Page 14 of 27

Table 2. A-EDS simulation parameters and results for the equilibrium sampling simulations. For comparison, resulting free-energy differences obtained with BAR from reference simulations are given. Energy units are kJ∙mol-1. A-EDS Equilibrium Simulation Input Parameters

BAR

Resulting Data

Resulting Data

System

𝑬𝒎𝒂𝒙

𝑬𝒎𝒊𝒏

∆𝑭𝑹𝒊≠𝟏

∆𝑭𝟏,𝒊

Unique State Visits (states 𝒊 ≠ 𝟏)

∆𝑭𝟏,𝒊

SPC-MeOH

-20

-80

10

5.6 ± 0.4

58

6.5 ± 0.1

FDI

-20

-210

0

0.1 ± 0.1

15136

n.a. (0.0)

1x SPC

-28

-38

26

26.0 ± 0.2

319

26.1 ± 0.1

0

0.0 ± 0.3

476

0

-0.4 ± 0.3

486

0

-0.3 ± 0.3

465

0

-0.4 ± 0.3

477

5x SPC

15

-80

n.a. (0.0)

CI

-80

-380

50

49.1 ± 0.4

2800

48.7 ± 0.1

FCI

-65

-440

-165

-171.5 ± 2.5

10

-159.3 ± 0.1

Na-Cl

-80

-740

20

27.4 ± 0.3

1191

24.8 ± 0.1

Phe-Ala

-18

-25

1

2.7 ± 1.8

5

-0.3 ± 0.1

-110

-110.9 ± 0.1

7261

20

22.5 ± 0.1

3475

10

8.7 ± 0.1

7522

-55

-54.2 ± 0.1

2870

-230

-230.6 ± 0.7

2548

55

41.4 ± 0.7

1854

20

21.2 ± 1.2

631

-100

-103.5 ± 0.9

51

1x GRA2-M (unbound)

2x GRA2-M (bound)

n.a.

150

n.a.

n.a.

-20

n.a.

In Table 3, the relative binding free energies of the five different GRA2 ligands computed with AEDS are presented and compared to computational reference values obtained using thermodynamic integration (TI)53 and OSP and experimental reference values obtained using isothermal titration calorimetry (ITC).14 As expected from our previous work, no significant conformational changes were observed. The binding free energy values calculated with A-EDS


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


are in good agreement with the experimentally determined values with a root-mean-square error (RMSE) of 2.5 kJ∙mol-1. Moreover, the accuracy of the automatically calculated A-EDS binding free energies was comparable to the accuracy of the OSP method, which requires the manual design of a reference state Hamiltonian. Table 3. Free-energy differences between the five different GRA2 positive allosteric modulators of the BTD family in solution and when bound to the protein and resulting relative binding free energies computed using A-EDS. Note that the bound state involves both binding sites of the dimeric GRA2. Reference values computed with TI and OSP and experimentally determined using ITC are taken from reference 14. Energy units are kJ∙mol-1. Compound

BPAM429

BPAM408

BPAM344

BPAM442

BPAM411

RMSE to ITC

(R5 fluorinated)

(R6 fluorinated)

(R7 fluorinated)

(R8 fluorinated)

(not fluorinated)

(fluorinated compounds)

𝟐 ∆𝑭𝒖𝒏𝒃𝒐𝒖𝒏𝒅 𝟏,𝒊≠𝟏

-108.4 ± 0.2

17.4 ± 0.2

-221.8 ± 0.2

45.0 ± 0.2

0.0

∆𝑭𝒃𝒐𝒖𝒏𝒅 𝟏,𝒊≠𝟏

-103.5 ± 0.9

21.2 ± 1.2

-230.6 ± 0.7

41.4 ± 0.7

0.0

A-EDS ∆∆𝑭𝒃𝒊𝒏𝒅

4.9 ± 0.9

3.8 ± 1.2

-8.8 ± 0.7

-3.6 ± 0.7

0.0

2.5

TI14 ∆∆𝑭𝒃𝒊𝒏𝒅

-2.9

3.4

-2.5

4.7

0.0

4.5

OSP14 ∆∆𝑭𝒃𝒊𝒏𝒅

7.0

0.1

-4.9

1.4

0.0

2.5

ITC14 ∆∆𝑮𝒃𝒊𝒏𝒅

2.8

1.4

-6.0

-0.9

0.0

Notably, the free energy offset parameters ∆𝐹𝑖𝑅 obtained from the A-EDS non-equilibrium parameter search simulations already match the equilibrium free-energy differences from A-EDS and BAR quite reasonably. Remaining discrepancies are attributed to the non-equilibrium character of the search simulation or, in case of the protein system GRA2-M, to incomplete sampling of the large conformational space of the protein during the comparatively short parameter search simulation. All parameter search simulations converged within 5-10 ns simulation time, except for system Phe-Ala, for which convergence was reached after 15-20 ns. Example energy trajectories during a parameter search simulation are given in Figure 2 for system SPC-MeOH. After the first 10 ns, the parameters seem converged and the system regularly switches between


15


Page 16 of 27

the two end-states. Energy trajectories for all performed A-EDS simulations are given in Figures S1-S11 in the Supporting Information. Similarly, the equilibrium simulations seem converged to within 3 kJ∙mol-1 of the final estimate in 1-15 ns of simulation time. This means that in an overall simulation time of 30 ns over two simulations, free-energy estimates are obtained for all systems that are very comparable to the BAR data obtained from 210 ns of simulations. More detailed analysis of the convergence characteristics of the method will be the subject of future work.

Figure 2. Energies and A-EDS parameters over time during a non-equilibrium A-EDS parameter search simulation (system SPC-MeOH). The energy of the system according to the Hamiltonian of state 1 (SPC, 𝐻1 (𝑟⃑)) is indicated with grey dots and the energy of the system according to state 2 (MeOH, 𝐻2 (𝑟⃑) − ∆𝐹2𝑅 ) is indicated with brown dots. The energy offset parameter of state 2 (∆𝐹2𝑅 ) is shown with a black solid line and A-EDS parameters 𝐸𝑚𝑖𝑛 and 𝐸𝑚𝑎𝑥 are shown with red solid lines. To illustrate the nature and the effectiveness of an A-EDS reference state Hamiltonian, the radial dipole orientation correlation function (ROCF) and radial distribution function (RDF) of the SPC water molecules around the A-EDS reference state of system Na-Cl are shown in Figure 3A. The Na-Cl reference state induces water orientation and occupancy both corresponding to Na+ and to Cl- in a single simulation. This becomes especially apparent if only the trajectory regions corresponding to one of the end-states are considered (red and blue dots). Figure 3B shows the probability distributions of the energy differences between the end-state Hamiltonians in the reference state simulation of system Na-Cl and the same probabilities reweighted54 to the end-state Hamiltonians. While the reference state samples the important energy regions for both end-states,


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


sampling is concentrated to intermediate energy regions not important for either end-state. We were able to fine-tune the A-EDS reference state Hamiltonian of system Na-Cl to sample more important energy regions simply by increasing the A-EDS parameter 𝐸𝑚𝑎𝑥 by 50 kJ∙mol-1. However, this also resulted in greatly reduced occurrences of state-transitions because of the increased energy barrier between the states (Figure S7 and S8 in the Supporting Information). In ∗ general, a too high value of ∆𝐸𝑚𝑎𝑥 will lead to insufficient sampling while a too low value of ∗ ∆𝐸𝑚𝑎𝑥 will lead to additional sampling of alternative, irrelevant phase-space. A

B

Figure 3. (A) Radial dipole orientation correlation functions (ROCF) and radial distribution functions (RDF) of the SPC water oxygen atoms around the simulated ions. Data from the A-EDS reference state of system Na-Cl are shown with thick black lines and reference data from simulations of Na+ and Cl- are shown with thin red and blue lines, respectively. Data from the AEDS reference state trajectory regions of the Na+ state (𝐻1 (𝑟⃑) ≤ 400 kJ∙mol-1) is shown with red dots and of the Cl- state (𝐻2 (𝑟⃑) − ∆𝐹2𝑅 ≤ 400 kJ∙mol-1) with blue dots. (B) Probability distributions of the energy differences 𝐻2 (𝑟⃑) − 𝐻1 (𝑟⃑) of the A-EDS reference state of system Na-Cl (black line), and of the energy differences reweighted to state 1 (Na+, red line) and state 2 (Cl-, blue line). The probability distribution of the energy difference obtained with a different A-EDS reference Hamiltonian, in which the parameter 𝐸𝑚𝑎𝑥 has been increased by 50 kJ∙mol-1, is shown with an orange line. We introduced a novel approach for the construction of an EDS reference state Hamiltonian which, in contrast to previously suggested approaches, conserves local energy minima of the combined end-states. Our approach allows for simple tuning of the EDS reference Hamiltonian parameters, without the need for complex search algorithms. The parameter search is intuitive, as only a single parameter has to be chosen prior to the parameter search simulation: the maximum


17


Page 18 of 27

∗ energy barrier height between the end-states (∆𝐸𝑚𝑎𝑥 ). The free-energy offset parameters which

can be obtained from a parameter search simulation ideally correspond to the free-energy differences between the combined end-states. Hence, the A-EDS approach allows for the automated calculation of multiple free-energy differences from a single MD simulation. Even if only applied on pairwise free-energy differences, this implies an enormous step forward for efficiency and user-friendliness of free-energy calculations in e.g. drug design. ASSOCIATED CONTENT Supporting Information. The following files are available free of charge. Energy trajectories for all performed non-equilibrium and equilibrium A-EDS simulations (Figures S1-S11) (PDF) AUTHOR INFORMATION Notes The authors declare no competing financial interests. ACKNOWLEDGMENT The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC).


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


REFERENCES 1.

Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.;

Dahlgren, M. K.; Greenwood, J.; et al., Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. J Am Chem Soc 2015, 137 (7), 2695-703. 2.

Tembe, B. L.; Mc Cammon, J. A., Ligand-receptor interactions. Comput Chem 1984, 8 (4),

281-283. 3.

Hansen, N.; van Gunsteren, W. F., Practical Aspects of Free-Energy Calculations: A

Review. J Chem Theory Comput 2014, 10 (7), 2632-47. 4.

Zwanzig, R. W., High‐Temperature Equation of State by a Perturbation Method. I.

Nonpolar Gases. J Chem Phys 1954, 22 (8), 1420-6. 5.

Mark, A. E.; Xu, Y.; Liu, H.; van Gunsteren, W. F., Rapid non-empirical approaches for

estimating relative binding free energies. Acta Biochim Pol 1995, 42 (4), 525-35. 6.

Liu, H.; Mark, A. E.; van Gunsteren, W. F., Estimating the Relative Free Energy of

Different Molecular States with Respect to a Single Reference State. J Phys Chem 1996, 100 (22), 9485-9494. 7.

Oostenbrink, B. C.; Pitera, J. W.; van Lipzig, M. M. H.; Meerman, J. H. N.; van Gunsteren,

W. F., Simulations of the Estrogen Receptor Ligand-Binding Domain: Affinity of Natural Ligands and Xenoestrogens. J Med Chem 2000, 43 (24), 4594-4605. 8.

Pitera, J. W.; van Gunsteren, W. F., One-Step Perturbation Methods for Solvation Free

Energies of Polar Solutes. J Phys Chem B 2001, 105 (45), 11264-11274.


19


9.

Page 20 of 27

Oostenbrink, C., Efficient free energy calculations on small molecule host-guest systems -

a combined linear interaction energy/one-step perturbation approach. J Comput Chem 2009, 30 (2), 212-21. 10. Oostenbrink, C., Free energy calculations from one-step perturbations. Methods Mol Biol 2012, 819, 487-99. 11. Raman, E. P.; Vanommeslaeghe, K.; Mackerell, A. D., Jr., Site-Specific Fragment Identification Guided by Single-Step Free Energy Perturbation Calculations. J Chem Theory Comput 2012, 8 (10), 3513-3525. 12. Graf, M. M.; Zhixiong, L.; Bren, U.; Haltrich, D.; van Gunsteren, W. F.; Oostenbrink, C., Pyranose dehydrogenase ligand promiscuity: a generalized approach to simulate monosaccharide solvation, binding, and product formation. PLoS Comput Biol 2014, 10 (12), e1003995. 13. Lai, B.; Nagy, G.; Garate, J. A.; Oostenbrink, C., Entropic and Enthalpic Contributions to Stereospecific Ligand Binding from Enhanced Sampling Methods. J Chem Inf Model 2014, 54 (1), 151-158. 14. Norholm, A. B.; Francotte, P.; Goffin, E.; Botez, I.; Danober, L.; Lestage, P.; Pirotte, B.; Kastrup, J. S.; Olsen, L.; Oostenbrink, C., Thermodynamic characterization of new positive allosteric modulators binding to the glutamate receptor A2 ligand-binding domain: combining experimental and computational methods unravels differences in driving forces. J Chem Inf Model 2014, 54 (12), 3404-16. 15. Jandova, Z.; Fast, D.; Setz, M.; Pechlaner, M.; Oostenbrink, C., Saturation mutagenesis by efficient free-energy calculation. J Chem Theory Comput 2017, 14 (2), 894-904.


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


16. Christ, C. D.; van Gunsteren, W. F., Enveloping distribution sampling: a method to calculate free energy differences from a single simulation. J Chem Phys 2007, 126 (18), 184110. 17. Han, K.-K., A new Monte Carlo method for estimating free energy and chemical potential. Phys Lett A 1992, 165 (1), 28-32. 18. Wu, D.; Kofke, D. A., Phase-space overlap measures. II. Design and implementation of staging methods for free-energy calculations. J Chem Phys 2005, 123 (8), 084109. 19. Chen, Y. G.; Hummer, G., Slow conformational dynamics and unfolding of the calmodulin C-terminal domain. J Am Chem Soc 2007, 129 (9), 2414-5. 20. Bennett, C. H., Efficient Estimation of Free Energy Differences from Monte Carlo Data. J Comput Phys 1976, 22 (2), 245-68. 21. Han, K. K., Multiensemble sampling: An alternative efficient Monte Carlo technique. Phys Rev E 1996, 54 (6), 6906-6910. 22. Lu, N.; Wu, D.; Woolf, T. B.; Kofke, D. A., Using overlap and funnel sampling to obtain accurate free energies from nonequilibrium work measurements. Phys Rev E 2004, 69 (5), 057702. 23. Christ, C. D.; van Gunsteren, W. F., Multiple free energies from a single simulation: extending enveloping distribution sampling to nonoverlapping phase-space distributions. J Chem Phys 2008, 128 (17), 174112. 24. Christ, C. D.; van Gunsteren, W. F., Simple, Efficient, and Reliable Computation of Multiple Free Energy Differences from a Single Simulation: A Reference Hamiltonian Parameter Update Scheme for Enveloping Distribution Sampling (EDS). J Chem Theory Comput 2009, 5 (2), 276-286.


21


Page 22 of 27

25. Christ, C. D.; Van Gunsteren, W. F., Comparison of three enveloping distribution sampling Hamiltonians for the estimation of multiple free energy differences from a single simulation. J Comput Chem 2009, 30 (11), 1664-79. 26. Riniker, S.; Christ, C. D.; Hansen, N.; Mark, A. E.; Nair, P. C.; van Gunsteren, W. F., Comparison of enveloping distribution sampling and thermodynamic integration to calculate binding free energies of phenylethanolamine N-methyltransferase inhibitors. J Chem Phys 2011, 135 (2), 024105. 27. Hansen, N.; Dolenc, J.; Knecht, M.; Riniker, S.; van Gunsteren, W. F., Assessment of enveloping distribution sampling to calculate relative free enthalpies of binding for eight netropsin-DNA duplex complexes in aqueous solution. J Comput Chem 2012, 33 (6), 640-51. 28. Lin, Z.; Liu, H.; Riniker, S.; van Gunsteren, W. F., On the Use of Enveloping Distribution Sampling (EDS) to Compute Free Enthalpy Differences between Different Conformational States of Molecules: Application to 310-, α-, and π-Helices. J Chem Theory Comput 2011, 7 (12), 38843897. 29. Lin, Z.; van Gunsteren, W. F., Enhanced conformational sampling using enveloping distribution sampling. J Chem Phys 2013, 139 (14), 144105. 30. Lin, Z.; van Gunsteren, W. F., Combination of Enveloping Distribution Sampling (EDS) of a Soft-Core Reference-State Hamiltonian with One-Step Perturbation to Predict the Effect of Side Chain Substitution on the Relative Stability of Right- and Left-Helical Folds of beta-Peptides. J Chem Theory Comput 2013, 9 (1), 126-34.


22

Page 23 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


31. Huang, W.; Lin, Z.; van Gunsteren, W. F., Use of enveloping distribution sampling to evaluate important characteristics of biomolecular force fields. J Phys Chem B 2014, 118 (24), 6424-30. 32. Lee, J.; Miller, B. T.; Damjanović, A.; Brooks, B. R., Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange. J Chem Theory Comput 2014, 10 (7), 2738-2750. 33. Lee, J.; Miller, B. T.; Damjanović, A.; Brooks, B. R., Enhancing Constant-pH Simulation in Explicit Solvent with a Two-Dimensional Replica Exchange Method. J Chem Theory Comput 2015, 11 (6), 2560-2574. 34. Lee, J.; Miller, B. T.; Brooks, B. R., Computational scheme for pH-dependent binding free energy calculation with explicit solvent. Protein Sci 2016, 25 (1), 231-43. 35. Sidler, D.; Schwaninger, A.; Riniker, S., Replica exchange enveloping distribution sampling (RE-EDS): A robust method to estimate multiple free-energy differences from a single simulation. J Chem Phys 2016, 145 (15), 154114. 36. Sidler, D.; Cristofol-Clough, M.; Riniker, S., Efficient Round-Trip Time Optimization for Replica-Exchange Enveloping Distribution Sampling (RE-EDS). J Chem Theory Comput 2017, 13 (6), 3020-3030. 37. Mori, T.; Hamers, R. J.; Pedersen, J. A.; Cui, Q., Integrated Hamiltonian sampling: a simple and versatile method for free energy simulations and conformational sampling. J Phys Chem B 2014, 118 (28), 8210-20.


23


Page 24 of 27

38. Miao, Y.; Feher, V. A.; McCammon, J. A., Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J Chem Theory Comput 2015, 11 (8), 3584-3595. 39. Hamelberg, D.; Mongan, J.; McCammon, J. A., Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. J Chem Phys 2004, 120 (24), 1191929. 40. Schmid, N.; Christ, C. D.; Christen, M.; Eichenberger, A. P.; van Gunsteren, W. F., Architecture, implementation and parallelisation of the GROMOS software for biomolecular simulation. Comput Phys Commun 2012, 183 (4), 890-903. 41. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J., Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces, Pullman, B., Ed. Springer Netherlands: 1981; Vol. 14, pp 331-42. 42. Walser, R.; Mark, A. E.; van Gunsteren, W. F.; Lauterbach, M.; Wipff, G., The effect of force-field parameters on properties of liquids: Parametrization of a simple three-site model for methanol. J Chem Phys 2000, 112 (23), 10450-10459. 43. Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F., A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6. J Comput Chem 2004, 25 (13), 1656-76. 44. Krintel, C.; Frydenvang, K.; Olsen, L.; Kristensen, M. T.; de Barrios, O.; Naur, P.; Francotte, P.; Pirotte, B.; Gajhede, M.; Kastrup, J. S., Thermodynamics and structural analysis of


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


positive allosteric modulation of the ionotropic glutamate receptor GluA2. Biochem J 2012, 441 (1), 173-8. 45. Martyna, G. J.; Klein, M. L.; Tuckerman, M., Nosé–Hoover chains: The canonical ensemble via continuous dynamics. J Chem Phys 1992, 97 (4), 2635-2643. 46. Hockney, R. W., The Potential Calculation and some Applications. Academic Press Inc.: New York, 1970; p 77 p. 47. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C., Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys 1977, 23 (3), 327-341. 48. Miyamoto, S.; Kollman, P. A., Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J Comput Chem 1992, 13 (8), 952-962. 49. Tironi, I.; Sperb, R.; Smith, P.; van Gunsteren, W. F., A Generalized Reaction Field Method for Molecular Dynamics Simulation. J Chem Phys 1995, 102, 5451. 50. Heinz, T. N.; van Gunsteren, W. F.; Hünenberger, P. H., Comparison of Four Methods to Compute the Dielectric Permittivity of Liquids from Molecular Dynamics Simulations. J Chem Phys 2001, 115 (3), 1125-1136. 51. Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F., Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem Phys Lett 1994, 222 (6), 529-539.


25


Page 26 of 27

52. Chodera, J. D.; Swope, W. C.; Pitera, J. W.; Seok, C.; Dill, K. A., Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations. J Chem Theory Comput 2007, 3 (1), 26-41. 53. Kirkwood, J. G., Statistical Mechanics of Fluid Mixtures. J Chem Phys 1935, 3 (5), 30013. 54. Torrie, G. M.; Valleau, J. P., Nonphysical Sampling Distributions in Monte Carlo FreeEnergy Estimation: Umbrella sampling. J Comput Phys 1977, 23 (2), 187-199.


26

Page 27 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


TOC GRAPHICS


27

Accelerated Enveloping Distribution Sampling: Enabling Sampling of

Recommend Documents