Linear Interaction Energy (LIE) Models for Ligand Binding in Implicit

256

J. Chem. Theory Comput. 2007, 3, 256-277

Linear Interaction Energy (LIE) Models for Ligand Binding in Implicit Solvent: Theory and Application to the Binding of NNRTIs to HIV-1 Reverse Transcriptase Yang Su,† Emilio Gallicchio,*,† Kalyan Das,‡ Eddy Arnold,‡ and Ronald M. Levy*,† BioMaPS Institute of QuantitatiVe Biology, Department of Chemistry and Chemical Biology, and Center for AdVanced Biotechnology and Medicine, Rutgers UniVersity, Piscataway, New Jersey 08854 Received August 8, 2006

Abstract: Expressions for Linear Interaction Energy (LIE) estimators for the binding of ligands to a protein receptor in implicit solvent are derived based on linear response theory and the cumulant expansion expression for the free energy. Using physical arguments, values of the LIE linear response proportionality coefficients are predicted for the explicit and implicit solvent electrostatic and van der Waals terms. Motivated by the fact that the receptor and solution media may respond differently to the introduction of the ligand, a novel form of the LIE regression equation is proposed to model independently the processes of insertion of the ligand in the receptor and in solution. We apply these models to the problem of estimating the binding free energy of two non-nucleoside classes of inhibitors of HIV-1 RT (HEPT and TIBO analogues). We develop novel regression models with greater predictive ability than more standard LIE formulations. The values of the regression coefficients generally conform to linear response predictions, and we use this fact to develop a LIE regression equation with only one adjustable parameter (excluding the intercept parameter) which is superior to the other models we tested and to previous results in terms of predictive accuracy for the HEPT and TIBO compounds individually. The new models indicate that, due to the different effects of induced steric strain of the receptor, an increase of ligand size alone opposes binding for ligands of the HEPT class, whereas it favors binding for ligands of the TIBO class.

1. Introduction The binding free energy of a ligand to a receptor is given by the difference of the free energies of inserting the ligand in the receptor and in solution. In principle the free energy for each process can be calculated exactly for a given force field using the free energy perturbation (FEP) or thermodynamic integration (TI) methods. In practice, however, the complexities involved in setting up suitable mutation paths * Corresponding author e-mail: [email protected] (E.G.) and [email protected] (R.M.L.). † BioMaPS Institute of Quantitative Biology and Department of Chemistry and Chemical Biology. ‡ Center for Advanced Biotechnology and Medicine and Department of Chemistry and Chemical Biology.

10.1021/ct600258e CCC: $37.00

and the long simulation times needed to reach convergence have limited the applicability of FEP and TI methods to the investigation of the variations of the binding free energy for small ligand modifications in the final stages of lead optimization.1 Linear Interaction Energy (LIE) models2,3 offer attractive approximate alternatives to the full FEP methodology because they require only the computation of average interaction energies at the end points of the mutation. LIE models can be described as empirical Quantitative Structure-Activity Relationships (QSAR) which employ physically motivated energetic estimators. As opposed to methods that predict the binding free energy on the basis of the structure of the ligand alone, LIE estimators also reflect properties of the ligand-receptor complex. LIE methods are expected to perform better than methods based on the ligand © 2007 American Chemical Society

Published on Web 11/10/2006

Linear Interaction Energy Models for Ligand Binding

J. Chem. Theory Comput., Vol. 3, No. 1, 2007 257

alone because they are more intimately related to the structure of the complex. One of the most popular LIE formulations employs the following regression expression for the binding free energy ∆Fb3,4 ∆Fb ) R∆V h vdw + β∆V h el + γ∆A h +δ

(1)

where ∆V h vdw, ∆V h el, and ∆A h are differences between quantities measured for the ligand complexed with the receptor and the ligand free in solution. V h vdw is the average van der Waals interaction energy between the ligand and its environment (the solvent or the receptor and the surrounding solvent). Similarly, V h el is the average electrostatic interaction energy between the ligand and its environment. A h is the average solvent accessible surface area of the ligand. These quantities are typically obtained from Molecular Dynamics (MD) or Monte Carlo (MC) simulations started from known or modeled conformations of the ligand and the ligand-receptor complex. R, β, γ, and δ are empirical adjustable parameters whose values are obtained by fitting the model over a training set of ligands of known binding affinity. A trained LIE model can then be used to predict the binding free energy of ligands of unknown affinity provided that the LIE estimators for these ligands can be reliably calculated. Each LIE estimator in eq 1 reflects physical forces that affect binding. ∆V h vdw and ∆V h el measure the balance between the desolvation penalty caused by the loss of ligand-solvent interactions and the gain of ligand-receptor interactions which favor binding, whereas the surface area estimator measures the hydrophobic driving force toward complexation. These LIE estimators do not take into account explicitly thermodynamic forces related to the receptor that affect the binding affinity, such as the desolvation of receptor atoms and reorganization free energy of the receptor for accommodating the ligand. Nevertheless it can be shown (see below) that, under the assumption of linear response, these effects are, in fact, included in the model and are encoded in the values of the LIE regression coefficients. LIE models have their origin in physical theories of solvation based on the linear response approximation to the free energy,5-7 applied to the problem of binding free energy estimation.2,8-11 The introduction of the ligand in either the solution or receptor environments can be regarded as a perturbation applied to the system. If the system responds perfectly linearly (as formally defined below) to the perturbation, the free energy of introducing the ligand can be shown to be exactly proportional to the interaction energy between the ligand and its environment. Given that the ligand constitutes a large perturbation to the system, it is unlikely that linear response applies to the entire processes of introducing the ligand in solution and in the receptor. The LIE regression equation (eq 1) assumes instead that linear response applies to the individual processes of introducing hydrophobic, van der Waals, and electrostatic interactions, albeit with different linear response proportionality coefficients. Even so, nonlinearities in practice limit the applicability of LIE models to within a related class of ligands. This is reflected in the fact that in practice LIE relationships are used for estimating relatiVe binding free energies of similar ligands rather than absolute binding free energies.

The accuracy of a LIE model therefore hinges on whether the perturbation corresponding to the mutation of each ligand into another is small enough so that linear response applies to the relative binding free energy. The effect of absolute binding free energies is collectively absorbed by the intercept parameter δ; deviations from linear response are expected to be reflected in the limited range of applicability of a LIE model. In this paper we investigate a series of outstanding issues with regard to LIE models and their applications to ligand binding in structural biology. The first question is to what extent linear response applies to a given ligand set and the consequences of deviations from linear response in terms of the accuracy of the LIE model. Although in principle addressing this question requires comparing LIE predictions with relative free energies evaluated using the rigorous FEP and TI methods, we take a first step in this direction by comparing the values of the LIE adjustable parameters obtained by fitting training sets of ligand binding data with those expected based on linear response theory. The second question concerns the form of the LIE regression equation. Equation 1 assumes that the response of the solution and receptor environments, as measured by the LIE coefficients R, β, and γ, is the same. To our knowledge this assumption is ubiquitous in LIE applications. We develop and validate an alternative formulation in which the processes of insertion of the ligand in the solution and in the receptor are decoupled so that each is allowed to have different linear response proportionality coefficients. Finally, based on linear response techniques we analyze the appropriate form of the LIE estimators when the solvent is treated implicitly. We show that the form for the electrostatic LIE estimator we derive based on linear response theory agrees with the corresponding estimator recently proposed by Carlsson et al.12 and differs from the more empirical expression proposed earlier by Zhou et al.4 We develop, following the linear response formalism, electrostatic, van der Waals, and cavity implicit solvent estimators that best represent the corresponding LIE estimators in explicit solvent and show that these lead to improved accuracy. We apply these ideas to ligand-protein complexes that have been studied previously using the LIE method: the binding of the HEPT and TIBO classes of Non-Nucleoside Reverse Transcriptase Inhibitors (NNRTI) to the HIV-1 Reverse Transcriptase (RT) enzyme.13 HIV-1 RT is essential for the life cycle of the virus. It converts the single stranded genomic RNA into double stranded DNA which is subsequently integrated into the host chromosome and passed on to all progeny cells.14 Computer-aided structure-based drug discovery technologies have made a significant contribution to the development of medicinally active NNRTI anti-AIDS compounds.15-17 Three of these compounds, dapivirine, etravirine, and rilpivirine, are currently in clinical trials. The goal is the design of NNRTIs of greater potency and resilience with respect to common drug-resistance mutations.18 However the mode of binding of NNRTIs, which includes extensive receptor conformational reorganization and mainly nonspecific hydrophobic ligand-receptor contacts, does not offer obvious chemical modification leads for

