Cluster Hydration Model for Binding Energy Calculations of Protein

J. Phys. Chem. B 2009, 113, 809–817

809

Cluster Hydration Model for Binding Energy Calculations of Protein-Ligand Complexes Katsumi Murata,*,† Dmitri G. Fedorov,‡ Isao Nakanishi,†,§ and Kazuo Kitaura†,‡ Department of Theoretical Drug Design, Graduate School of Pharmaceutical Sciences, Kyoto UniVersity, Sakyo-ku, Kyoto, 606-8501 Japan, and National Institute of AdVanced Industrial Science and Technology (AIST), Tsukuba, Ibaraki, 305-8568 Japan ReceiVed: June 6, 2008; ReVised Manuscript ReceiVed: October 5, 2008

We present a new model for predicting the binding affinity in protein-ligand complexes based on an explicit solvent treatment. The model is referred to as the cluster hydration model. Test calculations were performed on complexes of FK506-binding protein (FKBP) and its six ligands. The calculated binding energies had a good correlation with experimental binding affinities; the correlation coefficient r2 was 0.93. Moreover, we examined the competition between the residue-ligand interactions and the desolvation energies of residues. The results suggested that, for the stabilization of the complexes in solution, the contributions of nonpolar residues at the binding site are larger than those of charged and polar ones because the electrostatic interaction energies between the residues and the ligand were mostly canceled by the desolvation penalties. Introduction The theoretical prediction of the protein-ligand binding affinities is essential in a wide range of biophysical issues, from molecular recognition to rational drug design. A great deal of effort has been invested in this area, and several approaches have been developed, with a different tradeoff between accuracy and computer time. The most rigorous approaches are the freeenergy perturbation (FEP) and the thermodynamic integration (TI).1,2 However, their applications are limited to the calculation of free-energy differences between variants with a little structural difference of ligands and still require an enormous CPU time. Extremely fast QSAR-type approaches, which are well suited for high-throughput in silico screening of compounds, usually include many empirical parameters, and their reliability varies on the case-by-case basis.3-6 The linear interaction energy (LIE) method and the combined molecular mechanics and Poisson-Boltzmann surface area solvation model (MM/PBSA) are a compromise between accuracy and speed. LIE has been successfully applied to a number of protein-ligands systems.7-11 In the method, the binding free energy is evaluated as a sum of the electrostatic and van der Waals (vdW) interaction energy terms. On the basis of the linear response approximation, the coefficient of 0.5 is used for the electrostatic term, while the coefficient of the van der Waals term is an adjustable parameter, determined to reproduce experimental values. MM/PBSA is also widely used for the calculations of the binding free energy, in which the average interaction energy obtained from MD simulations is used, and the solvation energy is estimated by the PBSA solvation model.12 Despite its attractive simplicity, this method ignores many essential features of the solvent because of the absence of an explicit molecular description. In this work, we present a new model, the cluster hydration model, for the calculation of binding energies of protein-ligand * To whom correspondence should be addressed. E-mail: murata@ pharm.kyoto-u.ac.jp. † Kyoto University. ‡ National Institute of Advanced Industrial Science and Technology (AIST). § Present address: Kinki University.

complexes in solution. In this model, we focus on the conglomerate (cluster) of solvent molecules, which replaces the protein (ligand) in bulk solvent during the protein-ligand complex formation. To test this solvation model, we applied it to FK506 binding protein (FKBP), which is a small and fairly rigid protein. FKBP is an immunophilin that blocks early T-cell activation via calcineurin inhibition upon binding to the immunosuppressive drug FK506. The calculated binding energies had a good correlation with experimental binding affinities. Moreover, we examined the competition between the protein-ligand interaction and partial desolvation accompanying the complex formation. Method Section Cluster Hydration Model. The scheme of the cluster hydration model (CHM) is given in Figure 1. During the complexation of a protein (P) and a ligand (L) in aqueous solution, water molecules occupying their corresponding binding regions are released into bulk solvent. The water cluster in the protein binding region, which is replaced by the ligand in the (P) complex (P · L), is defined as the hydration water cluster WL (P) and is considered to form the complex P · WL with the protein. Similarly, the water cluster around the binding site of the ligand, (L) which is replaced by the protein in P · L, is defined as WP and (L) is considered to form the complex WP · L with the ligand. The (L) (P) (B) (B) released WP and WL become WP andWL , respectively, in (B) (B) (L) (B) bulk water and form the complex WP · WL . WP and WP have a similar shape but a little different water molecule configuration because the former is somewhat deformed in the protein binding (P) (B) pocket. WL and WL have also a little different conformation of solvent molecules. In addition, the sum of the solvation (P) (L) energies of the isolated state (P · WL and WP · L) is assumed to be exactly equal to the sum of the solvation energies of the (B) (B) complex state (P · L and WP · WL ), so that their difference is 0. This assumption is justified because their solvent-accessible surface areas and surface characteristics are similar by construction. Thus, our model avoids the numerical calculations of the entire solvation energies of solutes, which are known to be very difficult to estimate by MD simulations. In CHM, the association process of the protein and ligand in aqueous solution is described as

10.1021/jp805007f CCC: $40.75  2009 American Chemical Society Published on Web 12/31/2008

810 J. Phys. Chem. B, Vol. 113, No. 3, 2009

Murata et al.

P · WL(P) + WP(L) · L f P · L + WP(B) · WL(B) (L)

(P)

(P)

(1) (B)

The hydration water clusters, WP , WL , WL , and WL , are defined based on the hydrated structures of protein and ligand and the structure of bulk water obtained from MD simulations described in the next section. All systems appearing in eq 1 are immersed in bulk water. The binding energy, ∆Ebind, between the protein and the ligand is calculated as follows int def int j P,L jW j int j int ∆Ebind ) E +E (B),W(B) - EP,W(P) - EW(L),L + ∆EL P L L P (2) int j X,Y means the average interaction energy between X where E and Y. The partial desolvation accompanying the complex formation is estimated from the dissociation of the complexes composed of the solute (P or L) and water clusters followed by the recombination of these water clusters. The last term of eq 2 is the deformation energy of the ligand

internal j L(P j internal ∆ELdef ) E · L) - EL(free)

every 50 ps in a total 2 ns run described above, and 80 (40 × 2) separate simulation runs for each structure of the protein and ligand were performed for 200 ps after 200 ps equilibration runs, where the rigid protein or ligand was placed in an 81 Å × 81 Å × 81 Å water box with the same conditions as those described above. The complex, from which the protein and ligand were extracted, was superimposed on the protein (or the ligand) in the equilibrated water boxes, and the clusters of water molecules (P) covered by the ligand (or the protein) were defined asWL (L) (orWP ). The criterion to include a water molecule was that the oxygen atom of the water molecule was within the van der Waals radius of any atom of the solute molecule. Therefore, (L) (P) the volumes of water clusters, WP and WL , are nearly equal to those of P and L, respectively, as well as their solventaccessible surface areas (data not shown). (B) (B) For the water clusters WP and WL , 5000 snapshots were extracted from MD trajectory at every 400 fs for each complex. The extracted protein-ligand complexes were superimposed on the equilibrated pure water box, and water molecules covered (B) (B) by the protein and the ligand were defined as WP and WL , respectively, where the same criterion as that described above was applied to pick up water molecules.

