Statistical Mechanics of LigandâReceptor Noncovalent ... - FloRe - UniFI

Subscriber access provided by University of Newcastle, Australia

Article

The statistical mechanics of ligand-receptor non covalent association, revisited: binding site and standard state volumes in modern alchemical theories Piero Procacci, and Riccardo Chelli J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01192 • Publication Date (Web): 07 Apr 2017 Downloaded from http://pubs.acs.org on April 9, 2017

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 31

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

The statistical mechanics of ligand-receptor non covalent association, revisited: binding site and standard state volumes in modern alchemical theories Piero Procacci∗ and Riccardo Chelli Department of Chemistry, University of Florence, Via Lastruccia n. 3, Sesto Fiorentino, I-50019 Italy E-mail: [email protected]

Abstract The present paper is intended to be a comprehensive assessment and rationalization, from a statistical mechanics perspective, of existing alchemical theories for binding free energy calculations of ligand-receptor systems. More in detail, the statistical mechanics foundation of non covalent interactions in ligand-receptor systems is revisited, providing a unifying treatment that encompasses the most important variants in the alchemical approaches, from the seminal Double Annihilation Method [W. L. Jorgensen and J. K. Buckner and S. Boudon and J. TiradoRives, J. Chem. Phys. 89,3742, 1988], to the Double Decoupling Method [M. K. Gilson and J. A. Given and B. L. Bush and J. A. McCammon, Biophys. J. 72, 1047 1997] and the Deng and Roux alchemical theory [Y. Deng and B. Roux, J. Chem. Theory Comput. 2, 1255 2006]. Connections and differences between the various alchemical approaches are highlighted and discussed. ∗ To

whom correspondence should be addressed

1 ACS Paragon Plus Environment


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 The determination of the binding free energy in ligand-receptor systems is the cornerstone of drug discovery. In the last decades, traditional molecular docking techniques in computer assisted drug design have been modified, integrated or superseded using methodologies relying on a more realistic description of the drug-receptor system. It has becoming increasing clear that a microscopic description of the solvent is a crucial ingredient to reliably rank the affinity of putative ligands for a given target. 1,2 In the framework of atomistic molecular dynamics (MD) simulations with explicit solvent, several methods for rigorously determining the absolute binding free energy in noncovalently bonded systems have been devised. Most of these methodologies are based on the so-called alchemical route 3–5 (see Refs.

6–9

for recent reviews). In this approach, proposed for the

first time by Jorgensen et al. 10,11 and named Double Annihilation Method (DAM), the absolute binding free energy of the ligand-receptor system is obtained by setting up a thermodynamic cycle as indicated in Figure 1 in which the basic quantities to be estimated are the annihilation or decoupling free energies of the ligand in the bound state and in bulk water, indicated hereinafter with ∆Gb and ∆Gu , respectively. These decoupling free energies correspond to the two closing branches of the cycle and are obtained by discretizing the alchemical path connecting the fully interacting and fully decoupled ligand in a number of intermediate nonphysical states, and then running MD simulations for each of these states. Alchemical states are defined by a λ coupling parameter entering the Hamiltonian, varying between 1 and 0, such that at λ = 1 and at λ = 0 one has the fully interacting and gas-phase ligand, respectively. ∆Gb and ∆Gu are usually recovered as a sum of the contributions from each of the λ windows by applying the free energy perturbation method 12 (FEP). Alternatively, and equivalently, one can compute the canonical average of the derivative of the Hamiltonian at the discrete λ points, obtaining the decoupling free energy via numerical thermodynamic integration (TI). 13 Finally, the thermodynamic cycle is closed by computing the difference between the two decoupling free energies along the alchemical path, ∆Gb and ∆Gu , obtaining the dissociation free energy in solution. 2 ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

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


Rsol--- Lgas

∆Ggas=0

Rsol + Lgas

-∆Gu

∆Gb

Rsol--- Lsol

∆G0

Rsol + Lsol

Figure 1: The alchemical thermodynamic cycle for computing the absolute dissociation free energy, ∆G0 , in ligand-receptor systems. The subscript “sol” and “gas” indicates solvated and gasphase (decoupled) species, respectively. The ligand-environment interactions must be turned off in the solvated complex and in bulk solvent obtaining ∆G0 = ∆Gb − ∆Gu . Gilson et al. 14 criticized Jorgensen’s theory 10 by pointing out that the resulting dissociation free energy does not depend upon the choice of standard concentration. In order to define a reference chemical potential for the decoupling ligand when bound to the receptor, Gilson introduced a “restraint” that keeps the ligand in a volume Vr around the binding place. This restraint is shown to yield 14–16 an additive standard state dependent correction to the dissociation free energy of kB T ln(Vr /V 0 ), interpreted as a chemical potential difference of the ligand at concentration 1/Vr and at the standard concentration 1/V 0 . Gilson et al. named this variant of DAM, Double Decoupling Method (DDM). According to Gilson, the effect of progressively strengthening the restraint, 17 leading to a more negative correction, should be balanced by a larger work integral so that “errors will occur only when the integration region defined by the restraint volume becomes so small that conformations that ought to make important contributions to the work integral are missed.” DDM was further developed by Karplus and coworkers 15 who introduced a versatile set of harmonic restraints to be used in MD simulations, restricting both the position and the orientation



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

of the ligand. Hamelberg and McCammon 18 posed the “questions as to how sensitive are the calculated free energies to the strength of the restraining potentials and whether an arbitrary choice of the restraint potential would yield the correct results”. The recipe they proposed is to set the force constant in the harmonic restraining potential according to the mean fluctuation of the restrained coordinate in a preliminary unrestrained simulation. Subsequently, Deng and Roux 19 elaborated a DDM variant whereby the standard state correction, “in the strong restraint limit”, is no longer dependent on the imposed restraint volume Vr hence resolving the alluded inconsistency 17,20 in the Gilson theory. However, it does apparently require the estimate of an unknown translational, rotational and conformational binding site volume Vsite in the complex, 21 via an independent unrestrained simulation of the bound state or via an auxiliary FEP for turning on the restraint when the fully coupled ligand is in the binding site. 19,22,23 In a series of recent papers, Fujitani and coworkers, 24–27 in their “Massively Parallel Computation of Absolute Binding Free Energy (MP-CAFEE)”, successfully applied the unrestrained DAM approach to several drug-receptor systems, in many cases accurately predicting the experimental dissociation free energy via FEP. Most importantly, these authors directly compared their DAM/FEP values, ∆Gbind = ∆Gu − ∆Gb , to the experimental value ∆G0 , openly criticizing the DDM standard state correction 25 : “as far as we know there is no theoretical or experimental proof that [the standard state corrected] ∆G0 meets the definition of the absolute binding energy.[..] Therefore, we directly compare ∆Gbind with ∆G0 .” The standard state correction issue can be bypassed altogether by computing relative binding free energies, 3,28 due to the transmutation of an unrestrained ligand into another in the same binding site and in the solvent. Relative binding free energy calculations neglect altogether the possibility of a change of binding site volume due to the transmutation or even the possibility that the transmuting ligand may leave the binding site at some alchemical state. This approach is hence limited to the assessment of binding affinities in strictly congeneric series of ligands, with the tacit assumption of a constant binding site volume upon transmutation.


Page 4 of 31

Page 5 of 31

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 conclusion, the question of the standard state correction, or, equivalently, the issue of the binding site volume in drug-receptor dissociation free energy calculations is either ignored, as in calculations of relative free energies, or treated using methodologies relying on the definition of arbitrary sets of constraints whose effects on the resulting free energy has never been convincingly assessed. In any case, the standard state issue and the strictly connected binding site volume definition, that is indeed crucial for a reliable MD-based tool in drug discovery, is still far from being settled. In this contribution, we try to address these issues, revisiting the the DAM and DDM theory with a spotlight on the relationship between the binding site volume, restraint volume and standard state volume in noncovalent association, providing a unifying treatment encompassing Jorgensen, Gilson and Deng and Roux theories.