258 J. Chem. Theory Comput., Vol. 3, No. 1, 2007

binding free energy optimization. In this context, LIE models have been applied to the estimation of the binding affinities of NNRTI inhibitors with encouraging results.19,20 The renewed interest of our laboratories in the LIE modeling of HIV-1 RT inhibitors is motivated in part by recent progress in obtaining a crystal structure of HIV-1 RT in complex with an inhibitor of the N-acyl hydrazone class.21 The crystal structure identifies a novel non-nucleoside binding site adjacent to but distinct from the NNRTI binding site. The NAH inhibitors targeting this site are expected to suffer from little or no cross-resistance from existing drug resistance mutations, thus providing new options for novel therapeutic strategies in the treatment of AIDS. This crystal structure provides the initial framework for the application of LIE methodologies for lead optimization and the development of a new class of inhibitors of HIV-1 RT. The work presented in this paper on NNRTIs provides the theoretical and computational basis for the application of LIE modeling in implicit solvent for binding free energy estimation of the new NAH class of HIV-1 RT inhibitors which is currently being investigated in our lab. In the following section we review the theoretical foundations of the LIE method and discuss some of the approximations. We then derive a LIE formalism appropriate for situations where the solvent is modeled implicitly. We propose novel LIE regression equations which emerge naturally from the statistical mechanics derivation. We then apply these regression equations to analyze the binding of a series of 20 HEPT inhibitors and 37 TIBO inhibitors of HIV-1 RT. The models are validated using jack-knife prediction tests. The predictive ability of different forms and parametrizations of the LIE equations are compared. We compare features of the HEPT and TIBO binding modes and the corresponding LIE estimators. We conclude the paper with a discussion of the accuracy and physical interpretation of LIE models.

Su et al.

charging parameter λ, the interactions of the ligand with the surrounding environment (the receptor and/or the solvent). It is convenient to think about these processes in stages. First the ligand cavity is formed, then the van der Waals interactions between the ligand and the environment are turned on, and finally ligand-environment electrostatic interactions are established. At each stage the λ-dependent potential energy of the system is U(x;λ) ) U0(x) + λV(x)

where V represents the ligand-environment interactions which are being added, and U0, the reference potential energy, contains receptor-receptor, solvent-solvent, and ligandintramolecular interactions as well as the ligand-environment interactions established in the previous stages. Starting from the expression of the configurational partition function Z(λ) ) e-F(λ)/kT )

∫dx e- λV(x)/kTe-U (x)/kT 0

(3)

it is straightforward to show that the first and second derivatives of the free energy F(λ) ) - kT lnZ(λ) with respect to λ correspond, respectively, to the first and second moments of the probability distribution of V: ∂F ) 〈V〉λ ∂λ

(4)

∂2F ) -〈(δV)2〉λ/kT ∂λ2

(5)

If the second moment is approximately constant along the thermodynamic path from λ ) 0 to λ ) 1, the third and higher order derivatives of F can be neglected, and it is possible to express F(λ) as a Taylor series (known as the cumulant expansion of the free energy) centered at λ ) 0 and truncated at the second order

2. Theory and Methods 2.1. Linear Response Approximation. The idea of adopting a linear response approximation expression for binding free energy estimation was first stated by Lee et al.,8 who suggested the use of a “two-point” Linear Response Approximation (LRA) estimation formula previously derived for electrostatic solvation;6 see, for example, discussion in ref 10. This formula was later simplified by Åqvist et al.2 who introduced the Linear Interaction Energy (LIE) model which, unlike the more accurate “two-point” LRA formula,22 does not require the evaluation of estimators at more than one state of the ligand. The LIE method was used by Åqvist and collaborators to estimate relative binding affinities of endothiapepsin and HIV-protease inhibitors.2,9 We review here some basic concepts related to the linear response approximation and derive a linear response expressionseq 12 belowsfor the LIE proportionality coefficients, which will be used in the following to interpret the values of the LIE fitting coefficients obtained from the analysis of experimental binding affinities. The insertion of the ligand into either the receptor or the solvent can be thought of as turning on, by means of a

(2)

c F(λ) - F(0) ) 〈V〉0λ - λ2 2

(6)

c ) 〈(δV)2〉/kT

(7)

where

is assumed constant for 0 e λ e 1. It is of interest to note that when the fluctuations of the interaction potential V are Gaussian distributed this assumption is verified and eq 6 is exact.6 According to eq 6 and under these assumptions, the free energy is quadratic with respect to the charging parameter. This is a manifestation of linear response behavior defined as when the solute-environment average interaction energy 〈V〉λ is linearly related to the charging parameter. Indeed, using eqs 4 and 5, and the same assumptions that have lead to eq 6, we obtain 〈V〉λ ) 〈V〉0 - cλ

(8)

which confirms that 〈V〉λ is linearly related to λ. By evaluating eqs 6 and 8 at λ ) 1 we obtain the free energy change and average interaction energy for adding solute-environment interactions under the assumption of linear response


∆F ) F(1) - F(0) ) 〈V〉0 -

c 2


(9)

and 〈V〉1 ) 〈V〉0 - c

(10)

Linear interaction energy models are based on the assumption that ∆F is proportional to 〈V〉1:2 ∆F ) R〈V〉1

(11)

It is therefore of interest to compute, under the assumption of linear response, the proportionality coefficient R given by the ratio of ∆F to 〈V〉1. From eqs 9 and 10 ∆F 〈V〉0 - c/2 ) 〈V〉1 〈V〉0 - c

(12)

The limiting values for this ratio are

[

1/2 if |〈V〉0| , c ∆F ) 1 if |〈V〉0| . c 〈V〉1

]

(13)

Thus, under the assumption of linear response, the limiting values of the ratios between the free energy of adding soluteenvironment interactions and the corresponding average interaction energy are 1/2 and 1. The free energy change is half the interaction energy when the fluctuation of V, measured by c ) 〈(δV)2〉/kT, is much larger than 〈V〉0, the mean ligand-environment interaction energy calculated within the ensemble of conformations obtained in the absence of ligand-environment interactions. In the opposite limit at which the fluctuations of V are much smaller than 〈V〉0, eq 12 gives ∆F/〈V〉1 ) 1, that is the free energy change is equal to the average solute-environment interaction energy. In the following we apply linear response to the problem of binding free energy estimation and use eq 12 and physical arguments to derive values of the linear response coefficients. We will first review the derivation in explicit solvent and then examine the case in which the solvent is treated implicitly. 2.2. LIE Models in Explicit Solvent. 2.2.1. Hydration Free Energy - Explicit SolVent. To successfully apply linear response ideas to the hydration free energy estimation, the process of hydration is described as occurring in stages. First the solute cavity is formed in the solvent, then solute-solvent van der Waals interactions are turned on, and finally solutesolvent electrostatic interactions are established. The process of cavity formation is dominated by excluded volume effects and solvent reorganization and does not conform well to the linear response formalism. When using a hard-sphere cavity interaction potential, the solute-solvent interaction energy is zero when the solute and the solvent do not overlap, and it is infinite when overlaps occur. In this limit the average solute-solvent interaction energy is identically zero because conformations in which solute cavity-solvent overlaps exist do not appear in the ensemble. Computational studies of cavity formation have generally been conducted using a continuous but sharp solute-solvent repulsive interaction potential.23 In these cases, however, due to strong nonlin-