(3)

internal j L(P j internal where E · L) and EL(free) are the average internal energies of the ligand in the complex P · L and that in the isolated state in solution, respectively, and they are calculated from MD simulations for the corresponding states. Note that in this work we assume that the deformation energy of the protein is negligible because FKBP is a fairly rigid protein. In practice, the neglect of the deformation energy of the protein cancels the noise that would result from sampling insufficiency and the error inherent in force field energies. Simulation Protocol and Definition of Hydration Water Clusters. FKBP complexes with six ligands L1-L6 (Figure 2) were selected to test CHM, and we performed MD simulations to define the hydration water clusters used in the model. Initial geometries for the simulations were taken from the Protein Data Bank (PDB code, L1: 1FKB; L2: 1FKF; L3: 1FKG; L4: 1FKI).13-15 The geometries of the L5 and L6 complexes were not reported; therefore, they were modeled using the L3 complex as a template. The geometries of these ligands were optimized at the Hartree-Fock level with the 6-31G* basis set using Gaussian03,16 and RESP charges were calculated and assigned with the ANTECHAMBER module from AMBER8.17 First, rough geometry optimizations (a gradient norm less than 0.1 kcal mol-1 Å-1) were performed on the six complexes with the SANDER module from AMBER8. Then, the complexes were placed in the 81 Å × 81 Å × 81 Å water box, whose density was set to 1 g/cm3. In all MD simulations, TIP3P for water,18 the ff99 force field for the protein,19 and the gaff force field for the ligand20 were employed. For the six complexes, we separately performed MD simulations of 2 × 106 time steps (or 2 ns total) under constant NVT conditions and, before collecting data, also performed MD simulations of 400 ps for equilibration. The temperature during MD simulations was controlled by using the Berendsen algorithm21 and was set to 298 K. The SHAKE algorithm22 was applied to all bonds involving hydrogen atoms, and the unit time step was set to 1 fs. In addition, MD simulations were performed for pure water under the same conditions to define (B) (B) WP and WL . (L) (P) To define WP and WL , the structures of the protein and ligand were extracted from the complex in the 40 snapshots at

Results and Discussion Hydration Water Clusters. In Table 1, the average number (L) (P) of water molecules of the hydration water clusters WP , WL , (B) (B) WP , and WL are listed for all complexes (L1-L6). The values (L) of WP are much the same, which is reasonable because of the (P) same protein in all of the complexes. The values of WL follow the order of the ligand size (L1 > L2 > L3 ≈ L4 > L5 ≈ L6). (B) (B) The values of WP and WL are essentially the same as those (L) (P) (L) (B) of WP and WL , respectively. WP is smaller than WP by two to three molecules. This is because the ligand distorts (destabilizes) the water conglomerate bound to it, reducing the density of water molecules (since the volume is chosen to be nearly the same). (L) (P) We note that the average number of atoms of WP and WL is much less than the number of atoms in the protein (P) and (L) ligands (L), respectively, 894 in WP versus 1666 in P and 66 (P) in WL versus 126 in L (in the case of L2). The surface density of atoms in the water clusters is much less than those in the corresponding molecules (P or L). This is because covalently bonded atoms in a molecule are closer than nonbonded atoms in a molecular (water) cluster. This fact might be important to account for why the protein-ligand vdW interaction energy can compete against the desolvation vdW energy, which is the sum of the protein-solvent and ligand-solvent vdW interaction energies (vide infra). The probability distributions of the number of water molecules in the clusters corresponding to the protein and ligands are shown in Figure 3a and b, respectively, where only L2, L4, and L6 complexes are presented as examples. The probability (B) (B) distributions of WP and WL are essentially the same as those (L) (P) of WP and WL , respectively, and are not shown. The (L) (P) probability distributions of both WP and WL appear to be canonical, indicating that the simulations were properly per(L) formed. The probability distribution of WP for all ligands is very similar (Figure 3a) because the protein is the same. The (P) peak of the probability distribution of WL (Figure 3b) shifts to a larger number, and the distribution becomes broader as the ligand size or ligand volume increases from L6 to L2. This is reasonable because more configurations of water molecules are included in the larger volume clusters.

Binding Energy Calculations of Protein-Ligand Complexes

J. Phys. Chem. B, Vol. 113, No. 3, 2009 811

(P)

(L)

(B)

(B)

Figure 1. Schematic illustration of the cluster hydration model (CHM). Water clusters enclosed by dotted lines (WL , WP , WL , and WP ) are (P) (P) (L) (L) considered as units forming complexes. WL in the “P · WL complex” (a) and WP in the “WP · L complex” (b) dissociate when the complex (P) (B) formation of P and L occurs (c), and they form the “WL · WP complex” (d).

Figure 3. Probability distributions of the number of water molecules (L) (P) in (a) WP and (b) WL . The solid, dashed, and dotted curves show the results for L2, L4, and L6, respectively.

Figure 2. Structures of six FKBP ligands.

TABLE 1: Average Number of Molecules in Hydration (L) (B) (P) (B) Water Clusters WP , WP , WL , and WL (L)

(B)

(P)

(B)

ligand

WP

WP

WL

WL

L1 L2 L3 L4 L5 L6

297.7 297.9 297.8 297.5 297.9 297.5

300.9 300.0 300.2 299.6 299.0 299.4

25.2 22.1 12.1 12.3 7.9 8.1

25.3 22.9 13.0 13.8 8.5 9.1

Binding Energy. The calculated binding energies and each term of eq 2 are presented in Table 2 for all of the ligands. Throughout this paper, the negative energy means stabilization, and the positive one denotes destabilization. Considering enthalpies and omitting the entropy contribution to the energetics, we aim at the correlation with the experimental binding free energy and not so much at the absolute value reproduction. A good correlation (r2 ) 0.93) was obtained between the calculated energies and experimental binding free energies (Figure 4). We found that the conformation of ligand L3 is considerably distorted in the complex. The two phenyl groups are stacked in

