Computation of the contribution from the cavity effect to protein-ligand binding free energy

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 DeltaG(cav) 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 DeltaG(cav) are negative and favorable for binding. The main contribution to DeltaG(cav) 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.     

In silico prediction of drug solubility: 2. Free energy of solvation in pure melts

The solubility of drugs in water is investigated in a series of papers and in the current work. The free energy of solvation, DeltaG*(vl), of a drug molecule in its pure drug melt at 673.15 K (400 degrees C) has been obtained for 46 drug molecules using the free energy perturbation method. The simulations were performed in two steps where first the Coulomb and then the Lennard-Jones interactions were scaled down from full to no interaction. The results have been interpreted using a theory assuming that DeltaG*(vl) = DeltaG(cav) + E(LJ) + E(C)/2 where the free energy of cavity formation, DeltaG(cav), in these pure drug systems was obtained using hard body theories, and E(LJ) and E(C) are the Lennard-Jones and Coulomb interaction energies, respectively, of one molecule with the other ones. Since the main parameter in hard body theories is the volume fraction, an equation of state approach was used to estimate the molecular volume. Promising results were obtained using a theory for hard oblates, in which the oblate axial ratio was calculated from the molecular surface area and volume obtained from simulations. The Coulomb term, E(C)/2, is half of the Coulomb energy in accord with linear response, which showed good agreement with our simulation results. In comparison with our previous results on free energy of hydration, the Coulomb interactions in pure drug systems are weaker, and the van der Waals interactions play a more important role.     

Calculations of the Gibbs Energy of Formation of Cavities in Water

The Gibbs energy Δ_cav G of formation of cavities in water was calculated by thermodynamic integration. The cavities corresponded to organic molecules of various volumes and shapes with characteristic radii of 3–7 Å and spheres of radii 3–6 Å. Statistical integrals were calculated by the Monte Carlo simulation of an ensemble of water molecules with periodic boundary conditions at 25°C and pressure 1 atm. Interaction between water molecules was described by the TIP4P four-point nonpolarizable model. It was found that the proportionality of Δ_cav G to cavity volume V obtained earlier for spheres with radii not exceeding 5 Å remained valid for cavities corresponding to organic molecules of different volumes and shapes, including cavities with the characteristic radii exceeding 5 Å. A possible explanation of the retention of the Δ_cav G ∞ V dependence at large characteristic radii was suggested in terms of the two-peak binomial model of the cavitation effect. The conclusion was drawn that the Δ_cav G = αV dependence could be used for calculating the nonpolar part of the Gibbs energy of solvation in implicit water models.     

Computation of hydration free energies of organic solutes with an implicit water model

A new approach for computing hydration free energies DeltaG(solv) of organic solutes is formulated and parameterized. The method combines a conventional PCM (polarizable continuum model) computation for the electrostatic component DeltaG(el) of DeltaG(solv) and a specially detailed algorithm for treating the complementary nonelectrostatic contributions (DeltaG(nel)). The novel features include the following: (a) two different cavities are used for treating DeltaG(el) and DeltaG(nel). For the latter case the cavity is larger and based on thermal atomic radii (i.e., slightly reduced van der Waals radii). (b) The cavitation component of DeltaG(nel) is taken to be proportional to the volume of the large cavity. (c) In the treatment of van der Waals interactions, all solute atoms are counted explicitly. The corresponding interaction energies are computed as integrals over the surface of the larger cavity; they are based on Lennard Jones (LJ) type potentials for individual solute atoms. The weighting coefficients of these LJ terms are considered as fitting parameters. Testing this method on a collection of 278 uncharged organic solutes gave satisfactory results. The average error (RMSD) between calculated and experimental free energy values varies between 0.15 and 0.5 kcal/mol for different classes of solutes. The larger deviations found for the case of oxygen compounds are probably due to a poor approximation of H-bonding in terms of LJ potentials. For the seven compounds with poorest fit to experiment, the error exceeds 1.5 kcal/mol; these outlier points were not included in the parameterization procedure. Several possible origins of these errors are discussed.     

On the cavitation energy of water