Corpora non agunt nisi fixata: 29 Alchemical theory with restraints Molecular recognition in host-guest or drug-receptor non covalent interactions is based on a highly specific molecular complementary, 30 translating in the existence of a single overwhelmingly prevalent binding “pose”. The latter can be defined using an appropriate set of coordinates that are functions of the ligand and receptor Cartesian coordinates x. A natural coordinate in ligand-receptor binding is represented by the vector R connecting the center of mass (COM) of the ligand to the COM of the receptor. The vector R defines the precise location of the ligand COM on the receptor surface in the bound state. The orientation of the ligand in the binding site can be specified using the three Euler angles Ω between the receptor frame and the ligand frame or using alternative angular coordinates based on a subset of ligand and receptor atoms. 15,16 The conformational state in the binding site can be further specified by providing a set of internal coordinates χ of the ligand and the receptor. On the overall, a set of translational, orientational and conformational coordinates, Ω(x),χ χ (x)}, may be used to structurally define the complex. From a theoretical Y(x) ≡ {R(x),Ω standpoint, this is formally achieved, as we shall see further on, by defining an Heaviside step function Θ(Y) that is equal to 1 in a domain marking the relevant ligand-receptor coordinates Y for the associated state. 31 5 ACS Paragon Plus Environment


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 6 of 31

In most applications for the determination of the dissociation free energy using the alchemical method, a set of harmonic restraints are introduced on the Y set of coordinates in order to keep the ligand in the binding pose, while the decoupling process proceeds. The restraint potential in the Y coordinates has the general form 1 Vr (x) = (Y(x) − Yc )t K(Y(x) − Yc ) 2

(1)

where K is the diagonal matrix of the harmonic force constants. Note that the function e−βVr may be interpreted as a non normalized multivariate Gaussian distribution in the space defined by the Ω(x),χ χ (x)}, i.e. coordinates Y = {R(x),Ω T Σ −1 (Y(x)−Y ) c r

1

e−βVr (x) = e− 2 (Y(x)−Yc )

(2)

where the diagonal covariance matrix Σ r is defined as Σ r = kB T K−1

(3)

The parameters Yc in the restraint potential should be chosen so as to match the true mean values of the vector Y in the unrestrained complex. In Refs., 32,33 in the context of single molecule pulling experiments, a simple relation was derived between the free energy of the driven system (i.e. with Hamiltonian including the harmonic potential of an external device coupled to a specific molecular distance R) and the free energy of the system with unperturbed Hamiltonian along the driven coordinate (i.e. the potential of mean force along R). The relation proposed by Marsili (Eq. 7 in Ref. 32 ) can be straightforwardly applied to the Y coordinate for any of the restrained λ alchemical state in DDM as:

ρr (Y|Yc , λ ) Gr (Yc , λ ) = G(Y, λ ) +Vr (Y − Yc ) + kB T ln ρ(Y∗ )


(4)

Page 7 of 31

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 we have used the subscript r when referring to quantities of the restrained system and where Z −β [H(x,λ )+Vr (Y(x)−Yc )] Gr (Yc , λ ) = −kB T ln C dxe # "R ρ(Y) dxδ (Y − Y(x))e−β H(x,λ ) = −k T ln G(Y, λ ) = −kB T ln R B ρ(Y∗ ) dxδ (Y∗ − Y(x))e−β H(x,λ )

(5) (6)

Here, Gr (Yc , λ ) is the free energy of the restrained system (C is a constant that makes argument of the logarithm adimensional) and G(Y, λ ) is the free energy of the unrestrained system at Y with respect to some immaterial reference state at Y∗ . In Eqs. 5 and 6, H(x, λ ) is the Hamiltonian at the alchemical state λ , with x encompassing all solvent, ligand and receptor coordinates. ρr (Y|Yc , λ ) ≡ hδ (Y − Y(x))iYc ,λ , finally, is the canonical probability density evaluated at Y for the system restrained around Yc with free energy at the alchemical state λ given by Eq. 5. In the alchemical decoupling of the complex (left branch of the cycle in Figure 1), one computes, via FEP or TI, the free energy difference between the states at λ = 1 (interacting ligand) and λ = 0 (non interacting ligand), subject to the restraint potential Vr (x), Eq. 1. We can express this difference, using the general Eq. 4, therefore obtaining the Y∗ independent relation Ωc ,χ χ c , 1) Ωc ,χ χ c , 0) − Grb (Rc ,Ω Ωc ,χ χ c ) = Grb (Rc ,Ω ∆Grb (Rc ,Ω Ω,χ χ ) + kB T ln = ∆G(R,Ω

Ω,χ χ |Rc ,Ω Ωc ,χ χ c , 0) ρbr (R,Ω r Ω,χ χ |Rc ,Ω Ωc ,χ χ c , 1) ρb (R,Ω

(7)

Ω,χ χ } denoting the set of transwhere, for the time being, have used the expanded notation Y = {R,Ω Ωc ,χ χ c) lational, rotational and internal ligand-receptor coordinates. In Eq. 7, the quantity ∆Grb (Rc ,Ω corresponds to decoupling free energy of the restrained complex, while "R

Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,0) dxδ (R − R(x))δ (Ω Ω,χ χ ) = −kB T ln R ∆G(R,Ω Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,1) dxδ (R − R(x))δ (Ω

# (8)

Ω,χ χ }. Note that, is the decoupling free energy of the unrestrained system evaluated at Y = {R,Ω Ωc ,χ χ c } in going from the initial (coupled) since there is no change in the parameters Yc = {Rc ,Ω



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 31

to the final (decoupled) state, there can’t be correspondingly no change in the harmonic potential Ω,χ χ } due to the restraint. energy at Y = {R,Ω Similarly, it can be shown that, if we restraint the dissociated state of the ligand-receptor system χ c using a harmonic potential, the decoupling (R = ∞) at the rotational/conformational state Ω c ,χ free energy of the ligand in this restrained dissociated state (i.e. in the right branch of the cycle of Figure 1) is given by

Ωc ,χ χ c ) = ∆G(R∞ ,Ω Ω,χ χ ) + kB T ln ∆Gru (Ω

Ω,χ χ |Ω Ωc ,χ χ c , 0) ρur (Ω r Ω,χ χ |Ω Ωc ,χ χ c , 1) ρu (Ω

(9)

where, "R

Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,0) dxδ (R∞ − R(x))δ (Ω Ω,χ χ ) = −kB T ln R ∆G(R∞ ,Ω Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,1) dxδ (R∞ − R(x))δ (Ω

# (10)

Ω,χ χ. is the decoupling free energy of the unrestrained ligand in the unbound state evaluated at R∞ ,Ω Here R∞ represents a ligand-receptor COM distance that is large enough to allow the ligand and the Ω,χ χ |Ω Ωc ,χ χ c , 1) and receptor to interact only with the solvent when λ 6= 0. Similarly to Eq. 7, ρur (Ω Ω,χ χ |Ω Ωc ,χ χ c , 0) are the probability densities of the Ω ,χ χ coordinates for the dissociated (free) ρur (Ω ligand, orientationally and conformationally restrained, evaluated at λ = 1 and at λ = 0, respecΩ,χ χ ) in Eq. 10 represents the reversible work to bring the ligand in the dissociated tively. ∆G(R∞ ,Ω state from the bulk into the gas-phase when the system is in the ro-vibrational states defined by the χ . All relative ligand-receptor rotational states at R∞ (i.e. for the separated species) have vector Ω ,χ Ω,χ χ ) ≡ ∆G(R∞ ,χ χ ) depends only on the conformational state equal weights 1/8π 2 so that ∆G(R∞ ,Ω χ. In this regard, we notice that the restrained coordinate χ in Eqs. 7, 9 may in principle include also conformational coordinates of the receptor. In such a case, since only the ligand is annihilated in right branch of the cycle of Figure 1, the alchemical contribution to the decoupling free energy comes only from the ligand with no contribution from the χ -restrained unbound protein, Ωc ,χ χ c ) neither in the probability density ratio ρur (Ω Ω,χ χ |Ω Ωc ,χ χ c , 1)/ρur (Ω Ω,χ χ |Ω Ωc ,χ χ c , 0). nor in ∆Gru (Ω 8 ACS Paragon Plus Environment

Page 9 of 31

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


Nonetheless, as we shall see further on, restraining the conformational state of the protein does indeed have an impact on the alchemically computed dissociation free energy. By subtracting Eq. 10 and Eq. 8, we obtain "R

Ω,χ χ ) − ∆G(R∞ ,Ω Ω,χ χ) = ∆G(R,Ω + = =

# Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,0) dxδ (R − R(x))δ (Ω −kB T ln R + Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,1) dxδ (R − R(x))δ (Ω # "R Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,0) dxδ (R∞ − R(x))δ (Ω kB T ln R Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,1) dxδ (R∞ − R(x))δ (Ω "R # Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,1) dxδ (R∞ − R(x))δ (Ω −kB T ln R Ω −Ω Ω(x))δ (χ χ − χ (x))e−β H(x,1) dxδ (R − R(x))δ (Ω Ω,χ χ) −W (R,Ω (11)

where we have exploited the fact that the probability densities of the decoupled ligand and receptor Ω,χ χ ) represents the reversible work for bringing the separated (λ = 0) are uniform over R. W (R,Ω χ arrangement into the corresponding bound conformation at R. (R∞ ) ligand and receptor in the Ω ,χ Ω,χ χ ) for R∞ does not depend on Ω but may depend on χ , it is convenient As the function W (R,Ω Ω,χ χ ) in terms of potentials of mean force, evaluated with to express the reversible work W (R,Ω respect to an arbitrary reference conformational state χ ∗ for the separated species, i.e. Ω,χ χ ) ≡ w(R,Ω Ω,χ χ ) − w(R∞ ,χ χ) W (R,Ω

(12)

with " R

# −β H(x,1) Ω Ω χ χ dxδ (R − R(x))δ (Ω −Ω (x))δ (χ − (x))e Ω,χ χ ) = −kB T ln R w(R,Ω Ω −Ω Ω(x))δ (χ χ ∗ − χ (x))e−β H(x,1) dxδ (R∞ − R(x))δ (Ω "R # −β H(x,1) Ω Ω χ χ dxδ (R∞ − R(x))δ (Ω −Ω (x))δ (χ − (x))e χ )) = −kB T ln R w(R∞ ,χ Ω −Ω Ω(x))δ (χ χ ∗ − χ (x))e−β H(x,1) dxδ (R∞ − R(x))δ (Ω

(13)

χ ) represents the reversible work to bring the separated system from the χ ∗ In Eq. 13, w(R∞ ,χ Ω,χ χ ) is the reversible conformational arrangement to the new conformational state χ , and w(R,Ω work to bring the separated system at χ ∗ and with relative orientation expressed by Ω to the bound



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 31

Ω,χ χ. state defined by the coordinates R,Ω In order to close the alchemical thermodynamic cycle, we subtract Eq. 9 from Eq. 7 and use Eq. 11, by setting Y = Yc : Ωc ,χ χ c ) − ∆Gru (R∞ ,Ω Ωc ,χ χ c ) = −[w(Rc ,Ω Ωc ,χ χ c ) − w(R∞ ,χ χ c )] + ∆Grb (Rc ,Ω Ωc ,χ χ c |0) ρbr (Rc ,Ω Ωc ,χ χ c |0) ρur (Ω + kB T ln r − ln r (14) Ωc ,χ χ c |1) Ωc ,χ χ c |1) ρb (Rc ,Ω ρu (Ω where, in order to to remove the excess baggage from the formalism, we have redefined the probability density for the restrained system as

Ωc ,χ χ c |Rc ,Ω Ωc ,χ χ c, λ ) Ωc ,χ χ c |λ ) ≡ ρbr (Rc ,Ω ρbr (Rc ,Ω Ωcχ c |λ ) ≡ ρur (Ω Ωc ,χ χ c |Ω Ωc ,χ χ c, λ ) ρur (Ω (15)

with λ = 0 or λ = 1. Equation 14 expresses the fact that the dissociation free energy with a set of harmonic restraints of the kind of Eq. 1 and computed in alchemical simulations via FEP or TI, namely the quantity

Ωc ,χ χ c ) = ∆Grb (Rc ,Ω Ωc ,χ χ c ) − ∆Gru (R∞ ,Ω Ωc ,χ χ c) ∆Grd (Rc ,Ω

(16)

is equal to minus the reversible work for forming the complex plus a correction related to the logarithm of the ratio of the canonical probability densities for the restrained decoupled and coupled bound and separated states. It should be stressed that for Eq. 14 to be valid, the canonical probΩ,χ χ |0, 1) and ρur (Ω Ω,χ χ |0, 1), must be both evaluated with ability densities at the end states, ρbr (R,Ω the restraint in place. The strength of the harmonic potentials, restraining the bound state (left Ωc ,χ χ c or restraining the rotational/conformational free ligand (right branch) branch) around Rc ,Ω χ c can be set independently. around Ω c ,χ Ωc ,χ χ c ) of How does then the FEP or TI computed DDM dissociation free energy ∆Grd (Rc ,Ω 10 ACS Paragon Plus Environment

Page 11 of 31

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 restrained system relate to the standard dissociation free energy ∆G0d ? Or, equivalently, how Ωc ,χ χ c ) at {Rc ,Ω Ωc ,χ χ } relate to the dissociation constant does the potential of mean force w(Rc ,Ω 0

Kd /C0 = e−β ∆Gd ? We have seen that the rotational states Ω , marking the orientation of the ligand in the binding site have all equal probability of 1/8π 2 when the molecules are separated in the bulk (no matter what the conformational states of the partners are). The χ conformational states of the separated species are in general not uniformly distributed and can be thus rationalized in terms of conformational basins with uneven weights. To this end, following Ref., 21 we rewrite the Ω,χ χ ) as dissociation constant in terms of the potential of mean force w(R,Ω 1 Kd

Ω,χ χ) Ω,χ χ )J(R,Ω Ω)dRdΩ Ω]dχ χ V [ e−β w(R,Ω Θ(R,Ω R R Ω,χ χ ) J(R,Ω Ω)dRdΩ Ω]dχ χ [ e−β w(R,Ω R R −β w(R,Ω Ω,χ χ) Ω,χ χ )J(R,Ω Ω)dRdΩ Ω]dχ χ 1 [ e Θ(R,Ω R 2 χ −β w(R ,χ ) ∞ 8π χ e dχ R −β w(R ,χχ ) R −β W (R,Ω Ω,χ χ) ∞ Ω,χ χ )J(R,Ω Ω)dRdΩ Ω]dχ χ 1 e [ e Θ(R,Ω R 2 χ −β w(R ,χ ) ∞ 8π χ e dχ

R R

= limV →∞ = =

(17)

Ω) is the Jacobian of the transformation from Cartesian to the relative translational where J(R,Ω Ω and where Θ(R,Ω Ω,χ χ ) is a Heaviside step function (sometimes reorientational coordinates R,Ω ferred as the indicator function 14,21,34 ) defining the region of existence of the complex between Ω,χ χ space. In the second line of Eq. 17, in the denominator the receptor and the ligand in the R,Ω Ω,χ χ ) can be taken to be we have exploited the fact, that, as V → ∞ (infinite dilution limit), w(R,Ω χ ) everywhere so that we can perform the integral on the R,Ω Ω coordinates obtainequal to w(R∞ ,χ ing the factor 8π 2V . In the last equality we have used the definition Eq. 13 for the reversible work Ω,χ χ ) in terms of the potential of mean force (PMF) with respect to an arbitrary reference W (R,Ω conformational state. The integral in brackets in the numerator of Eq. 17 (that has the dimension of a volume and χ ) times the exponential of square radiants) can be expressed in terms of an effective volume Vb (χ Ωχ ,χ χ , corresponding to the minimum value of the the potential of mean force at the point Rχ ,Ω ligand receptor reversible work function W , given a χ conformation. This effective binding site χ )) depends in the general on the conformational state χ of the ligand and the receptor. volume Vb (χ 11 ACS Paragon Plus Environment


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 31

Ω space of the reversible work function W The same holds true for the minimum value in the R,Ω at fixed χ . Conformational states (e.g. extended vs compact configurations of a flexible ramified ligand) can be characterized by different inertia moments with an impact on the size and shape of χ )). In particular, Vb (χ χ ) may be equal to zero (or equivalently the indicator function the volume Vb (χ Ω,χ χ ) is everywhere zero) for conformational states that are incompatible with the complex. Θ(R,Ω Ωχ and of the corresponding In summary, given the conformation χ , the choices of the point Rχ ,Ω χ ) are intertwined so as to satisfy Eq. 17. effective volume Vb (χ χ ), In order to define more precisely this χ dependent effective roto-translational volume Vb (χ Ω we rewrite the integral in brackets in Eq. 17 with respect to the six external coordinates X = R,Ω χ ), i.e. using a multivariate Gaussian distribution of appropriate covariance Σ b (χ Z

Ω,χ χ) −β W (R,Ω

e

Z

Ω,χ χ )J(R,Ω Ω)dRdΩ Ω = Θ(R,Ω

Ω,χ χ) Ω)dRdΩ Ω e−β W (R,Ω J(R,Ω

χ) Vb (χ Ωχ ,χ χ) −β W (Rχ ,Ω

= e

Z

J(Rχ Ωχ )

1

T Σ −1 (χ b χ )(X−Xχ )

e− 2 (X−Xχ )

dX

(18)

so that q χ ) = J(Rχ Ωχ ) (2π)6 |Σ Σb (χ χ )| Vb (χ

(19)

If we choose a diagonal covariance matrix satisfying Eq. 19, and given that the Jacobian has Ωχ ) = JT (Rχ )JΩ (Ω Ωχ ), then the volume Vb (χ χ ) may be factored in a the factorized form J(Rχ ,Ω translational and rotational binding site volume:

χ ) = VT (χ χ )VΩ (χ χ) Vb (χ

(20)

q χ ) = JT (Rχ ) (2π)3 |Σ ΣTb (χ χ )| VT (χ q ΣΩ χ )| χ ) = JΩ (Ω Ωχ ) (2π)3 |Σ VΩ (χ b (χ

(21)

with


Page 13 of 31

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


χ ) and Σ Ω χ ) are 3 × 3 diagonal covariance matrices corresponding to the translational where Σ Tb (χ b (χ and orientational block. Using Eqs. 18 and 19, Eq. 17 may be rearranged as 1 = Kd

Z

χ ρ(χ χ) dχ

χ ) −β W (Rχ ,Ω Vb (χ Ωχ ,χ χ) e 8π 2

(22)

e−β w(R∞ ,χχ ) χ e−β w(R∞ ,χχ ) dχ

(23)

where χ) = R ρ(χ

expresses the canonical probability density for the χ conformational state of the separated species.

χ ) as a function of the conformational coordiFigure 2: Free energy of the separated state, w(R∞ ,χ χ i ) 6= 0) are indicated in green color nates. Binding conformational basins (for which Vb (χ Ω,χ χ ) is non zero only for Suppose now, as schematically depicted in Figure 2, that function Θ(R,Ω χ i of the χ coordinates, i.e. that binding can occur in Nc conformationNc disconnected domains ∆χ χ i ) around the local minimum of the ally distinct poses, each with a roto-translational volume Vb (χ Ωi ,χ χ i . Then, by appropriately selecting the volumes Vb (χ χ i ) and the points PMF at the point Ri ,Ω



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 31

Ωi ,χ χ i in the Nc conformational basins, we can rewrite the integral Eq. 22 as Ri ,Ω Nc χ) 1 V (χ Ωi ,χ χ i) χ i ) b 2i e−β W (Ri ,Ω = ∑ W (χ Kd 8π i

(24)

where we have defined the conformational weights of the Nc poses as Z

χ i) = W (χ

χ )dχ χ, ρ(χ

i = 1, 2....Nc

(25)

χi ∆χ c χ i ) 6= 1. This occurs, as depicted in Figure 2, if Note that in general we may have that ∑N i W (χ

only some of the all possible conformational basins in the domain spanned by the χ coordinates are involved in binding. Equation 24 was originally derived by Gallicchio 34 see Eqs. 19 and 20 of Ref. 34 ). If the binding involves only one conformational state or basin, i.e. if the indicator function Ω,χ χ ) is non zero around the point Rc ,Ω Ωc ,χ χ c (or, equivalently, the volume Vb (χ χ ) is non zero Θ(R,Ω χ c around the point χ c ) then we can write only within an interval ∆χ 1 Kd χ c) = where W (χ

R

χ )dχ χ χ c ρ(χ ∆χ

χ c) = W (χ

χ c ) −β [w(Rc ,Ω Vb (χ Ωc ,χ χ c )−w(R∞ ,χ χ c )] e 8π 2

(26)

χ c ) is the roto-translational binding site volume in the binding and Vb (χ

χ c ) can be therefore identified, in first instance, with the canonical conformational basin at χ c . W (χ weight in dilute solution of the binding ligand/receptor conformation for the separated species. If such single binding conformation has a low weight for the separated species, then the drug and/or the receptor experiences substantial conformational changes upon binding and the free energy χ ) or gain in the association process comes either from the volume (or entropy, vide infra) term Vb (χ Ωc ,χ χ c )−w(R∞ ,χ χ c )] from the solvent mediated enthalpic gain due to the e−β [w(Rc ,Ω term. The free energy

difference that is referred to the separated species at χ c and at χ u (where the latter is the most favorable conformational state of the separated species) is called binding reorganization energy, measuring energetic strain and entropic factors that oppose binding. 35,36


Page 15 of 31

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


Taking into account that ∆G0d = −kB T ln(Kd V 0 ), Eq. 26 can be equivalently written in terms of dissociation free energy as

∆G0d

χ c) Vb (χ χ c) Ωc ,χ χ c ) − w(R∞ ,χ χ c )] + kB T ln + kB T lnW (χ = −[w(Rc ,Ω 8π 2V 0

(27)

Going back to Eq. 14, Eq. 27 provides the sought relationship between the reversible work Ωc ,χ χ c ) and the standard dissociation free energy ∆G0d in the context of alchemical theW (Rc ,Ω ory with restraints. If we use Eq. 27 into Eq. 14 and using the definition Eq. 16, we finally find

∆G0d

Ωc ,χ χ c ) + kB T = ∆Grd (Rc ,Ω

Ωc ,χ χ c |0) ρbr (Rc ,Ω χ c) Ωc ,χ χ c |0) Vb (χ ρur (Ω ln − ln r − kB T ln r Ωc ,χ χ c |1) Ωc ,χ χ c |1) 8π 2V 0 ρb (Rc ,Ω ρu (Ω (28)

where we have defined the overall binding site volume as χ c ) = VT (χ χ c )VΩ (χ χ c )W (χ χ c) Vb (χ

(29)

Equation 28 is the central result of this paper and provides a general relation embracing (as we shall see further on) all current alchemical theories from the DAM approach, 10,11 with (apparently) no restraints, to the Deng and Roux method 19 with strong restraints. Equation 28 says that the standard or absolute dissociation free energy, ∆G0d , is given by the simulation dissociation free energy on the restrained systems, ∆Grd (Yc ), plus a composite standard state correction due to the allowance volume and probability density ratios that is strictly depending on the imposed translational and (possibly) orientational/conformational harmonic restraints. When restraints on χ c ) and the “rotational volume”, VΩ (χ χ c ), R, Ω and χ are imposed, the “translational volume” VT (χ defined in Eq. 21 are both a function of the restrained conformational coordinates. When only translational and orientational restraints are imposed, and the canonical sampling of the accessible conformational space has been attained for all λ states in either branch of the cycle, then the volume Vb = VT VΩ has, correspondingly, a translational and orientational component that are inde15 ACS Paragon Plus Environment


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 31

pendent of the conformational state in the bound or free state. By the same token, the probability density ratios in Eq. 28 become independent of the (unrestrained) conformational coordinate χ . If the harmonic restraints are imposed only on translations, i.e., Y = R, provided that the orientational sampling of the translationally restrained bound state is canonical for all λ states in either branch of the cycle, then only the volume VT and the probability density ratio for bound state alone are needed in the correction. χc Finally, we remark that while the restraining orientational and conformational parameters Ω c ,χ (if present) must be identical when decoupling the bound and free ligand, the force constant K can be selected independently. By imposing an infinitely stiff force constant when decoupling Ωc ,χ χ c |0) = the orientationally and conformationally restrained ligand in bulk, we have that ρur (Ω Ωc ,χ χ c |1). Equation 28 can hence be compactly re-written as ρur (Ω ∆G0d

= ∆Grd (Yc ) + kB T

ρbr (Yc |0) Vb (Yc ) − k T ln ln B 8π 2V 0 ρbr (Yc |1)

(30)

with the composite correction term depending only on the bound state properties and the volume Vb (Yc ) depending on the vector Yc possibly via its conformational components χ c . Deng’s and Roux’s theory: stiff restraint regime When K → ∞ also in the bound state, i.e. in the so-called stiff-spring regime, 32,37 the last logarithmic terms on the rhs of Eq. 28 are zero since the probability densities for the restrained bound and free ligand in the λ = 1 and λ = 0 states becomes identical. According to Eq. 14, the alchemically Ωc ,χ χ c ), can be thus taken to be equal to determined dissociation free energy (Eq. 16), ∆Grd (Rc ,Ω Ω,χ χ ) evaluated at the restraint value {Rc ,Ω Ωc ,χ χ c }, i.e. minus the reversible work W (R,Ω Ωc ,χ χ c ) = −W (Rc ,Ω Ωc ,χ χ c) lim ∆Grd (Rc ,Ω

Σr |→0 |Σ


(31)

Page 17 of 31

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


Using Eqs. 29 and 31, Eq. 28 can be hence rearranged as

∆G0d

χ c) χ c) VΩ (χ VT (χ χ c) Ωc ,χ χ c ) + kB T ln + kB T ln + kB T lnW (χ = −W (Rc ,Ω V0 8π 2

(32)

χ c ) and VΩ (χ χ c ) and the factor W (χ χ c ) are related to the oscillations of the R,Ω Ω,χ χ The volumes VT (χ coordinates when the ligand and the receptor are in the conformational basin with minimum χ c . As suggested in Ref. 18 these volumes should be determined in a preliminary unrestrained simulation of the complex. 38 Ω,χ χ , it is tacitly assumed that the ligand can perform When imposing strong restraints on R,Ω Ωc ,χ χ c defining the binding site and that the binding only small librations around the minimum Rc ,Ω state is characterized by a single conformational pose around χ c (i.e. Eq. 26 holds). One can see the three logarithmic terms in Eq. 32 as a translational, rotational and conformational entropy loss of the bound state, producing a penalty in the binding affinity. Hence, the more tightly the χ ), VΩ (χ χ ) and W (χ χ ) will be and the ligand is bound in the pocket, the smaller the “volumes” VT (χ larger is the entropy loss due to association. We may hence say that Eq. 32 constitutes the staΩc ,χ χ c ), tistical mechanics foundation of a sophisticated Docking approach: the function W (Rc ,Ω Ωc ,χ χ c , and of is evaluated by the alchemical decoupling of the complex, tightly restrained at Rc ,Ω χ c . In molecular Docking, the solvent averaged energetic the free ligand, tightly restrained at Ω c ,χ Ωc ,χ χ c ), is evaluated on the end-points using molecular mechanics Poissoncontribution, ∆E (Rc ,Ω Boltzmann surface area (MM/PBSA)) 39,40 or the molecular mechanics generalized Born surface area (MM/GBSA) 40–42 models, while the elusive volume entropic contributions are either evaluated using MD methodologies 43,44 or by simplified analytic estimates. 45 χ c ) ≡ VT (χ χ c )VΩ (χ χ c )W (χ χ c) It is important to stress that the size and the units of the volume Vb (χ depends on the choice of the coordinates Y that are used to define the binding site via the restraint χ c ) can be somehow estimated in independent unrestrained simulapotential. Provided that Vb (χ χ ) and VΩ (χ χ c )) and of the free ligand (needed tions of the complex 38 (needed form measuring VT (χ χ c )), Eq. 32 allows to compute the absolute dissociation free energy from the for measuring W (χ



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

difference of the decoupling free energies of the free ligand and of complex obtained by FEP or TI, Ωc ,χ χ c ligand-receptor position by a set of strong restraints where the latter is tightly kept at the Rc ,Ω of the form Eq. 1. While R and Ω are relative ligand receptor external coordinates involving the atoms of both species, the χ coordinates in general involve only atoms of the ligand, typically expressed using dihedral angles. 15,19,23 This is so, since when the ligand is annihilated in the left branch of the alchemical cycle, the free protein must not experience any restraint. If χ is chosen so as to include (or affect) also conformational coordinates of the protein, then the restraint on χ affects the motion (and the sampling) of the protein when the ligand is annihilated. In such case, an unrestrained simulation of the free protein is in principle necessary to determine the overall χ c ) for the separated species entering Eq. 32. conformational weight W (χ Eq. 32 was previously derived using a different route by Boresch et al. 15 and by Deng and Roux. 19 Eq. 32 does not explicitly depend on the imposed restraint potential, 38 hence resolving the inconsistency when one lets the volume Vr (i.e. in the stiff restraint limit) go to zero in the Gilson’s standard state correction term kB T ln(Vr /V0 ). In the strong restraint approach, however, the estimate of the dissociation free energy crucially depends on the estimate of the unknown χ c ) that can vary by several kcal mol−1 , 19 hence spanning more than binding site volume Vb (χ three orders of magnitude in the inhibition constant. Moreover, the parameters Yc in the restraint potential, Eq. 1, should be preferentially chosen such that they coincides with the corresponding mean values of the unrestrained bound state Yc = hYib , where the subscript b indicate that the mean must be taken over the canonical configurations of bound state. If any of the Yic components of the Yc vector differs from the corresponding equilibrium value hYi ib , then the alchemically decoupling system is subject to a systematic strain potential that will be reflected in the computed reversible work W (Y). Probably, as noted in Ref., 16 the major weakness in DDM with strong restraints lies in the χ c ). First choice of the restrained coordinates themselves, that impact on the size and units of Vb (χ of all, the number and the nature of the ligand and receptor conformational coordinates participating to binding is not known from the start. Secondly, whatever their choice, due to the inherent


Page 18 of 31

Page 19 of 31

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


conformational flexibility of many ligands and of virtually all receptors, these coordinates will be coupled to other ligand and receptor coordinates, so that restraining them may make hard or even prevent the sampling of configurational states that are relevant for binding. In some sense, DDM theory with strong restraints appears essentially to be based on the traditional picture of “lock and key” model, 46 with a systematic underestimation of receptor and ligand conformational reshaping upon binding (“induced fit” model 47 ).

Gilson, Bush, Given and McCammon theory: Intermediate restraint regime We now assume that we impose only translational and orientational restraints and that these restraints, when applied to the bound state (left branch of the cycle), are weak enough to allow the ligand-receptor system to canonically sample all χ conformational states that are important for binding. We then set an infinitely stiff orientational restraint in the right branch of the cycle so that Eq. 30 can be used, i.e.

∆G0d

Ωc ) + kB T = ∆Grd (Rc ,Ω

Vb ln 8π 2V 0

Ωc |0) ρbr (Rc ,Ω − kB T ln r Ωc |1) ρb (Rc ,Ω

(33)

where the binding site volume Vb = VT VΩ accounts for translational and rotational contribution. Ωc |0) ≡ The probability density of the restrained fully decoupled ligand in the binding site, ρbr (Rc ,Ω ρbr (Yc |0) is given by ρbr (Yc |0) = R

1 − 12 (Y−Yc )T Σ −1 r (Y−Yc )

e

= J(Y)dY

1 Vr

(34)

Σr ) defines the temperature dependent allowance where we have used Eqs. 1 and 3 and Vr ≡ Vr (Σ restraint volume such that Vr contains the unknown roto-translational binding site volume Vb . Ωc |1) ≡ ρbr (Yc |1), The probability density of the fully coupled restrained bound state, ρbr (Rc ,Ω can be written in terms of a product of two multivariate Gaussian distributions 48 with covariance



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 31

Σ−1 Ω} space, i.e. matrix Σ −1 = Σ −1 r defined in the Y = {R,Ω b +Σ ρbr (Yc |1) = R

e−β w(Yc ) e−β [w(Y)+V (Y−Yc )] J(Y)dY e−β w(Yc )

=

T −1

1

T Σ −1 (Y−Y ) c r

1

e−β w(Yc ) J(Yc ) e− 2 (Y−Yc ) Σ b (Y−Yc ) e− 2 (Y−Yc ) q p Σ−1 |1 +Σ r Σb| Σr +Σ Σb | (2π)6 |Σ p = = Vb Σr | Vb (2π)6 |Σ R

dY (35)

1

T −1

where we have set Σ b so that Vb e−β w(Y) J(Y)dY = e−β w(Yc ) J(Yc ) e− 2 (Y−Yc ) Σ b (Y−Yc ) with p Σb |. The effective covariance Σ b defining the volume Vb no longer depends Vb = J(Yc ) 2π 6 |Σ R

R

on the conformational states, whose contribution is supposed to be implicitly integrated away in Ω). Inserting Eqs. 35 and 34 into Eq. 33, we find the PMF w(R,Ω

∆G0d

Ωc ) + kB T = ∆Grd (Rc ,Ω

Vr ln 8π 2V 0

q Σ−1 + kB T ln |1 +Σ r Σb|

(36)

In the assumption that the last term is small and can be neglected (i.e. Vr Vb ), and factoring the restraint volume Vr in translational and orientational parts VI , ξI , then Eq. 36 is identical to the Equation proposed by Gilson. 49 Note also that, when Vr Vb , the probability density at Yc of the fully coupled bound state, ρ(Yc |1) (Eq. 35) becomes identical to that of the fully decoupled bound state, ρ(Yc |0), thus recovering the stiff regime result involving Vb only (i.e. obtaining Eq. 32 lacking the conformational term). As stated shown in Ref. 18 “[the choice of the harmonic force constant in the restraint potential] raises important questions as to how sensitive are the calculated free energies to the strength of the restraining potentials and whether an arbitrary choice of k [ the force constant of the restraint ] would yield the correct results.” With this regard, Eq. 36 establishes a “precise statistical mechanics definition of the restraining potential” 18 in terms a ratio between the known restraint volume Vr and the underlying binding site volume Vb and of its impact in the alchemically computed standard dissociation free energy. In particular, as conjectured in Ref. 14 and computationally verified in


Page 21 of 31

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


Ref., 34 if we increase Vr (or, equivalently, if we decrease the ligand concentration) by weakening the restraint, the corresponding gain in the dissociation energy due to the logarithmic term in Eq. Ωc ), provided that 36 should be exactly compensated by a decrease in the free energy ∆Grd (Rc ,Ω in the simulation of the restrained complex all states that are made accessible by the weakened restrained potential have been canonically sampled for all λ states along the alchemical path. Because of this compensation, Eq. 36 becomes in principles independent of the arbitrary choice of Vr when Vr Vb , allowing to compute the absolute dissociation free energy without knowing the underlying binding site volume Vb . On the other side, DDM with weak restraint potentials should be handled with due care by practitioners. In case of highly symmetric ligands like the case of benzene bound to T-lysozime, 19 for example, orientational restraints may prevent the sampling of the bound conformations that are defined by a mere exchange of the atom labels due to rotational operations of the symmetry group of the ligand (say σ ), underestimating the orientational volume in the bound state and hence the dissociation free energy. If the orientational restraint prevents the sampling of any of the equivalent σ = 12 states of benzene, then the free energy should be corrected by an additive volume term kB T ln 12 apparently due to “symmetry”. 50,51 If instead the restraints are engineered so that they allow the sampling of the bound states that are generated by rotations around the six-fold axis of the benzene molecule but not of those that can be generated by rotation around the 2-fold symmetry axis, then the correction factor reduces to kB T ln 2. Symmetry numbers therefore enters in the Gilson’s theory as volume correction factors. In other words, the imposed restraint may prevent the sampling of orientationally equivalent configuration, even if only a translational restraint is imposed. In the case of the benzene ligand, if the restraint on R is too tight (or the simulation too short), then the ligand may not be allowed to partially leave the binding pocket so that it can rotate around the symmetry axis in order to sample different configuration of the atom labels, i.e. in order to visit the underlying orientational volume VΩ . This kind of “symmetry” corrections should also apply to the Deng and Roux DDM variant Ωχ ) was estimated from the librations of the unrestrained ligand in with strong restraints, if Vb (Ω



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 31

just one of the equivalent orientational poses, therefore underestimating the actual orientational volume.

Jorgensen’s theory: unrestrained (DAM) regime What happens when instead we let the force constant K → 0 in Eqs. 28 and 30? The ligand in unrestrained simulations is mistakenly considered as “free to wander” while actually its concentration (or restraint volume Vbox ) is externally imposed by the periodic boundary conditions (PBC) used in the simulation. As noted in Refs., 15,23 the resulting Vbox dependent free energy should hence be related to the standard dissociation free energy as

∆G0d = ∆G(DAM) + kB T ln

Vbox V0

(37)

This result can be obtained by applying the K → 0 limit of the general Equation 30 and assuming that only a restraint on R is imposed in the left branch of the cycle, i.e.

∆G0d where VT = eβ w(Rc )

ρbr (Rc |0) VT r − kB T ln r = lim ∆Gd (Rc ) + kB T ln V0 ρb (Rc |1) Σ|→0 |Σ

−β w(R) J(R)dR Vb e

R

(38)

is the allowance oscillation volume of the COM vector

distance R in the complex irrespective of the ligand-receptor orientational and conformational Σ| → 0 , the restraint the probability density of the decoupled ligand in coordinates. In the limit |Σ bound state is uniform and is given by

lim ρbr (Rc |0) =

Σ|→0 |Σ

1 Vbox

.


(39)

Page 23 of 31

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 probability density of the coupled system at R = Rc , ρbr (Rc |1) , is simply given by e−β w(Rc ) −β w(R) J(R)dR Vbox e 1 1 1 i= h = ' T ) β w(Rc ) VT (1 + cr ) VT e VT 1 + (VboxV−V T

lim ρbr (Rc |1) = R

Σ|→0 |Σ

where we have used the fact that the constant cr =

−β w(R) J(R)dR+V −β w(R) J(R)dR = box −VT VT e Vbox e

R

R

(Vbox −VT ) β w(Rc ) e VT

'

Vbox −β ∆G0 e V0

(40)

and where

can be neglected as long as eβ ∆G0 Vbox /V 0 .

Plugging Eqs. 40 and 39 into Eq. 38, and defining ∆Gd (DAM) = lim|ΣΣ|→0 ∆Grb (Rc ) − ∆Gu , we finally recover Eq. 37. DAM theory has indeed a “restraint translational volume” that sets the ligand concentration to 1/Vbox , so that, when supplemented with the standard state correction kB T ln VVr0 , becomes a particular case of the Gilson theory, Eq. 36, with Vr = Vbox . On the other hand, Eq. 37 holds 0

if: i) the simulations are converged for all λ states; ii) the MD box volume is such that eβ ∆G Vbox /V 0 . While the second proviso is usually met in alchemical simulations for tight binding ligands (Kd < 1µM), as pointed out in Refs., 15,23 convergence of an unrestrained ligand-receptor subject only to PBC is simply out of the reach for the present computer power. 52 Nonetheless, as discussed in the introduction, the unrestrained DAM approach is still used in ligand-receptor systems, 25,27 providing in general good estimates of the dissociation free energies. In this regard, it has been recently proposed, 9,53,54 that these unrestrained DAM simulations actually express a non equilibrium measure of the dissociation free energy that is close to the real equilibrium value because of well grounded theoretical reasons stemming from nonequilibrium thermodynamics.