the parallel-displaced conformation in the isolated state, but inside of the protein-ligand complex, one phenyl group of L3 is stacked with the aromatic ring of Tyr82 in the T-shaped edge facing the protein (Figure 5a and b). The energy difference between the two conformations of L3 was estimated to be ∼2 kcal/mol. The geometric distortion of the other ligands in the complexes is small, and their deformation energies are within 0.4 kcal/mol. Thus, the difference of deformation energies among the ligands cannot be ignored for estimating the relative binding energy. It is encouraging that an excellent linear correlation was obtained between the binding energies calculated by CHM and experimental binding affinities. However, our model does not consider the entropy contribution, which is indispensable for an accurate prediction of the binding affinities. In MM/PBSA,24-26 the vibrational entropy of the solute is estimated by the harmonic oscillator approximation or the quasi-harmonic approximation. Yet, the validity and accuracy of these methods remain obscure. The addition of the entropic contribution can improve the results, especially for flexible proteins and ligands. Protein-Ligand Interaction Including the Desolvation Penalty. An advantage of CHM over implicit solvent models is in its ability to calculate the partial desolvation energy, which allows one to survey the competition between the protein-ligand interaction and the desolvation penalty of the protein. This may


Murata et al.

TABLE 2: The Calculated Binding Energies ∆Ebind between FKBP and Its Ligandsa ligand L1 L2 L3 L4 L5 L6

int j P,L E

int (B) jW (B) E P ,WL

int (P) j P,W E L

int (L) jW E P ,L

-91.7 (0.3) -76.7 (0.2) -62.3 (0.1) -55.1 (0.1) -42.2 (0.1) -48.9 (0.1)

-52.4 (0.5) -42.0 (0.4) -33.8 (0.3) -32.2 (0.3) -24.1 (0.2) -24.7 (0.2)

-68.9 (0.5) -58.8 (0.4) -46.7 (0.3) -43.6 (0.3) -27.4 (0.2) -34.9 (0.2)

-56.3 (0.2) -41.7 (0.2) -30.6 (0.2) -26.5 (0.2) -21.6 (0.1) -22.4 (0.1)

∆ELdef

∆Ebind

∆Hexptb

∆Gexptc

0.4

-18.4 (0.6) -18.2 (0.5) -17.0 (0.4) -17.0 (0.4) -16.8 (0.3) -16.2 (0.3)

-16.3

-13.2

-14.7

-12.8

N/A

-10.9

N/A

-9.5

N/A

-8.4

N/A

-7.8

0.1 1.8 0.2 0.4 0.2

a All energies are given in kcal/mol. Error estimates are in parentheses. ∆Ebind is divided into components according to eq 2. b The enthalpies are calculated at 298 K from a constant heat capacity change and a reference enthalpy at 286 K in ref 23. c Derived from experimentally measured inhibition constants in ref 16.

TABLE 3: The vdW and Electrostatic (ele) Energies of Interaction Energy Components (eq 4) between FKBP (P), (P) Its Ligands (L), and Water Clusters (WL )a ligand L1 L2 L3 L4 Figure 4. Correlation between the calculated binding energies and the experimentally determined binding free energies (r2 ) 0.93).

L5 L6

vdWb j P,L E

-52.4 (0.1) -49.3 (0.1) -41.7 (0.1) -39.2 (0.1) -30.0 (0.1) -31.6 (0.1)

ele b j P,L E

vdW (P)b j ele (P)b vdWc ele d eff e j P,W E EP,WL ∆EP,L ∆EP,L ∆EP,L L

-39.3 -12.6 (0.2) (0.1) -27.4 -12.9 (0.1) (0.1) -20.6 -11.3 (0.1) (0.1) -15.8 -9.9 (0.1) (0.1) -12.2 -8.7 (0.1) (0.1) -17.2 -8.9 (0.1) (0.1)

-56.4 (0.5) -45.9 (0.4) -35.4 (0.3) -33.7 (0.3) -18.7 (0.2) -26.0 (0.2)

-39.9 (0.1) -36.5 (0.1) -30.4 (0.1) -29.3 (0.1) -21.3 (0.1) -22.7 (0.1)

17.1 (0.5) 18.6 (0.4) 14.8 (0.3) 17.8 (0.3) 6.6 (0.2) 8.7 (0.2)

-22.8 (0.5) -17.9 (0.4) -15.6 (0.3) -11.5 (0.3) -14.7 (0.2) -14.0 (0.2)

a All energies are given in kcal/mol. Error estimates are in (P) vdW ele int j X,Y j X,Y j X,Y parentheses. b E +E )E (X ) P and Y ) L or WL ). See vdW (P) d int vdW vdW ele ele j X,Y j P,L j P,W j P,L Table 2 for E . c ∆EP,L ) E - E . ∆E ) E P,L L eff vdW ele ele (P) e j EP,WL . ∆EP,L ) ∆EP,L + ∆EP,L.

(P)

Figure 5. The representative structures of L3 in (a) the protein-ligand complex and (b) the isolated state.

prove to be a useful viewpoint for designing higher affinity j int ligands. Since E P,WL(P) in eq 2 is regarded as the partial solvation energy of the protein, the effective interaction between the eff is defined as protein and ligand in solution ∆EP,L eff int int j P,L j P,W ∆EP,L )E -E (P) L

(4)

int (B) jW (B) Also, the partial desolvation energy of ligands (E P ,WL int (L) j EWP ,L) can be easily calculated from the data in Table 2. Because its contribution to the total binding energy is much eff , we focus on the latter in the smaller compared to that of ∆EP,L following discussion. The vdW and electrostatic energy components of the plain int vdW ele j X,Y j X,Y j X,Y )E +E ), desolvation energies, and the interaction (E effective interaction energy are listed in Table 3. The magnitude vdW vdW (P) j P,W j P,L is three- to four-fold larger than E . Since the vdW of E L interaction energy becomes larger as the number of contact j vdW atoms between molecules increases, a larger E P,L is a reflection of a denser surface population of atoms in ligand L than that in

WL , as discussed previously. We also note that in TIP3P, hydrogen atoms do not contribute to the vdW interaction, which is determined only by heavy atoms. Moreover, water molecules have the heavy-atom fraction of 1/3, whereas for the protein, it is about 1/2, and for ligands, it is 1/2 or larger. All of this makes ligands have a larger van der Waals interaction with proteins compared to water clusters. This tendency can be expected for most of protein-ligand interactions in complexes. In many cases, the vdW interaction energies between the protein and ligand in protein-ligand complexes increase with the size of the ligand molecules. The vdW interaction energy vdW j P,L in the FKBP-ligand complexes follows this propensity, E L1 > L2 > L3 ≈ L4 > L5 ≈ L6. The vdW contribution to the jvdW partial solvation energy of the protein E P,WL(P) also increases with the size of the corresponding ligand except L1 and L2; the interaction is larger for the latter (L2 > L1), despite the number of water molecules in L1 being larger than that in L2. The experimental geometries of the complexes (Figure 6) shows that both L1 and L2 ligands stick out from the binding pocket, and L1 appears to be more solvent-exposed than L2. This might be the reason why the larger L1 brings a smaller desolvation penalty than L2. The result suggests that ligand location in the binding site of the protein is one of the important factors to determine the desolvation energy. ele ele (P) j P,L j P,W is one- to two-fold larger than E . The magnitude of E L Because the water clusters are more polar than the ligands, the ele for all of the ligands is effective electrostatic interaction ∆EP,L



Figure 6. Superimposition of L1 (thick red sticks) and L2 (thin blue sticks) crystal structures, as well as a part of the protein (FKBP) shown as green lines.

positive; the electrostatic contribution to the partial solvation ele (P) j P,W ) is larger than the electrostatic energy of the protein (E L ele j P,L ). We note that interaction energy between P and L (E electrostatic contribution to the effective interaction depends on the polar and nonpolar characteristics of the binding site of the protein, as well as that of ligand, and no general tendency is expected. In the case of FKBP complexes studied in this work, ele is positive the effective electrostatic interaction energy ∆EP,L (destabilization), and it is the vdW interaction that stabilizes the protein-ligand complexes in solution (Table 3). Also, the total binding energy ∆Ebind (Table 2) is largely determined by eff (Table 3). the total protein-ligand effective interaction ∆EP,L Residue-Ligand Interaction Including the Desolvation Penalty. Next, we discuss the contribution of each residue to the ligand binding including the solvation penalty, using the FKBP-L2 complex as an example. To examine the desolvation int (P) j P,W is decomposed into residue penalty of each residue, E L contributions

int j P,W E (P) ) L