Free-energy-perturbation theory from molecular dynamics calculations has been used to obtain the DeltaG of adjoining cavities' formation in water. The DeltaGs for systems with three, five and seven cavities are compared with that of a single cavity of the same volume, and found to be in good agreement. The conditions under which the analytical formulation of the energy of cavity formation proposed by Pierotti holds are discussed. The data for a single cavity have been tabulated and can lend themselves to a simple numerical implementation in standard quantum chemical packages, which can be used when high accuracy for DeltaG(cav) is required.     

In silico prediction of drug solubility: 1. Free energy of hydration

As a first step in the computational prediction of drug solubility the free energy of hydration, DeltaG*(vw) in TIP4P water has been computed for a data set of 48 drug molecules using the free energy of perturbation method and the optimized potential for liquid simulations all-atom force field. The simulations were performed in two steps, where first the Coulomb and then the Lennard-Jones interactions between the solute and the water molecules were scaled down from full to zero strength to provide physical understanding and simpler predictive models. The results have been interpreted using a theory assuming DeltaG*(vw) = A(MS)gamma + E(LJ) + E(C)/2 where A(MS) is the molecular surface area, gamma is the water-vapor surface tension, and E(LJ) and E(C) are the solute-water Lennard-Jones and Coulomb interaction energies, respectively. It was found that by a proper definition of the molecular surface area our results as well as several results from the literature were found to be in quantitative agreement using the macroscopic surface tension of TIP4P water. This is in contrast to the surface tension for water around a spherical cavity that previously has been shown to be dependent on the size of the cavity up to a radius of approximately 1 nm. The step of scaling down the electrostatic interaction can be represented by linear response theory.     

Free energy of solvation from molecular dynamics simulation applying Voronoi-Delaunay triangulation to the cavity creation

The free energy of solvation for a large number of representative solutes in various solvents has been calculated from the polarizable continuum model coupled to molecular dynamics computer simulation. A new algorithm based on the Voronoi-Delaunay triangulation of atom-atom contact points between the solute and the solvent molecules is presented for the estimation of the solvent-accessible surface surrounding the solute. The volume of the inscribed cavity is used to rescale the cavitational contribution to the solvation free energy for each atom of the solute atom within scaled particle theory. The computation of the electrostatic free energy of solvation is performed using the Voronoi-Delaunay surface around the solute as the boundary for the polarizable continuum model. Additional short-range contributions to the solvation free energy are included directly from the solute-solvent force field for the van der Waals-type interactions. Calculated solvation free energies for neutral molecules dissolved in benzene, water, CCl4, and octanol are compared with experimental data. We found an excellent correlation between the experimental and computed free energies of solvation for all the solvents. In addition, the employed algorithm for the cavity creation by Voronoi-Delaunay triangulation is compared with the GEPOL algorithm and is shown to predict more accurate free energies of solvation, especially in solvents composed by molecules with nonspherical molecular shapes.     

Thermodynamic Studies on Amino Acid Solvation in Aqueous Urea

The present paper discusses the behaviour of transfer free energy of some amino acids from water to 4M, 6M and 8M aqueous urea. Dissection of transfer free energy into cavity term, interaction term and electrical term reveals that cavity forming free energy of transfer ΔG⁰_t (cav) plays an important role in dictating actual interaction of amino acids in aqueous urea. Cavity forming free energy of transfer has been estimated by using Scaled Particle Theory (SPT).     

Excluded volume effect for large and small solutes in water

The cavitation effect, i.e., the process of the creation of a void of excluded volume in bulk solvent (a cavity), is considered. The cavitation free energy is treated in terms of the information theory (IT) approach [Hummer, G.; Garde, S.; Garcia, A. E.; Paulaitis, M. E.; Pratt, L. R. J. Phys. Chem. B 1998, 102, 10469]. The binomial cell model suggested earlier is applied as the IT default distribution p(m) for the number m of solute (water) particles occupying a cavity of given size and shape. In the present work, this model is extended to cover the entire range of cavity size between small ordinary molecular solutes and bulky biomolecular structures. The resulting distribution consists of two binomial peaks responsible for producing the free energy contributions, which are proportional respectively to the volume and to the surface area of a cavity. The surface peak dominates in the large cavity limit, when the two peaks are well separated. The volume effects become decisive in the opposite limit of small cavities, when the two peaks reduce to a single-peak distribution as considered in our earlier work. With a proper interpolation procedure connecting these two regimes, the MC simulation results for model spherical solutes with radii increasing up to R = 10 A [Huang, D. H.; Geissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704] are well reproduced. The large cavity limit conforms to macroscopic properties of bulk water solvent, such as surface tension, isothermal compressibility and Tolman length. The computations are extended to include nonspherical solutes (hydrocarbons C1-C6).     

