Evaluation of Hydration Free Energy by Level-Set Variational Implicit-Solvent Model with Coulomb-Field Approximation.

In this article, we systematically apply a novel implicit-solvent model, the variational implicit-solvent model (VISM) together with the Coulomb-Field Approximation (CFA), to calculate the hydration free energy of a large set of small organic molecules. Because these molecules have been studied in detail by molecular dynamics simulations and other implicit-solvent models, they provide a good benchmark for evaluating the performance of VISM-CFA. With all-atom Amber force field parameters, VISM-CFA is able to reproduce well not only the experimental and MD simulated total hydration free energy but also the polar and nonpolar contributions individually. The correlation between VISM-CFA and experiments is R2 = 0.763 for the total hydration free energy, with a root-mean-square deviation (RMSD) of 1.83 kcal/mol, and the correlation to results from TIP3P explicit water MD simulations is R2 = 0.839 with a RMSD = 1.36 kcal/mol. In addition, we demonstrate that VISM captures dewetting phenomena in the p53/MDM2 complex and hydrophobic characteristics in the system. This work demonstrates that the level-set VISM-CFA can be used to study the energetic behavior of realistic molecular systems with complicated geometries in solvation, protein-ligand binding, protein-protein association, and protein folding processes.

Introduction

The recently developed variational implicit-solvent models (VISM),[1,2] coupled with the level-set numerical method,[3−7] provide a physically self-consistent description of molecular solvation. Here, we apply this method to a large set of organic molecules to understand its accuracy and properties. Central in the VISM is a mean-field free-energy functional of all possible solute–solvent interfaces, or dielectric boundaries, that separate the continuum (or implicit) solvent from solute atoms. Such a free-energy functional consists of surface energy, solute–solvent van der Waals interaction energy, and continuum electrostatic energy, all of which depend solely on a given solute–solvent interface. The surface energy includes local curvature correction characterized by a fitting parameter: the Tolman coefficient. The minimization of the free-energy functional determines the solvation free energy and stable equilibrium solute–solvent interfaces. In our previous work, we developed a level-set method to numerically relax the free-energy functional in the three-dimensional space.[3−6] Compared with molecular dynamics (MD) simulations, our extensive numerical results have demonstrated the initial success of this new approach to molecular solvation in efficiently capturing the hydrophobic interaction, multiple equilibrium states of hydration, and fluctuation between such states.[3,5,7]

Most existing implicit-solvent models,[11−14] particularly those used in MD simulations, are based on various predefined solute–solvent interfaces, such as the van der Waals surface (vdWS), solvent-excluded surface (SES), or solvent-accessible surface (SAS).[15−19] Such surfaces are used to compute the solvation free energy by the sum of the independently calculated surface energy and electrostatic energy. Recently, a Semi-Explicit Assembly (SEA) solvation model has been developed by Fennell et al.[20−23] SEA improves on traditional implicit-solvent models of solvation by accounting for local solute curvature, near-neighbor nonadditivities, and water dipoles. This method is efficient for solvation free energy calculation and describing the first shell water structure. However, the SEA solvation model may have difficulties capturing wetting and drying transitions for hydrophobic pockets, since the precomputation of the water solvation shell is based on simple convex spheres and pairwise additivity of effects is assumed.

In VISM, electrostatic and apolar (i.e., surface and van der Waals) contributions to the total solvation free energy are coupled in the free-energy functional. This coupling makes the free energy estimation self-consistent with physical processes, such as capillary evaporation,[3,5,7] many-body hydrophobic effects,[24] and hydrophobic–hydrophilic coupling effects.[1,2,25] Consequently, stable equilibrium solute–solvent interfaces determined by the level-set VISM can be quite different from vdWS, SES, or SAS, particularly for the description of hydrophobic interactions.[26−29] We believe that the most significant feature of VISM is its variational nature: It is based on the minimization of the free-energy functional that has a complex energy landscape with multiple local minima. An underlying system can fluctuate among these states and exhibit hysteresis[26,30] (i.e., a relaxation-pathway dependent ensemble of equilibrium states). All of these are difficult to capture by fixed-surface implicit-solvent models.