nres

int ∑ Ej i,W i

(P) L

(5)

where nres is the number of residues in the protein. Similally, int j P,L is also decompsed the protein-ligand interaction energy E nres j int into residue contributions, ∑i Ei,L . Then, the effective eff is defined as residue-ligand interaction energy ∆Ei,L eff int int j i,L j i,W ∆Ei,L )E -E (P) L

(6)

Using these quantities, the competition between the protein-ligand int int (P) j i,L j i,W interaction E and the desolvation penalty -E can be L disccussed for each individual residue i. int eff j i,L and effective ∆Ei,L residue-ligand interaction The plain E energies are shown in Figure 7a and b, respectively, for all residues. It is observed that the relative importance of residues is dramatically changed by considering the desolvation penalty. Hydrogen Bonding Residues. Table 4 shows the vdW and electrostatic energy components for the protein-ligand (L2) binding, listed for the top 12 residues ranked by the stabilization energy of the plain interaction and also for Asp41 and Arg42, which are important because they form a salt bridge with the

Figure 7. (a) Plain and (b) effective interaction energies between each residue of FKBP and L2.

top ranked residues (Lys35 and Asp37). It is observed that the rankings are dramatically changed in the effective interaction energies due to the desolvation penalty, especially the electrostatic desolvation energy. There are four residues, Asp37, Glu54, Ile56, and Tyr82, forming hydrogen bonds with the ligand (Figure 8), and they are all high (within top five) in the ranking of the plain interactions. However, Asp37, Glu54, and Tyr82 fall to the 107, 102, and 11th places, respectively, in the ranking of the effective interaction energies. Among the hydrogen bonding residues, only nonpolar Ile56 keeps its position in both the plain and effective interactions. As for non-hydrogen bonding residues, the ranking of Tyr26 is raised from sixth in the plain interactions to first in the effective interactions. The ranking change of Arg42 is remarkable, from 106th to 2nd. The desolvation penalty of negatively charged residues, Asp37 and Glu54, exceeds the hydrogen bond energies, and their effective interaction energy becomes positive (destabilization). Asp37 forms a hydrogen bond with the ligand by the negatively charged side-chain carboxyl group. The hydrogen bond (-10.2 kcal/mol) is rather weak for the COO- · · · HO type since the group is also connected with the positively charged side chain of Arg42 by the salt bridge. Then, the large partial desolvation energy (17.7 kcal/mol, mostly electrostatic) of the COO- group exceeds the hydrogen bond energy, resulting in a destabilizing contribution to the protein-ligand binding. On the other hand, Glu54, which also has a negatively charged side chain, forms a hydrogen bond with the ligand by the mainchain carbonyl group, and the partial desolvation occurs around the group with the solvation of the side-chain carboxyl group remaining intact. Thus, both the hydrogen bond (-7.3 kcal/ mol) and the partial desolvation energies (8.6 kcal/mol, mainly electrostatic) are smaller than those of Asp37. However, the hydrogen bond stabilization is also beaten by the desolvation penalty, and small destabilization (1.3 kcal/mol) is left for the effective interaction. Tyr82 forms a (Ph)OH · · · O type hydrogen bond with the ligand, and the stabilization energy is the largest (-11.8 kcal/ mol) among the four hydrogen bonding residues. The partial desolvation energy is 10.5 kcal/mol (mostly electrostatic), and


Murata et al.

TABLE 4: VdW and Electrostatic Components (kcal/mol) of the Plain and Effective (eff) Residue-Ligand (L2) Interaction Energies (the latter corrected by the desolvation (desolv) energy) vdWb

electrostaticb

total interaction

residue

typea

plain

desolv

eff

plain

desolv

eff

plain (rank)c

eff (rank)d

Tyr82 Asp37 Ile56 Val55 Glu54 Tyr26 Trp59 Phe46 Phe99 Phe36 His87 Lys35 Asp41 Arg42 Lys35 + Asp41e Asp37 + Arg42e

OH · · · O COO- · · · H NH · · · O vdw- ∼ eleCO · · · H vdw- > elevdw- > elevdwvdwvdwvdweleele+ vdw- ∼ ele+ elevdw- ∼ ele-

-5.0 -2.0 -5.2 -4.5 -1.7 -5.2 -3.8 -4.8 -2.6 -2.1 -1.4 -0.1 0.0 -2.0 -0.1 -4.0

0.5 -0.4 1.7 1.4 0.1 1.3 1.2 1.6 0.8 0.6 0.4 0.0 0.0 0.8 0.0 0.4

-4.5 -2.4 -3.5 -3.1 -1.6 -3.9 -2.6 -3.2 -1.8 -1.5 -1.0 -0.1 0.0 -1.2 -0.1 -3.6

-6.8 -8.2 -3.2 -3.9 -5.6 -1.7 -1.5 0.4 -1.1 -0.8 -0.5 -1.3 0.3 3.7 -1.0 -4.5

10.0 18.1 3.1 3.8 8.5 0.6 2.3 0.6 0.5 0.8 5.2 -1.1 1.4 -6.3 0.3 11.8

3.2 9.9 -0.1 -0.1 2.9 -1.1 0.8 1.0 -0.6 0.0 4.7 -2.4 1.7 -2.6 -0.7 7.3

-11.8 (1) -10.2 (2) -8.4 (3) -8.4 (4) -7.3 (5) -6.9 (6) -5.3 (7) -4.4 (8) -3.7 (9) -2.9 (10) -1.9 (11) -1.4 (12) 0.3 (93) 1.7 (106) -1.1 -8.5