Transfer Thermodynamics of Protein in Denaturing and Stabilizing Media

Free energies of transfer (ΔG_t) of RibonucleaseA (RNaseA) from water to aqueous solutions of urea (4 M, 6 M and 8 M), a protein denaturing solvent as well as ΔG_t of RibonucleaseA, β‐Lactoglobulin, α‐Chymotripsin and ChymotrypsinogenA from water to aqueous glycerol (10%, 20%, 30% and 40%), a protein stabilizing solvent has been dissected into cavity term [ΔG_t(cav)] and interaction term [ΔG_t(int)]. The interaction free energy includes all types of interactions like hard‐soft, hydrogen bonding, electrostatic, etc. The cavity forming free energies have been calculated using the standard version of scaled particle theory (SPT) with well‐reported SPT parameters. It has been found that transfer free energies of cavity terms ΔG_t(cav) for native protein from water to urea‐water and water to aqueous glycerol follow almost opposite trends. This primarily indicates there may be some correlation between cavity creation energies and protein denaturing and stabilizing ability of a solvent. The results are in agreement with those obtained from preferential binding coefficient studies in these media.     

On the molecular origins of volumetric data

We use a statistical thermodynamic approach and a simple thermodynamic model of hydration to examine the molecular origins of the volumetric properties of solutes. In this model, solute-solvent interactions are treated as a binding reaction. The free energy of hydration of the noninteracting solute species coincides with the free energy of cavity formation, while the free energy of solute-solvent interactions is given by the binding polynomial. By differentiating the relationship for the free energy of hydration with respect to temperature and pressure, one obtains the complete set of equations describing the thermodynamic profile of hydration, including enthalpy, entropy, volume, compressibility, expansibility, and so forth. The model enables one to rigorously define in thermodynamic terms the hydration number and the related concept of hydration shell, which are both widely used as operational definitions in experimental studies. Hydration number, nh, is the effective number of water molecules solvating the solute and represents the derivative of the free energy of hydration with respect to the logarithm of water activity. One traditional way of studying hydration relies on the use of volumetric measurements. However, microscopic interpretation of macroscopic volumetric data is complicated and currently relies on empirical models that are not backed by theory. We use our derived model to link the microscopic determinants of the volumetric properties of a solute and its statistical thermodynamic parameters. In this treatment, the partial molar volume, V degrees, of a solute depends on the cavity volume, hydration number, and the properties of waters of hydration. In contrast, the partial molar isothermal compressibility, K degrees T, and expansibility, E degrees, observables, in addition to the intrinsic compressibility and expansibility of the cavity enclosing the solute, hydration number, and the properties of waters of hydration, contain previously unappreciated relaxation terms that originate from pressure- and temperature-induced perturbation of the equilibrium between the solvated solute species. If significant, the relaxation terms may bring about a new level of nonadditivity to compressibility and expansibility group contributions that goes beyond the overlap of the hydration shells of adjacent groups. We apply our theoretical results to numerical analyses of the volume and compressibility responses to changes in the distribution of solvated species of polar compounds.     

Excess densities and equimolar surfaces for spherical cavities in water