In our recent work, we developed a level-set VISM with the Coulomb-field approximation (CFA) for electrostatic energy.[31−36] The CFA allows us to derive a simple analytical formula for the effective boundary force defined as the negative functional derivative of the total free energy with respect to the location change of the solute–solvent interface. This force is used as the “normal velocity” in our level-set numerical optimization.

The CFA of the electrostatic free energy is in the form of a volume integral occupied by solvent. However, we recognize that like any other sharp-interface model, the final solute–solvent surface of VISM and the dielectric solute–solvent boundary (DSSB) required by CFA do not necessarily coincide for proper description of polar solvation free energies. We discuss the subtleties in the DSSB definition and propose a heuristical fix based on a parallel shift of the VISM surface. This leads to overall satisfying results.

In this work, we apply the level-set VISM with CFA to a large set of molecules with different complexity to further evaluate the theory and method. Specifically, we consider a set of 504 small organic compounds whose experimental values of absolute hydration free energies are commonly used to benchmark the performance and optimize the parameters of various solvation models. These molecules have a wide variety of physical chemical properties commonly encountered in drug like molecules, including saturated and unsaturated hydrocarbons, aromatic and heterocyclic rings, halides, and molecules with polar functional groups. They are categorized into 17 groups.[9,10] Here, we also use these molecules to evaluate the performance of VISM-CFA. We emphasize here that there are only two adjustable parameters in our method: the Tolman coefficient τ and the DSSB shifting parameter ξ. This Tolman coefficient is the coefficient characterizing the curvature effect on solvent surface tension. It is set to be a constant for all molecular systems. We observe very good agreement between our results and the previous MD simulations and experimental results, in terms of both the total solvation energy and individual polar and apolar contributions. We notice that the results are sensitive to the 12–6 Lennard-Jones (LJ) parameters in the all-atom Amber force field. This sensitivity was also observed in MD simulations.[9]

Subsequently, we apply the VISM-CFA to study a more complicated biomolecular system: the p53/MDM2 complex. It has been shown that the tumor suppressor gene p53 is responsible for maintaining the integrity of the human genome and plays a vital role in DNA repairing machinery.[37] Loss of p53 tumor suppressor activity is observed in about 50% of human cancers. Understanding the interactions between p53 and MDM2 has important implications for the design of new anticancer agents. The wild type p53/MDM2 (PDB code 1YCR) complex involves a hydrophobic binding pocket. The highly hydrophobic binding cavity in MDM2 is an excellent example to show VISM’s ability to capture dewetting phenomena.[38]

We notice that several related issues, such as the coupling of the solvent boundary to the overall energy, the curvature effect to surface energy, and dewetting transition, have been discussed in the literature.[39−42] We also notice that some related models and methods have been proposed,[14,29,43−47] providing helpful resources to refine new methods.

The rest of the paper is arranged as follows: In section II, we briefly introduce the VSIM with the CFA, and the level-set method for minimizing numerically the VISM functional. In section III, we apply this VISM-CFA method to calculate the hydration free energies of the 504 organic compounds, and to study a biological macromolecular system in terms of the dewetting phenomena in the hydrophobic pockets. Section IV concludes our studies.

Theory and Methods

We describe here briefly the variational implicit-solvent model (VISM) with the Coulomb-field approximation (CFA). Details are described previously.[1,2,8,48]

VISM with CFA

Let Ω denote the region of a solvation system. It is divided by a solute–solvent interface Γ into the solute region Ωm and the solvent region Ωw. We assume that there are N solute atoms located at X1,...,X inside Ωm and with point charges Q1,...,Q, respectively. We consider the free energy of the solvation system as a functional of all possible solute–solvent interface Γ’s under the CFA.[1,2,8]Here, the first term Pvol(Ωm) is the volumetric part of energy for creating the solute cavity Ωm with P being the pressure difference between the solvent liquid and solute vapor. This term can often be neglected for systems at the nanometer scale, since the pressure difference P under normal conditions is very small.