Conclusions In this paper, the statistical mechanics of non covalent bonding in drug-receptor systems has been revisited. It has been shown that all existing alchemical theories in binding free energy calculations can be rationalized in term of a unifying treatment encompassing the original unrestrained 23 ACS Paragon Plus Environment


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

Figure 3: Lock and key model (left) and Induced fit model (right) for the PMF in the binding site domain in the Y coordinates for the associated state DAM, 10,11 the Gilson’s restrained DDM variant 14 and the sophisticated docking approach proposed by Deng and Roux. 19 According to the general Eq. 30 relating the alchemical determined free energy difference to the standard dissociation free energy, the cited alchemical theories differ in the definition (explicit or implicit) of the binding site volume through the enforcement of a set of appropriately selected restraint potentials. The strong restraint approach 19 evaluates the potential of mean force at the bottom of the well Yc , relying on a precise knowledge of the binding pose volume Vb entering in the standard state correction. As such, the Deng and Roux alchemical theory can be considered as a sophisticated docking approach where only a single pose is assessed in the context of the traditional picture of the lock and key model (see Figure 3). The strong restraint approach, by limiting the sampling space, is in general assumed to converge more easily. 19 On the other hand, we have seen that this methodology may underestimate the effective binding site volume by preventing the sampling of conformational states that are important for binding. In the intermediate regime, conformational (and possibly orientational as well) restraints are lifted and restraints on external coordinates are weakened so that the ligand can now move in