For hard spheres with a radius up to 10 A in TIP4P water under ambient conditions, the author studies how the excess number of molecules at the accessible surface depends on the radius of the cavity. Simulation results derived from excess volumes are discussed in terms of radial distribution functions (rdfs), which compare well with extended simple point charge and theoretical rdfs from the literature. The excess number of molecules at the accessible surface inserted in the expression which refers to an arbitrary dividing surface enables one to find the position of the equimolar surface. The surface tension corresponding to this dividing surface was obtained from values of the free energy of cavity formation. For radii in the range of the simulation data, its behavior with curvature is quite different from that usually shown in the literature. A model, which describes how the excess number of molecules at the accessible surface changes with the radius, is discussed in the large length limit by examining consistent rdfs described by a simple analytical form. The inclusion in the model of a logarithmic term has also been considered. Comparison with theoretical results from the literature shows a good agreement for a cavity with a radius of 20 A. For a radius of 100 A and beyond, the model predicts instead sharper density profiles. Such differences have a poor effect on the surface tension at the equimolar surface.     

Accurate and efficient generalized born model based on solvent accessibility: derivation and application for LogP octanol/water prediction and flexible peptide docking

A novel method for fast and accurate evaluation of the generalized Born radii in macromolecular solvation electrostatics calculations is proposed, based on the solvent accessibility of the first two solvation layers around an atom. The reverse generalized Born radii calculated by the method have correlation coefficient of 98.7% and RMSD of 0.031 A(-1) with the values obtained using a precise but significantly slower numerical boundary element solution. The method is applied to derive an estimate of the free solvation energy difference between octanol and water and to predict LogP octanol-water. A nine-parameter model is optimized on an 81 compound training set and applied to predict LogP(ow) for an external evaluation set of 19 drug molecules with RMSD of 0.9. The new GB approximation is also tested in Monte Carlo docking simulations of the fully flexible p53 peptide fragment to MDM2. The best energy solution found in the simulations has RMSD of 2.8 A to the X-ray structure.     

In silico prediction of drug solubility. 3. Free energy of solvation in pure amorphous matter

The solubility of drugs in water is investigated in a series of papers. In this work, we address the process of bringing a drug molecule from the vapor into a pure drug amorphous phase. This step enables us to actually calculate the solubility of amorphous drugs in water. In our general approach, we, on one hand, perform rigorous free energy simulations using a combination of the free energy perturbation and thermodynamic integration methods. On the other hand, we develop an approximate theory containing parameters that are easily accessible from conventional Monte Carlo simulations, thereby reducing the computation time significantly. In the theory for solvation, we assume that DeltaG* = DeltaGcav + ELJ + EC/2, where the free energy of cavity formation, DeltaGcav, in pure drug systems is obtained using a theory for hard-oblate spheroids, and ELJ and EC are the Lennard-Jones and Coulomb interaction energies between the chosen molecule and the others in the fluid. The theoretical predictions for the free energy of solvation in pure amorphous matter are in good agreement with free energy simulation data for 46 different drug molecules. These results together with our previous studies support our theoretical approach. By using our previous data for the free energy of hydration, we compute the total free energy change of bringing a molecule from the amorphous phase into water. We obtain good agreement between the theory and simulations. It should be noted that to obtain accurate results for the total process, high precision data are needed for the individual subprocesses. Finally, for eight different substances, we compare the experimental amorphous and crystalline solubility in water with the results obtained by the proposed theory with reasonable success.     

Charge-dependent cavity radii for an accurate dielectric continuum model of solvation with emphasis on ions: aqueous solutes with oxo, hydroxo, amino, methyl, chloro, bromo, and fluoro functionalities