-1.3 (11) 7.5 (107) -3.6 (3) -3.2 (4) 1.3(102) -5.0 (1) -1.8 (8) -2.2 (7) -2.4 (6) -1.5 (10) 3.7 (106) -2.5 (5) 1.7 (102) -3.8 (2) -0.8 3.7

a Hydrogen bonds are shown as A · · · B; other interactions are classified based on the vdW (vdw) and electrostatic (ele) components of the eff (P) plain interaction energy. The superscript denotes stabilization (-) or destabilization (+). b The effective value ∆Ei,W is the sum of the plain L int (P) int j i,L j i,W interaction E (in complex P · L) and the residue i’s desolvation penalty, -E ; see eq 6. c Ranking by the plain interaction energies. L d Ranking by the effective interaction energies. e Salt bridges formed between charged residues.

Figure 8. The binding mode of FK506 (L2) at the active site. The dashed lines show hydrogen bonds.

the effective interaction is small and negative. Ile56 makes a NH · · · O type hydrogen bond between the main-chain NH and a ligand carbonyl oxygen. The partial desolvation energy (4.8 kcal/ mol) is small compared to the hydrogen bond energy (-8.4 kcal/ mol), resulting in the effective interaction energy of -3.6 kcal/ mol being the largest stabilization among the four hydrogen bonds. Thus, the contribution of the nonpolar residue (Ile56) to the ligand binding is the largest among the hydrogen bonding residues. In the seed/lead optimization, to enhance the binding affinity of a ligand, one effective strategies is to modify the ligand to form additional hydrogen bonds with the target protein. Actually, the binding affinity is often enhanced an order of magnitude per hydrogen bond.27-29 On the contrary, the desolvation penalty is considered to entirely cancel the stabilization brought by proteinligand hydrogen bonds30,31 in some cases such as the binding of glucose to glycogen phosphorylase.32 Analyses of the desolvation penalty based on CHM will be helpful to understand the real role of hydrogen bonds in the protein-ligand binding. Non-hydrogen Bonding Residues. Tyr26 holds the top position in the effective interaction energy ranking. The side-chain OH group of the residue forms a hydrogen bond with the COOgroup of Asp37, and unlike Tyr82, this residue does not form a

Figure 9. (a) Plain and (b) effective interaction energies between each residue of FKBP and L1.

hydrogen bond with the ligand. The interaction energy is mainly brought by the vdW interaction between the aromatic ring and the ligand (-5.2 out of -6.9 kcal/mol in the plain interaction energy). The electrostatic desolvation energy is much smaller than that of Tyr82 (0.6 versus 10.0 kcal/mol) because the COO- group of Asp37 is near the OH group and only the limited part of the OH group is desolvated. The moderate desolvation penalty relative to the plain interaction energy results in the top ranking in the effective residue interaction energies. Another polar residue, Trp59, interacts with the ligand by its aromatic ring, and the main components of the plain interaction and desolvation energies are the vdW interactions. Neither is large, and the ranking in the effective interaction energies is not changed very much. The positively charged residues, Arg42, Lys35, and His87, are not directly in contact with the ligand. The plain interaction energies of Arg42 and Lys35 are not large, and their rankings are 106- and 12th, respectively. It should be noted that the positive plain interaction energy of Arg42 (+1.7 kcal/mol) is changed into negative (-3.8 kcal/mol) in the effective interaction energy. This



TABLE 5: VdW and Electrostatic Components (kcal/mol) of the Plain and Effective (eff) Residue-Ligand (L1) Interaction Energies (the latter corrected by the desolvation (desolv) energy) vdWb

electrostaticb

total interaction

residue

typea

plain

desolv

eff

plain

desolv

eff

plain (rank)c

eff (rank)d

Glu54 Tyr82 Val55 Asp37 Ile56 Gln53 Trp59 Tyr26 Lys52 Phe46 Phe36 Phe99 Arg42 Asp37 + Arg42e Glu54 + Lys52e

CO · · · H OH · · · O vdw- > eleCOO- · · · H NH · · · O CO · · · H vdw- > elevdw- > eleelevdwvdwvdwele+ ele- > vdwele- > vdw-

-6.3 -3.6 -6.6 -1.1 -5.1 -1.1 -4.0 -3.7 -0.5 -3.6 -2.5 -2.3 -0.5 -1.6 -6.8

1.0 -0.1 2.1 -0.4 1.5 -0.1 1.2 1.1 0.1 1.1 0.8 0.6 0.2 -0.2 1.1

-5.3 -3.7 -4.5 -1.5 -3.6 -1.2 -2.8 -2.6 -0.4 -2.5 -1.7 -1.7 -0.3 -1.8 -5.7

-6.3 -7.9 -3.8 -8.4 -3.6 -6.8 -1.5 -1.3 -3.7 -0.2 -1.0 -1.0 1.4 -7.0 -10.0

16.5 9.8 5.6 17.4 3.7 5.1 2.5 0.0 -2.0 0.6 0.8 0.6 -6.8 10.6 14.5

10.2 1.9 1.8 9.0 0.1 -1.7 1.0 -1.3 -5.7 0.4 -0.2 -0.4 -5.4 3.6 4.5

-12.6 (1) -11.5 (2) -10.4 (3) -9.5 (4) -8.7 (5) -7.9 (6) -5.5 (7) -5.0 (8) -4.2 (9) -3.8 (10) -3.6 (11) -3.3 (12) 0.9 (98) -8.6 -8.5

4.9 (105) -1.8(14) -2.7 (7) 7.5 (106) -3.5 (4) -2.9 (5) -1.8 (16) -3.9 (3) -6.1 (1) -2.1 (11) -1.9 (13) -2.1 (12) -5.7 (2) 1.8 -1.2

a Hydrogen bonds are shown as A · · · B, and other interactions are classified based on the vdW (vdw) and electrostatic (ele) components of eff (P) plain interaction energy. The superscript denotes stabilization (-) or destabilization (+). b The effective value ∆Ei,W is the sum of the plain L int (P) c int j i,L j i,W interaction E (in complex P · L) and the residue i’s desolvation penalty, -E ; see eq 6. Ranking by the plain interaction energies. L d Ranking by the effective interaction energies. e Salt bridges formed between charged residues.