The second term is the surface energy, where γ(x) is the surface tension given by γ(x) = γ0(1 – 2τH(x)), where γ0 is the constant macroscopic surface tension for a planar solvent liquid–vapor interface, τ is the first order correction coefficient termed here as the Tolman coefficient,[49] and H(x) is the mean curvature defined as the average of the two principal curvatures.

The third term is the energy of the van der Waals interaction between the solute atoms and the continuum solvent. The parameter ρw is the solvent number density. Each U(|x – x|) is the interaction potential density between the solute atom at x and solvent at x in the solvent region Ωw. We define U to be the Lennard-Jones (LJ) potentialwith ε = (εεww)1/2 and σ = (σ + σww)/2 being the potential parameters between solute atoms and water molecules taken directly from explicit water force fields, ε and σ are the depth of the van der Waals potential well and diameter of the ith solute atom, and εww and σww are the depth of the van der Waals potential and diamter of the water molecule. All these parameters are taken from the all atom force field directly without any modification.

CFA and DSSB Definition

The last term in eq 1 is the electrostatic contribution to the solvation free energy. It is defined by the Born cycle[50] as the difference of the energies of two states. The first is a reference state, and the second is the solvated state. A natural reference state is the charged solute molecules in a vacuum. In this case, the electric potential ψ1 is given bywhere ε0 is the vacuum permittivity and εm is the relative permittivity of the solute molecule. The corresponding electric field E1 = E1(x) and electric displacement field D1 = D1(x) are given by E1(X) = −∇ψ1(x) and D1(x) = εmε0E1(x), respectively. The electrostatic energy in this state is given byexcluding the self-interactions. In the second state, the solute molecules are immersed in the solvent, creating a solute–solvent interface Γ. The corresponding electric field E2 = E2(x) and the electric displacement field D2 = D2(x) are related by D2 = εε0E2, where the relative permittivity or dielectric coefficient ε = ε(x) is defined through the VISM solute–solvent interface Γ by ε(X) = εm if x ∈ Ωm and ε(x) = εw if x ∈ Ωw, and εw is the relative permittivity of the solvent. The electrostatic energy of this state is given byagain excluding the self-interactions. Note that G2 depends on Γ through the dielectric coefficient ε = ε(x). Now the electrostatic part of solvation free energy is G2[Γ] – G1. If we apply the CFA D2 ≈ D2 (cf. e.g., ref (51)), this part of the free energy is exactly given by the last term in eq 1.

In reality, the solute–solvent interface is not a sharp boundary.[14] Both the center-of-mass density distribution ρw(r) of water molecules and the local dielectric coefficient εw(r) are smooth functions. While in a sharp-boundary model, such as VISM, they have to be mapped onto an effective 2D-surface. Besides the value of water density, the VISM surface from optimizing eq 1 largely depends on the van der Waals term (cf., eq 2) for convex shapes in a molecule. This is because the short-range van der Waals repulsion is the only term to keep the VISM surface from collapsing in the convex part of a molecule. In the concave part of a molecule, solvent surface energy can also contribute to maintaining surface integrity together with van der Waals interactions. On the other side, it is known that the dielectric solute–solvent boundary (DSSB) may not overlap with the VISM surface in this implicit-solvent model.[2,14] Here, we propose an approximation that the DSSB could be obtained by parallel shifting the VISM surface similar to the size of a solvent molecule. This rests on the assumption that the DSSB is highly correlated with the VISM surface.[2] We found the best fitting value to be 1.4 Å for those organic molecules in water. This 1.4 Å shifting is an heuristic finding that works well for small organic molecules. Ideally, it is desirable to have a rigorous theoretical mapping to obtain the DSSB in a well-defined manner. We are working on the development and reporting our findings in the future. In the current study, the ad hoc shifting is applied for its simplicity and empirical accuracy.