Page 24 of 31

Page 25 of 31

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 binding pocket within a volume Vr without drifting off in the bulk. In this case, provided that Vr contains the underlying binding site volume Vb , Eq. 30 becomes equivalent to the Gilson formula Eq. 36 with Vb being replaced by the parameter Vr in the standard state correction. In the intermediate regime, the dissociation free energy estimate is in principle independent of the chosen parameter Vr due to the compensation effect on the alchemical free energy difference ∆Gb − ∆Gu when Vr is made larger, given that the simulations have converged at all λ states. When the force constant of the restraint potential is made infinitely small, only the translational restraints imposed by the PBC survives. The ligand must hence sample all states that are made accessible at the concentration 1/Vbox in all λ thermodynamic states, possibly including a fraction of unbound states at λ = 1. Alchemical decoupling becomes cumbersome and the DDM theory becomes equivalent to the DAM theory provided that the standard state correction kB T ln(Vbox /V 0 ) is used and that eβ G0 Vbox /V 0 . DAM, although it may be considered as a special case of the Gilson’s DDM theory with Vr = Vbox , is not viable in practice due to the difficulties in the sampling of ligand-receptor configurations in the entire simulation box for all λ thermodynamic states. If the restraint is weak and involves only the translational coordinate R, the ligand-receptor PMF and the underlying binding site volume Vb ≡ VT are modulated by the ro-vibrational coordinates of both the ligand and the receptor, as schematically depicted in Figure 3. As first remarked in Ref., 14 the integral defining the equilibrium constant in terms of the exponential of the PMF can be extended beyond the domain where the indicator function is equal to one with no appreciable change in Kd . It then follows that Kd depends only weakly on the specific definition of Vb , provided Ω regions of the binding site volume. 14,34,55 that such domain includes all of the important R,Ω The use of a weak restraint in the context of the alchemical theory implies weaker assumptions on the pose topology and nature, hence progressively shifting the approach towards a more realistic induced fit/conformational proofreading model in drug-receptor interaction. As remarked in the seminal papers, 14,19 the price to pay, is that the choice of a larger Vr expands the sampling of the accessible ligand-receptor configurational space in the left branch of the alchemical cycle to regions of low PMF yielding a negligible contribution to the dissociation constant. These regions,



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

