J. Phys. Chem. B 2007, 111, 11585-11591
11585
Three-Dimensional Distribution Function Theory for the Prediction of Protein-Ligand Binding Sites and Affinities: Application to the Binding of Noble Gases to Hen Egg-White Lysozyme in Aqueous Solution Takashi Imai,*,† Ryusuke Hiraoka,† Tomoyoshi Seto,‡ Andriy Kovalenko,§ and Fumio Hirata| Department of Bioscience and Bioinformatics, Ritsumeikan UniVersity, Kusatsu, Shiga 525-8577, Japan, Department of Anesthesiology, Shiga UniVersity of Medical Science, Otsu, Shiga 520-2192, Japan, National Institute for Nanotechnology, National Research Council of Canada, and Department of Mechanical Engineering, UniVersity of Alberta, 11421 Saskatchewan DriVe, Edmonton, Alberta T6G 2M9, Canada, and Department of Theoretical Studies, Institute for Molecular Science, Okazaki, Aichi 444-8585, Japan ReceiVed: June 22, 2007; In Final Form: July 25, 2007
The three-dimensional distribution function theory of molecular liquids is applied to lysozyme in mixtures of water and noble gases. The results indicate that the theory has the capability of predicting the protein-ligand binding sites and affinities. First, it is shown that the theory successfully reproduces the binding sites of xenon found by X-ray crystallography. Then, the ability of the theory to predict the size selectivity of noble gases is demonstrated. The effect of water on the selectivity is clarified by a theoretical analysis. Finally, it is demonstrated that the dose-response curve, which is employed in experiments for examining the binding affinity, is realized by the theory.
1. Introduction Molecular recognition by protein, or protein-ligand binding, is one of the most fundamental functions of proteins in biological processes. In addition to scientific interest, prediction of ligand binding sites and affinities is the starting point for drug discovery.1,2 Therefore, a large number of computational methods as well as experimental approaches have been proposed.1-6 The computational methodologies are divided into two categories, or stages. One is the prediction of ligand binding sites in the target protein. The binding sites are located, in the most common case, based on a purely geometric analysis of the protein structure, in which cavities or clefts in the protein are detected and regarded as potential binding sites.3 The binding sites can also be predicted by bioinformatics from multiple alignments of the amino acid sequences in the protein family.7 The other is the docking of a ligand molecule at the binding sites which are already known or predicted in advance. Possible docking structures are then evaluated on the basis of a force field or a scoring function.4-6 Although such docking programs are increasingly popular among the fields of bioscience and pharmacology,8 theoretical methodologies are not fully developed. One of the least developed questions is how to incorporate the effect of water into the binding affinity or free energy. Water participates in protein-ligand binding in the following two ways. Primarily, bulk water provides the reaction field acting on the binding. This effect includes the electrostatic screening and the hydrophobic interaction between protein and ligand molecules. * Author to whom correspondence should be addressed. E-mail: t-imai@ is.ritsumei.ac.jp. † Ritsumeikan University. ‡ Shiga University of Medical Science. § National Research Council of Canada, and Department of Mechanical Engineering, University of Alberta. | Institute for Molecular Science.
Moreover, individual water molecules can act as integral molecular components of the complex.9-11 In fact, water molecules are often found at the binding interface of proteinligand complexes, mediating the hydrogen bonds or simply filling void spaces. In spite of the evident significance of such water molecules, the effect of water is usually treated at the level of continuum solvation models,4-6 unless the interfacial water molecules are introduced in advance. Molecular simulation with explicit water molecules can, in principle, incorporate the effects of water stated above into the protein-ligand binding, if the simulation is performed correctly so as to sample the whole configuration space of water, including that inside the cavity and cleft of the protein. The binding free energy can be calculated accordingly, using advanced simulation techniques such as free energy perturbation and thermodynamic integration methods.12 However, such a free energy calculation requires very high computational expenses, even though the binding site and the complex structure are known in advance. If we have no spatial and structural information about the binding site, such a rigorous simulation is obviously no longer applicable to predict protein-ligand binding sites and affinities. An alternative approach is to employ distribution function theory of liquids based on statistical mechanics. Among several versions of this methodology, the three-dimensional reference interaction site model (3D-RISM) theory of molecular liquids13-16 is the most advanced and promising. In this theory, protein is immersed into the water-ligand mixture of a finite concentration and the 3D-spatial distribution function of ligand around the protein is obtained by solving the 3D-RISM statisticalmechanical equations. The ligand binding sites and affinities are obtained from the positions and heights of the peaks of the 3D distribution function, respectively. This calculation does not require any information in advance on the ligand binding sites and the water configurations around them. If a water-mediated
10.1021/jp074865b CCC: $37.00 © 2007 American Chemical Society Published on Web 09/08/2007
11586 J. Phys. Chem. B, Vol. 111, No. 39, 2007
Imai et al.
binding is stable, the corresponding peak of the distribution function is obtained and identified. It should be noted here that to obtain such distribution functions by molecular simulation is possible in theory, but not feasible in practice. We have recently applied the 3D-RISM theory to proteins in aqueous solutions and succeeded in obtaining their hydration structure and thermodynamic quantities.17,18 The most important finding of our interests in the present paper is that the theory reproduces the 3D-spatial distribution of water molecules isolated from bulk phase and confined in the protein interior as well as the ordinary hydration structure around the protein.19,20 Furthermore, the 3D-RISM theory was applied to a protein in aqueous electrolyte solutions, and the selective ion binding to a cleft of protein was successfully predicted.21,22 Those studies strongly imply that the theory is capable of predicting the watermediating protein-ligand binding sites and affinities described above. In this contribution, we present a preliminary study for the prediction of protein-ligand binding sites and affinities by the 3D-RISM theory. For this purpose, we choose the binding of noble gases to hen egg-white lysozyme as a simple model system, because the 3D structure of the lysozyme-xenon complex has been resolved by X-ray analysis. First, we demonstrate that two peaks of the distribution function of xenon are in agreement with the two binding sites of xenon found in the X-ray structure. Second, we discuss the size dependence of the coordination numbers of gases at the sites, which represents ligand-size selectivity in the binding. The effects of water on the binding affinities are also evaluated and discussed. Finally, the concentration dependence of the binding affinities of gases, which corresponds to the dose-response curve used in biochemistry, is examined. In this study, rather than to quantitatively compare the theoretical results with experimental data, our concern is to confirm the 3D-RISM theory is practically capable of predicting protein-ligand binding sites and affinities in the qualitative aspect. 2. Methods 2.1. 3D-RISM Theory. The detailed description of the 3DRISM theory is given in previous papers.13-16 Here, we provide only a brief outline of the theory. In this study, we consider lysozyme as solute and the mixture of water and noble gas as solvent. For the solute-solvent system at infinite dilution, the 3D-RISM integral equation
hγ(r) ) vv vv (|r′ - r|) + Fγ′hγ′γ (|r′ - r|))dr′ ∫ cγ′(r′)(wγ′γ ∑ γ′
(1)
is solved simultaneously with the 3D hypernetted chain (HNC) approximation
hγ(r) ) exp(-βuγ(r) + hγ(r) - cγ(r)) - 1
(2)
to obtain the total and direct correlation functions hγ(r) and cγ(r) of solvent site γ (water (oxygen and hydrogen), and noble gas) in the 3D space around the solute. In the above equations, uγ(r) is the interaction potential between solvent site γ and the vv vv (r) and hγ′γ (r) are whole solute specified on the 3D grid, wγ′γ respectively the site-site intramolecular and site-site total correlation functions of solvent, Fγ is the number density of molecular species to which solvent site γ belongs, and β ) 1/(kBT) is the inverse of the Boltzmann constant times temperature, all of which are the input to the 3D-RISM theory. The
TABLE 1: Lennard-Jones Potential Parameters for Noble Gases Nea G1b Ara a
σ/Å
/kcal mol-1
3.035 3.225 3.415
0.0369 0.1480 0.2484
From Guillot and Guissani.27
G2b Kra Xea b
σ/Å
/kcal mol-1
3.545 3.675 3.975
0.2963 0.3358 0.4267
Hypothetical gases (see text).
vv (r) are site-site total correlation functions of solvent hγ′γ calculated in advance from the ordinary RISM theory for radial correlations, applied to the solvent mixture. The 3D potential function uγ(r) is calculated on the supercell grid using the minimum image convention and the Ewald summation methods, from the conventional site-site pair potential which consists of the Lennard-Jones (LJ) and Coulombic terms. The convolution integral in eq 1 is calculated by the 3D fast Fourier transform technique. The long-range electrostatic asymptotics of all the correlation functions in eqs 1 and 2 are separated out and treated analytically, including the special corrections for the finiteness of the 3D supercell.13,16 2.2. Models and Parameters. The 3D atomic coordinates of hen egg-white lysozyme were taken from Protein Data Bank, code 1C10, which was determined by X-ray crystallographic analysis for a lysozyme-xenon complex.23 We employed Amber99 parameters24 for lysozyme and the SPC/E model25 for water, with a correction concerning the LJ parameters of the hydrogen sites26 (σ ) 0.4 Å and ) 0.05 kcal mol-1). The parameters of Guillot and Guissani27 were used for noble gases, Ne, Ar, Kr, and Xe. In this study, we additionally defined two hypothetical gases in order to investigate the size dependence in detail. One is of the middle size (LJ-σ) between Ne and Ar. The other is between Ar and Kr. They are referred to as G1 and G2, respectively. The values of LJ- for G1 and G2 were determined from the regression curve of on σ for the other four gases. The LJ potential parameters for the noble gases are compiled in Table 1. The LJ potential parameters between different molecular species were obtained from the LorentzBerthelot mixing rules. We use the mixtures of five different concentrations, 0.001, 0.01, 0.1, 0.5, and 1 M, of each gas at a temperature of 298.15 K for investigating the concentration dependence of the binding affinities. Since the solubility of noble gases is in the order of 10-3 M (0.5 × 10-3 M, 1.4 × 10-3 M, 2.5 × 10-3 M, and 4.4 × 10-3 M for Ne, Ar, Kr, and Xe, respectively),28 the mixtures of the higher concentrations are hypothetical ones to investigate the concentration dependence in wide range. Therefore, we determined the number densities of water and noble gases on the basis of the assumption of ideal mixing wherein the total number density of a mixture is constant at the value for pure water: Fwater + Fgas ) 0.033329 Å-3. It is worthwhile to briefly mention the equilibrium structures of the aqueous noble-gas solutions at the high concentrations. A recent molecular simulation study pointed out that methane forms clusters in aqueous solution at a high concentration of 2.5 M.29 If this is the case, the concentration dependences of the binding affinities would be much more complicated than we expected. However, it is found that the pair distribution functions of each noble gas are almost the same and indicate no unusual solution structure in the range of the concentrations adopted in this study. The 3D-RISM/HNC equations were solved on a grid of 128 × 128 × 128 points in a cubic supercell of size 64 × 64 × 64
Prediction of Protein-Ligand Binding Sites and Affinities
J. Phys. Chem. B, Vol. 111, No. 39, 2007 11587
Figure 1. Isosurface representation of the three-dimensional distribution functions of xenon and water (oxygen and hydrogen) around hen eggwhite lysozyme. The yellow, red, and white surfaces (or spots) show the areas where the distribution functions g(r) are greater than 8 for xenon, oxygen, and hydrogen, respectively. The blue surface is the molecular surface of lysozyme. Part (a): The overall view. Part (b): The distributions at the substrate binding site. In the lower panel, the isosurface of g(r) > 80 instead of g(r) > 8 is presented for xenon in order to disclose the highest peak in the region, which is indicated by the yellow arrow. The orange sphere represents the xenon found in the X-ray structure. Part (c): The distributions at the internal site. In the lower panel, the isosurface of g(r) > 0.02 instead of g(r) > 8 is presented for xenon in order to show the minor peak in the region, which is indicated by the yellow arrow. The orange sphere represents the xenon found in the X-ray structure.
11588 J. Phys. Chem. B, Vol. 111, No. 39, 2007
Imai et al.
Å3, which is sufficiently large to accommodate the protein molecule with enough solvation space. 2.3. Calculation of the Coordination Number at Binding Site. The coordination number of noble gas (ligand) at each binding site is calculated by
N ) Flig
∫V∈binding site glig(r)dr
(3)
where the integral ranges over the space identified as the binding site. The integration range for the internal site (see Results and Discussion) is apparent because it is completely isolated from the bulk phase. On the other hand, it is difficult to define the integration range for the substrate binding site (see Results and Discussion) because the distribution of noble gas at the site overlaps with that in the bulk phase. Therefore, we calculated the coordination number at the substrate binding site in the following way. First, the 3D distribution function is converted to the radial distribution function g(r). Next, 4πr2g(r) is fitted by several Gauss functions. Then, the Gauss functions identified as the peaks within the binding site are extracted. Finally, the Gauss functions are integrated to obtain the coordination number. 3. Results and Discussion 3.1. Three-Dimensional Distribution of Noble Gases around Protein. Figure 1a shows the 3D distribution functions of xenon and water (oxygen and hydrogen) around lysozyme, calculated by the 3D-RISM theory for lysozyme in water-xenon mixture at the concentration of 0.001 M. The molecular surface of the protein is colored blue. The regions where g(r) > 8 is colored different colors for different species: yellow for xenon, red for water oxygen, and white for water hydrogen. Of course, the surface painted blue is covered by water molecules weakly bound to the protein, which are not shown. A number of welldefined peaks, yellow and red spots, are found for xenon and water oxygen at the surface of the protein, which are separated from each other. The result demonstrates the great capability of the 3D-RISM theory to predict the “preferential binding” of ligands, which is nontrivial for any theoretical methods including molecular simulation. The distributions of ligand and water are simultaneously found in this result, which means that the peak of either the ligand or water is found at each site, depending on the ratio of their affinities to the site. Actually, Figure 1a indicates that there are the water- and xenon-preferred sites on the protein surface. Similar results are obtained for the other noble gases and the other concentrations. It is interesting to compare the distribution of xenon obtained by the 3D-RISM theory with the xenon sites in the X-ray structure,23 even though their conditions are different; the former is an aqueous solution under atmospheric pressure, whereas the latter is crystal under xenon gas pressure of 12 bar. There are two binding sites of xenon in lysozyme: one corresponds to the binding pocket of native ligands, which is referred to as the substrate binding site, and the other is located in a cavity inside of the protein, which is referred to as the internal site.23 Figure 1b compares the theoretical result of the 3D distribution of xenon with the X-ray xenon site at the substrate binding site. The location of a high and sharp peak found by the theory is in complete agreement with the X-ray xenon site. Figure 1c shows the result at the internal site. The xenon peak found there is actually a minor one; nevertheless, the location is again consistent with the X-ray site. It is interesting to note that the peaks of water are shifted away from the xenon binding site.
Figure 2. Coordination numbers of noble gases at the binding sites of lysozyme in the 0.001 M aqueous noble-gas solutions, plotted against the atomic diameter (Lennard-Jones σ) of the gases. (a) The substrate binding site. (b) The internal site.
3.2. Size Dependence of the Coordination Number of Noble Gases at the Binding Sites. Figure 2 shows the size dependence of the coordination number of noble gases at the two binding sites, which is calculated at the concentration of 0.001 M. At the substrate binding site, the coordination number becomes exponentially larger as the size of noble gas increases. At the internal site, the coordination number becomes larger with increasing the gas size up to σ ≈ 3.4 Å, while it decreases in the region where σ > 3.4 Å. As a result, argon has the largest binding affinity to the internal site. These results demonstrate that the 3D-RISM theory has the ability to describe ligand size selectivity in the binding, or molecular recognition. Although there are no corresponding experimental data, the present results serve as a representative test case. We now discuss the effect of water on the binding affinities of noble gases, which will be the true touchstone of the theory. To evaluate the effect of water, we define a coefficient given by
γ)
Flig N ) Nid Flig
∫V∈binding site glig(r)dr ∫V∈binding site gidlig(r)dr
(4)
where Nid and gid lig(r) are respectively the coordination number and the distribution function of noble gas (ligand) in the ideal state where no water exists. Thus, the coefficient γ indicates how significant the effect of water is in increasing (γ > 1) or
Prediction of Protein-Ligand Binding Sites and Affinities
J. Phys. Chem. B, Vol. 111, No. 39, 2007 11589
Figure 3. The ideal coordination numbers of noble gases at the binding sites and the effects of water for the 0.001 M solutions plotted against the atomic diameter (Lennard-Jones σ) of the gases. (a) The coordination number in the ideal state and in water at the substrate binding site. (b) The water-effect coefficient at the substrate binding site. (c) The coordination number in the ideal state and in water at the internal site. (d) The watereffect coefficient at the internal site.
decreasing (γ < 1) the coordination number relative to the ideal state. In this regard, we name γ as “water-effect coefficient.” The ideal distribution function is defined by
gid lig(r) ) exp(-βulig(r))
(5)
where ulig(r) is the direct interaction potential between protein and noble gas (ligand). The distribution function can be written using the corresponding expression
glig(r) ) exp(-βulig(r) + w(r))
(6)
where w(r) is the water-mediated interaction. (It should be noted that the correlation between noble gases is also involved in w(r) under a high concentration condition.) Substituting eqs 5 and 6 into eq 4, the water-effect coefficient is rewritten as
γ)
∫V∈binding site exp(-βulig(r) + w(r))dr ) 〈exp(w(r))〉 ∫V∈binding site exp(-βulig(r))dr
(7)
where the bracket denotes the potential weighed average within the binding site. The size dependences of the ideal coordination number and the water-effect coefficient, γ, are shown in Figure 3. At the substrate binding site, the coordination number for all the noble gases is reduced by water (Figure 3a). The reduction
effect becomes larger as the gas size increases (Figure 3b). The result is explained qualitatively as follows. Since the substrate binding site is hydrophilic, water prefers binding to the site more than noble gases do. Water has to expel noble gases from the site to interact with the site. The effect can increase with increasing gas size because the larger the gas size becomes, the more the site is shielded from water. In contrast, the coordination number is increased by water at the internal site (Figure 3c). The effect is more enhanced by an increase of the gas size (Figure 3d). Since the internal site is located in a hydrophobic cavity, it is reasonable to consider this effect as the hydrophobic effect. The noble gases larger than G2 are destabilized by some van der Waals overlap with the protein in the ideal state. Water moderates the destabilization by the hydrophobic effect. As a result, the water effect shifts the peak of the binding affinity to the larger gas, increasing the overall binding affinity. 3.3. Concentration Dependence of the Coordination Number of Noble Gases at the Binding Sites. It is well-known that the activity of protein plotted against the logarithm of ligand concentration generally produces a sigmoidal curve, which is the so-called dose-response curve. Experimentalists use the sigmoidal dose-response curve to obtain the equilibrium constant of the protein-ligand binding and the binding free energy. Following the experimental procedure, we plotted the coordination number of each noble gas against the logarithm of the gas concentration. The curves for the two binding sites
11590 J. Phys. Chem. B, Vol. 111, No. 39, 2007
Imai et al. effect of water through the theoretical framework. By obtaining the binding affinities or the coordination numbers of a series of ligands at the binding sites, the theory describes the dependence of the ligand properties on the binding affinities, i.e., the molecular recognition by the protein. The effect of water upon the molecular recognition is extracted in the theoretical manner. The theory can also realize the dose-response curve, which is used for calculating the equilibrium constant and the binding free energy. Applying the procedure described in this report to other ligand molecules instead of noble gases and another protein instead of lysozyme, we can predict the ligand binding sites and affinities for a given system, including the effect of water. We have already succeeded in using the 3D-RISM theory for a protein in aqueous electrolyte solutions21,22 as well as mixtures of water with small organic compounds such as urea and glycerol.30 Of course, the description of a ligand as cosolvent in this theory can be applied to larger molecules of drug compounds, too. However, calculation for a mixture including such large molecules is not as straightforward in practice and requires further improvement of the site-site approach. Work in this direction is in progress in our group.
Figure 4. Coordination numbers of noble gases at the binding sites plotted against the logarithm of gas concentration. Similar relations are found for G1 and G2. (a) The substrate binding site. (b) The internal site.
are shown in Figure 4. In this study, the complete sigmoidal curves were not obtained because the affinities between the protein and noble gases are considerably weak. Each of the plots probably exhibits the low-concentration region of the putative sigmoidal curve. If we could draw the complete sigmoidal curves, we would have succeeded in obtaining the equilibrium constants and the binding free energies without calculating the free energy directly. It should be emphasized here that the production of the doseresponse curve can be achieved only if the employed method can treat a highly dilute mixture, because the typical equilibrium constant of ligand binding is in the order of µM. The ordinary molecular simulation would never cover such highly dilute conditions. (For example, we have to prepare about 55,000,000 water molecules per one ligand molecule in the simulation box for a mixture of 1 µM.) In the 3D-RISM theory, the calculation can be done at an arbitrary concentration, just by setting the value of component density F in the equation. 4. Concluding Remarks In this contribution, we have tested the 3D-RISM theory as a theoretical tool for prediction of protein-ligand binding sites and affinities. The present study for the binding of noble gases to lysozyme has demonstrated the following capabilities of the 3D-RISM theory. The theory detects the ligand binding sites as well-defined peaks of the 3D distribution function of ligand. It should be emphasized that the result naturally includes the
Acknowledgment. This work was supported in part by Grant-in-Aid for Scientific Research on Priority Areas (No. 15076212) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan and by the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan. T.I. thanks the support by “High-Tech Research Center” Project for Private Universities, MEXT, Japan. A.K. acknowledges the support from the National Research Council of Canada. The isosurface illustrations in Figure 1 were produced with the visualization program VMD.31 References and Notes (1) Kitchen, D. B.; Decornez, H.; Furr, J. R.; Bajorath, J. Nat. ReV. Drug DiscoVery 2004, 3, 935. (2) Klebe, G. Drug DiscoVery Today 2006, 11, 580. (3) Sotriffer, C.; Klebe, G. Farmaco 2002, 57, 243. (4) Gohlke, H.; Klebe, G. Angew. Chem., Int. Ed. 2002, 41, 2644. (5) Halperin, I.; Ma, B.; Wolfson, H.; Nussinov, R. Proteins: Struct., Funct., Genet. 2002, 47, 409. (6) Brooijmans, N.; Kuntz, I. D. Annu. ReV. Biophys. Biomol. Struct. 2003, 32, 335. (7) Lichtrage, O.; Sowa, M. E. Curr. Opin. Struct. Biol. 2002, 12, 21. (8) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. Proteins: Struct., Funct., Genet. 2006, 65, 15. (9) Ladbury, J. E. Chem. Biol. 1996, 3, 973. (10) Levy, Y.; Onuchic, J. N. Annu. ReV. Biophys. Biomol. Struct. 2006, 35, 389. (11) Li, Z.; Lazaridis, T. Phys. Chem. Chem. Phys. 2007, 9, 573. (12) Wang, W.; Donini, O.; Reyes, C. M.; Kollman, P. A. Annu. ReV. Biophys. Biomol. Struct. 2001, 30, 211. (13) Kovalenko, A. In Molecular Theory of SolVation; Hirata, F., Ed.; Kluwer: Dordrecht, 2003; p 169. (14) Beglov, D.; Roux, B. J. Phys. Chem. B 1997, 101, 7821. (15) Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 1998, 290, 237. (16) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 10391. (17) Imai, T.; Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 2004, 395, 1. (18) Imai, T.; Kovalenko, A.; Hirata, F. Mol. Simul. 2006, 32, 817. (19) Imai, T.; Hiraoka, R.; Kovalenko, A.; Hirata, F. J. Am. Chem. Soc. 2005, 127, 15334. (20) Imai, T.; Hiraoka, R.; Kovalenko, A.; Hirata, F. Proteins: Struct., Funct., Bioinf. 2007, 66, 804. (21) Yoshida, N.; Phongphanphanee, S.; Maruyama, Y.; Imai, T.; Hirata, F. J. Am. Chem. Soc. 2006, 128, 12042. (22) Yoshida, N.; Phongphanphanee, S.; Hirata, F. J. Phys. Chem. B 2007, 111, 4588. (23) Prange´, T.; Schiltz, M.; Permot, L.; Colloc’h, N.; Longhi, S.; Bourguet, W.; Fourme, R. Proteins: Struct., Funct., Genet. 1998, 30, 61.
Prediction of Protein-Ligand Binding Sites and Affinities (24) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049. (25) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (26) Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1983, 77, 1451. (27) Guillot, B.; Guissani, Y. J. Chem. Phys. 1993, 99, 8075. (28) Krause, D., Jr.; Benson, B. B. J. Solution Chem. 1989, 18, 823.
J. Phys. Chem. B, Vol. 111, No. 39, 2007 11591 (29) Oostenbrink, C.; van Gunsteren, W. F. Phys. Chem. Chem. Phys. 2005, 7, 53. (30) Yamazaki, T.; Imai, T.; Hirata, F.; Kovalenko, A. J. Phys. Chem. B 2007, 111, 1206. (31) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33.