Level-Set Numerics

Now the free energy G[Γ] determines the effective boundary force, −δΓG[Γ], acting on the VISM solute–solvent surface Γ, where δΓ is the variational derivative with respect to the location change of Γ. It is only the normal component of this force that can affect the motion of such a solute–solvent surface. We denote by n = n(x) the unit normal vector at a point x on the solute–solvent surface Γ, pointing from the solute region Ωm to the solvent region Ωw. Then, the normal component of the effective boundary force is given by[3,52,53]where K = K(x) is the Gaussian curvature, defined as the product of the two principal curvatures, at a point x on Γ. This force will be used as the “normal velocity” in our level-set numerical calculations.

To minimize the free-energy functional (eq 1), we begin with an initial surface that encloses all the solute atoms located at x1,...,x. The initial interface may have a very large value of the free energy. We then move the surface in the direction of steepest descent of the free energy by the level-set method until a minimum is reached.

The starting point of the level-set method is the representation of a surface Γ using the (zero) level set of a function ϕ = ϕ(x): Γ = {x: ϕ(x) = 0}.[54−56] The motion of a moving surface Γ = Γ(t) with t denoting the time is then tracked by the evolution of the level-set function ϕ = ϕ(x,t) whose zero level set is Γ(t) at each t. Such evolution is determined by the level-set equationwhere v = v(x,t) is the normal velocity of a point x on the surface at time t. To apply the level-set method to minimize our free-energy functional, we choose the “normal velocity” v to move our surface in the direction of steepest descent of the free energy. This means that the normal velocity v is proportional to the normal component of the effective boundary force, v = MF, where M is the mobility constant which we take to be the unity. Thus, mathmatically v = F, and it is given by eq 5.[3,52,53]

With such a choice of the normal velocity, our level-set method is in fact an optimization method of the approximate steepest descent type. The word “approximate” here is due to the shift of VISM surface for electrostatic energy evaluation. The “time” here is the optimization step. The VISM free-energy functional has a complex energy landscape due to the capillary evaporation or “dewetting” energy barriers existing in an underlying molecular system. Different initial surfaces can lead to different local minima that are physically meaningful. In order to capture multiple local minima, we design three types of initial solute–solvent interfaces. The first one is a tight wrap: a surface that is close to the van der Waals surface of the atoms. The second one is a loose wrap: a surface that loosely encloses all the solute atoms. An example of such a loose wrap is a sphere of large radius. The third one is a combination of tight and loose wraps. We refer readers to the previous work[3−6,8] for all of the details of numerical solutions of the level-set eq 6.

Force Fields, Molecular Configurations, and Parameters

The coordinates and partial charges for all compounds in the database are obtained directly from the Supporting Information provided by Mobley et al.,[9] and the 12–6 LJ potential parameters of solute atoms were obtained from the Amber all-atom force field.[57] Extensive replica exchange molecular dynamics simulation (REMD) and Monte Carlo simulations (MC) of n-alkanes by Ferguson et al.[58] indicate that the dominating conformations for n-alkanes up to n-eicosane (C20) in length are extended in both the gas and solvated phases. The longest n-alkane in our data set is n-decane (C10). Hence, the solvation free energy obtained from the optimized gas phase conformation can be compared with ensemble average results from MD simulations or experiments. We also carry out conformational sampling and ensemble averaging for the n-decane and confirmed that the ensemble averaged solvation energy is very close to the one obtained from a single optimized gas phase conformation.

To compare our calculations with the results of TIP3P water MD simulations, we use the macroscopic planar surface tension γ0 = 0.076 kcal/mol·Å2 at 300 K obtained from the TIP3P water simulation.[59] The Tolman coefficient τ is chosen to be 1 Å.[6,8] The well depth between water molecules is εww = 0.152 kcal/mol, and the solvent molecular diameter is σww = 3.15 Å.[60] The number density of water is ρw = 0.0333/Å3, and the dielectric constant in a vacuum εm and TIP3P water εw are 1 and 92, respectively.[61]