remarkable change is accomplished by the negative desolvation energy of the residue (-5.5 kcal/mol). Arg42 is salt-bridged to Asp37. The desolvation penalty of the partner Asp37 is very large (+17.7 kcal/mol), suggesting that solvent water molecules around the salt bridge are arranged to stabilize Asp37 at the cost of destabilizing Arg42. This may commonly take place in the solvation of salt bridges. The desolvation energy of Lys35, which is saltbridged to Asp41, is also negative (-1.0 kcal/mol). Due to the negative desolvation energies, both Arg42 and Lys35 have raised rankings in the effective interaction energies. His87 is a standalone positive residue interacting weakly with the ligand. The residue is located at the lid of the binding pocket, and one side of the aromatic ring faces the ligand while the other side is exposed to solvent. Then, the desolvation occurs at a broad area, resulting in a rather large desolvation energy (+5.6 kcal/mol, mostly electrostatic). The nonpolar residues Phe36, Phe46, and Phe99 bind to ligands mainly by the vdW interaction. Their desolvation energies follow the propensity described previously. Val55 forms a CH · · · O type bond between the main-chain CRH and O of the lactone carbonyl group of the ligand. Both the vdW and electrostatic components contribute to the plain interaction and the desolvation energies. Thus, the contribution of the Val55 to the ligand binding is the second largest among the nonpolar residues next to the Ile56. We note that the plain interaction energies of nonpolar residues as well as those of polar and hydrogen bonding residues are consistent with the results of quantum chemical calculations, which include the electron correlation effects.33 Salt-Bridged Residues. It is useful to analyze the desolvation of salt-bridged charged residues. As discussed above regarding saltbridged Asp37 and Arg42, the solvation of a cation is destabilized due to the strong solvation of the partner anion. The solvation of the salt-bridged ion pairs is deeply connected to each other. The plain interaction of the two ions with a ligand is also coupled, so that the stabilizing electrostatic interaction of one ion is often accompanied by a destabilizing interaction of the other ion because of their opposite charges. The plain interaction energies and the desolvation energies of the individual charged residues in salt bridges are very large, but the sum of them is moderate. We propose to consider the plain and effective interactions of salt-bridged residues together. From this point of view, the effective interaction energies of the salt-bridged residues become insignifi-

TABLE 6: The Calculated Binding Energies ∆Ebind between FKBP and L2 for TIP3P and TIP4P Water Modelsa water model TIP3P TIP4P

int j P,L E

int (B) jW (B) E P ,WL

-76.7 (0.2) -77.2 (0.3)

-42.0 (0.4) -44.0 (0.4)

int (P) j P,W E L

int (L) jW E ∆ELdef P ,L

-58.8 -41.7 (0.4) (0.2) -59.6 -43.0 (0.4) (0.2)

0.1 0.1

∆Ebind -18.2 (0.5) -18.5 (0.5)

a All energies are given in kcal/mol. Error estimates are in parentheses. ∆Ebind is divided into components according to eq 2.

cant, as shown in Table 4; the sum of the effective interaction energies of salt-bridged Asp37 and Arg42 becomes very small (-0.8 kcal/mol), and that of another salt-bridged Lys35 and Asp41 pair becomes positive (+3.7 kcal/mol). Thus, as far as salt-bridged chargedresiduesareconcerned,theircontributiontotheprotein-ligand binding becomes less important, despite the fact that the apparent plain interaction energies are large. The Other Ligands. Although almost the same trends are obtained for the other ligands (L1, L3-L6), we show another example using the FKBP-L1 complex for comparison and also to demonstrate that the trends that we discussed are of a general eff j int nature. The plain E i,L and effective ∆Ei,L residue-ligand interaction energies are shown in Figure 9a and b, respectively, for all residues. It is also observed that the relative importance of residues is dramatically changed by considering the desolvation penalty. Table 5 shows the vdW and electrostatic energy components for the protein-ligand (L1) binding, listed for the top 12 residues ranked by the stabilization energy of the plain interaction and also for Arg42, which is important because this residue forms a salt bridge with the top ranked residue (Asp37). The top-ranking residues listed in Tables 4 and 5 are almost the same (10 out of 12). The new ranked residues in Table 5 are Gln53 and Lys52. Gln53 forms a hydrogen bond with the ligand by the main-chain carbonyl group. The partial desolvation energy (5.0 kcal/mol) cancels the hydrogen bond energy (-7.9 kcal/mol), resulting in the effective interaction energy of -2.9 kcal/mol. Thus, the contribution of the polar residue (Gln53) to the ligand binding is the second largest among the hydrogen bonding residues (the largest residue is Ile56). The plain interaction energy of Lys52 (-4.2 kcal/mol) is smaller than the effective interaction energy (-6.1 kcal/mol). However, Lys52 is salt-bridged to Glu54, and the desolvation penalty of the


Murata et al.

TABLE 7: The vdW and Electrostatic (ele) Energies of Interaction Energy Components (eq 4) between FKBP (P), L2 (L), and (P) Water Clusters (WL ) for TIP3P and TIP4P Water Modelsa water model TIP3P TIP4P

c

vdWb j P,L E

ele b j P,L E

vdW (P)b j P,W E L

ele (P)b j P,W E L

vdWc ∆EP,L

ele d ∆EP,L

eff e ∆EP,L

-49.3 (0.1) -49.6 (0.2)

-27.4 (0.1) -27.7 (0.2)

-12.9 (0.1) -11.9 (0.1)

-45.9 (0.4) -47.7 (0.4)

-36.5 (0.1) -37.6 (0.1)

18.6 (0.4) 20.0 (0.5)

-17.9 (0.4) -17.6 (0.5)

(P) a vdW ele int int j X,Y j X,Y j X,Y j X,Y All energies are given in kcal/mol. Error estimates are in parentheses. b E +E )E (X ) P and Y ) L or WL ). See Table 2 for E . vdW (P) d ele (P) e vdW vdW ele ele eff vdW ele j P,L j P,W j P,L j P,W ∆EP,L )E -E . ∆E ) E E . ∆E ) ∆E + ∆E . P,L P,L P,L P,L L L

TABLE 8: The vdW and Electrostatic Components (kcal/mol) of the Plain and Effective (eff) Residue-Ligand (L2) Interaction Energies (the latter corrected by the desolvation (desolv) energy) for the TIP4P Water Model vdWb

electrostaticb

total interaction

residue

typea

plain

desolv

eff

plain

desolv

eff

plain (rank)c

eff (rank)d

Tyr82 Glu54 Val55 Asp37 Ile56 Tyr26 Trp59 Phe46 Phe99 Phe36 Glu107 Lys35 Asp41 Arg42 Lys35 +Asp41e Asp37 +Arg42e

OH · · · O CO · · · H vdw- ∼ ele+ COO- · · · H NH · · · O vdw- > elevdw- > elevdwvdw- > elevdw- > eleeleeleele+ vdw- ∼ ele+ eleele- > vdw-

-4.8 -1.8 -4.7 -1.8 -5.7 -5.4 -3.8 -4.8 -2.9 -1.9 0.0 -0.1 0.0 -0.9 -0.1 -2.7

0.6 0.1 1.4 -1.1 2.0 1.2 1.0 1.7 0.6 0.6 0.0 0.0 0.0 0.4 0.0 -0.7

-4.2 -1.7 -3.3 -2.9 -3.7 -4.2 -2.8 -3.1 -2.3 -1.3 0.0 -0.1 0.0 -0.5 -0.1 -3.4

-5.7 -7.6 -4.2 -6.8 -2.8 -1.2 -1.4 0.4 -1.2 -0.9 -1.7 -1.6 0.5 1.8 -1.1 -5.0