earities in the response of the solvent to the introduction of the solute cavity, the cavity hydration free energy is poorly correlated with the average repulsive cavity interaction potential. For cavity hydration free energy estimation a term proportional to the solute surface has been shown in some cases to be a reasonably good estimator for the free energy of cavity formation in water.24,25 The proportionality coefficient γ between the cavity hydration free energy and the surface area can be interpreted as a surface tension coefficient. Indeed molecular simulations have obtained values of this proportionality coefficient similar to the value of the experimental air-water surface tension coefficient.23,26 Assuming linear response for the processes of introducing van der Waals and electrostatic solute-solvent interactions, the free energy change at each stage is proportional to the appropriate solute-solvent interaction energy (the solutesolvent van der Waals interaction energy and the solutesolvent electrostatic interaction energy, respectively) averaged in the state corresponding to the end of each stage (the uncharged solute and the fully interacting solute, respectively). The linear response coefficients are given by eq 12 and are, in general, different for each stage. Finally, making the approximation that all averages (denoted by 〈‚‚‚〉w) are calculated when the solute fully interacts with the solvent, a LIE model for the hydration free energy is obtained27 ∆Fh = R〈Vvdw〉w + β〈Vel〉w + γ〈A〉w

(14)

where R and β are the linear response proportionality coefficients for the van der Waals and electrostatic stages of the hydration process, Vvdw and Vel are the solute-solvent van der Waals and electrostatic interaction energies, respectively, γ is an empirical surface tension coefficient, and A is the solvent accessible surface area of the solute. As we show below, the form of the cavity hydration term can be justified in terms of linear response when the solvent is modeled implicitly. Regression equations based on eq 14 have been parametrized against known experimental hydration free energies of small molecules.27 However, linear response theory provides a way to estimate some of these parameters from first principles. It has been observed that the charging free energy of ionic and polar solutes in water is approximately proportional to the average solute-solvent electrostatic interaction energy with a linear response proportionality coefficient β of 1/2. According to eq 13 this occurs when the solvent reaction field in the absence of solute charges is much smaller than the fluctuations of the solvent reaction field. These conditions have been indeed verified by numerical studies, confirming that in general water behaves as a good linear dielectric medium.6,28,29 These observations have constituted the basis for electrostatic linear response free energy models for solutions2,6,7 as well as for the success of continuum dielectric models of water.28,30-35 It is well-known that the process of adding solute-water van der Waals interactions has different characteristics than the process of adding electrostatic interactions. It has been shown that the free energy change for adding attractive van der Waals interactions can be well approximated by the average solute-water van der Waals interaction energy.23,36,37


This implies that in this case ∆F/〈Vvdw〉1 = 1, which, under the assumption of linear response occurs (see eq 13) when the mean solute-solvent attractive van der Waals interaction energy in the absence of solute-solvent van der Waals interactions, is larger than the variance of the same quantity divided by kT. This is due to the fact that, although the solvent responds linearly to the solute perturbation, this response is smaller than the attractive van der Waals solvent field that exists at the solute location in the absence of solute-solvent van der Waals interactions. Contrary to dipolar fields in water that tend to cancel each other in the absence of a polarization source, van der Waals interactions are always additive. Based on this analysis we conclude that the R coefficient in the LIE regression equation should assume a value near 1. Explicit solvent simulations have generally confirmed this theoretical prediction.23 Carlson and Jorgensen27 have instead reported that a value of R significantly smaller than 1 is obtained by fitting the LIE equation to experimental hydration free energies of small molecules. We believe that the small value of R obtained by Carlson and Jorgensen is caused by compensation between the van der Waals and surface area fitting coefficients. The correct relative magnitude of these two coefficients is difficult to pinpoint because they correspond to highly correlated descriptors (the van der Waals solute-solvent interaction energy and the solute surface area). Indeed both the R and γ coefficients obtained by Carlson and Jorgensen are smaller than well established theoretical predictions.24,36,38 We have reanalyzed the data from Tables 2 and 3 of ref 27. By setting R ) 1 and allowing for a nonzero intercept (to fit relative hydration free energies rather than the absolute ones) we achieved a nearly equivalent fit to the experimental hydration free energies. Specifically, using the Coulombic+van der Waals+solvent-accessible surface area model of Carlson and Jorgensen we reproduce the parameters reported previously,27 R ) 0.49, β ) 0.42, and γ ) 20 cal/mol Å2, with a crossvalidated correlation of Rpred2 ) 0.83, whereas our model with R set to 1 gives β ) 0.49 and γ ) 62 cal/mol Å2 with Rpred2 ) 0.81. The RMSD from the experimental free energies of hydration of the two models are also similar, 0.88 and 0.98 kcal/mol, respectively. Furthermore, the value of γ (62 cal/mol Å2) we obtained from the analysis of the data of Carlson and Jorgensen is closer to the experimental value of the macroscopic vacuum-water surface tension (104 cal/ mol Å2)39 and is similar to the value of the microscopic surface tension parameter obtained from explicit solvent estimates of the work of cavity formation in water (73 cal/ mol Å2).23 These observations indicate that the simulation data obtained by Carlson and Jorgensen is consistent with the linear response behavior for these solutes. In conclusion, this analysis shows that eq 14 should provide a good approximation of hydration free energies with the following choice of LIE coefficients: R = 1, β = 1/2, and γ = 73 cal/mol Å2, the previously reported explicit solvent estimate.23 2.2.2. Binding Free Energy Estimation - Explicit SolVent. The binding free energy ∆Fb of a ligand to a receptor is taken as the difference of the work ∆Fc of creating the ligand in the receptor and the work ∆Fh of creating the ligand in

Su et al.

