Perturbation Potentials to Overcome Order/Disorder Transitions in Alchemical Binding Free Energy Calculations

TL;DR

This paper explores how order/disorder transitions affect alchemical free energy calculations in protein-ligand systems, proposing analytical models and strategies to improve sampling efficiency and accuracy.

Contribution

It introduces an analytical model combined with existing formalism to analyze and mitigate order/disorder transition effects in alchemical binding free energy simulations.

Findings

Order/disorder transitions hinder sampling efficiency.

Analytical models can predict and analyze these transitions.

Strategies to design better alchemical schedules are proposed.

Abstract

We investigate the role of order/disorder transitions in alchemical simulations of protein-ligand absolute binding free energies. We show, in the context of a potential of mean force description, that for a benchmarking system (the complex between the L99A mutant of T4 lysozyme and 3-iodotoluene) and for a more challenging system relevant for medicinal applications (the complex of the farnesoid X receptor and inhibitor 26 from a recent D3R challenge) that order/disorder transitions can significantly hamper Hamiltonian replica exchange sampling efficiency and slow down the rate of equilibration of binding free energy estimates. We further show that our analytical model of alchemical binding combined with the formalism developed by Straub et al. for the treatment of order/disorder transitions of molecular systems can be successfully employed to analyze the transitions and help design alchemical schedules and soft-core functions that avoid or reduce the adverse effects of rare binding/unbinding transitions. The results of this work pave the way for the application of these techniques to the alchemical estimation with explicit solvation of hydration free energies and absolute binding free energies of systems undergoing order/disorder transitions.

Tables

Table 1. Table 1: Summary of the alchemical calculations reported as part of this work

Simulation id	soft-core function	$u_{\max }$	$a$	biasing potential	number of replicas	simulation length
T4L99A/3-iodotoluene
1	TanhSC	$50$	–	Linear	$16$	$6$ ns
2	RatSC	$50$	$1/16$	Linear	$16$	$6$ ns
3	RatSC	$50$	$1/16$	ILog	$16$	$6$ ns
FXR/26
1	TanhSC	$50$	–	Linear	$24$	$6$ ns
2	RatSC	$50$	$1/16$	Linear	$24$	$6$ ns
3	RatSC	$50$	$1/16$	ILog	$24$	$6$ ns
4	TanhSC	$100$	–	Linear	$24$	$6$ ns
5	RatSC	$100$	$1/14$	Linear	$24$	$6$ ns
6	RatSC	$100$	$1/14$	ILog	$24$	$6$ ns

Table 2. Table 2: Optimized parameters for the analytical model of binding for the three complexes studied in this work

	weight	$p_b$	${\overline{u}}_b$111In kcal/mol	${\sigma }_b^a$	${ϵ}^a$	${\tilde{u}}^a$	$n_l$
T4L99A/3-iodotoluene
mode 1	$7.5\times {10}^{-3}$	$1\times {10}^{-5}$	$-11.0$	$1.95$	$20$	$-20$	$5.5$
mode 2	$0.10$	$1\times {10}^{-6}$	$-4$	$2.8$	$20$	$-4$	$5.5$
mode 3	$0.8925$	$0$	$100$	$10$	$100$	$100$	$8$
FXR/FXR-26
mode 1	$1\times {10}^{-5}$	$1.5\times {10}^{-10}$	$-28$	$2.9$	$6$	$-6$	$19$
mode 2	$\sim 1$b	$0$	$1\times {10}^5$	$1\times {10}^4$	$100$	$1\times {10}^5$	$30$

Table 3. Table 3: Number of binding and unbinding transitions

Protocol	$n_{\text{bind}}$	$n_{\text{unbind}}$
T4L99A/3-iodotoluene222$u_{\text{lower}}=-10$, $u_{\text{upper}}=25$ kcal/mol, $t>50$ ps
TanhSC	$2$	$2$
RatSC	$52$	$55$
ILog	$29$	$29$
FXR/26 ($u_{\max }=50$ kcal/mol)333$u_{\text{lower}}=-10$, $u_{\text{upper}}=25$ kcal/mol, $t>50$ ps
TanhSC	$0$	$1$
RatSC	$0$	$1$
ILog	$41$	$40$
FXR/26 ($u_{\max }=100$ kcal/mol)444$u_{\text{lower}}=-20$, $u_{\text{upper}}=50$ kcal/mol, $t>50$ ps
TanhSC	$0$	$1$
RatSC	$2$	$4$
ILog	$22$	$22$

Table 4. Table 4: Standard binding free energy estimates

Protocol	$\Delta G_{\text{site}}^{\circ }$	$t_{\text{eq}}$
T4L/3-iodotoluene
TanhSC	$-2.53\pm 0.140$	$2.0$
RatSC	$-2.26\pm 0.086$	$0.9$
ILog	$-2.32\pm 0.097$	$1.1$
FXR/26 ($u_{\max }=50$ kcal/mol)
TanhSC	$-17.4\pm 0.3$555from last 1 ns of 6 ns per replica	–
RatSC	$-12.4\pm 0.3$a	–
ILog	$-13.2\pm 0.2$	$4.0$
FXR/26 ($u_{\max }=100$ kcal/mol)
TanhSC	$-17.8\pm 0.3$a	–
RatSC	$-14.4\pm 0.3$	$4.0$
ILog	$-13.7\pm 0.2$	$3.2$

Table 5. Table 5: Alchemical schedule of the integrated logistic perturbation function for the T4L/3-iodotoluene complex

$\lambda$	${\lambda }_1$	${\lambda }_2$	$\alpha$666kcal/mol-1	$u_0$777kcal/mol	$w_0$888kcal/mol
0.000	0.000	0.000	0.010	4.000	0.000
0.067	0.000	0.067	0.167	4.000	-0.544
0.133	0.000	0.133	0.333	4.000	-0.811
0.200	0.000	0.200	0.500	4.000	-1.077
0.267	0.000	0.267	0.667	4.000	-1.344
0.333	0.000	0.333	0.833	2.504	-1.112
0.400	0.000	0.400	1.000	1.000	-0.677
0.467	0.111	0.400	1.000	1.000	-0.289
0.533	0.222	0.400	1.000	1.000	-0.178
0.600	0.333	0.400	1.000	1.000	-0.067
0.667	0.444	0.444	1.000	1.000	0.000
0.733	0.556	0.556	1.000	1.000	-0.000
0.800	0.667	0.667	1.000	1.000	0.000
0.867	0.778	0.778	1.000	1.000	0.000
0.933	0.889	0.889	1.000	1.000	0.000
1.000	1.000	1.000	1.000	1.000	0.000

Table 6. Table 6: Alchemical schedule of the integrated logistic perturbation function for the FXR/26 complex with u max = 50 subscript 𝑢 max 50 u_{\rm max}=50 kcal/mol

$\lambda$	${\lambda }_1$	${\lambda }_2$	$\alpha$999kcal/mol-1	$u_0$101010kcal/mol	$w_0$111111kcal/mol
0.000	0.000	0.000	0.350	2.000	0.000
0.043	0.000	0.056	0.350	2.000	-0.223
0.087	0.000	0.111	0.350	2.000	-0.442
0.130	0.000	0.167	0.350	2.000	-0.665
0.174	0.000	0.222	0.350	2.000	-0.884
0.217	0.000	0.278	0.350	2.000	-1.107
0.261	0.000	0.334	0.350	2.000	-1.329
0.304	0.000	0.389	0.350	2.000	-1.548
0.348	0.000	0.445	0.350	2.000	-1.771
0.391	0.000	0.500	0.350	2.000	-1.990
0.435	0.072	0.500	0.350	2.000	-1.848
0.478	0.143	0.500	0.350	2.000	-1.707
0.522	0.215	0.500	0.350	1.500	-1.314
0.565	0.286	0.500	0.350	1.000	-0.924
0.609	0.357	0.500	0.350	0.500	-0.533
0.652	0.429	0.500	0.350	0.000	-0.141
0.696	0.500	0.500	0.350	0.000	0.000
0.739	0.572	0.572	0.350	0.000	0.000
0.783	0.643	0.643	0.350	0.000	0.000
0.826	0.714	0.714	0.350	0.000	0.000
0.870	0.786	0.786	0.350	0.000	0.000
0.913	0.857	0.857	0.350	0.000	0.000
0.957	0.929	0.929	0.350	0.000	0.000
1.000	1.000	1.000	0.350	0.000	0.000

Table 7. Table 7: Alchemical schedule of the integrated logistic perturbation function for the FXR/26 complex with u max = 100 subscript 𝑢 max 100 u_{\rm max}=100 kcal/mol