however, may be the only access to the sampling alternate low free energy poses contributing to binding. This sampling expansion in DDM with weak restraints may be more relevant for low value of the alchemical coupling parameter λ (when the ligand starts to be loosely bound to the binding site), making the overall convergence of the FEP or TI simulation harder.

References (1) Muddana, H. S.; Fenley, A. T.; Mobley, D. L.; Gilson, M. K. J. Comput. Aided Mol. Des. 2014, 28, 305–317. (2) Yin, J.; Henriksen, N. M.; Slochower, D. R.; Shirts, M. R.; Chiu, M. W.; Mobley, D. L.; Gilson, M. K. J. Comput. Aided Mol. Des. 2016, 1–19. (3) Tembe, B. L.; McCammon, J. A. Comput. Chem. 1984, 8, 281–283. (4) Jorgensen, W.; Ravimohan, C. J. Chem. Phys. 1985, 83, 3050–3054. (5) Gao, J.; Kuczera, K.; Tidor, B.; Karplus, M. Science 1989, 244, 1069–1072. (6) Mobley, D. L.; Klimovich, P. V. J. Chem. Phys. 2012, 137, 230901. (7) Chodera, J.; Mobley, D.; Shirts, M.; Dixon, R.; K.Branson,; Pande, V. Curr. Opin. Struct. Biol 2011, 21, 150–160. (8) Gumbart, J. C.; Roux, B.; Chipot, C. J. Chem. Theory Comput. 2013, 9, 974–802. (9) Procacci, P. J. Mol. Graph. Model. 2017, 71, 233 – 241. (10) Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; TiradoRives, J. J. Chem. Phys. 1988, 89, 3742– 3746. (11) Jorgensen, W. L. Acc. Chem. Res. 1989, 22, 184–189. (12) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420–1426.