solution. For the process of creating the ligand in the receptor an expression similar to eq 14 has been proposed3,9 where the interaction energies include ligand-water interactions as well as ligand-receptor interactions, and A is taken as the solvent accessible surface area of the ligand in the receptor-ligand complex. By taking the difference between the LIE estimates of the insertion free energies in the receptor and in solution and assuming that the same values of the LIE coefficients R, β, and γ are appropriate for both, the following binding free energy LIE model is obtained ∆Fb ) ∆Fc - ∆Fh = R(〈Vcvdw〉c - 〈Vfvdw〉w) + β(〈Vcel〉c - 〈Vfel 〉w) + γ(〈A〉c - 〈A〉w) (15) where V c represents an interaction energy between the ligand and the environment (receptor and solvent), and V f is the corresponding ligand-solvent interaction energy in the absence of the receptor (free ligand), A is the solvent accessible surface area of the ligand, 〈‚‚‚〉c represents an ensemble average with the ligand in the receptor pocket, and 〈‚‚‚〉w represents the corresponding average in solution. Equation 15 and its variations are widely used in binding free energy prediction applications.40-42 In these studies the LIE coefficients R, β, and γ are obtained by fitting eq 15 to known experimental binding free energies. In principle, linear response theory arguments could be applied, as for the case of hydration free energy estimation, to gain insights into the expected values of these coefficients. It is less clear, however, to what extent the receptor/solvent environment can be considered an ideal linear dielectric medium10,43 and what value of the proportionality coefficient to use to estimate the electrostatic charging free energy from the ligandenvironment electrostatic interaction energy. The LIE equation eq 15 implicitly assumes that the same electrostatic proportionality coefficient, β, applies to both the charging process in solution and in the protein receptor. However, contrary to the water environment, many proteins produce strongly anisotropic electrostatic fields. It is therefore reasonable to assume that a non-negligible electrostatic field exists at the binding site even in the absence of ligand charges. Under the assumption of linear response, eq 12 indicates that the presence of a residual average electrostatic potential at the ligand charge sites in the absence of ligand charges (that is 〈Vel〉0 * 0) will cause the ratio ∆F/〈Vel〉1 to deviate from the ideal solution value of 1/2. This analysis predicts that, although assuming β ) 1/2 is a reasonable first guess, in general it would be advantageous to adopt a LIE regression equation in which the electrostatic estimator is split into a receptor environment component and a solution environment component each multiplied by an independent LIE coefficient.44 Similarly, the LIE eq 15 implicitly assumes that the work of cavity formation in the protein receptor can be also estimated by the ligand surface area using a single surface tension coefficient applicable to both the water and receptor environments. However, due to the complex reorganization of the receptor binding pocket induced by ligand binding, the solute surface area is likely a poor descriptor for the free energy of ligand cavity formation in protein receptors. In this work we explore the alternative approach



of modeling the work of inserting the ligand in the receptor independently from the work of inserting the ligand in solution. Most of the limitations of the LIE regression equation underlined here are mitigated in practice by cancellation of errors. The main goal of LIE studies is to obtain relative binding free energies within a family of related ligands, rather than absolute binding free energies. The overall deviation of estimated absolute binding free energies from the experimental affinities is accounted for by an adjustable intercept parameter that is often added to the LIE regression equation (see eq 15).44 This intercept parameter obviously does not affect relative binding free energies estimated from the LIE equation. For example, the differences of cavity formation free energies between ligands of similar shape could be well represented by a surface area descriptor even though the individual absolute values are not. In this context the limits of the LIE regression equation are manifested in the limited range of applicability of a particular parametrization rather than in the accuracy of the LIE parametrization for a particular ligand set. A better understanding of the properties of the LIE equation could therefore lead to improved LIE ligand coverage and to rules that, based on properties of the ligand, associate a particular parametrization to a particular class of ligands. 2.3. The AGBNP Implicit Solvent Model. The Analytical Generalized Born plus Non-Polar (AGBNP) implicit solvent model35 is based on an analytical pairwise descreening implementation of the Generalized Born model and a nonpolar hydration free energy model consisting of an estimator for the solute-solvent van der Waals dispersion energy and a surface area term corresponding to the work of cavity formation. In the Generalized Born (GB) model33 the electrostatic component of the hydration free energy is estimated as

( )

1 1 1 Gel = GGB ) 2 in w

qiqj

∑ij f

(16) ij

where in is the dielectric constant of the interior of the solute, w is the dielectric constant of the solvent (in this work in ) 1 and w ) 80), qi and qj are the charges of atom i and j, and fij ) xrij2 + BiBj exp(-rij2/4BiBj)

(17)

where rij is the distance between atoms i and j, and Bi and Bj are the Born radii of atoms i and j defined below. The summation in eq 16 runs for all atom pairs (i, j) including i ) j. The diagonal i ) j terms can be separated from offdiagonal terms i * j yielding the equivalent expression

( )(∑

1 1

1

GGB ) 2 in w

i

qi2 Bi

+2

ij

)

1 1 1 ) Bi Ri 4π

∫d3r (r -1 r )4 Ωi

(18)

The first summation at the right-hand side of eq 18 is the sum of the GB self-energies of the atoms of the molecule, and the second term is the sum of the GB pair-energies. The self-energy of atom i corresponds to the solvation energy of

(19)

i

where Ωi is the bounded region corresponding to the solute volume excluding the atomic sphere corresponding to atom i, and Ri is the van der Waals radius of atom i. 1/Ri is the inverse Born radius of atom i in the absence of all the other solute atoms. The second term on the right-hand side of eq 19 takes into account the displacement of the solvent dielectric due to the other solute atoms. In pairwise solute descreening schemes this term is approximated by a pairwise sum48,49 over the volumes of the neighboring atoms, which are traditionally empirically adjusted to account for atomic overlaps. The AGBNP model instead makes use of a parameter-free geometrical algorithm to calculate the volume scaling coefficients used in the pairwise descreening scheme. The same algorithm is also used to calculate atomic surface areas. This feature is particularly advantageous in ligand binding applications when parametrizations of volume scaling coefficients are not available for some chemical groups. It has been shown that AGBNP gives excellent agreement for the GB self-energies and surface areas in comparison to accurate, but much more expensive, numerical evaluations.35 The nonpolar model adopted in this work differs from most other implicit hydration free energy models in that the nonpolar component Gnp of the hydration free energy is subdivided into cavity and solute-solvent van der Waals interaction terms Gnp ) Gcav + Gvdw

qiqj

∑ i 0) between the C descriptor, which is always negative (see Tables 3 and 4), and the binding free energy. This is indeed the case for MDL2 and MDL3 for the TIBO compounds. However the same models applied to the HEPT compounds yield a negative γ coefficient (see Table 5). For both ligand sets the value of γ for MDL2 and MDL3 is significantly smaller than for MDL1. It appears therefore that in these models the γ coefficient absorbs physical effects correlated to the size of the binding site which oppose binding and are not directly related to the hydrophobic effect. A likely candidate in this role is the reorganization free energy of the receptor, that is the work required to deform the binding site region to accommodate the ligand. It is reasonable to assume that the larger the binding site the larger the work required to deform the receptor structure. This contribution effectively reduces the benefit of having a large ligand-receptor contact surface area and could result in the smaller values for the γ LIE coefficients obtained for MDL2 and MDL3 as compared to MDL1, which adopts estimators that combine solution and receptor environments descriptors and is therefore less able to capture these effects. We hypothesize that for the HEPT compounds the conformational strain contribution overcomes the favorable hydrophobic desolvation effect, resulting in a negative value of γ. The prediction that the receptor reorganization free energy plays a more important role for the HEPT compounds than the TIBO compounds could also be a consequence of the wider distribution of ligand sizes of the HEPT compounds as measured by the ligand cavity formation descriptor which is proportional to the solute surface area. The range of variation of the CL descriptor of the HEPT compounds is 13.0 kcal/mol as compared 5.8 kcal/mol of the TIBO compounds (see Tables 3 and 4). The larger variation of receptor reorganization free energies embedded in the experimental binding free energies of the HEPT compounds makes the statistical fit of this term via linear regression more significant. Although deviations from ideality of LIE parameters can also be caused by the residual average electrostatic potential created by the receptor in absence of solute charges,42 we do not believe that this is the case in the present system. As eq 12 shows, linear response theory predicts that a nonzero residual electrostatic interaction energy (〈V〉0 in eq 12) tends to increase the value of the corresponding electrostatic LIE coefficient from its ideal value of 1/2. If the electrostatic LIE descriptor EC+2EL is negative (as for most HEPT compounds), this effect would appear, to an LIE model that

Su et al.