$\lambda$	${\lambda }_1$	${\lambda }_2$	$\alpha$121212kcal/mol-1	$u_0$131313kcal/mol	$w_0$141414kcal/mol
0.000	0.000	0.000	0.250	5.000	0.000
0.043	0.000	0.037	0.250	5.000	-0.288
0.087	0.000	0.075	0.250	5.000	-0.583
0.130	0.000	0.112	0.250	5.000	-0.871
0.174	0.000	0.149	0.250	5.000	-1.158
0.217	0.000	0.187	0.250	5.000	-1.453
0.261	0.000	0.224	0.250	5.000	-1.741
0.304	0.000	0.261	0.250	5.000	-2.029
0.348	0.063	0.261	0.250	4.000	-1.593
0.391	0.125	0.261	0.250	3.000	-1.160
0.435	0.188	0.261	0.250	1.000	-0.463
0.478	0.250	0.261	0.250	0.000	-0.030
0.522	0.313	0.313	0.250	0.000	0.000
0.565	0.375	0.375	0.250	0.000	0.000
0.609	0.438	0.438	0.250	0.000	0.000
0.652	0.500	0.500	0.250	0.000	0.000
0.696	0.563	0.563	0.250	0.000	0.000
0.739	0.625	0.625	0.250	0.000	0.000
0.783	0.688	0.688	0.250	0.000	0.000
0.826	0.750	0.750	0.250	0.000	0.000
0.870	0.813	0.813	0.250	0.000	0.000
0.913	0.875	0.875	0.250	0.000	0.000
0.957	0.938	0.938	0.250	0.000	0.000
1.000	1.000	1.000	0.250	0.000	0.000

Equations

$$\Delta G_b^{\circ }=\Delta G_{\text{site}}^{\circ }+\Delta G_b$$

$$\Delta G_{\text{site}}^{\circ }=-\frac{1}{\beta }\ln C^{\circ }{V}_{\text{site}}$$

$$U_{\lambda }(x)=U_0(x)+W_{\lambda }(u)$$

$$u(x)=U_1(x)-U_0(x)$$

$$p_0(u)=\sum _i{c}_i{p}_{0i}(u)$$

$$p_{0i}(u)=p_{bi}g(u;{\bar{u}}_{bi},{\sigma }_{bi})$$

$$+(1-p_{bi}){\int }_0^{+\infty }{p}_{WCA}(u^{\prime };n_{li},{ϵ}_i,{\tilde{u}}_i)g(u-u^{\prime };{\bar{u}}_{bi},{\sigma }_{bi})d{u}^{\prime }$$

$$p_{WCA}(u;n_l,ϵ,\tilde{u})=n_l{\left[1-\frac{(1+x_C{)}^{1/2}}{(1+x{)}^{1/2}}\right]}^{n_l-1}\frac{H(u)}{4{ϵ}_{LJ}}\frac{(1+x_C{)}^{1/2}}{x(1+x{)}^{3/2}},$$

$$p_{\lambda }(u)=\frac{1}{K(\lambda )}{p}_0(u)\exp \left[-\beta W_{\lambda }(u)\right]$$

$$K(\lambda )={\int }_{-\infty }^{+\infty }{p}_0(u)\exp \left[-\beta W_{\lambda }(u)\right]=⟨\exp \left[-\beta W_{\lambda }(u)\right]{⟩}_{\lambda =0}$$

$$\Delta G_b(\lambda )=-\frac{1}{\beta }\ln K(\lambda )$$

$${\lambda }_0(u)-\frac{\partial W_{\lambda }(u)}{\partial u}=0$$

$${\lambda }_0(u)\equiv \frac{1}{\beta }\frac{\partial \ln p_0(u)}{\partial u}$$

$$\frac{d\Delta G_b(\lambda )}{d\lambda }=⟨u{⟩}_{\lambda }.$$