Dielectric continuum solvation models are widely used because they are a computationally efficacious way to simulate equilibrium properties of solutes. With advances that allow for molecular-shaped cavities, they have reached a high level of accuracy, in particular for neutral solutes. However, benchmark tests show that existing schemes for defining cavities are unable to consistently predict accurately the effects of solvation on ions, especially anions. This work involves the further development of a protocol put forth earlier for defining the cavities of aqueous solutes, with resulting advances that are most striking for anions. Molecular cavities are defined as interlocked spheres around atoms or groups of atoms in the solute, but the sphere radii are determined by simple empirically based expressions involving the effective atomic charges of the solute atoms (derived from molecular electrostatic potential) and base radii. Both of these terms are optimized for the different types of atoms or functional groups in a training set of neutral and charged solutes. Parameters in these expressions for radii were fitted by minimizing residuals between calculated and measured standard free energies of solvation (DeltaG(s)*), weighted by the uncertainty in the measured value. The calculations were performed using density functional theory with the B3LYP functional and the 6-311+G** basis set and the COnductor-like Screening MOdel (COSMO). The optimized radii definitions reproduce DeltaG(s)* of neutral solutes and singly charged ions in the training set to within experimental uncertainty and, more importantly, accurately predict DeltaG(s)* of compounds outside the training set, in particular anions (J. Phys. Chem. A 2003, 107, 5778). Inherent to this approach, the cavity definitions reflect the strength of specific solute-water interactions. We surmise that this feature underlies the success of the model, referred to as the CD-COSMO model for Charge-Dependent (also Camaioni-Dupuis) COSMO model. These findings offer encouragement that we can keep extending this scheme to other functional groups and obtain better accuracy in using continuum solvation models to predict equilibrium properties of aqueous ionic solutes. The approach is illustrated for a number of test cases, including the determination of acidities of an amine base, a study of the tautomerization equilibrium of a zwitterionic molecule (glycine), and calculating solvation energies of transition states toward a full characterization of reaction pathways in aqueous phase, here in S(N)2 exchange reactions. The calculated reaction barriers in aqueous solution are in excellent agreement with experimental values.     

N-(1-piperidinepropionyl)amphotericin B methyl ester (PAME)--a new derivative of the antifungal antibiotic amphotericin B: searching for the mechanism of its reduced toxicity

N-(1-piperidinepropionyl)amphotericin B methyl ester (in short, PAME), a low-toxicity amphotericin B derivative, has been investigated in Langmuir monolayers at the air/water interface alone and in mixtures with cellular membrane sterols (a mammalian sterol, cholesterol, and a fungal sterol, ergosterol) and a model phospholipid (DPPC). The analysis of the strength of interaction between PAME and both sterols as well as DPPC was based, on surface pressure measurements and analysis of the isothermal compressibility (C(s)(-1)), the mean area per molecule (A(12)), the excess free energy of mixing (DeltaG(Exc)) and the total free energy of mixing (DeltaG(M)). It has been found that the interactions between PAME and sterols are attractive; however, their strength is significantly weaker for mixtures of PAME with cholesterol than with ergosterol. This casts light on the improved selectivity of PAME toward fungal cells. The strongest interactions, found for PAME/DPPC mixtures, proved an important role of DPPC in the mechanism of reduced toxicity of PAME as compared to amphotericin B. Due to stable complex formation between PAME and DPPC the antibiotic is immobilized with DPPC molecules, which reduces the concentration of free antibiotic, which is capable of interacting with membrane sterols.     

Continuum polarizable force field within the Poisson-Boltzmann framework

We have developed and tested a complete set of nonbonded parameters for a continuum polarizable force field. Our analysis shows that the new continuum polarizable model is consistent with B3LYP/cc-pVTZ in modeling electronic response upon variation of dielectric environment. Comparison with experiment also shows that the new continuum polarizable model is reasonable, with accuracy similar to that of B3LYP/cc-pVTZ in reproduction of dipole moments of selected organic molecules in the gas phase. We have further tested the validity to interchange the Amber van der Waals parameters between the explicit and continuum polarizable force fields with a series of dimers. It can be found that the continuum polarizable model agrees well with MP2/cc-pVTZ, with deviations in dimer binding energies less than 0.9 kcal/mol in the aqueous dielectric environment. Finally, we have optimized atomic cavity radii with respect to experimental solvation free energies of 177 training molecules. To validate the optimized cavity radii, we have tested these parameters against 176 test molecules. It is found that the optimized Poisson-Boltzmann atomic cavity radii transfer well from the training set to the test set, with an overall root-mean-square deviation of 1.30 kcal/mol, an unsigned average error of 1.07 kcal/mol, and a correlation coefficient of 92% for all 353 molecules in both the training and test sets. Given the development documented here, the next natural step is the construction of a full protein/nucleic acid force field within the new continuum polarization framework.     

The importance of excluded solvent volume effects in computing hydration free energies