adopts an ideal LIE electrostatic coefficient of 1/2, as an unaccounted term faVorable to binding. It seems therefore unlikely that the residual electrostatic field of the receptor is the origin of the observed negative value of γ resulting from an unaccounted effect unfaVorable to binding, such as the reorganization free energy of the receptor. The results obtained with LIE models that do not include the factor of 2 for the implicit solvent electrostatic estimators are significantly inferior than the other models with the same number of adjustable parameters. It is predicted based on linear response theory that β should be half the value of β′. This is indeed approximately the case for models in which β (which corresponds to the ligand-receptor electrostatic interaction energy) is allowed to vary independently from β′ (the coefficient corresponding to the implicit solvent electrostatic estimator). Contrary to linear response predictions however β is always smaller than the theoretical value of 1/2, and β′ is always smaller than the theoretical value of 1. This is the case for the calculated R and ω parameters as well. It appears that the calculated parameters differ from their theoretical values by a constant proportionality factor. This would occur if a linear relationship with a proportionality coefficient different from 1 exists between the effective binding free energies calculated from the IC50’s as kT lnIC50 and the actual binding free energies. We plan to further investigate this issue by studying complexes for which direct measurements of the binding free energies have been reported. The IC50’s used in this work have been obtained from cell survival assays,53-55 rather than enzymatic rate inhibition assays which are more directly related to the binding affinity. In cell-based assays the IC50 is calculated from the number of cells in a colony that remain viable after infection with the HIV-1 virus in the presence of the inhibitor. A number of environmental factors such as cell absorption and metabolism and molecular factors62 such as stoichiometry of binding and induced RT dimerization could affect the observed effective binding free energies. The results for HIV-1 RT NNRTI systems similar to those studied here5,55 are in general agreement with our finding that, even though the optimal values of β and R are smaller than expected, their relative magnitude is consistent with linear response predictions. To support the hypothesis that the binding free energies calculated from the measured IC50’s are proportional to the actual binding free energies, we show in Table 6 the ratios between the β, β′, and ω parameters and the value of the R parameter obtained from a series of models (including those in Table 5). The theoretical values of these ratios based on linear response are β/R ) 1/2, β′/R ) 1, and ω/R ) 1. It can be clearly seen from Table 6 that the relative magnitudes of the LIE coefficients are in good agreement with the theoretical predictions even though their absolute values do not. Although this result indicates that the effective free energies of binding derived from the IC50’s deviate from the actual binding free energies by a constant proportionality factor, further investigations are needed to confirm this hypothesis. The constant proportionality factor between the theoretical and calculated LIE coefficients appears to be close to 1/2.



Table 6: Ratios between the LIE Parameter β, β′, and ω and the van der Waals LIE Parameter R Extracted from the Fits to the Experimental HEPT and TIBO Binding Data for the Model Listed Compared to Linear Response Theoretical Values HEPT set

TIBO set

model

β/R

β′/R

ω/R

β/R

β′/R

ω/R

theory MDL1 MDL1aa MDL2 MDL2ab MDL2bc

0.50 0.55 0.79 0.47 0.45 0.68

1

1

0.50 0.67 0.67 0.51 0.55 0.64

1

1

a

1.18

1.06

1.11 1.07 1.18

1.37

1.18

1.11 1.24 1.24

Eq 41. b Eq 42 with β′ ) 2β. c Eq 42.

The values of the R, β′, and ω coefficients are consistently near 0.5 (see Table 5), whereas the theoretical value for these parameters based on linear response is 1. Models in which ω is set to 1 perform significantly worse (data not shown) than MDL2 in which ω is set to 1/2 (Table 5). Similarly, the calculated value of the β coefficients is close to 1/4 as compared to the theoretical value of 1/2. Based on these observations we have attempted to construct a minimally parameterized LIE model for both the HEPT and TIBO compounds by setting the LIE coefficients to their theoretical linear response values reduced by a factor of 1/2. The remaining two parameters for which we do not have a firm theoretical prediction (the parameter γ corresponding to the cavity descriptor C, and the intercept δ) were allowed to vary to fit to the experimental binding free energies. This model (MDL3 in Table 5) is the most successful in predicting the binding free energies of the HEPT and TIBO compounds separately as measured by the jack-knife correlation coefficient Rpred2 and root-mean-square deviation RMSDpred. A noticeable difference between the LIE coefficients of model MDL3 for the HEPT and TIBO compounds is the value of γ, which is negative (γ ) -0.18) for the HEPT compounds and positive (γ ) 0.19) for the TIBO compounds. The high statistical significance of this difference is highlighted by the poor results obtained when fitting MDL3 to the combined HEPT and TIBO set. In this fit the γ coefficient is close to zero, an intermediate value between the value appropriate for the HEPT compounds and that appropriate for the TIBO compounds, and, as a result, the model is unable to fit accurately the experimental binding free energies. As discussed above we believe that this difference is due to the reorganization free energy of the receptor which opposes binding of the HEPT analogues more than the TIBO analogues, whereas hydrophobicity favors association relatively equally. Based on this result we predict that a hypothetical increase of ligand size, which leaves electrostatic and van der Waals receptor-ligand interactions unchanged, would have opposite effects on the binding affinities of the HEPT and TIBO analogues. For the TIBO analogues the hydrophobic gain due to the increase in the amount of receptor surface area buried by a larger ligand overcomes the opposing effect due the increase of receptor reorganization free energy, leading to stronger binding. For the HEPT analogues, instead, the hydrophobic gain is more

than offset by the increase of receptor reorganization free energy, leading to weaker binding. In practice, however, because it is not possible to modify the size of the ligand without also affecting electrostatic and van der Waals receptor-ligand interaction energies, all of the energetic components of the LIE model need to be considered to accurately predict the effect of ligand modifications. The success we obtained with MDL3 in predicting the binding free energies of the HEPT and TIBO compounds based on theoretically derived parameters and only one adjustable LIE coefficient (excluding the intercept) raises the question of whether this result is a direct consequence of the receptor and solution environments obeying linear response or simply a coincidental occurrence for the receptor system and ligand sets we analyzed. The results of previous LIE studies on a variety of ligand binding systems are mainly inconclusive on this issue, partly due to difficulty of interpreting results obtained with different LIE regression expressions. A LIE explicit solvent study of P450cam complexes by Paulsen and Ornstein63 has shown that values of LIE coefficients near to their theoretical values of R ) 1 and β ) 1/2 reproduced well the experimental binding free energies, whereas in a related study Almlo¨f et al.44 reproduced the binding free energies of a similar set of inhibitors of P450cam with a significantly smaller value of the R LIE coefficient. As noted by Amlo¨f et al.44 the difference between the values of R between the two studies is due to the intercept parameter, considered only in the latter study, whose optimal value was found to depend on the hydrophobicity of the receptor site.44 Wang et al.64 have investigated different values of the van der Waals parameter R in conjunction with β set to 1/2 without and intercept parameter. They found that for various ligand-protein complexes the optimal value of R varied from 1 to e 0.1 depending on the hydrophobicity of the binding pocket. In a related study Wang et al.65 have compared explicit solvent LIE predictions to rigorous free energy thermodynamic integration (TI) calculations for complexes with streptavidin and established that for these systems values of R and β near their theoretical values reproduced well the experimental binding free energies for ligands for which LIE and TI produce consistent results. Cases in which TI produced results in better agreement with the experiments could be rationalized by the inadequacy of the LIE model to properly take into account the free energy cost for reorganizing the receptor pocket.65 Earlier LIE studies2,9 reported smaller optimal values for the R LIE parameter when setting β to 1/2. Clearly, much remains to be done to better understand the factors that influence the values of the LIE coefficients obtained on different systems with different LIE regression equations. It is conceivable that LIE regression models, such as those adopted in the present work, which include an estimator for the work of cavity formation3,4 are more likely to yield electrostatic and van der Waals LIE coefficients consistent with linear response predictions. The values of electrostatic and especially van der Waals LIE coefficients obtained using LIE models without an explicit cavity formation estimator,44 and especially those without an intercept parameter,63,65 are


