J. Phys. Chem. B 2008, 112, 15355–15360
15355
Computation of the Contribution from the Cavity Effect to Protein-Ligand Binding Free Energy F. V. Grigoriev,* S. N. Gabin, A. N. Romanov, and V. B. Sulimov Research Computing Center, M.V. LomonosoV Moscow State UniVersity, Russia, 119992 Moscow, VorobjoVi Gory, MGU NIVC ReceiVed: May 10, 2008; ReVised Manuscript ReceiVed: September 22, 2008
We present results of the investigation of the cavity creation/annihilation effect in view of formation of the protein-ligand (PL) complexes. The protein and ligand were considered as rigid structures. The change of the cavity creation/annihilation free energy ∆Gcav was calculated for three PL complexes using the thermodynamic integration procedure with the original algorithm for growing the interaction potential between the cavity and the water molecules. The thermodynamic cycle consists of two stages, annihilation of the cavity of the ligand for the unbound state and its creation at the active site of the protein (bound state). It was revealed that for all complexes under investigation, the values of ∆Gcav are negative and favorable for binding. The main contribution to ∆Gcav appears due to the annihilation of the cavity of the ligand. All computations were made using the parallel version of CAVE code, elaborated in our preceding work. Introduction The prediction of ligand affinity to target a protein is an important task of drug design. A lot of approaches were developed, including empirical scoring functions,1-3 continuum solvent methods,4-7 quantum mechanical approaches,6 molecular dynamic simulations,8-11 and special methods for the effective investigation of the configurational space of the protein-ligand (PL) complexes.12 The effect of partial desolvation of both the ligand and protein as a result of formation of the PL complex plays an important role in all methodologies predicting the PL binding free energy. The essential part of the desolvation effect appears due to the cavity formation/annihilation, where the “cavity” is defined as a void of excluded volume in water prepared in order to accommodate the solute particle. The free energy Gcav of this effect is defined as a change of the system free energy resulting from the cavity creation. Because the cavity creation is the energetically unfavorable process, the value of Gcav is positive. The cavity creation/annihilation contribution to the free energy of the PL complex formation can achieve 15 kcal/mol.4,10 The quantitative theory of the cavity creation/annihilation effect is an active area of the recent research.13-26 One of the approaches for the computation of the corresponding free-energy change is based on Monte Carlo (MC) simulations providing the mean force potential, which represents the interaction of the dissolved molecule with the ensemble of water molecules.16,18,25,26 In our preceding work,25 the original two-step algorithm for growing the interaction potential between the hard cavity wall and the water molecules was developed. On the basis of this algorithm, the values of Gcav were calculated for a large set of organic molecules employing the thermodynamic integration procedure with MC sampling of the microstates representing the ensemble of water molecules. In the present work, we used this algorithm to calculate the change of the cavity creation/ annihilation free energy ∆Gcav arising due to the formation of PL complexes. * To whom correspondence should be addressed. E-mail: fedor.grigoriev@ gmail.com. Phone and Fax: (095)9393253.
The thermodynamic cycle determining ∆Gcav was organized as follows. (i) Annihilation of the cavity corresponding to a ligand. This is an energetically favorable process; its contribution to ∆Gcav is negative and equal to -Gcav(L), where Gcav(L) means the free energy of cavity creation. (ii) Creation of the new cavity corresponding to the ligand at the active site of the protein. It could be expected that this process is energetically unfavorable because it includes removing of water molecules from the active site of the protein. In principle, such a contribution can be estimated as Gcav(PL) Gcav(P), (where Gcav(PL) and Gcav(P) denote free energies of cavity formation for the PL complex and protein, respectively) by applying the approaches developed for calculation of the cavity creation/annihilation energy for macromolecules.10,20 However, in the present work, we calculated Gcav(PL) - Gcav(P) directly from the MC simulation where the initial and final states of the system represented the dissolved protein and the PL complex respectively. The advantage of this method is that instead of treating Gcav(PL) and Gcav(P) separately, we obtain the difference Gcav(PL) - Gcav(P) directly. It is worth reminding that an accurate calculation of Gcav(PL) and Gcav(P) is a complicated task owing to large sizes inherent to both particles. We applied this two-stage scheme to three PL complexes. Their initial structures were taken from the Protein Data Bank27 and specially adapted for modeling as described below. The MC simulations were carried out using the CAVE (CAVity free Energy) program specially developed for this goal.28 Methods Section The quantity∆Gcav is calculated as a difference of the corresponding values for PL and the ligand
∆Gcav ) ∆Gcav(PL) - Gcav(L)
(1)
Here, designations L and PL correspond to the ligand and the PL complex, respectively. The corresponding thermodynamic cycle is schematically shown in Figure 1.
10.1021/jp8041439 CCC: $40.75 2008 American Chemical Society Published on Web 11/08/2008
15356 J. Phys. Chem. B, Vol. 112, No. 48, 2008
Grigoriev et al.
Figure 1. The thermodynamic cycle for calculation of the ∆Gcav.
Every term in eq 1 is calculated using the thermodynamic integration (TI) procedure for a computation of the free-energy change ∆G resulting from a quasi-equilibrium creation of the external potential (Ucav_L) according to29
∆G(L) )
∫01 〈∂Ucav_L(λ)/∂λ〉λdλ
(2)
where
〈∂Ucav_L(λ)/∂λ〉λ )
∫Γ (∂Ucav_L(λ)/∂λ)exp(-β(Ucav_L(λ) + Uww))dΓ (3) ∫Γ exp(-β(Ucav_L(λ) + Uww))dΓ
and λ is a dimensionless parameter which governs the evolution of Ucav_L; Uww is the potential energy of the intermolecular interaction of water molecules, and β ) 1/(kBT). The integration is performed over the total configurational space Γ. The cavity potential Ucav_L is organized as follows. For the individual sphere corresponding to every atom, the cavity potential Ucav is defined as25
( )
Ucav(R) ) U01(Rcav - R)θ θ(x) ) 1
if 0 e x e 1
( )
R R + U02θ Rcav Rcav θ(x) ) 0 if x g 1
where Rcav is the cavity radius, and R represents the distance between the cavity center and the oxygen atom of a water molecule. Therefore, for the ligand composed of Nat atoms interacting with N water molecules, the Ucav_L is defined as25 N
Ucav_L )
Nat
∑∑ i)1 j)1
(
( )
U01(Rjcav - Rij)θ
( ))
Rij Rij + U02θ Rjcav Rjcav
(4) Here, Rij is the distance between the oxygen atom of the ith water molecule and the jth atom of the solute, and Rjcav denotes the radius of the excluded volume sphere corresponding to atom j. It is defined as Rjcav ) Rsol + Rj, where Rsol is the effective radius of the water molecule and Rj is the radius of the atomic sphere of the jth atom. The radii of atomic spheres were calculated as recommended in ref 30 when applied to electrostatic free-energy computations in terms of the generalized Born procedure31 based on force field MMFF94.32 The potential of the interaction of protein with water was organized similar eq 4.
The difference between evaluations of ∆Gcav(L) and ∆Gcav(PL) is that for the ∆Gcav(L), the growing procedure of the potential in eq 4 is performed in bulk water, whereas for ∆Gcav(PL), this procedure is performed at the active site of the protein. The averages in eq 3 were found via MC modeling using the CAVE (CAVity free Energy) program. The NPT ensemble with periodic boundary conditions at temperature T ) 25 °C and pressure P ) 1 atm was applied. Sampling microstates of the ensemble was performed according to the Metropolis scheme.33-35 Rigid four-site nonpolarizable model TIP4P36 was used for water molecules. To take into account the statistical error of the MC simulation, the averaging of the ∆Gcav over the set of the eight independent runs was used for computations of both ∆Gcav(PL) and Gcav(L). These runs differed from each other by random number arrays and by initial microstates of the ensembles. The total number of the sampled microstates during the MC simulation was close to 3 × 108. Empirical expressions for Gcav were also considered (see the Discussion). In order to test their validity, it was necessary to calculate the values of the solvent-accessible surface area (SASA) and volume V by this boundary. The SASA was calculated as follows. In agreement with the definition of the cavity potential (eq 4), both the ligand and protein were described as a set of spheres, each of them being characterized by the radius Rjcav (for the jth atom). Every spherical surface was partitioned using the spherical coordinates into elements (tesserae) along the parallels and meridians.37 The total SASA value was calculated as a sum of the tesserae areas, excluding hidden ones appearing inside of spheres, corresponding to any atom of the ligand or protein. Volume V surrounded by the SASA was calculated as V ) b. (1/3)∫S bdS r b following the Stokes theorem, ∫V diVa bdV ) ∫S b adS All computations were made by means of facilities of the MSU Computational Center38 and Joint SuperComputer Center (JSCC) of the Russian Academy of Sciences (RAS).39 For more details on the TI procedure and MC simulation, see ref 25. Preparation of Structures for Modeling. The initial structures of the PL complexes were borrowed from the Protein Data Bank (PDB).27 Hydrogen atoms were added using the Reduce program,40 which optimizes the positions of OH, SH, NH3+ fragments, hydrogen atoms of methyl groups, Asn, and Gln. Further, we applied a truncation procedure to the protein in order to extract only its active site. Truncated structures containing no less than ∼500 atoms (including hydrogens) were prepared using the program elaborated in our group; it keeps all protein residues that have at least one heavy atom at a distance no more than a cutoff radius of Rcut from any heavy atom of the ligand. The value Rcut was 8 Å, being chosen in order to allow for all
Cavity Effect to Protein-Ligand Binding
J. Phys. Chem. B, Vol. 112, No. 48, 2008 15357 of both the forward and the reverse paths were carried out for every PL complex and ligand. For ligands, we computed also the quantity GcavV(L) defined in eq 5 in order to verify our approximate scheme designed for treating cavity creation free energies of small- and mediumsized molecules using their molecular volumes V25
GVcav ) ξV + ηχns
(5)
where ξ ) 0.0361 (kcal/mol)Å-3, η ) 22.3 kcal/mol. The nonspherical index χns was introduced in ref 25 as
Figure 2. Structural formulas of the ligand from PL complexes under investigation.
protein fragments forming close contacts with a ligand to be kept as specific hydrophobic pockets. The dependence of the ∆Gcav on the Rcut was investigated (see the Discussion). All coordinates of atoms belonging to the protein and ligand were kept fixed under the simulation. The initial states of the water molecule ensemble with the cavity corresponding to a truncated protein were prepared as follows. The equilibrated at 300 K and 1 atm water box was taken, and water molecules inside of the cavity were removed. This structure was equilibrated until the temperature, average total potential energy of the ensemble, and the volume converged. Then, insertion of the cavity potential corresponding to a ligand at the active site of the protein was started in order to evaluate ∆Gcav(PL). The equilibrated water box was chosen as the initial state for a calculation of Gcav(L). Results Three PL complexes were investigated, benzamidine-trypsin (3ptb), neuraminidase-oseltamivir (2ht7), and streptavidin-biotin (1swt). The ligands are shown in Figure 2. There are several reasons for choosing these complexes. First, their protein ingredients are relatively small and suitable for modeling. Second, the active sites of those proteins are relatively rigid, justifying the rigid protein approach for this case. At last, a significant variety of active site shapes is observed for these complexes. Therefore, the trypsin has a deep and narrow active site, while the streptavidin has a more open one. Benzamidine is a relatively small ligand, while the other ones are medium-sized. For benzamidine, the main binding motive is the electrostatic interaction of its positively charged amidine group with negatively charged Asp189 and polar atoms of Gly219 and Ser190. As a result of this interaction, a salt bridge and two charged hydrogen bonds are formed. Biotin forms hydrogen bonds with Ser27, Tyr23, Ser88, and Ser45. Oseltamivir forms both neutral and charged hydrogen bonds with negatively charged Glu119 and positively charged Arg371, Agr292, and Arg152. The results of TI computations are shown in the table. There, 〈∆Gcav(PL)〉 and δ∆Gcav(PL) denote mean and root-mean-square deviations of ∆Gcav(PL), respectively, and ∆Gcav is the total change of cavity creation/annihilation free energy resulting from a formation of the PL complex. The TI was carried out along both the forward and the reverse paths, the forward direction being defined in agreement with the scheme in Figure 1. Corresponding results in the table are labeled as forw. and rev., respectively. To estimate the statistical error, eight computations
χns ) 1 -
RV Rs
(6)
where effective radii Rs and RV are defined via SASA (S), Rs ) (S/4π)1/2, and volume V, RV ) (3V/4π)1/3. This quantity characterizes quantitatively the cavity shape, measuring the degree of its deviation from a sphere (the value χns ) 0 corresponds to a sphere). The change of the cavity creation/annihilation free energy as a result of accommodation of the ligand into the active site of the protein was also empirically estimated using the SASA values as:
∆GcavS(PL) ) σ(S(PL) - S(P)) ) σ∆S(PL)
(7)
where S(PL) and S(P) are the values of SASA for the PL complex and protein, respectively, and σ is the surface tension, σ ) 0.1036 kcal/mol/Å2 for T ) 298 K41). It should be stressed that the eq 7 provides only part of the cavity effect, and to obtain the total ∆Gcav value, it is necessary to add the -Gcav(L) term. The proportionality of the cavity creation/annihilation free energy to the SASA for macromolecules (effective radius R > 12-15 Å) was found earlier.16 The values 〈∆Gcav〉 were calculated as a sum of values of -〈Gcav(L)〉 and 〈∆Gcav(PL)〉 averaged over the forward and reverse paths, and the error δ∆Gcav was calculated as the square root of the sum of squared errors of the components incorporated into 〈∆Gcav〉. To investigate the dependence of ∆Gcav on Rcut, three structures corresponding to Rcut ) 4, 6, and 8 Å were prepared for the trypsin-benzamidine (3ptb) complex. Also, modeling with all trypsin atoms was performed. We used the rectangular boxes with sizes SX,Y,Z defined as
SX,Y,Z ) S(PL)X,Y,Z + 16 Å where S(PL)X,Y,Z denotes the sizes of the truncated PL complex for every dimension and 16 Å is the doubled cutoff radius for the interaction between water molecules in the TIP4P water model. The box sizes varied from 40 Å for Rcut ) 4 Å to 60 Å when all atom structures were considered. We found 〈∆Gcav(PL)〉 ) 0.27, 0.5, 0.45, and 0.34 kcal/mol for three truncated structures and the full-atoms protein, respectively. For other complexes, similar results were obtained for truncated structures with Rcut ) 6 and 8 Å. The difference between 〈∆Gcav(PL)〉 values is less than 0.2 kcal/mol.
15358 J. Phys. Chem. B, Vol. 112, No. 48, 2008
Grigoriev et al.
Figure 3. (a) Dependence of water density n(r) on the distance from the surface of the unbound proteins for it is truncated structures with cutoff radii of Rcut ) 4, 6, and 8 Å (1, 2, 3) and for all atoms (4); (b) the average number (Nas) of water particles inside of the active site of the protein and in its vicinity.
Discussion As seen from the table, the values 〈∆Gcav(PL)〉 are relatively small for all PL complexes under investigation. There are no essential differences of the 〈∆Gcav(PL)〉 values for all values of the Rcut. To clarify this result, we calculated the density profiles for water around the proteins without ligands (see Figure 3a). To obtain the distance dependence of the density n(r), the space around the protein was partitioned into shells with thicknesses of 0.1 Å. To define n(r), we calculated the average number of water molecules dN(r), which are separated from the closest atom of the protein by the distance in between r and r + dr. The volume of the corresponding shell was calculated as dV(r) ) dr(S(r) + S(r + dr))/2, where S(r) and S(r + dr) are the areas of the surfaces, consisting of spheres with radii Rjcav + r and Rjcav + r + dr, respectively (see the Methods Section). The density was defined as n(r) ) (dN(r))/(dV(r)). To validate this methodology, the density profiles for a sphere with a radius of 5 Å were calculated. The positions and sizes of the maxima and minima coincide with earlier reported results.24 It is seen from Figure 3a that the water density falls off while approaching the protein surface. This effect becomes more appreciable with an increase in the size of a solute. Similar behavior of the density was obtained from the RDF modeling for spheres of varying radii.17,24 Its interpretation as water
depletion from the layer around a hydrophobic surface was suggested earlier.17 Such a general concept is useful for understanding the degree of occupation of the protein active site by water molecules in the absence of a ligand. Let us denote as Nas the average number of water particles inside of the active site of the protein and in its vicinity. Its dependence on the distance separating a water molecule from the boundary of the active site is illustrated by Figure 3b. At r ) 0, Nas means the number of waters occupying the ligand volume in the active site of the protein. Practically no water molecules (Nas < 0.5 for all cases) are observed when the ligand is withdrawn. This observation continues in the close vicinity of the active site. Only when r > 1 Å does the Nas value reaches unity. On the other hand, the number of water molecules NL which can be located inside of the volume occupied by a ligand at normal water density can be estimated as NL ) VL/30 Å3, where VL is the ligand volume and 30 Å3 is the volume per one water molecule. Let us notice that the VL value should be calculated as the volume inside of the surface which is constructed on Rj radii (not on Rsol + Rj; see the Methods Section). We obtained NL ) 4.6, 7.5, and 8.1 for ligands belonging to 3ptb, 2ht7, and 1swt complexes, respectively. Thus, the water density in the active site of the protein when the ligand
Cavity Effect to Protein-Ligand Binding
J. Phys. Chem. B, Vol. 112, No. 48, 2008 15359
TABLE 1: The Calculated Components of Cavity Creation/Annihilation Free Energies, Their Statistical Errors (δ) in kcal/mol, Change of the SASA ∆S(PL) (in Å2), Excluded Volume V(L) (in Å3), and the Nonspherical Index χns(for comments, see the text) ∆Gcav(PL) comp. 1swt
3ptb
2ht7
〈∆Gcav(PL)〉 ( δ∆Gcav(PL)
∆GcavS(PL), ∆S(PL)
forw.
rev.
forw.
rev.
0.93 0.52 1.42 0.94 0.81 0.74 1.41 1.28 0.16 0.18 0.20 0.17 0.66 0.29 0.26 0.24 0.19 0.22 0.18 0.30 0.20 0.13 0.22 0.17
1.10 1.21 0.92 0.92 1.05 0.83 1.12 1.20 0.30 0.02 0.25 0.40 0.40 0.12 0.05 0.01 0.19 0.02 0.40 0.41 0.21 0.02 0.32 0.01
1.01 ( 0.31
1.04 ( 0.13
-3.9, -39.0
0.27 ( 0.15
0.19 ( 0.15
-0.7, -7.1
0.20 ( 0.05
0.20 ( 0.02
2.9, 29.2
is withdrawn is essentially less in comparison with the density of the bulk water. As a result, absolute values of the free energy change (i.e., 〈∆Gcav(PL)〉) accompanying the corresponding void formation are small. Indeed, the free energy grows up owing to the depletion of water from the cavity volume. Provided protein/ water interactions are treated in terms of the potential in eq 4, the water occupation of the active site remains extremely weak. This is why the energetic cost of its total elimination is negligible. Therefore, according to eq 1, the main contribution to 〈∆Gcav〉 arises due to the annihilation of the ligand’s cavity in water (-Gcav(L)). It is interesting to compare the cavity free-energy changes calculated using MC simulations with their empirical estimates based on the excluded molecular volume V (eq 5) and the SASA computation (eq 7). As shown in the table, the changes in SASA (and in ∆GcavS(PL)) due to the insertion of the ligand into the active site of the protein ∆S(PL) are relatively small, and they can be either negative or positive. The values of ∆Gcav(PL) obtained in MC simulations are always positive, and their absolute values are less than ∆GcavS(PL) for all PL complexes. The negative sign of the ∆S(PL) can be explained by reducing the curvature of a surface of the active site of the protein caused by the ligand if its shape corresponds precisely to the shape of the pocket of the active site. On the other hand, the insertion process for the ligand at the active site is energetically unfavorable, being independent of the conformity of protein and ligand shapes, because the excluded volume for solvent molecules is always created. It seems that negative values of ∆GcavS(PL) appear as the artifact of the calculation method (eq 7) for the change in the cavity free energy. Besides, for the case of the unbound protein, the solvent-accessible surface has a negative value of the curvature radius near the active site (as shown in Figure 1). For such a case, the empirical approach based on the calculation of SASA is invalid. The excess of absolute values of ∆GcavS(PL) over ∆Gcav(PL) can be explained by introducing the surface tension value σ as a proportionality factor between ∆GcavS(PL) and ∆S(PL). Indeed, applying σ in eq 7 is justified by relatively large sizes of both
〈Gcav(L)〉 ( δGcav(L)
Gcav(L)
∆GcavV(L), V(L), χns
〈∆Gcav〉 ( δ∆Gcav
forw.
rev.
forw.
rev.
28.2 28.3 29.8 30.1 28.4 29.3 30.4 29.1 18.8 19.2 19.3 19.1 19.8 19.9 20.5 18.8 36.2 34.4 32.3 35.1 31.8 36.4 34.6 36.7
28.5 28.2 28.8 29.1 29.0 29.0 29.2 29.5 19.2 19.0 19.1 19.2 19.4 19.5 18.8 19.4 34.1 36.2 33.2 32.9 35.1 31.2 33.7 35.4
29.2 ( 0.8
28.9 ( 0.4
28.1, 730.2, 0.077
-28.0 ( 0.5
19.4 ( 0.6
19.2 ( 0.2
17.8, 465.1, 0.045
-18.9 ( 0.3
34.7 ( 1.7
34.0 ( 1.5
34.5, 905.1, 0.082
-34.2 ( 1.1
the PL complex and the protein; therefore, we can invoke the macroscopic continuum theory for treating the cavity creation/ annihilation effect.16 On the other hand, the change in the surface shape due to insertion of the ligand is relatively small, and it is probable that the corresponding process cavity effect should be described in terms of microscopic approach, implying a reduction of the empirical proportionality factor σ. As is shown in the table, the values of 〈Gcav(L)〉 and ∆GcavV(L) for ligands differ slightly; therefore, the approximation in eq 5 works well for the set of considered ligands. It should be stressed that in the present work, only the case of the rigid ligand and protein was studied. We hope that this simplification can be used for small ligands forming complexes with proteins with nonflexible active sites (for instance, the trypsin-benzamidine complex). In opposite cases, it is necessary to take into account the contribution to the cavity effect from both complex compounds. It seems that the relatively simple way to estimate the corresponding effect is to evaluate Gcav(L) and use eq 5 via the molecular volume and SASA for different conformations of the ligand (in order to avoid time-consuming MC simulations). The change in ∆Gcav(PL) can be estimated in a similar way by means of eq 7. Thus, as follows from our results, the values of ∆Gcav are negative and favorable for binding in all PL complexes under investigation. It is worth a reminder that the cavity effect can be a source of positive compensation and unfavorable for binding changes in translation/rotation fractions of free energy ∆Gtr. These changes can be calculated in the frame of the ideal gas approach4,42
(
∆Gtr ) RT 1 +
(
)
2πmkBT 3 1 - ln F + ln(πIAIBIC) + ln 2 2 h2 8π2kBT 3 ln 2 h2
( ))
where m and IA,B,C are the mass and principal inertia moments of the ligand, h is Plank’s constant, and F is the number density at 1mol/L. The obtained vales constitute 18.2, 20.3, and 20.6
15360 J. Phys. Chem. B, Vol. 112, No. 48, 2008 for benzamidine, biotin, and oseltamivir, respectively (in kcal/ mol). Caused by weak logarithmic dependence of ∆Gtr on the mass and principal inertia moments, the differences in the ∆Gtr for ligands prove to be small. While ∆Gcav is roughly proportional to the ligand volume, with increasing size of the ligand, the absolute value of the ∆Gcav will increase faster than the ∆Gtr. Therefore, the compensation effect will be more essential for ligands with relatively large size. Conclusions We investigated the cavity creation/annihilation effect in the context of the formation of the protein-ligand (PL) complexes. The corresponding contributions to the free energy of the PL complex formation ∆Gcav are calculated for three PL complexes, 1swt (streptavidin-biotin), 3ptb (benzamidine-trypsin), and 2ht7 (neuraminidase-oseltamivir). The protein and ligand were considered as rigid structures. Calculation of ∆Gcav was performed using the thermodynamic integration procedure with the original two-step algorithm for growing the interaction potential between the cavity and the water molecules. The Monte Carlo sampling was used for the search in the configurational space of the system. For all complexes under investigation, the calculated values of ∆Gcav are negative and favorable for binding. The change in the cavity free energy due to the annihilation of the ligand’s cavity at the unbound state (denoted as -Gcav(L)) contributes most significantly to the total change of the cavity free energy of the protein-ligand complex, while the contribution from the insertion of the ligand at the active site of the protein is relatively small. The value of -Gcav(L) can easily be estimated empirically in terms of the cavity volume and nonspherical index. The cavity effect can play an important role in compensation of the unfavorable binding changes in the translation/rotation free energy ∆Gtr. The absolute values of ∆Gcav for small ligands like benzamidine are expected to be close to ∆Gtr, whereas for ligands with larger sizes, the ∆Gcavcan exceed ∆Gtr significantly. Therefore, the compensation effect will be more essential for the medium-sized and large ligands. Acknowledgment. This work was supported by the Russian Foundation of Fundamental Research (Project No. RFBR 0603-33171). References and Notes (1) Schapira, M.; Totrov, M.; Abagyan, R. J. Mol. Recognit. 1999, 12, 177. (2) So, S.-S.; Karplus, M. J. Comput. Aided. Mol. Des. 1999, 13, 243. (3) Huang, S.-Y.; Zou, X. J. Comput. Chem. 2006, 27, 1866. (4) Schwarzl, S. M.; Tschopp, T. B.; Smith, J. C.; Fischer, S. J. Comput. Chem. 2002, 23, 1143.
Grigoriev et al. (5) Nikitina, E.; Sulimov, V.; Grigoriev, F.; Kondakova, O.; Luschekina, S. Int. J. Quantum Chem. 2006, 106, 1943. (6) Romanov, A. N.; Jabin, S. N.; Martynov, Y. B.; Sulimov, A. V.; Grigoriev, F. V.; Sulimov, V. B. J. Phys. Chem. A 2004, 108, 9324. (7) Carlsson, J.; Ander, M.; Nervall, M.; Aqvist, J. J. Phys. Chem. B. 2006, 110, 12034. (8) Luo, H.; Sharp, K. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 10399. (9) Li, Z.; Lazaridis, T. J. Phys. Chem. B 2006, 110, 1464. (10) Grigoriev, F. V.; Luschekina, S. V.; Romanov, A. N.; Sulimov, V. B.; Nikitina, E. A. Biochemistry (Moscow) 2007, 72, 785. (11) Carlsson, J.; Aqvist, J. Phys. Chem. Chem. Phys. 2006, 8, 5385. (12) Ruvinsky, A. J. Comput. Chem. 2007, 28, 1364. (13) Hummer, G.; Garde, S.; Garcia, A. E.; Pratt, L. R. Chem. Phys. B 2000, 258, 349. (14) Sun, S. X. Phys. ReV. E 2001, 64, 021512. (15) Wolde, P. R.; Sun, S. X.; Chandler, D. Phys. ReV. E 2001, 65, 011201. (16) Huang, D. M.; Geissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704. (17) Huang, D. M.; Chandler, D. J. Phys. Chem. B 2002, 106, 2047. (18) Gallicchio, E.; Kubo, M. M.; Levy, R. J. Phys. Chem. B 2000, 104, 6271. (19) Alexandrovsky, V. V.; Basilevsky, M. V.; Leontyev, I. V.; Mazo, M. A.; Sulimov, V. B. J. Phys. Chem. B 2004, 108, 15830. (20) Basilevsky, M. V.; Grigoriev, F. V.; Leontyev, I. V.; Sulimov, V. B. J. Phys. Chem. A 2005, 109, 6939. (21) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (22) Hyang, D. M.; Gaissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704. (23) Floris, F. M. J. Phys. Chem. B 2005, 109, 24061. (24) Floris, F. M. J. Chem. Phys. 2007, 126, 74505. (25) Grigoriev, F. V.; Basilevsky, M. V.; Gabin, S. N.; Romanov, A. N.; Sulimov, V. B. J. Phys. Chem. B 2007, 111, 13748. (26) Floris, F. M.; Selmi, M.; Tani, A.; Tomasi, J. J. Chem. Phys. 1997, 107, 6353. (27) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235. (28) Grigoriev, F. V.; Romanov, A. N.; Sulimov, V. B. Numerical Methods and Programming 2007, 8, 334 In Russian. (29) Kirkwood, J. G.; Theory of Liquids; Alder. B. J., Ed.; Gordon and Breach: New York, 1968. (30) Bordner, A. J.; Cavasotto, C. N.; Abagyan, R. A. J. Phys. Chem. B 2002, 106, 11009. (31) McCarrick, M. A.; Kollman, P. A. J Comput. Aided Mol. Des. 1999, 13, 109. (32) Halgren, T. A. J. Comput. Chem. 1996, 17, 490. (33) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1989. (34) Metropolis, N.; Ulam, S. J. Am. Stat. Assoc. 1949, 44, 335. (35) Wood, W. W. J. Chem. Phys. 1968, 48, 415. (36) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (37) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027. (38) High performance Internet accessible Computing Cluster for research and teaching. http://parallel.ru/cluster/. (39) The supercomputer centre RAS http://www.jscc.ru/. (40) Word, J. M. J. Mol. Biol. 1999, 285, 1735. (41) CRC Handbook of Chemistry and Physics, 78th ed.; CRC Press: Boca Raton, FL, 1997. (42) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976.
JP8041439