Page 26 of 31

Page 27 of 31

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


(13) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300–313. (14) Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. Biophys. J. 1997, 72, 1047–1069. (15) Boresch, S.; Tettinger, F.; Leitgeb, M.; Karplus, M. J. Phys. Chem. B 2003, 107, 9535–9551. (16) Mobley, D. L.; Chodera, J. D.; Dill, K. A. J. Chem. Phys. 2006, 125, 084902. (17) In the Gilson et. paper 14 there is a note added to the proofs alluding to a still unpublished work by Wang and Hermans. 20 In Ref., 20 Table 5, (where the Gilson’s correction is computed for the Lysozime-benzene complex at increasing restraint strength) raises an important question on the validity of the Gilson’s standard state correction when the restraint strength is made infinitely large. (18) Hamelberg, D.; McCammon, J. A. J. Am. Chem. Soc. 2004, 126, 7683–7689. (19) Deng, Y.; ; Roux, B. J. Chem. Theory Comput. 2006, 2, 1255–1273. (20) Hermans, J.; Wang, L. Journal of the American Chemical Society 1997, 119, 2707–2714. (21) Luo, H.; Sharp, K. Proc. Natnl. Acad. Sci. USA 2002, 99, 10399–10404. (22) Wang, J.; Deng, Y.; B.,; Roux, Biophys. J. 2006, 91, 2798–2814. (23) Deng, Y.; Roux, B. J. Phys. Chem. B 2009, 113, 2234–2246. (24) Fujitani, H.; Tanida, Y.; Ito, M.; Jayachandran, G.; Snow, C. D.; Shirts, M. R.; Sorin, E. J.; Pande, V. S. J. Chem. Phys. 2005, 123, 084108. (25) Fujitani, H.; Tanida, Y.; Matsuura, A. Phys. Rev. E 2009, 79, 021914. (26) Yamashita, T.; Ueda, A.; Mitsui, T.; Tomonaga, A.; Matsumoto, S.; Kodama, T.; Fujitani, H. Chem. Pharm. Bull. 2014, 62, 661–667. (27) Yamashita, T.; Ueda, A.; Mitsui, T.; Tomonaga, A.; Matsumoto, S.; Kodama, T.; Fujitani, H. Chem. Pharm. Bull. 2015, 63, 147–155. 27 ACS Paragon Plus Environment


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