8.1 9.2 4.5 17.6 2.7 0.3 1.8 0.5 0.8 0.9 1.9 0.3 0.7 -5.3 1.0 12.3

2.4 1.6 0.3 10.8 -0.1 -0.9 0.4 0.9 -0.4 0.0 0.2 -1.3 1.2 -3.5 -0.1 7.3

-10.6 (1) -9.5 (2) -8.9 (3) -8.7 (4) -8.5 (5) -6.6 (6) -5.2 (7) -4.4 (8) -4.1 (9) -2.8 (10) -1.7 (11) -1.7 (12) 0.5 (95) 0.9 (101) -1.1 -7.7

-1.8 (8) -0.1 (27) -3.0 (4) 7.9 (107) -3.8 (3) -5.1 (1) -2.4 (6) -2.2 (7) -2.7 (5) -1.3 (10) 0.2 (96) -1.3 (9) 1.2 (102) -4.0 (2) -0.2 3.9

a Hydrogen bonds are shown as A · · · B, and other interactions are classified based on the vdW (vdw) and electrostatic (ele) components of eff (P) is the sum of the plain plain interaction energy. The superscript denotes stabilization (-) or destabilization (+). b The effective value ∆Ei,W L int (P) int j i,L j i,W interaction E (in complex P · L) and the residue i’s desolvation penalty, -E ; see eq 6. c Ranking by the plain interaction energies. L d Ranking by the effective interaction energies. e Salt bridges formed between charged residues.

partner Glu54 is very large (+17.5 kcal/mol). This suggests that solvent water molecules around the salt bridge are arranged to stabilize Glu54 at the cost of destabilizing Lys52, which is also confirmed for other salt-bridged residues in the FKBP-L2 complex, as discussed above. Comparison of Two Water Models. In the present study, we used the TIP3P18 water model for explicit water. However, could we obtain the same results for other water models? Water models can be approximately classified into three groups (three-, four-, and five-site models). While TIP3P and SPC34 are classified as three-site models, water models which have the extra midpoint site such as TIP4P18 and COS/G335 are classified as four-site models. TIP5P18 and ST236 exemplify five-site models which have the two lone pair sites. In order to examine the effect of other water models, we performed MD simulations of 1 × 106 time steps (or 1 ns total) for the FKBP-L2 complex using the TIP4P water model under the same conditions as those for the TIP3P water model. In Tables 6 and 7, the calculated binding energy and each component for TIP4P are compared with those for TIP3P. Although the interaction energies for TIP4P are slightly higher than those for TIP3P, the results are in good agreement with those for TIP3P. The reason for these slightly higher energies for TIP4P might be that the negative charge is removed from the O atom and is placed at an extra point on the bisector of the HOH angle. Overall, the differences are within ∼1 kcal/mol, which is roughly the order of magnitude of the statistical errors of calculations except for int (B) int (B) jW jW (B). The difference of E (B) (2.0 kcal/mol) is slightly E P ,WL P ,WL higher than the other differences. However, this difference is j int j int canceled by the sum of E P,WL(P) (0.8 kcal/mol) and EWP(L),L(1.3 kcal/ mol), and thus, ∆Ebind for TIP4P is in excellent agreement with that for TIP3P.

Table 8 shows the vdW and electrostatic energy components for the protein-ligand (L2) binding for TIP4P. The top 12 residues ranked by the plain interaction and also that for Arg42, which forms a salt bridge with the top ranked residue (Asp37), are listed. The top-ranking residues in Tables 4 and 8 are almost the same (11 out of 12), and the top eight residues ranked by the effective interaction energies are exactly the same. This remarkable change of the ranking in the effective interaction energies for TIP4P is also accomplished by the desolvation penalty. We conclude that the desolvation penalty is the main reason for stabilization of the protein-ligand complex in TIP4P as well as TIP3P. Conclusions In this work, we have proposed a new model for calculating the binding energy between the protein and ligand in solution, which is referred to as the cluster hydration model. We applied this model to six FKBP complexes and obtained a good linear correlation (r2 ) 0.93) between the calculated binding energies and experimental binding affinities. This model can also be used to study other processes, such as protein-DNA binding, or chemical reactions in solution. One of the most useful features of CHM is the ability to estimate the desolvation penalty of each amino acid residue. The analysis of the partial desolvation energies for all FKBP complexes revealed that nonpolar and polar residues in the binding site dominate in stabilizing the FKBP complexes, while charged residues do not contribute to the stabilization. We also found that some hydrogen bonds between the protein and the ligand do not contribute to stabilize the complex due to large desolvation penalties.

Binding Energy Calculations of Protein-Ligand Complexes Our results suggest that the desolvation penalty should be taken into account in designing ligands with a higher affinity, and CHM can be used for this purpose in the seed/lead optimization process. Acknowledgment. We acknowledge the Research Center for Computational Science, Okazaki Research Facilities, National Institutes of Natural Sciences, and the Academic Center for Computing and Media Studies, Kyoto University for computer resources. This work was supported by Astellas Pharma Inc. (formerly Fujisawa Pharmaceutical Co.) and a Grant-in-Aid for Young Scientists (19790032), Ministry of Education, Culture, Sports, Science and Technology (Japan). References and Notes (1) Beveridge, D. L.; Dicapua, F. M. Free-energy via molecular simulation: Applications to chemical and biomolecular systems. Annu. ReV. Biophys. Bio. 1989, 18, 431–492. (2) Straatsma, T. P.; McCammon, J. A. Computational alchemy. Annu. ReV. Phys. Chem. 1992, 43, 407–435. (3) Tomioka, N.; Itai, A.; Iitaka, Y. J. A method for fast energy estimation and visualization of protein-ligand interaction. J. Comput.-Aided Mol. Des. 1987, 1, 197–210. (4) Bo¨hm, H. J. The development of a simple empirical scoring function to estimate the binding constant for a protein-ligand complex of known three-dimensional structure. J. Comput.-Aided Mol. Des. 1994, 8, 243–256. (5) Wallqvist, A.; Jernigan, R. L.; Covell, D. G. A preference-based free-energy parameterization of enzyme-inhibitor binding: Applications to HIV-1-protease inhibitor design. Protein Sci. 1995, 4, 1881–1903. (6) Verkhivker, G.; Appelt, K.; Freer, S. T.; Villafranca, J. E. Empirical free energy calculations of ligand-protein crystallographic complexes. I. Knowledge-based ligand-protein interaction potentials applied to the prediction of human immunodeficiency virus 1 protease binding affinity. Protein Eng. 1995, 8, 677–691. (7) Åqvist, J.; Medina, C.; Samuelsson, J.-E. A New Method for Predicting Binding Affinity in Computer-Aided Drug Design. Protein Eng. 1994, 7, 385–391. (8) Hanson, T.; Åqvist, J. Estimation of Binding Free Energies for HIV Proteinase Inhibitors by Molecular Dynamics Simulations. Protein Eng. 1995, 8, 1137–1144. (9) Åqvist, J.; Hanson, T. On the Validity of Electrostatic Linear Response in Polar Solvents. J. Phys. Chem. 1996, 100, 9512–9521. (10) Ersmark, K.; Nervall, M.; Hamelink, E.; Janka, L. K.; Clemente, J. C.; Dunn, B. M.; Blackman, M. J.; Samuelsson, B.; Åqvist, J.; Hallberg, A. Synthesis of Malarial Plasmepsin Inhibitors and Prediction of Binding Modes by Molecular Dynamics Simulations. J. Med. Chem. 2005, 48, 6090– 6106. (11) Bortolato, A.; Moro, S. In silico binding free energy predictability by using the linear interaction energy (LIE) Method: bromobenzimidazole CK2 Inhibitors as a case study. J. Chem. Inf. Model 2007, 47, 572–582. (12) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; Cheatham, T. E., III. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889–897. (13) Van Duyne, G. D.; Standaert, R. F.; Schreiber, S. L.; Clardy, J. Atomic structure of the rapamycin human immunophilin FKBP12 complex. J. Am. Chem. Soc. 1991, 113, 7374–7375. (14) Van Duyne, G. D.; Standaert, R. F.; Karplus, P. A.; Schreiber, S. L.; Clardy, J. Atomic structure of FKBP-FK506, an immunophilin-immunosuppressant complex. Science 1991, 252, 839. (15) Holt, D. A.; Luengo, J. I.; Yamashita, D. S.; Oh, H. J.; Konialian, A. L.; Yen, H. K.; Rozamus, L. W.; Brandt, M.; Bossard, M. J.; Levy, M. A.; Eggleston, D. S.; Liang, J.; Schultz, L. W.; Stout, T. J.; Clardy, J. Design, synthesis, and kinetic evaluation of high-affinity FKBP ligands and the X-ray crystal structures of their complexes with FKBP12. J. Am. Chem. Soc. 1993, 115, 9925–9938. (16) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.;