Results and Discussion

504 Small Organic Molecules

In Figure 1a,b, we compare the total ΔGsolv from VISM with the experimental solvation free energies and the explicit TIP3P water MD simulation results, respectively. The correlation coefficient (R2) to experimental data is 0.763 with RMSD = 1.83 kcal/mol. This is comparable to the best implicit-solvent models where the correlation coefficients R2 range from 0.66 to 0.81.[10] The correlation coefficient is 0.839 with RMSD = 1.36 kcal/mol in comparison with the explicit TIP3P water MD simulations results. This is also comparable to other implicit-solvent models with R2 ranging from 0.794 to 0.911.[10] Note that other implicit-solvent models use a larger set of atomic radii as fitting parameters, while VISM-CFA directly uses the explicit all-atom force-field parameters without any modifications. The only adjustable parameters introduced in VISM-CFA are a constant Tolman coefficient for all molecules and a dielectric boundary shifting parameter related to solvent molecular size. On the basis of previous studies,[5,6] the optimal Tolman length is around 1 Å, and the shifting parameter can be taken from the size of a solvent molecule. This means that there could be effectively no fitting parameters in VISM formulation beyond all-atom force fields. For cases with extremely strong electrostatic potential, such as single ions, the shifting parameter might be adjusted based on the electrostatic potential.[2]

Figure 1

The correlations between the solvation free energies calculated by VISM-CFA and (a) the experimental total solvation free energy and (b) the ones calculated by explicit TIP3P water MD simulations.

In addition to the total hydration free energy, it is important to ensure that the individual components of the solvation free energy are captured correctly. In level-set VISM-CFA calculations, the range of nonpolar energies for the entire set of small molecules is less than 4 kcal/mol, much smaller than the corresponding polar solvation free energies in those molecules. Therefore, the correlation of the total solvation free energies is largely determined by the polar contribution. First, we investigate the correlation of electrostatic contributions between VISM-CFA and the explicit TIP3P water MD simulation values. In Figure 2, the polar components of explicit solvent MD simulations are compared to those obtained by the VISM-CFA. The corresponding correlation coefficient R2 is 0.899 with a RMSD = 0.92 kcal/mol. The correlation coefficient for polar interactions between TIP3P water MD simulations and other implicit-solvent models ranges from 0.79 to 0.93.[10] From these comparisons, it is clear that the VISM-CFA is a good approximation for the polar contributions to total solvation free energies. Currently, we are developing more accurate methods such as the Poisson–Boltzmann (PB) equation and Yukawa-Field Approximations (YFA)[52,53] to describe electrostatic interactions. We will report our related results in separate publications.

Figure 2

The correlation of electrostatic contributions to solvation free energy between the explicit solvent (TIP3P) MD simulation results and the VISM-CFA calculations.

In Figure 3, we compare the nonpolar solvation free energy between VISM-CFA and MD simulation results. The correlation coefficient is 0.415. The nonpolar part of VISM solvation free energies includes contribution from the van der Waals interactions and the surface energy. Both of them range from 2 to 22 kcal/mol. The net nonpolar energies ranging from −1 to 4 kcal/mol are the results of the small differences between these two relatively large values. This is generally a challenging issue that requires the physical model to be very accurate for both terms. Any bias toward one or the other would result in systematic errors. This issue is largely avoided in most other implicit-solvent models, where the nonpolar term is collectively represented by single surface energy with a fitting surface tension. In those models, the surface tension is drastically smaller than and has no physical connection to realistic surface tension obtained from MD simulations or experiments. In VISM, the surface tension is taken directly from MD simulations or experiments. For the nonpolar solvation free energies in VISM, it is very sensitive to the vdW potential parameters. In the MD simulations by Mobley et al., they also found that the systematic error could arise from wrong carbon well-depths.[9] In our case, the vdW interaction energy and surface energy are highly correlated with each other as shown in Figure 3. In Figure 3b, the correlation coefficient (R2) between VISM vdW solvation free energy and surface component is 0.835. The surface energy favors a small solute–solvent interface while the vdW interactions prefer more contact areas; the delicate balance between these distinct components is the driving force for the wet and dry transitions.[5]