(28) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. J. Am. Chem. Soc. 2015, 137, 2695–2703. (29) This quote (due to Paul Ehrlich, Lancet 1913, 182, 445-451 ) opened the introduction in Ref. 19 and translates as "bodies do not interact unless they are bound". (30) Lehn, J.-M. Supramolecular Chemistry: Concepts and Perspectives; Wiley-VCH: Weinheim (Bundesrepublik Deutschland), 1995. (31) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press, 1987. (32) Marsili, S.; Procacci, P. J. Phys. Chem. B 2010, 114, 2509–2516. (33) Marsili, S.; Procacci, P. J. Phys. Chem. B 2016, 120, 7037–7037. (34) Gallicchio, E.; Lapelosa, M.; Levy, R. M. J. Chem. Theory Comput. 2010, 6, 2961–2977. (35) Yang, C.-Y.; Sun, H.; Chen, J.; Nikolovska-Coleska, Z.; Wang, S. J Am Chem Soc 2009, 131, 13709–13721. (36) Gallicchio, E. Comput Mol Biosci 2012, 2, 7–22. (37) Park, S.; Schulten, K. J. Chem. Phys. 2004, 120, 5946–5961. (38) As done in Refs., 19,22,23 one could compute the volumes VT and VΩ by first turning on the corresponding restraint potentials in the fully coupled bound state by means of an additional FEP procedure, and then computing analytically the cost for turning them off in the decoupled state. As repeatedly stated in Refs., 19,23 due care must be taken that the ligand, when fully coupled, never leaves the binding site in the early stages of the additional FEP process on the restraints. 28 ACS Paragon Plus Environment