instead more likely to absorb effects that are not directly related to electrostatic and van der Waals interactions and are thus not as amenable to rationalization and generalization in terms of linear response theory. Furthermore, because solute-solvent van der Waals interaction energies are not perfectly correlated with the solute surface area,23,66,67 surface area-based cavity estimators provide nonredundant information content in addition to the van der Waals estimator, potentially leading to better descriptive accuracy and transferability. However, although it has been generally recognized that the free energy of cavity formation in water approximately scales as the surface area of the solute,23,38 the validity of using a surface area estimator to describe the work of ligand cavity formation in the receptor pocket remains to be fully addressed.

7. Conclusions In this paper we have addressed a number of outstanding issues in regard to the theory and practice of LIE modeling for binding free energy prediction. We have reviewed and clarified the linear response theory on which LIE methods are based. Following linear response formalism, we derived expressions for the LIE estimators when the solvent is treated implicitly. We showed that these estimators include descriptors related to the desolvation of receptor atoms which were not considered in a previously reported implicit solvent LIE model.4 The form of the estimators we derived are consistent with those proposed previously in the context of the LRA method and the PDLD implicit solvent models10 and for the LIE electrostatic estimator with a GB model.12 We have also developed a novel class of LIE models that, contrary to current practice, attempt to model independently the processes of insertion of the ligand in the receptor and in solution. These models are motivated by the fact that, potentially, the receptor and the solvent respond differently to the introduction of the ligand. We have applied these ideas to the problem of the binding free energy estimation of a series of NNRTI inhibitors of HIV-1 RT. LIE descriptors were collected from 57 molecular dynamics simulations of HIV-1 RT complexed with 20 HEPT inhibitors and 37 TIBO inhibitors. Based on the measured binding affinities and the calculated descriptors we developed a series of LIE models. We presented results for three of these models and tested several others. The first model, MDL1, is comparable to previously reported LIE studies for NNRTI binding in both explicit and implicit solvent. Relative to these studies, the predictive accuracy of our MDL1 model is generally superior when applied to the same ligand sets, suggesting that the AGBNP implicit solvation model provides sufficient accuracy for LIE modeling. This result also indicates that the LIE estimators we derived are more appropriate to describe the energetics of the binding process in implicit solvent than previously reported alternatives. The second and third models, MDL2 and MDL3, are novel models designed to treat the hydration free energy of the ligand and the work of inserting the ligand in the receptor independently. MDL2, which has the same number of parameters as MDL1, is found to be superior to MDL1 and to LIE models developed by others for predicting

Su et al.

the binding affinities of the HEPT and TIBO analogues. This result may indicate that LIE models that treat the two insertion processes independently can lead to better prediction accuracy. MDL3, a minimally parametrized version of MDL2 in which some parameters are set to their values predicted by linear response, is found to be superior to MDL1 and MDL2 in terms of predictive ability (jack-knife tests) for the HEPT and TIBO sets individually but not for the two sets combined. We hypothesize that this is caused by the larger sizes of some of the HEPT compounds which induce steric strain of the receptor; an effect which is not taken explicitly into account by the LIE models. We examined the values of the LIE coefficients obtained from the regression analysis and established that, although they are smaller than expected, their relative magnitudes generally conform to linear response predictions. We are planning calculations to test the applicability of linear response to other protein-ligand binding systems by computing directly the relative free energies of inserting ligands in solution and the protein environment and comparing the results to the corresponding first order cumulant descriptor. The discovery that linear response behavior is generally applicable to protein-ligand binding could provide the basis for the development of minimally parametrized LIE models. Given the substantial reduction of adjustable parameters achievable when assuming linear response for some of the interaction energy contributions and the potential improvement in predictive ability, minimally parametrized models based on linear response, such as the ones presented here, offer a productive route to binding free energy predictions using sparse experimental binding assay data for lead optimization in structure-based drug design. Acknowledgment. This work has been supported in part by a grant from the National Institutes of Health (GM30580). We thank Dr. Anthony Felts and Dr. Zhiyong Zhou for helpful discussions.

Appendix In this Appendix we derive eq 30. Let us first consider the process of turning on the ligand charges in the receptor environment. Conceptually we will divide this process into two steps. First the electrostatic ligand-receptor interactions are turned on and then interactions between the ligand and the implicit solvent continuum are turned on. If we assume that the receptor responds to the perturbation as a perfect linear dielectric, the free energy change, ∆Fel(1), corresponding to the first step is approximately ∆Fel(1) = β〈Vel〉1

(46)

where 〈Vel〉1 is the average of the ligand-receptor Coulomb interaction energy in the absence of ligand-continuum solvent electrostatic interactions, and β ) 1/2 based on eq 13 under the present assumptions. Based on linear response the free energy corresponding to turning on the ligand-implicit solvent continuum electrostatic interactions is ∆Fel(2) ) β′〈Gcel - Gc′el〉c

(47)



where Gcel is the Generalized Born energy of the receptorligand complex and Gc′el is the Generalized Born energy of complex when the ligand charges are turned off, and the average is taken in the ensemble in which the ligand is fully charged. To derive eq 47 we write the λ-dependent potential for this process as U(λ) ) U0 + λV, where U0 ) Uexpl + Gcav + Gvdw + Gc′el is the reference potential and V ) Gcel Gc′el is the perturbation. When the ligand and the receptor are both assumed rigid, the free energy change for adding electrostatic solute-continuum solvent interactions reduces to ∆F(2)el ) Gcel - Gc′el, which is the same as eq 47 with β′ ) 1 and with the average replaced by the constant value of the argument. Hence, based on the discussion presented above, in the general case when the ligand and the receptor are flexible and linear response applies, we expect the optimal value of the β′ coefficient in eq 47 to assume values near 1. Thus, assuming that in eq 46 the mean ligand-receptor electrostatic interaction energy 〈Vel〉 is the same whether the ligand-continuum solvent electrostatic interactions are present or not and assuming that β′ ) 2β, eqs 46 and 47 can be combined to give the linear response expression for the free energy of turning on the ligand charges in the receptor environment ∆Fel = β〈Vel + 2(Gcel - Gc′el)〉c

(48)

where the average is taken with the ligand fully interacting with the receptor and solution environments, and c′ corresponds to the state in which the ligand is uncharged. It is straightforward to show using eq 18 that the difference between Gcel and Gc′el for a given ligand-receptor complex conformation is the sum of the GB self-energies of the ligand atoms plus the sum of the GB pair energies between ligand atoms and between ligand atoms and receptor atoms. The expression for the LIE estimator for the free energy of adding the van der Waals ligand-receptor and ligandwater interactions to the ligand cavity involves the attractive part of the Lennard-Jones potential and the implicit solvent solute-solvent van der Waals interaction energy Gvdw. Following the same derivation as for the electrostatic case and assuming R = 1 for both the explicit and implicit contributions, we obtain c′ c′′ ∆Fvdw ) R〈Vvdw LJ + (Gvdw - Gvdw)〉c′

(49)