$$u_{\mathrm{sc}}(u)=\{\begin{array}{ll}u& u\le 0\\ u_{\max }{f}_{\text{sc}}(y)& u>0\end{array}$$

$$f_{\text{sc}}(y)=\tanh (y)$$

$$f_{\text{sc}}(y)=\frac{z^a-1}{z^a+1},$$

$$p_0^{\mathrm{sc}}(u_{\mathrm{sc}})=p_0[u(u_{\mathrm{sc}})]u^{\prime }(u_{\mathrm{sc}})$$

$${\lambda }_0^{\mathrm{sc}}(u_{\mathrm{sc}})=\frac{1}{\beta }\frac{\partial \ln p_0^{\mathrm{sc}}(u_{\mathrm{sc}})}{\partial u_{\mathrm{sc}}}=\frac{1}{\beta }\left[{\lambda }_0(u)u^{\prime }(u_{\mathrm{sc}})+\frac{u^{\prime \prime }(u_{\mathrm{sc}})}{u^{\prime }(u_{\mathrm{sc}})}\right]$$

$$W_{\lambda }(u)=\frac{{\lambda }_2-{\lambda }_1}{\alpha }\ln \left[1+e^{-\alpha (u-u_0)}\right]+{\lambda }_{2}u+w_0$$

$$\frac{\partial W_{\lambda }(u)}{\partial u}=\frac{{\lambda }_2-{\lambda }_1}{1+e^{-\alpha (u-u_0)}}+{\lambda }_1$$

$$\frac{1}{T_S(E)}\equiv \frac{\partial S(E)}{\partial E}=k_B\frac{\partial W_T(E)}{\partial E},$$

$$E$$

$$\Omega (E)$$

$$\frac{1}{k_{B}T}$$

$$\frac{1}{k_B{T}_S(E)}$$

$$W_T(E)$$

$$u_C(r)=\{\begin{array}{ll}{u}_{\mathrm{LJ}}(r)-\tilde{u},& r<r_C\\ 0& r\ge r_C,\end{array}$$

$$u_B(r)=\{\begin{array}{ll}\tilde{u},& r<r_C\\ u_{\mathrm{LJ}}(r)& r\ge r_C,\end{array}$$

$$\tilde{u}=u_{LJ}(r_C)>-ϵ$$

$$P_C(u_C)=H(u_C)\frac{V_C-V(u_C)}{V_C}$$

$$r(u_C)=\frac{r_0}{(1+x_C{)}^{1/6}};u_C\ge 0$$

$$x_C=\sqrt{\frac{u_C}{ϵ}+\frac{\tilde{u}}{ϵ}+1}$$

Full text

I Introduction

The accurate modeling of molecular recognition remains one of the most challenging problem in computational molecular biophysics.GallicchioSAMPL4 ; gaieb2018d3r Molecular simulations of molecular binding are limited by the incomplete description of the chemical system,Mobley:Dill:review:2009 ; pal2016SAMPL5 ; henderson2018ligand ; cruzeiro2018redox by the quality of energy functionsharder2015opls3 ; Albaugh2016advanced ; wang2017building ; roos2019opls3e as well as, and probably to a less understood extent, by insufficient conformational sampling.Lapelosa2011 ; Mobley2012 ; bodnarchuk2014strategies ; procacci2019solvation

One of the primary applications of computational models of molecular binding is to provide an estimate of the standard free energy of binding, $\Delta G_{b}^{\circ}$, or, equivalently, the equilibrium constant, $K_{b}$, for the association equilibrium $R+L\leftrightharpoons RL$, between two molecules $R$ and $L$. For example, the binding of a drug molecule to a protein receptor.Jorgensen2009 While many strategies have been proposed,michel2010prediction ; Gallicchio2011adv ; limongelli2013funnel ; deng2018comparing alchemical relative binding free energy methodslybrand:mccammon1986 are emerging as the leading approaches in pharmaceuticals applications.Shirts:Mobley:Chodera:2007:review ; Wang2012 ; homeyer2014binding ; abel2017advancing ; lee2018gpu ; Zou2019chemarxiv In the latter context, in particular, there have been significant advances in controlling errors related to conformational sampling.wang2011replica ; hauser2018predicting

We focus here on strategies to improve the estimation of absolute binding free energies by alchemical methods.Jorgensen:Buckner:Boudon:Rives:88 ; Gilson:Given:Bush:McCammon:97 ; Deng2009 ; Chodera:Mobley:cosb2011 ; kilburg2018assessment While more challenging to obtain than relative binding free energies, absolute binding free energies provide a more stringent assessment of protocols and force fieldsmobley2017predicting and are better suited for ranking dissimilar compounds in virtual screening applications,GallicchioSAMPL4 ; bfzhang2016rnaseh as well as for investigating binding specificity.deng2017resolving

In this work, we illustrate how, analogously to conventional chemical systems as a function of temperature, alchemical systems can undergo order/disorder phase transitions along the alchemical path, and that these cause entropic bottlenecks which hinder the equilibration and convergence of binding free energy estimates. We analyze this phenomenon using the analytic theory of alchemical binding we recently proposedkilburg2018analytical and the formalism developed by Straub et al.kim2010generalized ; lu2013order for the modeling of conventional phase transitions. We then use this knowledge to design novel alchemical perturbation potentialsMcCammon:Straatsma:92 and soft-core functionspitera2002comparison which eliminate or reduce the adverse effects of the transitions. This advance leads to faster equilibration and to more robust binding free energy estimates for challenging systems hard to treat with conventional protocols.

In this context, an order/disorder transition occurs at a critical value of the alchemical progress parameter at which two conformational states are in equilibrium even though one of them has a much weaker ligand-receptor interaction energy than the other. The equilibrium is established because the high energy state is entropically favored relative to the low energy state to the same extent that it is disfavored energetically–again, in analogy with, for example, a gas being in equilibrium with its liquid even though the molar enthalpy of the gas is much less favorable than that of the liquid. Similarly, during binding free energy alchemical calculations, weakly coupled states are created in which ligand and receptor have much more conformational freedom and conformational entropy than more ordered coupled states. When it is not coupled to the receptor, the ligand can freely rotate and translate within the binding site region. Similarly, protein sidechains can experience a wide range of conformations when they are not forming interactions with the ligand. To transition from the unbound state to the bound state, the complex has to go through a tight entropic bottleneck related to the small likelihood to find a bound pose in the absence of ligand-receptor interactions; a frustrated process not unlike those of protein foldingAlexeiV.Finkelstein2002 ; Zheng2007 and crystallization.yu2014order

The equilibrium between ordered and disordered states is manifested in bimodal binding energy distributions with maxima separated by wide and poorly sampled binding energy gaps.Lapelosa2011 The binding free energy estimate is heavily influenced by the relative weight of competing modes.kilburg2018analytical However, relative populations of ordered and disordered states can not be easily established because of the rare crossing events between them. The novel protocols we develop in this work are designed to reduce the binding energy gap and increase the number of binding/unbinding crossing events.

The paper is organized as follows. We first review the alchemical methodology and the analytic theory we employ. We then discuss methods to monitor and characterize the order/disorder transitions and ways to tune soft-core functions and alchemical perturbation functions to avoid them. We then illustrate the methods on two molecular complexes which display order/disorder biphasic behavior. The first is the complex between the L99A mutant of T4 lysozyme and 3-iodotoluene, a well known benchmarking system.gill2018binding The second is the more challenging complex between the farnesoid X receptor and an inhibitor from a recent D3R grand challenge,gaieb2018d3r which is more representative of typical medicinal applications. We show that in each case, order/disorder can be addressed leading to more efficient calculations and more robust binding free energy estimates.

II Theory and Methods

II.1 Alchemical Transformations for Binding Free Energy Estimation

We adopt a well-known statistical mechanics formulation of biomolecular non-covalent binding in which the standard free energy of binding is written as:Gilson:Given:Bush:McCammon:97 ; Boresch:Karplus:2003 ; Gallicchio2011adv

[TABLE]

where

[TABLE]

is the reversible work for transferring a ligand molecule from an ideal solution at the standard concentration $C^{\circ}=1$M of ligand molecules, to a binding site region of volume $V_{\text{site}}$, and $\Delta G_{b}$ is the excess component corresponding to the desolvation process and the establishing of ligand-receptor interactions. In this work, the excess component is computed by means of an alchemical free energy perturbation schedule based on a potential energy function $U_{\lambda}(x)$, parametric on the alchemical progress parameter $\lambda$, which interpolates between the unbound state of the complex, described by the potential function $U_{0}(x)$, and that corresponding to the bound state, $U_{1}(x)$. When the solvent is modeled explicitly, a double-decouplng processGilson:Given:Bush:McCammon:97 ; Chodera:Mobley:cosb2011 is used involving two alchemical legs, the first in which the ligand is decoupled from the solvent to reach a “vacuum” intermediate state, and a second leg in which the ligand is placed into the binding site and then coupled to the receptor.

Here we use a solvent potential of mean force formulation which allows to transfer the ligand directly from solution to the receptor binding site in a single alchemical process.Gallicchio2010 ; Gallicchio2011adv ; kilburg2018assessment In this formulation, in which the degrees of freedom of the solvent are averaged out, the effective potential energy function includes a standard molecular mechanics component describing covalent and non-bonded interactions plus an implicit solvation component to model hydration effects.Gallicchio:Levy:2004 In our implementation, the alchemical potential energy is expressed as

[TABLE]

where $x$ represents the set of atomic coordinates of the molecular complex,

[TABLE]

is the binding energy function, $U_{0}(x)$ is the effective potential energy of the uncoupled state, $U_{1}(x)$ is the effective potential energy of the coupled state, and $W_{\lambda}(u)$ is the alchemical perturbation function, which varies parametrically with $\lambda$. As defined in Eq. (4), the binding energy function is defined as the potential energy difference between the coupled and uncoupled states of the complex in conformation $x$.Gallicchio2011adv

In order to reproduce the physical coupled and uncoupled states of the complex at the beginning and end states of the alchemical transformation, it is necessary that the alchemical perturbation function is defined such that $W_{0}(u)=0$ and $W_{1}(u)=u$ at $\lambda=0$ and $\lambda=1$, respectively. The linear function $W_{\lambda}(u)=\lambda u$ satisfies this requirement and is the standard choice for the alchemical perturbation function.Su:Gallicchio:Levy:2007 To obtain the binding free energy, set of samples of the binding energies, $u_{i}$, are collected during molecular dynamics simulations performed at a sequence of $\lambda$ values between [math] and $1$. The excess free energy profile as a function of $\lambda$, $\Delta G_{b}(\lambda)$, is obtained by multi-state reweightingShirts2008a using the UWHAM method.Tan2012 The excess free energy of binding in Eq. (1) is by definition the value of free energy profile at $\lambda=1$, $\Delta G_{b}=\Delta G_{b}(1)$.

The binding-energy representation of the perturbation energy used in this work has the advantage that an analytic theory is available to describe the statistics of the binding energy function and, by variable transformation, of any other quantity, such as the alchemical perturbation energy, which depends on it (see below).kilburg2018analytical In the following, we exploit the theory to analyze order/disorder transitions along the alchemical path and derive soft-core functions and alchemical perturbation functions with superior replica exchange efficiency the conventional linear perturbation function $W_{\lambda}(u)=\lambda u$.

II.2 Analytic Theory of Alchemical Molecular Binding

The theory of alchemical binding we recently developedkilburg2018analytical provides analytic expressions for the probability densities of the binding energy and of the binding free energy profile as a function of the alchemical progress parameter $\lambda$ for any binding energy-based alchemical perturbation potential $W_{\lambda}(u)$. The parameters of the model are obtained by fitting the prediction of the model to the binding energy probability distributions extracted from atomistic simulations of the complex.

Briefly (see reference kilburg2018analytical, for the full derivation), the central quantity of the model is $p_{0}(u)$, the probability density of the binding energy $u$ in the uncoupled ensemble ($\lambda=0$ in this context).Gallicchio2010 ; Gallicchio2011adv The model is based on the assumption that the statistics of the random variable $u$ is the superposition of two processes, one that describes the sum of many “soft” background ligand-receptor interactions and that follows central limit statistics, and another process that describes “hard” atomic collisions and that follows max statistics. The probability density $p_{0}(u)$ is expressed as the superposition of probability densities of a small number of binding modes

[TABLE]

where $c_{i}$ are adjustable weights summing to $1$ and $p_{0i}(u)$ is the probability density corresponding to mode $i$ described analytically as (see reference kilburg2018analytical, and appendix A.1 for the derivation):

[TABLE]

where $g(u;\bar{u},\sigma)$ is the normalized Gaussian density function of mean $\bar{u}$ and standard deviation $\sigma$ and

[TABLE]

where $x=\sqrt{u/\epsilon+\tilde{u}/\epsilon+1}$ and $x_{C}=\sqrt{\tilde{u}/\epsilon+1}$. The model for each mode $i$ depends on a number of adjustable parameters corresponding to the following physical quantitieskilburg2018analytical :

• •

$c_{i}$: relative population of binding mode $i$

• •

$p_{bi}$: probability that no atomic clashes occur while in binding mode $i$

• •

$\bar{u}_{bi}$: the average background interaction energy of binding mode $i$

• •

$\sigma_{bi}$: the standard deviation of background interaction energy of binding mode $i$

• •

$n_{li}$: the effective number of statistical uncorrelated atoms of the ligand in binding mode $i$

• •

$\epsilon_{i}$: the effective $\epsilon$ parameter of an hypothetical Lennard-Jones interaction energy potential describing the receptor-ligand interaction in binding mode $i$

• •

$\tilde{u}_{i}$: the binding energy value above which the collisional energy is not zero in binding mode $i$

The parameters above, together with the weights $c_{i}$, are varied to fit the binding energy distributions obtained from numerical simulationskilburg2018analytical (see Fig. 4 and Table 2 for examples).

All other quantities of the alchemical transformation can be obtained from $p_{0}(u)$.Gallicchio2010 ; Gallicchio2011adv In particular, given $p_{0}(u)$, the probability density for the binding energy $u$ for the state with perturbation potential $W_{\lambda}(u)$ is

[TABLE]

where $\beta=1/k_{B}T$,

[TABLE]

is the excess component of the equilibrium constant for binding and

[TABLE]

is the corresponding excess binding free energy profile. Note that for a linear perturbation potential, $W_{\lambda}(u)=\lambda u$, Eqs. (9) and (10) state that the binding free energy profile is related to the double-sided Laplace transform of $p_{\lambda}(u)$.

II.3 Order/Disorder Transitions along Alchemical Transformation

In this context, an order/disorder phase transition occurs whenever two conformational states with significantly different average energies and entropies are in equilibrium. In these circumstances, transfer of population from one state to another takes place in a narrow range of the controlling thermodynamic parameter, typically temperature. For example, for the two-level system whereby one state is energetically favored, and the other is entropically favored, the width of the temperature transition $\delta T$ becomes increasingly narrow as the energy gap $\Delta E$ increases ($\delta T=4k_{B}T^{2}_{m}/\Delta E$ as measured by the derivative of the population of the upper level with respect to temperature at the midpoint temperature $T_{m}$).zuckerman2010statistical In general, the hallmark of an order/disorder transition is the presence of a bimodal energy distribution with a sudden transfer of population from the low energy phase to the higher energy phase as the temperature is increased. The nature of temperature-activated order/disorder transitions and the conformational sampling bottlenecks they cause has been thoroughly investigated by John Straub and collaborators in a series of publications.Kim2010b ; kim2010generalized ; lu2013order ; lu2014investigating

In this work, we extend the analysis conducted by Straub et al. to the case of order/disorder transitions occurring along the alchemical path for binding. Because there are many more ways that ligand atoms clash with receptor atoms than configurations free of clashes, at small $\lambda$ weakly coupled configurations are entropically favored relative to coupled states. Conversely, fully coupled configurations, which are free of severe clashes and are stabilized by favorable intermolecular interactions, predominate at large values of $\lambda$. Analogously to the temperature-driven transitions, coupled and weakly coupled states separated by a sizeable energy gap, are in equilibrium at a critical value of the alchemical parameter $\lambda$. The occurrence of an order/disorder equilibrium is evident from the presence of distinct modes of the binding energy distribution separated by a large binding energy gap. Hence, as in the case of variations in temperature, we identify order/disorder transitions by looking for $\lambda$-states with binding energy probability densities $p_{\lambda}(u)$ with two or more modes.

By setting to zero the derivative of the logarithm of Eq. (8), we find that the stationary points of $p_{\lambda}(u)$ obey the relationship

[TABLE]

where we have introduced the $\lambda$-function

[TABLE]

of the system. As defined, the $\lambda$-function is analogous to the energy-dependent statistical temperature function $T_{S}(E)$ in the analysis of Straub et al.Kim2010b for the canonical ensemble. Like the density of states $\Omega(E)$ and the statistical temperature $T_{S}(E)$ in the canonical ensemble, and $p_{0}(u)$ in the alchemical ensemble considered here, $\lambda_{0}(u)$ is a fundamental property of the molecular system and is independent of $\lambda$ and of the specific alchemical perturbation. Since in this model $p_{0}(u)$ is expressed analytically, the $\lambda$-function is also obtained in analytical form using Eq. (12).

Eq. (11) leads to the graphical construction illustrated in Fig. 1 to locate the maxima and minima of $p_{\lambda}(u)$.Kim2010b In the case of a linear perturbation ($W_{\lambda}(u)=\lambda u$) the stationary points of the distribution of binding energies $u$ at $\lambda$ occur when $\lambda_{0}(u)=\lambda$, that is when the $\lambda$-function intersects a horizontal line corresponding to the set value of $\lambda$ (Figure 1). For a general perturbation function $W_{\lambda}(u)$, the stationary points occur when $\lambda_{0}(u)$ intersects the function $\partial W_{\lambda}(u)/\partial u$ as in Eq. (11). As illustrated in Figure 1, a order/disorder transition occurring with the linear perturbation function can be avoided by using a suitable, non-linear, perturbation function.Kim2010b

The graphical construction in Figure 1 also easily indicates the regions of binding energies where $p_{\lambda}(u)$ in increasing or decreasing. When $\lambda_{0}(u)>\partial W_{\lambda}(u)/\partial u>0$, that is when the curve representing the derivative of the perturbation function is below the $\lambda$-function, $p_{\lambda}(u)$ has a positive derivative and increases with increasing $u$. Conversely, when $\partial W_{\lambda}(u)/\partial u$ is above the $\lambda$-function, $p_{\lambda}(u)$ is a decreasing. These behaviors are confirmed by our results as shown in Figure 5.

II.4 Soft-Core Binding Energy Functions to Reduce Order/Disorder Transitions

The alchemical potential of Eq. (3) with the energy function $u(x)$ defined by Eq. (4) leads to an unstable free energy estimation near the decoupled state.Steinbrecher2011 This issue, which is generally known as the “end-point catastrophe”,pitera2002comparison is due to the singularity of the derivative of the free energy profile near the decoupled state.simonson1993free For a linear bias, $W_{\lambda}(u)=\lambda u$, it is straightforward to show from Eq. (3) that the derivative of the binding free energy profile is the average binding energy

[TABLE]

Near the decoupled state, when ligand and receptor atoms interact weakly, it is very likely to encounter configurations in which ligand atoms and receptor atoms clash, leading to very unfavorable binding energies. The average binding energy and first derivative of the binding free energy profile diverge near the decoupled state at $\lambda=0$. (This behavior can also be confirmed by considering the first moment of $p_{0}(u)$ using Eqs. (5) and (6).)

To address the end-point singularity, we introduce a soft-core binding energy functions that caps the binding energy to some maximum value $u_{{\rm max}}$:Buelens2012

[TABLE]

where $u$ is the binding energy function, $y=u/u_{{\rm max}}$, and $f_{\text{sc}}(y)$ is a function that smoothly goes from zero at $y=0$ to one as $y$ goes to infinity (see below). The soft-core binding energy function $u_{{\rm sc}}=u_{{\rm sc}}(u)$ can be interpreted as a map from the unbound domain $[-\infty,+\infty]$ of $u$ to the $[-\infty,u_{{\rm max}}]$ domain of $u_{{\rm sc}}$, which is bounded from above at $u_{\max}$. The end states and their free energy difference are virtually unaffected by the replacement of $u$ with $u_{{\rm sc}}$. At the coupled state the binding energy distribution lays almost exclusively at negative binding energy values where $u$ and $u_{{\rm sc}}$ are the same [see Eq. (14)]. At the uncoupled state ($\lambda=0$) on the other hand, the biasing potential is zero regardless of the form of the binding energy function.

We consider two choices for $f_{\text{sc}}(y)$: a hyperbolic tangent

[TABLE]

and a rational function designed to lessen the strength of order/disorder transitions

[TABLE]

where $z=1+2y/a+2(y/a)^{2}$ and $a$ is an adjustable dimensionless exponent. Both forms of $f_{\text{sc}}(y)$ above are invertible and lead to a $C(2)$-smooth soft-core binding energy function suitable for molecular dynamics applications when included in Eq. (14).

The soft-core binding energy function of Eqs. (14)–(16) replaces the binding energy function $u$ in the binding energy-based alchemical potential of Eq. (3). As an illustration, Fig. (2) shows the effect of the soft-core function on the binding energy between two particles interacting by a Lennard-Jones interaction. While both soft-core functions cap the binding energy function to $u_{\rm max}$, the rational soft-core function leads to a smaller high-energy plateau region which, as shown below, reduces the strength of order/disorder transitions near the decoupled states (see below).

II.4.1 Probability Density and $\lambda$-Function of Soft-Core Binding Energy

The soft-core binding energy function (14) is interpreted as a redefinition of the binding energy of the system. Application of the alchemical theory for the soft-core definition of the binding energy hinges on knowledge the probability density, $p_{0}^{{\rm sc}}(u_{{\rm sc}})$, of $u_{{\rm sc}}$ at $\lambda=0$. This can be obtained from the probability density, $p_{0}(u)$, of $u$ using the variable transformation formula:

[TABLE]

where $u=u(u_{{\rm sc}})$ is the inverse of the soft-core function [Eq. (14)]. Since $p_{0}(u)$ is available in analytical form [Eq. (5)] and Eqs. (15) and (16) are analytically invertible and differentiable, Eq. (17) provides an analytical expression of the probability density of the soft-core binding energy at $\lambda=0$.

To $\lambda$-function for the soft-core binding energy is obtained by differentiating the logarithm of Eq. (17) with respect to $u_{sc}$, yielding

[TABLE]

where $\lambda_{0}(u)$ is the $\lambda$-function of the binding energy without soft-core [Eq. (12)] evaluated at $u=u(u_{{\rm sc}})$.

It can be shown from Eq. (18), and Eqs. (14), (15), and (16), that the $\lambda$-functions with both of the soft-core functions we are considering diverges to $+\infty$ at $u_{\rm sc}=u_{\rm max}$ where the first derivative of the soft-core function is zero. Interestingly, the characteristics of the singularity are system-independent (the term $\lambda_{0}(u)u^{\prime}(u_{{\rm sc}})$ is zero at the singularity) and they depend only on the nature of the soft-core function. As illustrated below, the divergence of $\lambda^{\rm sc}_{0}(u_{{\rm sc}})$, implies that all distribution functions of the soft-core binding energy, regardless of the form of the perturbation potential and of the alchemical schedule, will present a stationary point (a minimum, in fact) near $u_{\rm sc}=u_{\rm max}$. Because, $\partial w_{\lambda}(u)/\partial u$ is finite for any well-behaved perturbation function, near $p_{\lambda}(u)$ is necessarily an increasing function near $u_{\rm sc}=u_{\rm max}$. However, in practice the upward trend of $p_{\lambda}(u)$ may not be observed in alchemical states, such as bound states, for which the population near $u_{\rm max}$ is negligibly small.

II.5 Alchemical Perturbations Functions to Reduce Order/Disorder Transitions

As shown below, while suitable soft-core functions can reduce the binding energy gap across an order/disorder transition, they can also create order/disorder sampling bottlenecks elsewhere along the alchemical path (see for example Figure 6). We found that combining a soft-core function with new kinds of perturbation potentials to be particularly beneficial. One of the most effective perturbation potentials that we identified is

[TABLE]

which we named the integrated logistic biasing function. The parameters $\lambda_{2}$, $\lambda_{1}$, $\alpha$, $u_{0}$, and $w_{0}$ are functions of $\lambda$. The name of this function comes from the fact that its derivative is the logistic function (also know as Fermi’s function)

[TABLE]

which is sketched out in Figures 1 and 5(A).

The parameter $\lambda_{1}$ is the height of the horizontal branch at the low binding energy end, and $\lambda_{2}$ is the height at the high binding energy end. The parameter $u_{0}$ controls the position of the switch from $\lambda_{1}$ to $\lambda_{2}$, and $\alpha$ controls the range of the switch. The parameter $w_{0}$ is an overall energy offset. As illustrated below, bimodal behavior can be avoided by properly tuning the parameters of the integrated logistic function in regions of the alchemical path affected by order/disorder transitions. Conversely, the integrated logistic biasing function behaves as a linear biasing function away from the transition region.

II.6 Hamiltonian Replica Exchange Conformational Sampling

In this work, we employ the Hamiltonian Replica Exchange algorithmSugita2000 ; Felts:Harano:Gallicchio:Levy:2004 ; Ravindranathan:Gallicchio:Levy:2006 ; Okumura2010 in alchemical space to accelerate conformational sampling.Woods2003 ; Rick2006 ; Gallicchio2010 In the context of the alchemical potential energy function (3), Hamiltonian Replica Exchange consists of performing MD of multiple replicas of the system each assigned a value of the $\lambda$ parameter. Periodically, the assignment of replicas to $\lambda$ states is varied in such a way so as to preserve a canonical distribution of conformations at each $\lambda$. The algorithm allows each replica to explore a wide range of $\lambda$-states, from coupled to uncoupled, thereby accelerating the sampling of receptor-ligand intermolecular degrees of freedom. For the calculations reported here, we have employed the asynchronous implementation of Replica Exchange (ASyncRE)gallicchio2015asynchronous with the Gibbs Independence Sampling algorithmchodera2011replica for state reassignments (see Appendix A.2).

The occurrence of an order/disorder phase transition has been shown to limit the rate at which replicas diffuse in thermodynamic space,Kim2010b and hinder the ability of replica exchange to accelerate conformational sampling. This is because, at a phase transition, nearby thermodynamic states can be separated by a large energy gap. In the alchemical case, the probability that two replicas exchange their $\lambda$-states is large when the two replicas have similar binding energies, and it decreases rapidly when they are separated by a large binding energy gap. This effect is best appreciated in the case of a linear perturbation potential $W_{\lambda}(u)=\lambda u$ for two replicas on opposite sides of the order/disorder transition. The replica in the ordered state with more favorable binding energy is preferentially at the larger value of $\lambda$ whereas the replica in the disordered state is more likely to have higher binding energy at a smaller value of $\lambda$. In this case the probability of exchange, which is proportional to the factor $\exp(\Delta\lambda\Delta u)$, where $\Delta u>0$ is the difference in binding energies and $\Delta\lambda<0$ the difference in $\lambda$ values between the two replicas, decreases exponentially with increasing binding energy separation. In this work we show that Hamiltonian RE efficiency can be improved by using an alchemical perturbation potential, such as the one proposed above, that removes or softens order/disorder transitions.

II.7 Computational Details

The molecular structures of the L99A T4 lysozyme receptor (PDB ID 4W53) and of the Farnesoid X Receptor (FXR), see Figure 3, were prepared from their crystal structuresgaieb2018d3r ; merski2015homologous using the protein preparation wizard of Maestro (Schrodinger Inc) with default settings. Residues 1 through 71 of the T4 lysozyme receptor, which do not participate in ligand binding, were removed. The positions of the C$\alpha$ atoms of the receptors were loosely restrained using a flat-bottom harmonic potential with a tolerance of $1.5$ Å. 3-iodotoluene was placed in the binding site of T4L by superimposing it on the structure of bound toluene. The FXR-26 inhibitor bound to FXR was placed based on the available crystal structure.gaieb2018d3r Both ligands were prepared using the LigPrep facility of Maestro at pH 7. The single-decoupling binding free energy simulations were set up using the Single Decoupling Method (SDM) workflow (github.com/egallicc/openmm_sdm_workflow). OPLS-AA 2005 force field parametersKaminski:2001 ; Banks:Gallicchio:Levy:2005 were assigned using Desmond. The binding site volume for the T4L system is defined as any conformation in which the ligand center of mass is within $2.5$ Å of the center of mass of the $\mathrm{C_{\alpha}}$ atoms of residues 79, 84, 88, 91, 96, 104, 112, 113, 122, 133, and 150 of the T4L receptor. The binding site volume for the FXR system is similarly defined using the center of mass of the $\mathrm{C_{\alpha}}$ atoms of residues 273, 277, 291, 333, 336, 340, 356, and 369 of the FXR receptor. In both systems, the ligand was sequestered within the binding site by means of a flat-bottom harmonic potential with a force constant of $25$ kcal/mol Å2 applied to atoms with distances greater than $2.5$ Å. With these settings, the value of $\Delta G^{\circ}_{\rm site}$ in Eq. (2) is $1.97$ kcal/mol.

SDM alchemical calculations employed the OpenMMeastman2017openmm MD engine with the AGBNP (github.com/egallicc/openmm_agbnp_plugin) and SDM integrator plugins (github.com/rajatkrpal/openmm_sdm_plugin.git) using the OpenCL platform. The ASyncRE software,gallicchio2015asynchronous customized for OpenMM and SDM (github.com/baofzhang/async_re-openmm.git), was used for the Hamiltonian Replica Exchange in $\lambda$ space. We used 16 and 24 replicas, respectively, for the T4L and FXR complexes, set at equally spaced values of $\lambda$ between [math] and $1$. The $\lambda$-dependent parameters, listed in Tables 5, 6, and 7, of the integrated logistic schedule were chosen so as to avoid the order/disorder transitions near $u=0$, as described in the Results section. Molecular dynamics runs were conducted for a minimum of 6 ns per replica with a 1 fs time-step at $300$ K, exchanging approximately every 10 to 20 ps. A Langevin thermostat at $300$ K with a relaxation time constant of $2$ ps was used. Binding energy samples and trajectory frames were recorded every 5 ps. The calculations were performed on a farm of GPU servers at Brooklyn College and on the XSEDE Comet GPU HPC cluster using a mix of 780 Ti, 2080, Titan Xp, K80, and P100 NVIDIA GPUs.

III Results

The adverse effects of order/disorder transitions in alchemical calculations, and the ability of suitably crafted soft-core functions and perturbation potentials to avoid or ameliorate them, are illustrated here for two systems: the complex of L99A T4 lysozyme with 3-iodotoluene (T4L99A/3-iodotoluene) and the complex between the farnesoid X receptor with inhibitor 26 (FXR/26)gaieb2018d3r (Figure 3). The T4 lysozyme receptor is often used to test alchemical and conformational sampling methods.Boyce2009 ; Gallicchio2010 ; gill2018binding The FXR system is a much more challenging system and more representative of those encountered in medicinal applications. Both of these systems display order/disorder transitions along the alchemical path. The T4L/3-iodotoluene complex was simulated with the hyperbolic tangent (TanhSC), and the rational function soft-core (RatSC) functions with a linear perturbation schedule as well as with the integrated logistic perturbation potential (ILog) with the RatSC soft-core function. The FXR/26 complex was simulated with the same settings as T4L99A/3-iodotoluene and with an additional soft-core parameterization. In total, we report the results of 3 alchemical calculations for T4L/3-iodotoluene and 6 alchemical calculations for FXR/26 (see Table 1 for a summary).

The results are organized as follows. We first present the parameterization of the analytical model for the two systems (Figure 4 and Table 2) and the corresponding analytical predictions of the effects of different choices of soft-core functions and alchemical perturbations potentials (Figures 6 and 7). We then present the design of the alchemical schedule with the integrated logistic perturbation potential to address order/disorder transitions (Figure 8). The results of the alchemical calculations are shown next, with a particular focus on the effect of the various settings on the conformational sampling and replica exchange efficiency (Figures 9, 10, and Table 3), and the convergence of the binding free energies (Table 4 and Figure 11).

III.1 Parameterization of the Analytical Model of Alchemical Binding

The parameters for the analytical description of $p_{0}(u)$ [Eq. (5)]kilburg2018analytical for the complexes studied in this work are listed in Table 2. The good level of agreement between the analytical binding energy distributions, from Eq. (8), and the histograms obtained from the alchemical calculations are illustrated in Figure 4.

The model predicts the presence of three modes for the T4L/3-iodotoluene system and two modes for the FXR system (Table 2). In the case of T4L/3-iodotoluene, two of the modes (modes 1 and 2) correspond to alternative binding poses of 3-iodotoluene, one with strong interactions with the receptor (mode 1, with $\bar{u}_{b}=-11$ kcal/mol) and another (mode 2), more weakly bound, but approximately 13 times more likely to occur than the first in the uncoupled ensemble at $\lambda=0$ (see the corresponding values of the statistical weights in the second column of Table 2)). While they formally describe binding poses, the analytical model predicts that in the absence of receptor-ligand interactions the probability of occurrence of clash-free configurations is very small for these states (from $10^{-6}$ to $10^{-5}$, third column of Table 2). The third mode of T4L/3-iodotoluene (mode 3, with $89$% weight) corresponds to conformations in which the ligand is nearly freely rotating and translating within the binding site. We refer to modes such as this as “unbound”, keeping in mind however that in the alchemical approach the ligand is not allowed to leave the binding site. As reflected by the negligible $p_{b}$ parameter, the highly unfavorable average binding energy parameter, and the larger $n_{l}$ parameter (columns 3, 4, and 8), in this mode clashes with receptor atoms are more severe and occur with an overwhelmingly large probability.

The analytical model predicts a single binding pose for the FXR/26 complex (mode 1 in Table 2). The bound mode (mode 1) is favored relatively to the unbound mode (mode 2) by a more favorable interaction energy ($\bar{u}_{b}=-28$ vs. $10^{5}$ kcal/mol). At the same time, as measured by its statistical weight (second column in Table 2), the bound mode is predicted to be one hundred thousand times less likely than the unbound mode (mode 2) at $\lambda=0$. In addition, the probability of occurrence of configurations free of atomic clashes while in mode 1, is predicted to be very small ($1.5\times 10^{-10}$, third column of Table 2). Taking into account the small population of mode 1, the probability of observing conformations free of clashes at $\lambda=0$ is predicted to be as small as $1.5\times 10^{-15}$. The logarithm of this number is a measure of the entropy loss for binding.

Inspection of the molecular dynamics trajectories largely confirms the results of statistical analysis above. Indeed, 3-iodotolune is observed to visit predominantly two binding poses related by a swap of the positions of the methyl and iodo substituents within the binding site. The pose with weaker interactions (mode 2) is entropically favored and occurs more often at smaller values of $\lambda$. Conversely, the pose with stronger interactions occurs more frequently at values of $\lambda$ closer to 1. Conversely, and in agreement with the predictions of the analytical model, the inhibitor is observed to visit only one binding pose. In all complexes examined and with any of the simulation settings, we observed that, while in one of the bound modes, the ligand oscillates around the stable binding pose and that the range of the oscillations is greater near the uncoupled state at $\lambda=0$. Infrequently, and only if $\lambda$ is smaller than a critical value, the ligand transitions to a disordered state where it explores a wide range of positions and orientations within the binding site. Replicas mirror this behavior in the disordered unbound mode which, infrequently and only if $\lambda$ is larger than a critical value, transition to one of the binding modes where the ligand is ordered. As shown below, the frequency of order/disorder transitions such as these is influenced by the choice of the soft-core function and of the alchemical perturbation schedule.

III.2 $\lambda$-Functions

The $\lambda$ functions obtained from the analytical model [Eqs. (12) and (18)] and the parameters in Table 2 are shown in Figures 5, 6, and 7. As discussed above, the $\lambda$-function depends only on the chemical system and the choice of the soft-core function. Once parameterized, the analytical model yields, using Eq. (18), the $\lambda$ function for any choice of binding-energy based soft core function.

A non-monotonic behavior of the $\lambda$-function corresponds to multi-modal distributions which may be indicative of order/disorder transitions during the alchemical transformation. The data in Figure 6 for the T4L/3-iodotoluene system shows that the linear schedule and the TanhSC soft-core function is expected to yield an order/disorder transition between an ordered state with negative average binding energy and a disordered state with average binding energy close to $u_{\rm max}$. These states correspond, for example, to the intersections in Figure 6 of the horizontal line at $\lambda=0.217$ with the $\lambda$-function (as discussed above, the intersection near the upper limit of $u$ corresponds to a minimum followed by a maximum at $u=u_{\rm max}$). In this particular case, we expect that the ordered state and the disordered state are separated by a large energy gap of more than $50$ kcal/mol which leads, as presented below, to rare transitions and poor replica exchange efficiency. With the RatSC soft-core function (black curve in Figure 6), in contrast, the energy gap between low and high binding energy states is significantly reduced. For instance, at $\lambda=0.217$ the high energy state of $p_{\lambda}(u)$ is predicted to peak at $u\simeq 17$ kcal/mol rather than at $50$ kcal/mol with the TanhSC soft-core. As a result, the binding energy gap of the order/disorder transition is predicted to be reduced by roughly half.

As shown in Figure 7, the FXR/26 system displays order/disorder transitions of similar nature but of greater strength and complexity than T4L/3-iodotoluene. As indicated by the presence with the RatSC soft-core function of an additional, strong, peak of the $\lambda$-function at high binding energies, we expect that the alchemical coupling process for FXR/26 to display multi-modal behavior. Specifically, we expect to observe, starting at $\lambda=0$, sharp transfers of population from a highly disordered state with binding energies close to $u_{\rm max}$ to a less disordered state with binding energies in the $10$–$40$ kcal/mol range and then to coupled states at negative binding energies. With the same RatSC soft-core function as T4L/3-iodotoluene capped at $u_{\rm max}=50$ kcal/mol (solid black curve in Figure 7) we expect entropic bottlenecks against binding to occur at binding energies at around [math] and $35$ kcal/mol, respectively. These bottlenecks can be identified graphically by locating the intersections between horizontal lines at $\lambda$ and the $\lambda$-function corresponding to minima of the binding energy distribution function. The strength of the order/disorder transitions can be reduced using a less aggressive soft-core function with $u_{\rm max}=100$ kcal/mol (dashed black curve in Figure 7), at the expense of wider energy gaps between ordered and disordered states. With any of the soft-core settings, the $\lambda$-functions for the FXR/26 system we obtained indicate that the RatSC soft-core function to be significantly superior to the TanhSC soft-core functions in reducing the binding energy gaps between ordered and disordered states and lead to improved RE efficiency.

III.3 Design of the Integrated Logistic Alchemical Schedules

The $\lambda$-functions obtained above are used to design alchemical schedules based on the integrated logistic perturbation function from Eq. (19) to attempt to avoid order/disorder transitions or at least reduce their effects.Kim2010b The main design principle (see Figure 5) is to vary the $\lambda$ dependence of the parameters of the integrated logistic function so as that its derivative with respect to $u$ has a single intersection with the $\lambda$-function–yielding a binding energy distribution with a single maximum–or, when this is not easily achievable, at least at nearby points, thereby yielding maxima and minima separated by small energy gaps. The parameters of the optimized integrated logistic schedule for the systems studied in this work are listed in Tables 5, 6, and 7.

The general design strategy we followed in this work to parameterize the integrated logistic schedule is illustrated in Figure 8. The integrated logistic potential essentially allows using different $\lambda$ values in different ranges of the binding energy. The coupling transformation, which with the linear alchemical potential involves the progressive increase of the single coupling parameter $\lambda$, is divided into three phases. In the first phase (Figure 8A), the coupling parameter $\lambda_{2}$, for high binding energies is increased up to a critical value large enough to clear the maximum of the $\lambda$ function while $\lambda_{1}$, the coupling parameter for low binding energies, is left at zero. In the second phase (Figure 8B) $\lambda_{1}$ is now increased while $\lambda_{2}$ is left unchanged at the critical value. The second phase ends when $\lambda_{1}$ reaches the critical value. Finally, in the third phase (Figure 8B) $\lambda_{1}$ and $\lambda_{2}$ are increased in unison thereby acting as a linear perturbation up to $\lambda=1$. During the whole alchemical process $\partial W_{\lambda}(u)/\partial u$ intersects the $\lambda$-function at a single point, thereby yielding binding energy distributions with a single mode which progressively shifts to low binding energies without undergoing strong order/disorder transitions (Figure 9).

We followed the process outlined above to derive from the analytical $\lambda$-functions of the T4L/3-iodotoluene and FXR/26 systems (Figures 6 and 7) optimized integrated logistic schedules capable of avoiding the phase transition at $u\simeq 0$ (Tables 5, 6, 7). We have not attempted to design integrated logistic schedules aimed at resolving the phase transition of the FXR/26 system at large binding energies.

The binding energy probability distributions obtained from the molecular dynamics replica exchange calculations reported in Figures 9, largely confirm the predictions of the analytical model. With the TanhSC soft-core potential (panels A, B, and C), the binding energies of disordered states at $u\simeq u_{\rm max}$ abruptly shifts to negative values as $\lambda$ is lowered below a critical value ($\lambda=0.2$ approximately for both systems and settings considered). The distributions of disordered and ordered states on either side of this critical value are separated by a large binding energy gap which increases from $50$ to $100$ kcal/mol as $u_{\rm max}$ is increased.

Consistent with the predicted $\lambda$-functions (Figures 6 and 7)), the binding energy gap is reduced when using the RatSC soft-core potential (Figure 9, panels D, E, and F), mainly by shifting the distributions at $u=u_{\rm max}$ to lower values. In the case of T4L/3-iodotoluene the effect is very substantial to the point that some of the distributions of the disordered states overlap, albeit weakly, with those of ordered states. As shown in Figure 5 the distributions near the critical $\lambda$ are bimodal. The benefit of the RatSC function is not as significant for the FXR/26 system, which is characterized by a stronger order/disorder transition.

The integrated logistic alchemical perturbation potential is very effective at canceling the transition near $u=0$ (Figure 9, panels G, H, and I), especially for the T4L/3-iodotoluene system where it produces nearly homogeneous sampling of the whole binding energy range. The binding energy gap across the order/disorder transition For FXR/26 at large binding energies is reduced by a factor of five when using $u_{\rm max}=50$ kcal/mol (panel F) and by a factor of two with $u_{\rm max}=100$ kcal/mol (panel I). This result is due to the shifting of the distributions of the ordered state into the “no man’s land” region of binding energies between ordered and disordered states which are very poorly sampled with the linear potential.

III.4 Replica Exchange Efficiency

Replica exchange efficiency has been monitored in terms of the extent of diffusion of replicas in $\lambda$-space and binding energy space. The time trajectories of the binding energies sampled by the replicas are shown in Figure 10. They indicate that, while replicas rapidly equilibrate within the ordered low energy states and disordered high energy states of the complex, they more rarely interconvert between these states. We refer to these rare interconversions as binding or unbinding events. With the TanhSC soft-core potential replicas never, or only very rarely, undergo binding or unbinding transitions. The RatSC soft-core reduces the binding energy gap between ordered and disordered states and allowing the T4L/3-iodotoluene system to undergo numerous binding/unbinding transitions. However, only when also using the integrated logistic potential it is possible to observe binding/unbinding transitions for the FXR/26 system, which has a stronger order/disorder phase separation.

The number of binding and unbinding transitions, defined as the number of times that a replica goes from a disordered state to an ordered state or viceversa, are reported in Table 3. Here we identify ordered and disordered states based on their binding energy values. If $u>u_{\rm upper}$ we label the replica as disordered (unbound) and if $u<u_{\rm lower}$ we label it as ordered (bound). The values for $u_{\rm upper}$ and $u_{\rm lower}$ (see Table 3) were chosen based on the inspection of molecular dynamics trajectories. The trends are clear, confirm the qualitative analysis above, and do not depend on the specific choice of these parameters. In all cases, the TanhSC soft-core potential and the linear alchemical schedule lead to very few binding and unbinding transitions. As further discussed below, because of the lack of binding/unbinding events, the convergence of the calculations with the TanhSC soft-core potential is likely unreliable. The RatSC soft-core alone is sufficient to enable many interconversions for the T4L/3-iodotoluene system, one third more, in fact, than with the integrated logistic potential. The integrated logistic alchemical potential (ILog) produces many binding/unbinding events even for the more challenging FXR/26 system where the linear potential fails (Table 3). Based on this analysis, the combination of the RatSC soft-core and the integrated logistic potential emerges as the most consistent scheme to accelerate binding/unbinding transitions.

III.5 Binding Free Energy Estimates

The estimates of the standard free energy of binding for the systems and settings studied in this work are listed in Table 4. As the free energy is a thermodynamic state function, binding free energies should not depend on simulation settings, such as the choice of the soft-core function and the alchemical schedule. The data in Table 4 shows that all the simulation protocols yield consistent estimates for the T4L/3-iodotoluene complex, thereby validating the correctness of our implementation. This is particularly so for the simulations with the linear (RatSC) and integrated logistic (ILog) alchemical schedules (Figures 11B and 11C) that achieve rapid equilibration between coupled and uncoupled states (Figures 10B and 10C).

Equilibration analysis (Figure 11) indicates that it takes approximately 1 ns of simulation per replica to achieve equilibration for the T4L/3-iodotoluene system with the RatSC and ILog protocols (see below). Consistent with the smaller number of binding and unbinding events (Table 3), it takes about twice as long to achieve equilibration with the TanhSC protocol, although the steady downward drift of the binding free energy is a source of concern in this case (Figure 11A).

In this work, equilibration analysis is based on reverse cumulative profiles of the binding free energy with respect to the equilibration time, that is as a function the amount of data discarded from the beginning of the simulation.yang2004free ; chodera2016simple ; kilburg2018assessment Based on these (Figure 11), it is qualitatively clear that binding free energy estimates for T4L/3-iodotoluene are relatively robust.

It is less obvious to derive quantitatively converged and minimum variance estimates from reverse cumulative profiles. While the accuracy of the binding free energy improves as more unequilibrated samples are discarded, the precision of the estimate worsens as fewer samples are available. The trade-off between accuracy and precision is reflected in the steady increase in the size of the error bars in Figure 11 as the equilibration time increases. Several strategies have been proposed to extract optimal equilibration times and estimates from reverse cumulative profiles.yang2004free ; chodera2016simple ; kilburg2018assessment Here we take the simple approach of choosing the smallest equilibration time that gives a free energy estimate statistically indistinguishable from those at all subsequent equilibration times. Equilibration times chosen using this strategy are indicated by vertical lines in Figure 11 and are reported in the third column of Table 4.

Equilibration analysis based on the reverse cumulative plots indicates that the ILog protocol achieves equilibration of the binding free energy of the FXR/26 complex after 3 to 4 ns of simulation per replica, and the bias of the binding free energy estimate is less than half a kilocalorie per mole throughout the run (Figures 11F and 11I, and Table 4). In addition, the two variations of the ILog protocol with two different $u_{\rm max}$ parameters yield similar results ($-13.2$ and $-13.7$ kcal/mol, within nearly statistical uncertainty). These results, together with the relatively large number of binding/unbinding transitions (Table 3), confer high confidence in the binding free energy estimates for the FXR/26 complex with the ILog protocol.

In contrast, there appears to be an almost complete lack of equilibration of the binding free energy of FXR/26 with the linear perturbation and the TanhSC soft-core settings. The corresponding estimates are significantly more negative than with the other protocols and display little change as the equilibration time is increased (Figures 11D and 11G). This behavior is caused by the lack of unbinding transitions (Figures 11D and 11G) and it is of particular concern because conventional analysis would erroneously conclude in this case that the binding free energy estimate is well converged.

Similar, but less severe, issues occur for FXR/26 with the linear perturbation potential and the RatSC soft-core. With $u_{\rm max}=50$, despite the lack of binding/unbinding transitions, the binding free energy appears to equilibrate to the correct value on a time-scale longer than 5 ns (Figure 11E) per replica. With $u_{\rm max}=100$ the binding free energy is near the correct value at most equilibration times, but it also equilibrates slowly, and towards the end appears to start to steer in the wrong direction (Figure 11H). In both cases equilibration of the binding free energy estimate is uncertain. This is in contrast to the calculations with the integrated logistic potential which, from all accounts, appear equilibrated and converged (Figures 11F and 11I).

IV Discussion

In this work, we show that order/disorder phase transition can occur along the alchemical transformations and that these, by impeding rapid interconversions between conformational states across the thermodynamic phase boundary, are responsible for slow equilibration of binding free energy estimates. The bottlenecks caused by order/disorder transitions are particularly noticeable in Hamiltonian Replica Exchange sampling approaches, where the perturbation energy gap causes reduces the rates of acceptance of $\lambda$-exchanges among replicas.

Many methods are available to activate interconversions between low energy basins separated by high energy barriers,Okumura2010 ; Pierce2012 ; abel2017advancing ; xia2018improving especially when the slow collective variable is known.ciccotti2004blue ; barducci2008well ; bussi2018metadynamics Strategies to accelerate conformational sampling across entropic barriers, such as those considered here, are, however, not as established and understood. Probably the best known and still largely unsolved biophysical problem in this category is the process of protein folding, which is characterized by an entropic bottleneck to go from the disordered random coil state to the more ordered molten globule state.AlexeiV.Finkelstein2002 ; Zheng2007 Similarly, the process of protein-ligand binding involves loss of rotational and translational motion and the loss of conformational freedom of both the ligand and the receptor. In alchemical binding, this is reflected in the slow random search process of the rare ordered binding pose in the uncoupled ensemble, where the ligand and receptor are free to explore a wide range of configurations.Baron2008a ; pan2017quantitative

Here, we extend the formalism and non-Boltzmann conformational sampling simulation techniques developed by Straub at al. for the study of temperature-dependent phase transitions in spin systems, atomic clusters, and molecular liquids, to alchemical binding processes. There is a close analogy between the treatment of Straub et al.Kim2010b ; kim2010generalized ; lu2013order ; lu2014investigating for the canonical ensemble as a function of temperature and the alchemical ensemble as a function $\lambda$ considered here. Straub et al. considered non-Boltzmann sampling with a generalized reduced potential $W_{T}(E)$ where the modes of the energy distribution correspond to energies at which the inverse of the microcanonical statistical temperature function $T_{S}(E)$ equals the inverse of $k_{B}\partial W_{T}(E)/\partial E$:

[TABLE]

where $S(E)=k_{B}\ln\Omega(E)$ is the microcanonical entropy and $\Omega(E)$ is the density of states. In the case of canonical Boltzmann weighting, $W_{T}(E)=E/k_{B}T$, the condition Eq. (21) reduces to $1/T_{S}(E)=1/T={\rm constant}$. Near a first order phase transition the function $1/T_{S}(E)$ varies non-monotonically and intersects the constant line $1/T$ at three energy values corresponding to the phases in equilibrium and the unstable intermediate phase. In close analogy with Figure (1), the phase transition can be avoided by using a suitable non-linear sampling potential $W_{T}(E)$ whose derivative intersects $1/T_{S}(E)$ at only one value of $E$.Kim2010b

One of the main results of this work is the finding that the formalism of Straub et al. for temperature-driven order/disorder transitions extends nearly seamlessly to the alchemical ensemble for bindingGallicchio2010 ; Gallicchio2011adv by applying the following equivalences:

[TABLE]

We exploited these equivalences to design alchemical perturbation functions that avoid or reduce order/disorder phase transitions along the alchemical path. Using this strategies, we identified a particularly effective alchemical perturbation function in this respect, named the integrated logistic function [Eq. (19)], which we successfully employed to accelerate the sampling of two protein-ligand complexes affected by order/disorder transitions.

A critical step in achieving this result has been the description in analytical form of the $p_{0}(u)$ probability distribution function and the corresponding $\lambda$-function [Eq. (12)].kilburg2018analytical Estimation of the $\lambda$-function by numerical histogramming and finite difference techniques, in particular, leads to noisy results which are hard to interpret and utilize for the design of optimal perturbation functions. In contrast, analytical differentiation of $p_{0}(u)$ has proven a suitable strategy to obtain $\lambda_{0}(u)$ and formulate reliable predictions on the strength of order/disorder transitions, and on how to avoid them in numerical simulations.

The analytical route relies on a statistical theory of alchemical binding developed earlier.kilburg2018analytical The theory is based on a small number of physically-motivated parameters for each binding pose. The parameters are optimized by fitting the analytical predictions to binding energy distributions extracted from numerical simulations. Once the parameters are determined with some confidence, the theory is used to explore the effect of varying simulation conditions, such as alternative alchemical perturbation potentials and soft-core functions, as in the present work. At the core of the theory is the assumption that the binding energy is a random variable with two components, one that follows the central limit statistics and the other following max statisticsGumbelBook describing, respectively, favorable “soft” interatomic interactions and atomic clashes. Here the theory is used for the first time to describe order/disorder transitions in alchemical calculations.

V Conclusions

In this work, we investigate the role of order/disorder transitions in alchemical simulations of protein-ligand binding free energies. We show for a benchmarking system (the complex between the L99A mutant of T4 lysozyme and 3-iodotoluene) and for a more challenging system relevant for medicinal applications (the complex of the farnesoid X receptor and inhibitor 26 from a recent D3R challenge) that order/disorder transitions can significantly hamper Hamiltonian replica exchange sampling efficiency and slow down the rate of equilibration of binding free energy estimates. We further show that our analytical model of alchemical bindingkilburg2018analytical combined with the formalism developed by Straub et al. for the treatment of order/disorder transitions of molecular systemsKim2010b can be successfully employed to analyze the transitions and help design alchemical schedules and soft-core functions that avoid or reduce the asdverse effects of rare binding/unbinding transitions. Future work from our laboratory will explore the use of these techniques for alchemical free energy estimation with explicit solvation, specifically hydration free energies and double-decoupling binding free energies, for systems undergoing order/disorder transitions.

Acknowledgements.

We acknowledge support from the National Science Foundation (NSF CAREER 1750511 to E.G.). Molecular simulations were conducted on the Comet GPU supercomputer cluster at the San Diego Supercomputing Center supported by NSF XSEDE award TG-MCB150001 and on a GPU cluster at the Laboratory of Computational Molecular Biophysics at Brooklyn College of CUNY.

Appendix A Appendix

A.1 Derivation of the Pair Collisional Distribution Function

Consider the decomposition of the Lennard-Jones pair potential $u_{{\rm LJ}}(r)=4\epsilon[(\sigma/r)^{12}-(\sigma/r)^{6}]$ into repulsive and background components based on the cutoff distance $r_{C}$ (Fig. 12).

[TABLE]

[TABLE]

where

[TABLE]

is the value of the LJ potential at the cutoff distance $r_{C}$. The standard WCA decomposition of the LJ potential is a particular case of Eqs. (27) and (28) with a cutoff distance corresponding to the minimum of the LJ potential, $r_{C}=2^{1/6}\sigma$.

In this work, we model the collisional and background components of the binding energy between the receptor and the ligand in the uncoupled state are due to the repulsive and background component of the LJ potential, Eqs. (27) and (28) respectively, for some choice of the LJ cutoff energy $\tilde{u}$ [Eq. (29)]. This choice is justified as follows. The collisional component is defined as the contribution to the binding energy described by max statistics.kilburg2018analytical Max statistics applies when the interaction energy is dominated by the single pair of receptor and ligand atoms with the largest repulsive interaction energy. Thus, max statistics can set in only if the repulsive interaction energy potential, as in Eq. (27), grows rapidly for small interatomic distances, which are less likely. Conversely, the collisional energy should become negligibly small for large interatomic distances, which are more likely and tend to occur for many pairs of atoms at once. In contrast, the background component is defined as the contribution to the binding energy described by central limit statistics. Because central limit statistics assumes contributions from many interatomic interaction energies of similar magnitude and likelihood of occurrence, it cannot apply to a potential which grows indefinitively with decreasing interatomic distance. A sensible choice for the background component is therefore an interatomic potential that is constant below a certain interatomic distance, as in Eq. (28), thereby leaving its complement, Eq. (27), as the necessary description for the repulsive component of the LJ potential.

The choice of Eq. (27) as a viable description of the repulsive component of the interatomic potential is further supported by the requirement that it should grow sufficiently rapidly with decreasing interatomic distance such that the average and variance of the repulsive pairwise energy tend to infinity, so as to exclude the applicability of central-limit statistics to the collisional component. The probability distribution of the collisional component derived from Eq. (27) satisfy this requirement.

Finally, our collisional interatomic potential energy model does not include electrostatic interactions. This is partly because of the difficulty of representing the statistics of contributions of alternating sign coming from like and unlike partial charges. Nevertheless, we expect our LJ-based statistical model to be accurate at short interatomic distances where the $1/r^{12}$ repulsion dominates over the much weaker $1/r$ behavior of the electrostatic interactions. The contributions of electrostatic interactions to the background component of the binding energy are fully captured by the average background interaction energy parameter $\bar{u}_{B}$ (see above).

Now consider two particles interacting by the repulsive LJ potential (27) in which one particle (representing the receptor) is fixed at the origin and the other (representing the ligand) is uniformly distributed in a sphere of radius $r_{C}$ centered at the origin (Fig. 12). We will derive the probability density $p_{C}(u_{C})$ of the repulsive interaction energy $u_{{\rm C}}(r)$, where $r$ is the distance between the two particles, by differentiating the cumulative probability function $P_{C}(u_{C})$ defined as the probability that, given that the ligand particle is uniformly distributed in the sphere, the interaction energy $u_{{\rm C}}(r)$ is greater than the given value $u_{C}$.

The probability that the pair interaction energy is smaller than $u_{C}$ is given by:

[TABLE]

where the Heaviside function imposes the requirement that $u_{C}$ be larger than zero, $V_{C}$ is the volume of the sphere of radius $r_{c}$ and $V(u_{C})$ is the volume of the sphere of radius $r(u_{C})$, where $r(u_{C})$ is inter-particle distance at which the repulsive LJ potential has value $u_{C}$. From Eq. (27) we have

[TABLE]

where $r_{0}=2^{1/6}\sigma_{LJ}$ is the minimum of the Lennard-Jones pair potential and

[TABLE]

Inserting Eq. (31) into Eq. (30) gives

[TABLE]

Differentiating 33 with respect to $u_{C}$ and setting

[TABLE]

finally yields

[TABLE]

which expresses a normalized distribution as it can be verified by direct integration using the fact that

[TABLE]

Eq. (35) is the basis of the derivation of the collisional energy distribution function in Eq. (7) for a polyatomic ligand interacting with a polyatomic receptor as discussed in reference kilburg2018analytical, .

A.2 Gibbs Independence Sampling Algorithm

While it not central to the present analysis of alchemically-induced phase transitions, for completeness, here we describe the implementation of the Gibbs Independence Sampling algorithm we used in this work in conjunction with asynchronous replica exchange. The scheme has been used previously,pal2016SAMPL5 ; bfzhang2016rnaseh but has not yet been reported explicitly in the literature.

The space of assignments of $\lambda$-states to replicas is composed of $M!$ discrete members each corresponding to a permutation of $M$ replicas on $M$ $\lambda$-states. A permutation $r$ is defined by the set of indexes $\{k_{i}(r)\}$ specifying the state $k$ assigned to each replica $i$. In our specific implementation of Gibbs Independence Sampling, starting with current state assignment $r$, we consider transitioning to a set of $M-1$ new state permutations $s$ differing from $r$ by a single swap of $\lambda$-states between a chosen pivot replica $i$ and any other replica $j$. The discrete Metropolis transition matrixAllen:Tildesley of this random walk in state permutations is

[TABLE]

where the proposal probability $\alpha_{rs}$ is $1/(M-1)$ if the permutations differ by a single replica swap between replica $i$ and some other replica $j$ as described above and zero otherwise, and $v_{r}$ is the reduced generalized bias potential energy of the replica exchange ensemble corresponding to the state assignment $r$:

[TABLE]

where $W_{k}(u_{i})$ is the alchemical perturbation energy [Eq. (3)] of replica $i$ at $\lambda_{k}$.

Given a pivot replica $i$, the algorithm selects a replica $j$ for exchange with probability Eq. (37). The selected replica could include the same replica $i$, in which case no exchange occurs. The process is repeated sequentially for each replica $i$ so that at the end of the exchange cycle each replica has had the opportunity to exchange with any other replica. It is straightforward to show the random walk produced by this algorithm satisfy microscopic reversibility. To reduce computational cost, the set of energies $W_{k}(u_{i})$ is pre-computed at the beginning of each exchange cycle. We prefer this algorithm over alternatives,chodera2011replica because it automatically selects likely exchanges in an unbiased manner without having to specify neighboring states. Moreover, it naturally scales appropriately the number of exchange moves in proportion to the number of replicas.