Continuum dielectric methods such as the Born equation have been widely used to compute the electrostatic component of the solvation free energy, DeltaG(solv)(elec), because they do not need to include solvent molecules explicitly and are thus far less costly compared to molecular simulations. All of these methods can be derived from Gauss Law of Maxwell's equations, which yields an analytical solution for the solvation free energy, DeltaG(Born), when the solute is spherical. However, in Maxwell's equations, the solvent is assumed to be a structureless continuum, whereas in reality, the near-solute solvent molecules are highly structured unlike far-solute bulk solvent. Since we have recently reformulated Gauss Law of Maxwell's equations to incorporate the near-solute solvent structure by considering excluded solvent volume effects, we have used it in this work to derive an analytical solution for the hydration free energy of an ion. In contrast to continuum solvent models, which assume that the normalized induced solvent electric dipole density P(n) is constant, P(n) mimics that observed from simulations. The analytical formula for the ionic hydration free energy shows that the Born radius, which has been used as an adjustable parameter to fit experimental hydration free energies, is no longer ill defined but is related to the radius and polarizability of the water molecule, the hydration number, and the first peak position of the solute-solvent radial distribution function. The resulting DeltaG(solv)(elec) values are shown to be close to the respective experimental numbers.     

Water coordination structures and the excess free energy of the liquid

We assess the contribution of each coordination state to the hydration free energy of a distinguished water molecule, the solute water. We define a coordination sphere, the inner-shell, and separate the hydration free energy into packing, outer-shell, and local, solute-specific (chemical) contributions. The coordination state is defined by the number of solvent water molecules within the coordination sphere. The packing term accounts for the free energy of creating a solute-free coordination sphere in the liquid. The outer-shell contribution accounts for the interaction of the solute with the fluid outside the coordination sphere and it is accurately described by a Gaussian model of hydration for coordination radii greater than the minimum of the oxygen-oxygen pair-correlation function: theory helps identify the length scale to parse chemical contributions from bulk, nonspecific contributions. The chemical contribution is recast as a sum over coordination states. The nth term in this sum is given by the probability p(n) of observing n water molecules inside the coordination sphere in the absence of the solute water times a factor accounting for the free energy, W(n), of forming an n-water cluster around the solute. The p(n) factors thus reflect the intrinsic properties of the solvent while W(n) accounts for the interaction between the solute and inner-shell solvent ligands. We monitor the chemical contribution to the hydration free energy by progressively adding solvent ligands to the inner-shell and find that four-water molecules are needed to fully account for the chemical term. For a chemically meaningful coordination radius, we find that W(4) ≈ W(1) and thus the interaction contribution is principally accounted for by the free energy for forming a one-water cluster, and intrinsic occupancy factors alone account for over half of the chemical contribution. Our study emphasizes the need to acknowledge the intrinsic solvent properties in interpreting the hydration structure of any solute, with particular care in cases where the solute-solvent interaction strength is similar to that between the solvent molecules.     

An effective model for calculations of the hydrophobic component of the Gibbs energy of solution of small- and medium-sized molecules

A model for calculations of the hydrophobic component of the Gibbs energy G of solution of small- and medium-sized molecules in water was suggested. The model uses the decomposition G = G ^s + G ^ns , where G ^s is the contribution of the cavity volume characterized by the effective radius R _v calculated from the volume of the cavity replaced by a sphere equal in volume and G ^ns is the contribution caused by the difference of the shape of the cavity from spherical. The most substantial contribution to G for small- and medium-sized molecules is G ^s , which can be calculated once for spheres over the radius range under consideration and then used in calculations of G for cavities of arbitrary shapes. At the same time, the Gibbs energy part G ^ns specific with respect to the cavity shape is comparatively small and can be calculated in terms of a simple model with the use of nonspherical indexes, which quantitatively characterize the degree of cavity shape deviation from spherical. The model was verified by Monte Carlo calculations for a set of cavities with different volumes and shapes in an ensemble of water molecules with periodic boundary conditions under usual conditions. The error in G calculated using the model with respect to Monte Carlo results was ∼1 kcal/mol, which was close to the error of the Monte Carlo method itself for the problem under consideration. The time required to calculate G in terms of the model was a few seconds, whereas calculations of G by the Monte Carlo method took from several hours to several days.