where Vvdw LJ is the sum of the attractive components of the Lennard-Jones interactions68 between ligand and receptor atoms, Gc′vdw is the solute-solvent van der Waals interaction energy of the receptor-ligand complex, Gc′′ vdw is the solutesolvent van der Waals interaction energy of the complex in the absence of van der Waals interactions between the ligand and the solvent, and the average is taken in the state c′ with the uncharged ligand and with full ligand-receptor and ligand-solvent van der Waals interactions. We now consider the work of creating the ligand cavity within the receptor site. According to the scheme developed above we should consider as an LIE estimator for this process the average of the difference between the solvent potential of mean force in the presence of the solute cavity and in the absence of the solute cavity plus repulsive interactions

between the ligand atoms and receptor atoms. As above the average is assumed to be taken over the ensemble of conformations generated when the ligand cavity is present. The solvent potential of mean force difference includes two components: the change in Gcav in going from the receptor without the ligand cavity and the receptor with the ligand cavity, and, in addition, the change of the receptor-solvent van der Waals interaction energy and the change of Generalized Born energy of the receptor caused by the increase of the receptor atoms’ Born radii due to the introduction of the ligand cavity. The explicit ligand cavity-receptor interactions are modeled using the repulsive component of the LennardJones potential according to the WCA decomposition.68 Finally, assuming that linear response applies to processes involving Vrep LJ and the van der Waals and electrostatic implicit solvent components with the same proportionality coefficients as in eqs 48 and 49, we obtain the following expression for introducing the ligand cavity in the receptor c′′ p ∆Fcav ) R〈Vrep LJ + (Gvdw - Gvdw)〉c′′ + p c′′ p β〈2(Gc′′ el - Gel)〉c′′ + γ〈Gcav - Gcav〉c′′ (50)

where Vrep LJ is the sum of repulsive Lennard-Jones interacc′′ tions between the ligand and receptor atoms, Gc′′ vdw and Gel are, respectively, the solute-solvent van der Waals interaction energy and the Generalized Born energy of the complex with the ligand cavity with the ligand-environment van der Waals and electrostatic interactions turned off. Gpvdw and Gpel are, respectively, the solute-solvent van der Waals interaction energy and the Generalized Born energy of the receptor conformation obtained from the conformation of the complex after removal of the ligand, Gc′′ cav is the cavity free energy of the receptor-ligand complex, Gpcav is the cavity free energy of the complex after removal of the ligand, and 〈‚‚‚〉c′′ indicates averaging over complex conformations with the ligand cavity. When considering the three estimators from eqs 48-50 we now make the approximation to evaluate all of the averages in the final state in which the ligand fully interacts with the environment (state c). We choose to combine the vdw repulsive and attractive WCA components, Vrep LJ and VLJ , of the receptor-ligand Lennard-Jones interaction energies to form the total ligand-receptor Lennard-Jones interaction energy, VLJ. Also, when combining the cavity and electroc′′ static descriptors the Gc′′ vdw and Gel terms cancel out. We then obtain the following expression for the LIE regression equation for the free energy for creating the ligand in the receptor site ∆Fc = R〈VLJ + Gcvdw - Gpvdw〉c + β〈Vel + 2(Gcel - Gpel)〉c + γ〈Gccav - Gpcav〉c (51) which is eq 30. References (1) Free Energy Calculations in Rational Drug Design; Reddy, M. R., Erion, M. D., Eds.; Springer-Verlag: 2001. (2) A˙ qvist, J.; Medina, C.; Samuelsson, J.-E. Protein Eng. 1994, 7, 385-391.


Su et al.

(3) Jones-Hertzog, D. K.; Jorgensen, W. L. J. Med. Chem. 1997, 40, 1539-1549.

(23) Gallicchio, E.; Kubo, M. M.; Levy, R. M. J. Phys. Chem. B 2000, 104, 6271-6285.

(4) Zhou, R.; Friesner, R. A.; Ghosh, A.; Rizzo, R. C.; Jorgensen, W. L.; Levy, R. M. J. Phys. Chem. 2001, 105, 10388-10397.

(24) Pierotti, R. A. Chem. ReV. 1976, 76, 717-26.

(5) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265-322. (6) Levy, R. M.; Belhadj, M.; Kitchen, D. B. J. Chem. Phys. 1991, 95, 3627-3633. (7) Figueirido, F.; Del Buono, G.; Levy, R. M. Biophys. Chem. 1994, 51, 235-241. (8) Lee, F. S.; Chu, Z. T.; Bolger, M. B.; Warshel, A. Protein Eng. 1992, 5, 215-228.

(25) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570-4577. (26) Hummer, G.; Garde, S.; Garcı´a, A. E.; Paulaitis, M. E.; Pratt, L. R. Phys. Chem. B 1998, 102, 10469-82. (27) Carlson, H. A.; Jorgensen, W. L. J. Phys. Chem. 1995, 99, 10667-10673. (28) Warshel, A.; Russell, S. T. Q. ReV. Biophys. 1984, 17, 283422.

(9) Hansson, T.; A° qvist, J. Protein Eng. 1995, 8, 1137-1144.

(29) A° qvist, J.; Hansson, T. J. Phys. Chem. 1996, 100, 95129521.

(10) Sham, Y. Y.; Chu, Z. T.; Tao, H.; Warshel, A. Proteins: Struct., Funct., Genet. 2000, 39, 393-407.

(30) Cortis, C. M.; Friesner, R. A. J. Comput. Chem. 1997, 18, 1591-1608.

(11) Chen, X.; Tropsha, A. J. Chem. Theory Comput. 2006, 2, 1435-1443.

(31) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. J. Comput. Chem. 2002, 23, 128137.

(12) Carlsson, J.; Ande´r, M.; Nervall, M.; A° qvist, J. J. Phys. Chem. B 2006, 110, 12034-12041. (13) De Clercq, E. AntiViral Res. 1998, 38, 153-179. (14) Coffin, J. M.; Hughes, S. H.; Varmus, H. E. RetroViruses at the National Center for Biotechnology Information; Cold Spring Harbor: Cold Spring Harbor Laboratory: 1997. (15) Janssen, P. A. J.; Lewi, P. J.; Arnold, E.; Daeyaert, F.; de Jonge, M.; Heeres, J.; Koymans, L.; Vinkers, M.; Guillemont, J.; Pasquier, E.; Kukla, M.; Ludovici, D.; Andries, K.; de Be`thune, M.-P.; Pawels, R.; Das, K.; Clark Jr., A. D.; Frenkel, Y. V.; Hughes, S. H.; Medaer, B.; De Knaep, F.; Bohets, H.; De Clerck, F.; Lampo, A.; Williams, P.; Stoffels, P. J. Med. Chem. 2005, 48, 1901-1909. (16) Heeres, J.; de Jonge, M.; Koymans, L. M. H.; Daeyaert, F. F. D.; Vinkers, M.; Van Aken, K. J. A.; Arnold, E.; Das, K.; Kilonda, A.; Hoornaert, G. J.; Compernolle, F.; Cegla, M.; Azzam, R. A.; Andries, K.; de Be`thune, M.-P.; Azjin, H.; Pauwels, R.; Lewi, P. J.; Janssen, P. A. J. J. Med. Chem. 2005, 48, 1910-1918. (17) Das, K.; Clark, A. D., Jr.; Lewi, P. J.; Heeres, J.; De Jonge, M. R.; Koymans, L. M.; Vinkers, H. F.; Daeyaert, F.; Ludovici, D. W.; Kukla, M. J.; De Corte, B.; Kavash, R. W.; Ho, C. Y.; Ye, H.; Lichtenstein, M. A.; Andries, K.; Pauwels, R.; De Bethune, M. P.; Boyer, P. L.; Clark, P.; Hughes, S. H.; Janssen, P. A.; Arnold, E. J. Med. Chem. 2004, 47, 2550-2560. (18) Das, K.; Lewi, P. J.; Hughes, S. H.; Arnold, E. Prog. Biophys. Mol. Biol. 2005, 88, 209-231. (19) Rizzo, R. C.; Udier-Blagovic´, M.; Wang, D.-P.; Watkins, E. K.; Kroeger Smith, M. B.; Smith, R. H., Jr.; Tirado-Rives, J.; Jorgensen, W. L. J. Med. Chem. 2002, 45, 2970-2987. (20) Kroeger Smith, M. B.; Hose, B. M.; Hawkins, A.; Lipchock, J.; Farnsworth, D. W.; Rizzo, R. C.; Tirado-Rives, J.; Arnold, E.; Zhang, W.; Hughes, S. H.; Jorgensen, W. L.; Michejda, C. J.; Smith, R. H., Jr. J. Med. Chem. 2003, 46, 1940-1947. (21) Himmel, D. M.; Sarafianos, S. G.; Dharmasena, S.; Hossain, M. M.; McCoy-Simandle, K.; Ilina, T.; Clark Jr., A. D.; Knight, J. L.; Julias, J. G.; Clark, P. K.; Krogh-Jespersen, K.; Levy, R. M.; Hughes, S. H.; Parniak, M. A.; Arnold, E. 2006, submitted for publication. (22) Floria´n, J.; Goodman, M. F.; Warshel, A. J. Phys. Chem. B 2002, 106, 5739-5753.