Figure 3

The correlations between (a) VISM nonpolar solvation free energy and MD nonpolar solvation free energy and (b) the two nonpolar components of VISM-CFA: vdW energy and surface energy.

The total nonpolar contribution to the solvation free energy is typically assumed to be highly correlated with surface area in implicit-solvent models. The explicit TIP3P solvent simulations reveal that there is essentially no correlation between total nonpolar solvation free energy and the surface area even though the attractive and repulsive part correlated with surface area separately.[9] Through detailed analysis, our study reveals the underlying physics for these observations. Plots of VISM-CFA nonpolar components versus surface area are shown in Figure 4a. The correlation of the net nonpolar solvation free energies with surface area is R2 = 0.102 for the entire set. The net nonpolar solvation free energies in VISM-CFA are the summation of the surface energies and the vdW interactions between solute and solvent. In Figure 4b, it shows that the surface energy correlates with surface area R2 = 0.999, which is the result of the second term of eq 1. In VISM-CFA, since different components of solvation free energy are coupled with each other, the electrostatic contribution has strong influences on the nonpolar parts. Strong electrostatic interaction drives the solute–solvent surface closer to solute atoms until the vdW repulsion can counterbalance that. Under these circumstances, the change of vdW solvation energies depends less on the surface area due to the steepness of vdW potential. Therefore, nonpolar energy is not tightly related to the surface area. However, if the electrostatic part is relatively small, like in the cases of 25 alkanes (Figure 5), the surface energy strongly influences the nonpolar contribution, which makes the nonpolar contribution highly correlated with the surface area. Figure 5a shows that the correlation coefficient (R2) between nonpolar solvation free energies and the surface area is 0.878 for the alkanes. In this study, the VISM surface area is similar to the solvent accessible surface area (SAS), and they are highly correlated with each other with R2 = 0.922, as shown in Figure 4d.

Figure 4

Correlations between the solute–solvent surface area by VISM-CFA and solvation free energies of (a) the total nonpolar part, (b) the surface tension part, (c) the vdW part, and (d) the solvent accessible surface area calculated by Macromodel (GB) for the entire set.

Figure 5

The correlations between the solute–solvent surface area by VISM-CFA and solvation free energies of (a) the total nonpolar part, (b) the surface tension part, (c) the vdW part for alkane only set, and (d) the correlations between VISM nonpolar solvation free energy and MD nonpolar solvation free energy.

P53/MDM2 Complex

The results of the small molecule set demonstrate the accuracy of the level-set VISM-CFA method. It leads to a very good correlation with the experimental solvation free energies. Compared with conventional implicit-solvent models, the solute–solvent interface in VISM can adjust itself based on the interplay between polar and nonpolar interactions to reach the free energy minimum. As a result, it is capable of capturing the dewetting phenomena in protein association and dissociation processes. Conceptually, the hydrophobic pockets, critical for ligand binding[62] in biological systems, can be better described in VISM than traditional implicit-solvent models. These properties of VISM are not highly pronounced in small molecules as they are largely convex. The optimized VISM surfaces are very similar to SAS. However, in large biomolecules, complex geometry with deep pockets can significantly benefit from the coupling. In this section, we apply VISM-CFA to the p53-MDM2 protein complex.

To study the solute and solvent interface during the process of p53-MDM2 dissociation, we set up a series of configurations where the two components of p53-MDM2 are increasingly separated from d = 0 Å to d = 30 Å along the axis connecting their geometrical centers. This domain separation d is chosen to be the reaction coordinate. Note that d = 0 Å corresponds to their native configuration in the crystal structures (PDB code: 1YCR). In this study, we consider the protein complex with and without partial charges.