J. Phys. Chem. B, Vol. 113, No. 3, 2009 817 Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; Gaussian, Inc.: Wallingford, CT, 2004. (17) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Wang, B.; Pearlman, D. A.; Crowley, M.; Brozell, S.; Tsui, V.; Gohlke, H.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Schafmeister, C.; Caldwell, J. W.; Ross, W. D.; Kollman, P. A. AMBER 8; University of California: San Francisco, 2004. (18) Jorgensen, W. L.; Chandrasekhar, J.; Madurs, J.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. (19) Wang, J.; Cieplak, P.; Kollman, P. A. How does well a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules. J. Comput. Chem. 2000, 21, 1049–1074. (20) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174. (21) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684–3690. (22) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian eqations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327–341. (23) Connelly, P. R.; Thomson, J. A. Heat capacity changes and hydrophobic interactions in the binding of FK506 and rapamycin to the FK506 binding protein. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 4781–4785. (24) Swanson, J. M.; Henchman, R. H.; McCammon, J. A. Revisiting free energy calculations: a theoretical connection to MM/PB/SA and direct calculation of the association free energy. Biophys. J. 2004, 86, 67–74. (25) Xu, Y.; Wang, R. A computational analysis of the binding affinities of FKBP12 inhibitors using the MM-PB/SA method. Proteins: Struct., Funct., Bioinf. 2006, 64, 1058–1068. (26) Lee, M. S.; Olson, M. A. Calculation of absolute protein-ligand binding affinity using path and endpoint approaches. Biophys. J. 2006, 90, 864–877. (27) Phillips, G.; Guilford, W. J.; Buckman, B. O.; Davey, D. D.; Eagen, K. A.; Koovakkat, S.; Liang, A.; McCarrick, M.; Mohan, R.; Ng, H. P.; Pinkerton, M.; Subramanyam, B.; Ho, E.; Trinh, L.; Whitlow, M.; Wu, S.; Xu, W.; Morrissey, M. M. Design, synthesis, and activity of a novel series of factor Xa inhibitors: Optimization of arylamidine groups. J. Med. Chem. 2002, 45, 2484–2493. (28) Linus, S. L.; Sookhee, H.; Richard, G. B.; Nancy, N. T.; Laurie, A. C.; George, A. D.; Tung, M. F.; hun-Pyn, S.; Jing, C. X.; Mark, T. G.; William, K. H. Conformational analysis and receptor docking of N-[(1S,2S)3-(4-Chlorophenyl)-2-(3-cyanophenyl)-1-methylpropyl]-2-methyl-2-{[5-(trifluoromethyl)pyridin-2-yl]oxy}propanamide (Taranabant, MK-0364), a novel, acyclic cannabinoid-1 receptor inverse agonist. J. Med. Chem. 2008, 51, 2108–2114. (29) Michael, R. D.; Willy, B. T.; Akhilesh, B. B.; Scott Perrin Jr.; Jyotsna, T.; Alison, R.; Challa, V. K. Contributions of Hydroxyethyl Groups to the DNA Binding Affinities of Anthracene Probes. J. Phys. Chem. B 2006, 110, 20693–20701. (30) Ajay, A.; Murcko, M. A. Computational Methods to Predict Binding Free Energy in Ligand-Receptor Complexes. J. Med. Chem. 1995, 38, 2484– 2493. (31) Friesner, R. A.; Murphy, R. B.; Repasky, M. P.; Frye, L. L.; Greenwood, J. R.; Halgren, T. A.; Sanschagrin, P. C.; Mainz, D. T. Extra precision glide: Docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J. Med. Chem. 2006, 49, 6177–96. (32) Watson, K. A.; Mitchell, E. P.; Johnson, L. N.; Son, J. C.; Bichard, C. J. F.; Orchard, M. G.; Fleet, G. W. J.; Oikonomakos, N. G.; Leonidas, D. D.; Kontou, M.; Papageorgiouil, A. Design of inhibitors of glycogen phosphorylase: A study of R- and β-C-glucosides and 1-thio-β-D-glucose compounds. Biochemistry 1994, 33, 5745–5758. (33) Nakanishi, I.; Fedorov, D. G.; Kitaura, K. Molecular recognition mechanism of FK506 binding protein: an all-electron fragment molecular orbital study. Proteins: Struct., Funct., Bioinf. 2007, 68, 145–158. (34) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; p 331. (35) Yu, H.; van Gunsteren, W. F. Charge-on-spring polarizable water models revisited: from water clusters to liquid water to ice. J. Chem. Phys. 2004, 121, 9549–9564. (36) Stllinger, H. F.; Rahman, A. Improved simulation of liquid water by molecular dynamics. J. Chem. Phys. 1974, 60, 1545–1557.

JP805007F

Cluster Hydration Model for Binding Energy Calculations of Protein

Recommend Documents