(32) Lazaridis, T.; Karplus, M. Curr. Opin. Struct. Biol. 2000, 10, 139-145. (33) Bashford, D.; Case, D. A. Annu. ReV. Phys. Chem. 2000, 51, 129-152. (34) Gallicchio, E.; Zhang, L. Y.; Levy, R. M. J. Comput. Chem. 2002, 23, 517-529. (35) Gallicchio, E.; Levy, R. M. J. Comput. Chem. 2004, 25, 479499. (36) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 36833704. (37) Huang, D. M.; Chandler, D. J. Phys. Chem. B 2002, 106, 2047-2053. (38) Ashbaugh, H. S.; Kaler, E. W.; Paulaitis, M. E. J. Am. Chem. Soc. 1999, 121, 9243-9244. (39) CRC Handbook of Chemistry and Physics, 86th ed. ed.; Lide, D. E., Ed.; CRC Press: Boca Raton, FL, 2005. (40) Tounge, B. A.; Reynolds, C. H. J. Med. Chem. 2003, 46, 2074-2082. (41) Singh, P.; Mhaka, A. M.; Christensen, S. B.; Gray, J. J.; Denmeade, S. R.; Isaacs, J. T. J. Med. Chem. 2005, 48, 3005-3014. (42) Almlo¨f, M.; A° qvist, J.; Smålas, A. O.; Brandsdal, B. O. Biophys. J. 2006, 90, 433-442. (43) Zhang, L.; Gallicchio, E.; Friesner, R. A.; Levy, R. M. J. Comput. Chem. 2001, 22, 591-607. (44) Almlo¨f, M.; Brandsdal, B. O.; A° qvist, J. J. Comput. Chem. 2004, 25, 1242-1254. (45) Still, W. C. Tempczyk, A.; Hawley, R. C.; Hendrikson, T. J. Am. Chem. Soc. 1990, 112, 6127-6129. (46) Schaefer, M.; Karplus, M. J. Phys. Chem. 1996, 100, 15781599. (47) Onufriev, A.; Bashford, D.; Case, D. A. J. Phys. Chem. B 2000, 104, 3712-3720. (48) Schaefer, M.; Froemmel, C. J. Mol. Biol. 1990, 216, 10451066. (49) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 19824-19839.



(50) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225-11236.

L. Y.; Berne, B. J.; Friesner, R. A.; Gallicchio, E.; Levy, R. M. J. Comput. Chem. 2005, 26, 1752-1780.

(51) Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381.

(60) Smith, M. B. K.; Rouzer, C. A.; Taneyhill, L. A.; Smith, N. A.; Hughes, S. H.; Boyer, P. L.; Janssen, P. A. J. Moereels, H.; Koymans, L.; Arnold, E.; Ding, J.; Das, K.; Zhang, W.; Michejda, C. J.; Smith, R. H., Jr. Protein Sci. 1995, 4, 22032222.

(52) Simonson, T. Curr. Opin. Struct. Biol. 2001, 11, 243-252. (53) Rizzo, R. C.; Tirado-Rives, J.; Jorgensen, W. L. J. Med. Chem. 2001, 44, 145-154. (54) Ho, W.; Kukla, M. J.; Breslin, H. J.; Ludovici, D. W.; Grous, P. P.; Diamond, C. J.; Miranda, M.; Rodgers, J. D.; Ho, C. Y.; De Clercq, E.; Pauwels, R.; Andries, K.; Janssen, M. A. C.; Janssen, P. A. J. J. Med. Chem. 1995, 38, 794-802. (55) Smith Jr., R. H.; Jorgensen, W. L.; Tirado-Rives, J.; Lamb, M. L.; Janssen, P. A. J.; Michejda, C. J.; Kroeger Smith, M. B. J. Med. Chem. 1998, 41, 5272-5286. (56) Hopkins, A. L.; Ren, J.; Esnouf, R. M.; Willcox, B. E.; Jones, E. Y.; Ross, C.; Miyasaka, T.; Walker, R. T.; Tanaka, H.; Stammers, D. K.; Stuart, D. I. J. Med. Chem. 1996, 39, 1589-1600. (57) Das, K.; Ding, J.; Hsiou, Y.; Clark A. D., Jr.; Moereels, H.; Koymans, L.; Andries, K.; Pauwels, R.; Janssen, P. A.; Boyer, P. L.; Clark, P.; Smith, R. H., Jr.; Kroeger Smith, M. B.; Michejda, C. J.; Hughes, S. H.; Arnold, E. J. Mol. Biol. 1996, 264, 1085-1100. (58) Schro¨dinger, Inc., Portland, OR. (59) Banks, J. L.; Beard, J. S.; Cao, Y.; Cho, A. E.; Damm, W.; Farid, R.; Felts, A. K.; Halgren, T. A.; Mainz, D. T.; Maple, J. R.; Murphy, R.; Philipp, D. M.; Repasky, M. P.; Zhang,

(61) Ghosh, A.; Rapp, C. S.; Friesner, R. A. J. Phys. Chem. B 1998, 102, 10983-10990. (62) Cheng, Y.-C.; Prusoff, W. H. Biochem. Pharmacol. 1973, 22, 3099-3108. (63) Paulsen, M. D.; Ornstein, R. L. Protein Eng. 1996, 9, 567571. (64) Wang, W.; Wang, J.; Kollman, P. A. Proteins: Struct., Funct., Genet. 1999, 34, 395-402. (65) Wang, J.; Dixon, R.; Kollman, P. A. Proteins Struct. Funct. Genet. 1999, 34, 69-81. (66) Su, Y.; Gallicchio, E.; Levy, R. M. Biophys. Chem. 2004, 109, 251-260. (67) Levy, R. M.; Zhang, L. Y.; Gallicchio, E.; Felts, A. K. J. Am. Chem. Soc. 2003, 25, 9523-9530. (68) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237-47. CT600258E

Linear Interaction Energy (LIE) Models for Ligand Binding in Implicit

Recommend Documents