In Figure 6, we depict VISM surface solute–solvent interfaces at three different interdomain separations 10, 14, and 18 Å. It is clear that solvent is significantly excluded from the interdomain region even at d = 14 Å. The VISM optimized interfaces with atomic partial charges are tighter than those without them. This is due to the attractive nature of the electrostatic interactions between the solute and solvent. Indeed, it is clear from eq 5 (the last term) that the electrostatic force always points from the high dielectric solvent region into the low dielectric solute region. In Figure 7, we compare the molecular surface and the equilibrium VISM surface at d = 14 Å. In order to show the differences between the VISM surface and molecular surface, we slice the surface across the binding pocket. The green color represents the molecular interior enclosed by the traditional molecular surface, and the red color represents the enclosure of VISM. As expected, two surfaces are largely similar except for the hydrophobic region, especially inside the binding pocket and p53/MDM2 interdomain. For other parts, the VISM surface resembles the slightly expanded molecular surface. This is because it resembles the center of mass of the first layer water solvent while the molecular surface is generated by the contacting points of the probes. The expansion is about 1.4 Å. In contrast to a fixed-surface implicit-solvent model, this unique character of VISM, i.e., coupling the different contributions, enables it to obtain physically accurate descriptions of the hydration process.

Figure 6

The stable equilibrium solute–solvent interfaces of p53-MDM2 at three different domain separations obtained by the level-set VISM-CFA with loose initial surfaces. The top row is calculated without partial charges, and the bottom row is calculated with partial charges. The colors of the surface represent the values of curvature (blue for large positive and red for large negative curvature).

Figure 7

Slices of the VISM surface and the molecular surface across the hydrophobic packet at d = 14 Å. The green color represents the molecular surface enclosure, and the red color represents the VISM surface enclosure with partial charges.

The p53-MDM2 interface is hydrophobic (70% of the atoms at the interface are apolar).[63] Three hydrophobic residues (i.e., Phe19, Trp23, and Leu26) in p53 are located on the side facing the binding pocket of MDM2. Meanwhile, the MDM2 binding pocket contains 11 hydrophobic residues (Leu54, Leu57, Ile61, Met62, Tyr67, Val75, Val93, Phe86, Ile99, Phe91, and Ile103). The interaction between p53 and MDM2 is essentially hydrophobic.[64−66] This principle has been used to design inhibitors that block p53 and MDM2 binding,[63] and some of them are entering clinical trials.[67] Our ongoing MD simulations of this complex also show dewetting phenomena for this system. The detailed comparison and analysis will be presented in a forthcoming publication as it is outside the scope of this paper.

It should be noted that some other important factors, such as the hydrogen bonding between the protein and solvent, are not considered in the current model. In reality, it can add additional hydrophilicity to the energy functional. Another importance factor is that a uniform Tolman coefficient τ = 1 Å used in the current model could overestimate the energy penalty for concave interface curvature (H < 0), which favors dewetting. Higher order correction terms in the curvature expansion of the surface tension should be considered for the asymmetric reality between concavity and convexity on the small scale.[7] Without these factors, the dewetting phenomena in protein are overestimated quantitatively. Nonetheless, this still offers a sound physical picture of small scale dewetting in proteins.