Page 28 of 31

Page 29 of 31

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


(39) Srinivasan, J.; Cheatham, T. E.; Cieplak, P.; Kollman, P. A.; Case, D. A. J. Am. Chem. Soc. 1998, 120, 9401–9409. (40) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; ; Cheatham, T. E. Acc. Chem. Res. 2000, 33, 889–897. (41) Greenidge, P. A.; Kramer, C.; Mozziconacci, J.-C.; Wolf, R. M. J. Chem. Inf. Model. 2013, 53, 201–209. (42) Gallicchio, E.; Levy, R. M. J. Comput. Chem. 2004, 25, 479–499. (43) Gao, C.; Park, M.-S.; Stern, H. A. Biophys. J. 2010, 98, 901–910. (44) Hou, T.; Wang, J.; Li, Y.; Wang, W. J. Chem. Inf. Model. 2011, 51, 69–82. (45) Procacci, P. J. Comput. Chem. 2016, 37, 1819–1827. (46) Fischer, E. Ber. Dtsch. Chem. Ges. 1894, 27, 2984–2993. (47) Koshland, D. E. Proc. Natl. Acad. Sci. USA 1958, 44, 98–104. (48) Ahrendt, P. The Multivariate Gaussian Probability Distribution. (Accessed 10 February 2017); http://www2.imm.dtu.dk/pubdb/p.php?3312. (49) Equation 36 is identical to Eq. 28 of Ref. 14 in the assumption that i) the P∆V term is imΩc ) (simulations done in the NPT ensemble); ii) the alchemical plicitely included in ∆Grd (Rc ,Ω Ωc ) is corrected for symmetrically equivalent binding poses dissociation free energy ∆Grd (Rc ,Ω whose sampling is prevented by the restraint harmonic potential. (50) Gilson, M. K.; Irikura, K. K. J. Phys. Chem. B 2010, 114, 16304–16317. (51) Gilson, M. K.; Irikura, K. K. J. Phys. Chem. B 2013, 117, 3061–3061.



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 30 of 31

(52) In fact, the time scale of a random encounter in typical MD box containing a single drugreceptor pair can be straightforwardly estimated from the mean free path,

Vbox πd 2

(with d being

the mean radius of the receptor assumed to be much larger than that of the ligand), and the diffusion coefficient of a ligand in water, 56 typically obtaining collision rates of the order of 0.1:0.01 ns−1 , i.e. a random collision every 10 to 100 ns. (53) Procacci, P. Phys. Chem. Chem. Phys. 2016, 18, 14991–15004. (54) Sandberg, R. B.; Banchelli, M.; Guardiani, C.; Menichetti, S.; Caminati, G.; Procacci, P. J. Chem. Theory Comput. 2015, 11, 423–435. (55) Tan, Z.; Xia, J.; Zhang, B. W.; Levy, R. M. J. Chem. Phys. 2016, 144, 034107. (56) Wilke, C. R.; Chang, P. AlChE J. 1955, 1, 264–270.


Page 31 Journal of 31 of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

ACS Paragon Plus Environment

Statistical Mechanics of LigandâReceptor Noncovalent ... - FloRe - UniFI

Recommend Documents

Statistical Mechanics of LigandâReceptor Noncovalent ... - FloRe - UniFI