Figure 8 shows the differences of the solvation free energy ΔΔGsol of tight and loose initial surfaces and their individual components with and without partial charges. At d = 12 Å (the nearest atom distance between p53 and MDM2 is 4.94 Å), the level-set VISM indicates that the solvent is completely excluded from the interdomain region for the tight initial condition with or without partial charges. When the calculation starts from the loose initial condition, the water molecules are completely excluded from the interdomain region until d = 18 Å (the nearest atom distance between two components is 10.89 Å) with and without partial charges. The comparison between those optimized solvation free energies from loose and tight initial conditions reflects different equilibrium states. They are stable solutions of the free energy functional (eq 1). Generally, for solutes with significant apolar elements, VISM often reveals “wet and dry branches” in the solution of the equations. They are typically revealed by solving the free-energy functional with tight or loose initial surfaces, respectively.[7] The different branches appear to correspond to the wetting and drying fluctuations seen in the interface water in explicit solvent simulations, and during the actual binding process when a transition from the wet to the dry situation is expected. If thermal fluctuations are included in the VISM, the various energy branches could be sampled in a Boltzmann-weighted fashion.[7] The calculated VISM-CFA absolute hydration energy differences between the bound and unbound state is 237.29 kcal/mol. The explicit water free energy perturbation (FEP) calculation gives 306.7 ± 3.58 kcal/mol. They are in qualitative agreement considering that FEP for this system of 1688 atoms also has large errors.

Figure 8

The differences of solvation free energy ΔΔGsol of tight and loose initial surfaces and their relative components for both with and without partial charges.

Conclusions

In this work, we apply the recently developed level-set variational implicit-solvent model (VISM) with the Coulomb-field approximation (CFA)[8] to a large set of organic molecules. The calculated hydration free energies are compared with experimental and molecular dynamics (MD) simulation results. We find that our methods can provide very accurate hydration free energy estimation and are superior to other fixed surface implicit-solvent models in terms of the physical description and the number of fitting parameters. Moreover, the level-set numerical method is found to be robust in treating various kinds of molecules with complex geometry and charge distributions.

Through our studies we find that the free energy of a given molecular system is quite sensitive to the LJ parameters used in the vdW solute–solvent interaction. This is also observed in all atom force field MD simulations. Because of the steepness of the repulsive potential, small changes in those parameters can lead to large energy errors. Therefore, the accuracy of force field parameters has to be carefully evaluated before applying them in VISM calculations.

In reality, the solute–solvent interface is not a sharp interface regarding the density distribution of the water center-of-mass and the local dielectric coefficient. In our VISM, since sharp interfaces are used in this model, we introduce a shifting parameter ξ to uniquely relate the VISM surface (center-of-mass of water) to the DSSB surface. The final electrostatic contribution to the total hydration free energy is calculated using this DSSB. We are actively studying the implications of the shifting and ways to unify the VISM surface with DSSB. We will report our findings in a future publication.

In the VISM-CFA calculation, the solute–solvent interfaces are determined self-consistently through a functional variation. This is undoubtedly more expensive than fixed-surface models. On average, each calculation is about 2 to 3 orders of magnitude slower than traditional implicit-solvent models. On the other hand, VISM is several orders of magnitude more efficient than explicit MD simulations for capturing the wet–dry transition phenomena in biomolecular systems. The computational efficiency also depends on initial conditions, the grid resolution, error tolerance, etc. As we see from the comparison between the molecular surface and VISM surface, most parts of the surfaces are similar. The main differences occur primarily at the binding site where accurate description is critical. This is because the shape and hydrophobic properties of the binding sites are usually very different from the rest of the protein. We believe that a trade-off of speed for more accurate description of the binding site hydration is worthy. In practice, the VISM calculation can be significantly accelerated by restricting the calculation to the site of interests with the conventional molecular surface everywhere else. We would like to emphasize again that VISM is not meant to compete with traditional implicit-solvent models for speed. Instead, it is an excellent method for capturing the important wet and dry transitions in complex molecular systems, like the p53/MDM2 system. On the other hand, VISM is much more efficient than explicit water MD simulation which is the only robust method known to capture this behavior to date.

As the CFA is an approximation to the more accurate Poisson–Boltzmann (PB) continuum electrostatics, we are currently implementing PB in VISM. The results will be reported in a forthcoming publication. While PB provides a more accurate electrostatic description of the system, it adds an additional layer of numerical expense. As we have shown here, the CFA can capture the polar energetics with reasonable accuracy for both small and large molecular systems. Therefore, it can be combined with VISM-PB to accelerate the computations by only employing PB infrequently during the optimization.