Level-Set Variational Implicit-Solvent Modeling of Biomolecules with the Coulomb-Field Approximation.

Central in the variational implicit-solvent model (VISM) [Dzubiella, Swanson, and McCammon Phys. Rev. Lett.2006, 96, 087802 and J. Chem. Phys.2006, 124, 084905] of molecular solvation is a mean-field free-energy functional of all possible solute-solvent interfaces or dielectric boundaries. Such a functional can be minimized numerically by a level-set method to determine stable equilibrium conformations and solvation free energies. Applications to nonpolar systems have shown that the level-set VISM is efficient and leads to qualitatively and often quantitatively correct results. In particular, it is capable of capturing capillary evaporation in hydrophobic confinement and corresponding multiple equilibrium states as found in molecular dynamics (MD) simulations. In this work, we introduce into the VISM the Coulomb-field approximation of the electrostatic free energy. Such an approximation is a volume integral over an arbitrary shaped solvent region, requiring no solutions to any partial differential equations. With this approximation, we obtain the effective boundary force and use it as the "normal velocity" in the level-set relaxation. We test the new approach by calculating solvation free energies and potentials of mean force for small and large molecules, including the two-domain protein BphC. Our results reveal the importance of coupling polar and nonpolar interactions in the underlying molecular systems. In particular, dehydration near the domain interface of BphC subunits is found to be highly sensitive to local electrostatic potentials as seen in previous MD simulations. This is a first step toward capturing the complex protein dehydration process by an implicit-solvent approach.

Introduction

Recent years have seen the development of a new class of implicit-solvent models—the variational implicit-solvent model (VISM).[1,2] Coupled with the robust level-set numerical method, such models allow an efficient and quantitative description of molecular solvation.[3−7] Central in the VISM is a mean-field free-energy functional of all possible solute–solvent interfaces, or dielectric boundaries, that separate the continuum solvent from all solute atoms. In a simple setting, such a free-energy functional consists of surface energy, solute–solvent van der Waals interaction energy, and continuum electrostatic free energy, all depending solely on a given solute–solvent interface. The minimization of the functional determines the solvation free energy and stable equilibrium solute–solvent interfaces. Mathematically, such minimization leads to highly nonlinear geometrical partial differential equations that are hard to solve in general. In our previous work, we developed a level-set method to numerically relax the free-energy functional in the three-dimensional space.[3−5,7] Our extensive numerical results with comparison with molecular dynamics (MD) simulations demonstrated the success of this new approach to the solvation of nonpolar molecular systems in capturing the hydrophobic interaction, multiple equilibrium states of hydration, and fluctuation between such states.[3,5,6]

Most of the existing, surface type, implicit-solvent models[8−11] are based on various kinds of predefined solute–solvent interfaces, such as the van der Waals surface (vdWS), solvent-excluded surface (SES), or solvent-accessible surface (SAS).[12−16] Such a surface is used to compute the solvation free energy as the sum of the surface energy and electrostatic energy, with the two contributions being decoupled. In contrast, those contributions are coupled in the free-energy functional in the VISM, making free-energy estimation more consistent with physical processes, such as capillary evaporation,[3,5,6] many-body hydrophobic effects,[17] and hydrophobic–hydrophilic coupling effects.[1,2,18] Consequently, stable equilibrium solute–solvent interfaces determined by our level-set VISM can be quite different from vdWS, SES, or SAS, particularly when it comes to the description of hydrophobic interactions.[19−22] Perhaps the most significant feature of VISM is that its free-energy functional exhibits a complex energy landscape with multiple local minima corresponding to different equilibrium states. An underlying solvation system can fluctuate among these states and exhibit hysteresis (i.e., a relaxation-pathway-dependent ensemble of equilibrium states). All of these are difficult to capture with a fixed-surface type implicit-solvent model.

We notice that several related issues, such as coupling the solvent boundary to the optimization of overall energy, the curvature effect to surface energy, and dewetting transition, have been discussed in the literature.[23−26] Other related models and methods have also been proposed.[11,22,27−30]

Up to now, the development of level-set VISM has focused only on nonpolar systems. In this work, we propose and implement an efficient approach to treat the electrostatics in the framework of VISM. This approach is based on the Coulomb-field approximation (CFA) of the electric displacement field, without solving the Poisson or Poisson–Boltzmann equation as often done in an implicit-solvent model.[31−36] With such an approximation, we can express the electrostatic part of the solvation free energy as a volume integral over the solvent region that can be arbitrarily shaped. This is similar to the generalized Born model.[37,38] But, we do not compute generalized Born radii and hence introduce no additional parameters. Moreover, the CFA of electrostatics allows us to derive a simple analytical formula of the effective boundary force, defined as the negative functional derivative of the total free energy with respect to the location change of the dielectric boundary. This force is used as the “normal velocity” in our level-set numerical method. We emphasize that all of our level-set calculations are done in the three-dimensional space for arbitrarily shaped solute–solvent interfaces.

We apply our theory and method to several molecular systems of different complexity. The first is a one-particle system. We consider a one-atom system that is gradually charged. This simple radially symmetric system is used to test the convergence and accuracy of our numerical method. Such a system is also used to probe the effect of charge on the equilibrium solute–solvent interface and the minimum free energy. In addition, we consider the hydration of some single ions and compare our VISM calculations with experiments. The second system consists of two methane-like particles that are treated as atoms. For this system, we test the effect of charging on the potential of mean force along the center-to-center distance between the two atoms. The third system consists of two initially hydrophobic plates, where the plate charge is a free parameter and gradually increased. We study the influence of charging on the hydrophobic interactions, the system balance between dry and wet equilibrium states, and the resulting hysteresis. The last system is the protein BphC. In contrast to simple model systems, the complexity of such a protein in terms of size, geometry, and charge distribution, rigorously examines the applicability of our level-set VISM on realistic biological systems. We focus on the dehydration for this particular system and discuss our results in light of MD simulations.[39]

The rest of the paper is organized as follows. In section , we review the VISM functional, present the CFA of electrostatic free energy, and provide formulas of potentials of mean force. In section , we describe briefly our level-set numerical method. In section , we apply our level-set VISM to various model systems and the protein BphC. Finally, in section , we draw conclusions.

Theory

Free-Energy Functional

We consider the solvation of molecules (solutes) in a solvent (typically water). We divide geometrically the underlying system into three parts: the solute region Ωm, the solvent region Ωw, and the corresponding solute–solvent interface Γ, which is the dielectric boundary, cf. Figure 1. We assume that there are N atoms in the solute molecules that are located at x1, ..., x inside Ωm and carry point charges Q1, ..., Q, respectively.

Figure 1

The geometry of a solvation system with an implicit solvent. The solute region, solvent region, and solute–solvent interface are denoted by Ωm, Ωw, and Γ, respectively.

In VISM, one minimizes the solvation free-energy functional, proposed in refs (1) and (2), of all possible solute–solvent interfaces Γ:Here, the term P vol(Ωm), proportional to the volume of solute region Ωm, is the energy of creating the solute cavity, with P being the pressure difference between the solvent liquid and solute vapor. This term can often be neglected for systems on a nanometer scale, since the pressure difference P under normal conditions is very small. The surface integral term is the surface energy, where γ is the surface tension. It is known for systems of nanometer scale that the surface tension γ is no longer a constant. Corrections with a curvature effect are often needed. Here, we apply the correction:[40]where γ0 is the constant macroscopic surface tension for a planar solvent liquid–vapor interface, τ is the correction coefficient historically called the Tolman length,[40] and H is the mean curvature defined as the average of the two principal curvatures. We denote by Ggeom[Γ] the sum of the first two terms in eq II.1 and call it the geometrical part of the free energy.

For each i (1 ≤ i ≤ N), the volume integral in eq II.1 over the solvent region Ωw is the van der Waals (vdW) type interaction energy between the solute atom at x and the solvent molecules or ions at x that are coarse grained. The parameter ρw is the constant solvent density, and each U is a pairwise interaction potential. As typically used in MD simulations, we define U to be the Lennard-Jones (LJ) potentialThe parameters ε of energy and σ of length can vary with different solute atoms as in conventional force fields. We shall denote by GvdW[Γ] the third term (i.e., the summation term) in eq II.1 and call it the van der Waals (vdW) part of the free energy.

The last term Gelec[Γ] in eq II.1 is the electrostatic part of the solvation free energy. It is discussed in detail in the next subsection.

The Coulomb-Field Approximation

The electrostatic part of the solvation free energy Gelec[Γ] is defined by the Born cycle[41] 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 that our charged solute molecules are placed in a vacuum. In this case, the electric potential ψ1 in SI units is given bywhere ε0 is the vacuum permittivity and εm is the relative permittivity of the 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 bywhere the integral is taken over the entire space excluding small balls centered at singularities x with cutoff radii.

In the second state, the solute molecules are immersed in the solvent, creating a solute–solvent interface or dielectric boundary Γ. The electric potential ψ2 = ψ2(x) is the unique solution of the boundary-value problem of the Poisson equationwhere the relative permittivity or dielectric coefficient ε = ε(x) is defined through the dielectric boundary Γ byand εw is the relative permittivity of the solvent, cf. Figure 1. Since the dielectric boundary Γ can be arbitrarily shaped, there is in general no analytical formula for the potential ψ2(x). The corresponding electric field E2 = E2(x) and the electric displacement field D2 = D2(x) are defined by E2 = −▽ψ2 and D2 = εε0E2, respectively. The electrostatic energy of this state is given bywhere again the integral is taken over the entire space excluding small balls centered at singularities x with cutoff radii. Note that the dependence of G2 on Γ is through the dielectric coefficient ε = ε(x).

Now the electrostatic part of solvation free energy is Gelec[Γ] = G2[Γ] – G1. Invoking the Coulomb-field approximation (CFA) D2 ≈ D1 (cf. e.g., ref (38)), we obtain from eqs II.3 and II.6 thatHere, we assume that the respective cutoff small balls in the integrals in G1 and G2[Γ] are the same. As a result, the cutoff radii do not appear in the expression of Gelec[Γ].

Notice that for a single-atom, spherical solute, the solute region is the ball of some radius R > 0 centered at x1. Set x1 = 0 and write Q = Q1. Then the boundary-value problem of eqs II.4 and II.5 has a unique solution:Direct calculations verify that in this case D2(x) = D1(x) = Qx/(4π|x|3) for all x. Hence, the CFA is exact in this case.

In summary, our VISM free-energy functional with the CFA of electrostatics is given bywhere each interaction potential U is given by the LJ potential II.2.

The Effective Boundary Force

The solvation free-energy functional G[Γ] defined on a solute–solvent interface or dielectric boundary Γ gives rise to an effective force acting on the boundary Γ. We define this force to be −δΓG[Γ], the negative functional derivative of the free-energy functional G[Γ] with respect to the location change of the boundary Γ. It is only the normal component of this force that can affect the motion of such a boundary. We denote by n = n(x) the unit normal vector at a point x on the boundary Γ, pointing from the solute region Ωm to the solvent region Ωw, cf. Figure 1. Using the concept of shape derivatives, we can obtain the normal component of this boundary force as a function defined on the boundary Γ. This function is[3,42,43]where K = K(x) is the Gaussian curvature, defined as the product of the two principal curvaures, at a point x on Γ. This force will be used as the “normal velocity” in our level-set numerical calculations.

Potentials of Mean Force

One is often interested in the effective interaction between solutes that stems from direct solute–solute interactions and that is mediated by the solvent. The potential of mean force (PMF) is a general term for such effective interactions. It is usually defined with respect to a reaction coordinate as the difference between the free energy of solvated state at a given coordinate d and that at a fixed, reference coordinate dref. The choice of reaction coordinate is system dependent. For our applications, we consider a solvation system in which solute atoms are divided into two groups: x1, ..., x and x, ..., x. The relative positions of all atoms in the same group are fixed. We chose the reaction coordinate d to be the distance between the geometrical centers (∑x)/M and (∑x)/(N – M) of the two groups of solute atoms. For instance, d is the center-to-center distance between the two atoms in a two-atom system, and it is the plate-to-plate separation distance for a system of two parallel plates. For a two-domain protein, this is the domain–domain separation distance. The reference state is conveniently chosen with dref = ∞; i.e., the atoms in the first group are at infinite separation from those in the second group.

Fix now a finite coordinate d. Denote by Γ a corresponding VISM optimal surface, i.e., a stable equilibrium solute–solvent interface minimizing locally the VISM solvation free-energy functional. We define the (total) PMF as the sum of its separate contributionseach of them given byHere a quantity at ∞ is understand as the limit of that quantity at a coordinate d′ as d′ → ∞.

The term GvdW(d) is the sum of van der Waals interaction energies between all the solute atoms. The above definition of GvdWpmf (d) and some simple calculations lead towhere U is the LJ interaction potential between x and x. The last term in eq II.9 is the direct solute–solute LJ interaction. The electrostatic part Gelecpmf of the PMF is defined via the electrostatic free energies of the solvated states at d and ∞. Denote by G1(d) the vacuum electrostatic energy as defined in eq II.3. Since Gelec[Γ] = G2[Γ] – G1(d), we obtain by the above definition of Gelecpmf(d) that

As d becomes large, the VISM optimal solute–solvent interface Γ becomes the union of two separate VISM optimal solute–solvent interfaces ΓI and ΓII, both independent of d. They are obtained by minimizing the VISM free-energy functional for the corresponding groups of fixed, solute atoms. If we denote by G[ΓI] and G[ΓII] the corresponding minimum VISM free energies for these individual groups of atoms, thenSimilarly, each component of the VISM free energy is the sum of that for the two groups of solute atoms, i.e., G in the above equation can be replaced by Ggeom, or GvdW, or Gelec

For small d, the solute–solute vdW interaction, defined by the double summation term in eq II.9, can be very large and can therefore dominate over all other parts in the total PMF. In order to better understand the solvent influence in the PMF for small d, it is reasonable to look only at the PMF that excludes the solute–solute vdW interaction.

We remark that for a given reaction coordinate d there can be multiple stable equilibrium interfaces Γ that are local minimizers of the VISM free-energy functional. Different local minimizers Γ for the same coordinate d define multiple values G[Γ] of VISM local minimum free energies. Therefore, the PMF can have multiple branches along the reaction coordinate d and hence can lead to hysteresis. Our current level-set VISM has not yet included fluctuations that in principle should allow an underlying system to get out of such a local minimizer. Strictly speaking, therefore, our PMFs are different from those defined using a Boltzmann average over all possible minimizers. Rather, our PMFs reflect possible branches of the VISM free energy along the reaction coordinate d.

Numerical Methods

The Level-Set Optimization Method

To numerically minimize the free-energy functional II.7, we begin with an initial surface that encloses all of 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 steady state 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}.[44−46] With this representation of the surface, the unit normal n = n(x), the mean curvature H = H(x), and the Gaussian curvature K = K(x) of a point x at the surface are then given by n = ▽ϕ/|▽ϕ|, H = (1/2)▽·n, and K = n·adj(▽2ϕ)n, respectively, where ▽2ϕ is the 3 × 3 Hessian matrix of the function ϕ whose entries are all the second order partial derivatives ∂2ϕ of the level-set function ϕ, and adj(▽2ϕ) is the adjoint matrix of the Hessian ▽2ϕ. 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 solve the level-set eq III.1 numerically, one needs to extend the normal velocity v to the entire computational box or a band surrounding the surface Γ(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 steepest descent of the free energy. This means that the normal velocity v is the normal component of the effective dielectric boundary force, v = F, and is given by eq II.8.

With such a choice of the normal velocity, our level-set method is in fact an optimization method of the steepest descent type. The “time” here is the optimization step. As the VISM free-energy functional is quite nonconvex due to the capillary evaporation or “dewetting” energy barriers existing in an underlying molecular system,[21] different initial surfaces can then lead to different local minimizers that are of practical interest. In order to capture multiple local minimizers, 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.

Discretization

We now briefly describe how we solve numerically the level-set eq III.1 with our normal velocity v = F given in eq II.8. More details can be found in our previous publications.[3−5,7] To discretize the spatial variable in the level-set eq III.1 with the normal velocity F given in II.8, we rewrite this equation aswith A = A(x) (dropping the time dependence in ϕ) and B = B(x) given byWe discretize A = A(ϕ) by the central differencing with parameter correction. This means that we linearize A = A(ϕ) at the ϕ that is computed in the previous time step and change the parameter τ so that the linearized part is elliptic; that is, the time-dependent equation is parabolic. We discretize B|▽ϕ| using an upwinding scheme. In our implementation, we use a fifth-order WENO (weighted essential-no-oscillation) scheme.

We use the forward Euler method to discretize the time derivative in the level-set eq III.1:where ϕ((x) and v((x) are the approximations of v(x,t) and ϕ(x,t), respectively, at time t = kΔt (k = 1, 2, ...) and Δt is the time step. We chooseso that the CFL condition is satisfied, where C = C(ϕ) is the matrix defined by A(ϕ) = γ0C(ϕ):▽2ϕ and

Input all the parameters P, γ0, τ, ρw, ε, σ, x, and Q for all i = 1, ..., N; ε0; εm; and εw. Discretize uniformly a computational box containing x1, ..., x with the grid size h a third or half of 1 Å. Compute B(x) at all the grid points x that are not any points x1, ..., x. Generate an initial surface by defining a corresponding level-set function. We choose a level-set function to be negative inside and positive outside. Set k = 1 and start the time iteration.

Choose a narrow band that is centered at the surface and that has a width of about 12 to 16 grid points. At each grid point in the band, compute the gradient ▽ϕ, the Hessian ▽2ϕ, and the curvatures H and K using central differencing schemes. We employ for efficiency the local level-set method developed in ref (47) in which the zero boundary condition is used for the level-set function ϕ at the boundary of the band.

Compute the free energy II.7.

Calculate and extend the normal velocity F. The extension of the normal velocity is necessary, since the LJ potential changes rapidly near the surface, causing possibly numerical instabilities. In practice, we need only to extend the B part of the normal velocity.

Calculate Δt using III.3 and update the level-set function using the Euler scheme III.2.

Reinitialize the level-set function ϕ. To do so, we solve the equationwith the initial value ϕ = ϕ0 at t = 0 to obtain a steady-state solution. Here, ϕ0 is the level-set function before reinitialization, and the time t is different from that in the original level-set equation. The quantity sign(ϕ0) is the sign of ϕ0 and is approximated as ϕ0/(ϕ02 + h)1/2 with h being the spatial step size.[48]

Check if a steady state is reached. If not, locate the interface Γ by the level-set function obtained in the previous step, set k: = k + 1, and go back to step 2.

Applications

One Charged Particle

We first consider a solute consisting of a single atom located at the origin carrying a point charge Q > 0. A solute–solvent interface in this case is a sphere centered at the origin with a radius denoted R. For this system, we have an exact analytical formula of the free energy II.7:This simple one-dimensional function can be minimized numerically with very high accuracy. We use the following parameters valid for methane-like particles in water at normal conditions: P = 0, γ0 = 0.175 kBT/Å2 with T = 300 K, τ = 1 Å, ρw = 0.0333 Å–3, ε = ε1 = 0.3 kBT, σ = σ1 = 3.5 Å, εm = 1, and εw = 80. We also solve the level-set eq II.8 with the normal velocity F given in eq II.8, starting from a large sphere centered at the origin.

In Table 1, we display the results of our level-set VISM calculations on a three-dimensional grid (marked level-set) and those of one-dimensional numerical optimization based on the above analytical formula G[R] (marked analytical) that have basically the vanishing error. We compare the nonpolar and polar parts of the free energy (in units kBT) and the optimal (equilibrium) radius R (in units Ångstrom) for different values of the point charge Q (in units e, the elementary charge). We see that level-set VISM is indeed very accurate for this system. In fact, the relative error for both the minimum free energy and the optimal radius is within 0.5%. Note that the total free energy for Q = 1e is consistent with the typical magnitude of solvation free energies of monovalent ions, cf., e.g., Dzubiella et al.[2] Note also that the higher the charge, the smaller the resulting optimal radius R. This is because the dielectric medium (multipolar water) gains more and more electrostatic energy by penetrating deeper into the solute region. Consequently, the hydrophobic (geometrical) part of the energy decreases because the radius R decreases. The vdW part of the free energy on the other hand increases significantly as the water pushes deep into the repulsive LJ part. This clearly indicates that both the polar and nonpolar contributions are inherently coupled in the solvation free energy.

Table 1

Solvation Free Energy (in units kBT) and Optimal Radii (in units Å) of the Spherical Solute for Different Values of Point Charges Q > 0 (in units e)a

	nonpolar energy		polar energy		total energy		optimal radii
charge	level-set	analytical	level-set	analytical	level-set	analytical	level-set	analytical
0.0	4.3767	4.3135	0	0	4.3767	4.3135	3.1336	3.1249
0.5	4.7816	4.7178	–22.719	–22.792	–17.938	–18.074	3.0100	3.0099
1.0	8.9565	8.9158	–97.894	–98.223	–88.937	–89.307	2.8008	2.7937
1.5	20.307	20.308	–236.45	–237.28	–216.14	–216.97	2.6090	2.6020
2.0	39.569	40.125	–445.52	–447.70	–405.95	–407.58	2.4624	2.4517
2.5	67.318	68.899	–730.91	–735.07	–663.60	–666.17	2.3454	2.3332

The full 3D level-set minimization (level-set) is compared to the 1D numerical optimization of the free energy with the analytical formula G[R] (analytical).

We now apply our VISM one-particle framework to single ions K+, Na+, Cl–, and F–. For all these systems, we use the surface tension γ0 = 0.175kBT/Å2 with the temperature T = 298 K, the water density ρw = 0.0333 Å–3, the dielectric coefficients εm = 1 and εw = 80, and the Tolman length τ = 1 Å. In our VISM-CFA calculation of the hydration free energy of anions, Cl– and F–, we use the dielectric boundary that is obtained by shrinking the solute–solvent boundary by the water OH bond length of 1 Å to include the asymmetry effect of water as pointed out in Dzubiella et al.[2] In Table 2, we display our VISM solvation free-energy values and the experimental free energies[49] for the hydration of these ions. We see that a good agreement between our VISM and experimental results is reached.

Table 2

Hydration Free Energies (in kBT) Obtained by VISM and by Experiment[49] (converted from kJ/mol to kBT) for Single Ions K+, Na+, Cl–, and F–a

ions	ε (kBT)	σ (Å)	VISM	experiment
K+	0.008	3.85	–112.3	–117.5
Na+	0.008	3.49	–131.1	–145.4
Cl–	0.21	3.78	–126.7	–135.4
F–	0.219	3.3	–171.9	–185.2

The second and third columns are the LJ parameters, taken from Horinek et al.[50]

Two Charged Particles

We now consider a molecular system of two methane-like particles in water and treat these particles as atoms. We use the parameters P = 0, γ0 = 0.175 kBT/Å2 with T = 300 K, τ = 1 Å, ρw = 0.0333 Å–3, ε1 = ε2 = 0.3603 kBT, σ1 = σ2 = 3.4418 Å, εm = 1, and εw = 80. The point charges Q1 and Q2 are chosen to have the same magnitude, i.e., |Q1| = |Q2|. We denote by d the center-to-center distance between the two atoms.

To test the accuracy of the CFA, we compare the electrostatic solvation energy values computed by our CFA and that using the adaptive Poisson–Boltzmann solver (APBS).[33] We fix the charges Q1 = +0.2e and Q2 = −0.2e. We examine a set of d values and use the same vdW surfaces for both of the CFA and APBS calculations. We also calculate the CFA values of electrostatic solvation energies using our optimal level-set VISM surfaces. All of these energy values are listed in Table 3. We see from the table that the CFA agrees with APBS well, particularly for large d where the vdW surface breaks into two disjoint parts and vdW and VISM surfaces are very similar. For small d, the CFA shows about 15% deviation to the APBS result. We also notice that the energy values of CFA with VISM optimal surfaces are about 10% different from those obtained using CFA with vdW surfaces. This is a consequence of VISM optimal surfaces being different from the vdW surfaces.

Table 3

Comparison of the Electrostatic Solvation Energy Values (in kBT) Computed by the CFA and That by APBS, Where d (in Å) Is the Center-to-Center Distance of the Two Atoms, vdW Means van der Waals Surfaces, and VISM Means Optimal VISM Surfaces

d	APBS-vdW	CFA-vdW	CFA-VISM
4	–2.3062	–1.9552	–1.7740
6	–3.6343	–3.4594	–3.2810
8	–4.5973	–4.5224	–4.3041
10	–5.1607	–5.1068	–4.8745

We now study the solvent-mediated PMF of the system with the center-to-center distance d between the two atoms as the reaction coordinate. For each of a few selected coordinates d and charge combinations Q1 = Q2 = 0.1e and Q1 = −Q2 = 0.1, 0.2, and 0.5e, we relax an initial surface to a final, stable equilibrium solute–solvent surface and calculate different parts of the minimum solvation free energy. In Figure 2, we plot the different components of the solvent mediated PMF. The full vdW component of the PMF, including the solute–solute vdW interaction, is plotted in the inset of Figure 2b. The total PMFs excluding and including the solute–solute vdW interaction are plotted in the main frame and inset, respectively, of Figure 2d. For comparison, the Coulomb law of charge–charge interaction for the two charged atoms is also plotted in Figure 2c.

Figure 2

Different contributions to the PMF of the two-atom system vs the center-to-center distance d for different values of charges. (a) The geometrical part Ggeompmf. (b) The vdW part GvdWpmf. The solute–solute vdW interaction is excluded in the curves in the main frame but included in those in the inset. (c) The electrostatic part Gelecpmf. The Coulomb law of the charge–charge interaction is also plotted for comparison. (d) The total PMF Gtotpmf. The solute–solute vdW interaction is excluded in the curves in the main frame but included in those in the inset.

Remarkably the nonpolar parts (Figure 2a and b) display a desolvation barrier at d = 5–6 Å as observed in explicit water simulations.[51] In our implicit-solvent model, it appears due to increasingly concave surface parts upon merging of the isolated solute surfaces. Somehow, the magnitude of distortion of the overlapping first solvation shells is reflected in increasing concavity. The water-induced vdW interaction (Figure 2b) displays a maximum at around d ≃ 5 Å because the dispersion interaction of the solutes with the water is minimal here. The electrostatic part of the PMF follows roughly the anticipated Coulomb interaction as shown in Figure 2c with corrections stemming from the (image charge) repulsion between the low-dielectric atomic cavities in the electric field. For small d ≲ 4 Å, the electrostatic interaction is strongly enhanced as the solvent dielectric screening is diminished for the overlapping cavities. Charging influences the nonpolar (geometric and vdW) parts for d ≲ 7 Å where the solvation shells merge and are perturbed. At this fragile point, the system seems most sensitive to a change in electrostatics, even in an implicit-solvent model. As we see in Figure 2d, charging thus mostly affects the total PMF below d ≃ 7 Å, where the interface starts overlapping. This is partially due to the subtle influence of electrostatics on surface geometry and not directly due to the electrostatic solute–solute interaction. Note that this phenomenon would not be describable by a simple fixed-surface (e.g., vdWS, SES, or SAS) model, where the solvent surface is predefined. For small d ≲ 4 Å, the electrostatic contribution dominates due to weak solvent dielectric screening. Note again that the shown total PMF in the main frame is only the solvent-mediated part of the PMF, i.e., it excludes the solute–solute vdW interaction. This interaction is included in the PMF plotted in the inset.

Two Parallel Charged Plates

Here, we consider the strongly hydrophobic system of two parallel paraffin plates taken from the work of Koishi et al.[52] Each plate consists of 6 × 6 fixed CH2 atoms and has a square length of about 3 nm. The two plates are placed at a center-to-center distance d. For these hydrophobic plates, capillary evaporation takes place at distances d ≲ 15 Å.[3,52] In the following, we investigate how (a) the capillary evaporation, (b) the hydrophobic attraction, and (c) a possible hysteresis in the free energy are affected by charging up the plates. To this end, we assign central charges q1 and q2 to the first and second plates, respectively, with |q1| = |q2|. The total charges of these two plates are 36q1 and 36q2, respectively. We study like-charged and oppositely charged plates by choosing the values of (q1,q2) to be (0e,0e), (+0.1e,–0.1e), (+0.1e,+0.1e), (+0.2e,–0.2e), and (+ 0.2e,+0.2e).

In Figure 3, we show 2D surface cuts through the 3D equilibrium surface for d = 10 Å and different charges. For q1 = q2 = 0e, the final state indeed depends on the initial surface (tight vs loose). This indicates that two local minima are present corresponding to wet and dry states. The dry state has been observed before where we did not check for hysteresis.[3] The bimodal hydration seems to be a general effect due to bubble nucleation barriers.[21] Charging the plates at d = 10 Å has different consequences on the final state depending on the sign of the charge. If the plates are charged oppositely, capillary evaporation is surpressed, and the final state is wet. This is because that the strong electrostatic potential between the plates drags the polar water into the void. If the plates are equally like-charged, then a stable capillary bubble remains with only a slightly tighter surface when compared to the case q1 = q2 = 0e. This is because the oppositely directed electrostatic field cancels out in the void, and the water distribution is hardly affected. A few 3D snapshots of the VISM surfaces are shown in Figure 4 for d = 10 Å and different charges with loose initial configurations, where the mean curvature is color-coded. We observe clearly that charging displays tighter surfaces, diminished capillary evaporation (totally surpressed for opposite charges (q1,q2) = (+0.2e,–0.2e)), and more concave parts of the surface. These examples highlight the sensitive coupling between electrostatics and hydrophobicity in aqueous solvation, especially when the system is prone to capillary evaporation.[2] The surpression of the latter between plates by introducing hydrophilic patches has been observed in recent MD simulations.[18]

Figure 3

2D cuts through the center of the 3D stable equilibrium solute–solvent interfaces around the two plates at d = 10 Å with tight or loose initial surfaces for different charge combinations.

Figure 4

Stable 3D equilibrium solute–solvent surfaces of the two-plate system obtained by the level-set VISM calculations with loose initials at d = 10 Å. From left to right: atomic charges (q1,q2) = (0e,0e), (+0.2e,+0.2e), and (+0.2e,–0.2e). The color code represents the mean curvature being convex (red), flat (green), and concave (blue).

A more detailed, quantitative assessment is provided in Figures 5 and 6 where we plot the different components of the full PMF (including the solute–solute vdW interaction) with loose and tight initial surfaces, respectively. For the loose initials (Figure 5), the geometric part displays a strong attraction below a critical distance dc at which capillary evaporation begins. The crossover distance decreases from dc ≃ 14 Å down to 9 Å for (q1,q2) = (+0.2e,–0.2e). Note that the opposite charging has a much stronger effect than like-charging due to the electrostatic field distribution discussed above. Also, the solute–solvent vdW part of the interaction is strongly affected by electrostatics due to the very different surface geometries induced by charging. Both curves Ggeompmf and GvdWpmf demonstrate the strong sensitivity of nonpolar hydration to local electrostatics when capillary evaporation occurs and very “soft” surfaces are present. The total PMF is not just the independent sum of nonpolar and electrostatic interactions. For the surfaces resulting from the tight initials, cf. Figure 6, the situation is a bit less sensitive to electrostatics, as the final surface is closer to the vdW surface for dc ≳ 6 Å. Still, the nonpolar parts, in particular the vdW contributions, are affected by charging the plates as the surface geometry for distances d close to the desolvation barrier can be strongly perturbed by the electrostatics, see Figures 3 and 4.

Figure 5

Different components of the full PMF vs separation distance d between the two plates for different charge combinations (q1,q2) (see legend) obtained by the level-set VISM with loose initial surfaces.

Figure 6

Different components of the full PMF vs separation distance d between the two plates for different charge combinations (q1,q2) (see legend) obtained by the level-set VISM with tight initial surfaces.

Finally, Figure 7 graphically displays the bimodal behavior and hysteresis by plotting the two different PMF branches stemming from the equilibria of wet and dry states. For the neutral plates (cf. Figure 7a), a strong hysteresis is present for 6 ≲ d ≲ 14 Å. Adding charges influences the free-energy branches and hysteresis as shown in Figure 7b and c. However, only in the case of oppositely charged plates are the changes significant, as a strong electrostatic field develops in between the hydrophobic plates. Here, the water occupancy is most sensitive to local potentials.

Figure 7

Comparison of the two branches of PMF corresponding to the wet and dry states for the two plates carrying different charges (q1,q2) (see legend).

The Protein BphC

Our last example is biphenyl-2,3-diol-1,2-dioxygenase (BphC), a key enzyme of biphenyl biodegradation pathway in Pseudomonos sp. The functional unit of this protein is a homo-octamer, and each subunit consists of two domains. We set up a series of configurations where the two domains are increasingly separated from d = 0 to d = 20 Å apart, perpendicular to their interface. The domain separation d is chosen here to be the reaction coordinate. Note that the zero domain separation corresponds to the native configuration in the crystal structure (PDB code: 1dhy).

To ensure the VISM with CFA of electrostatics is suitable for this system, we compared electrostatic solvation energy computed by the CFA and that by APBS using the same hard sphere dielectric boundary. Across the entire range of the system configurations (i.e., from 0 to 20 Å domain separations), the CFA reproduced APBS solutions reasonably well with a mean deviation of 14%.

MD simulations reported in Zhou et al.[39] suggested that water density in the inter domain region is lower than in the bulk when two domains are within 6 Å. When electrostatic interactions between the water and protein are turned off, the strong dewetting transition occurs at as far as 9 Å domain separation. In comparison, we display in Figure 8 six surfaces of our level-set VISM at three different domain separations with and without atomic partial charges. At d = 8 Å, the level-set VISM identifies the interdomain region as partially solvent excluded when atomic charges are included, and as completely solvent excluded without any atomic partial charges. The interface wraps around the protein more tightly with charges than that without them due to the attractive nature of the polar interactions between the solute and solvent. The dark blue regions on the surfaces correspond to the deep trenches created by a strong polar interaction. At d = 14 Å, the uncharged and charged BphC molecules pose topologically distinct solute–solvent interfaces. With polar interactions, both domains are completely solvated. In contrast, the center of the domain interface still retains low water occupancy without electrostatic interactions. This is consistent with the results from atomistic simulations, where dewetting extends to a much greater region without polar interactions. Compared with traditional surfaces, such as vdWS, SES, or SAS, the VISM surfaces are topologically similar at small and large interdomain separations. However, a traditional surface would break into two independent surfaces at d > 4 Å regardless of charge distribution. Here, VISM is able to capture the entire transition process self-consistently.

Figure 8

The stable equilibrium solute–solvent interfaces of BphC at three different domain separations, obtained by the level-set VISM with loose initial surfaces. The top row is with atomic partial charges, and the bottom row is without partial charges.

Figure 9 shows different parts of the PMF of BphC with respect to the separation of two domains, with and without partial charges, obtained by our level-set calculations using tight and loose initial surfaces. These PMFs include the solute–solute vdW interactions. They indicate that both charging and initial surfaces can strongly affect the PMF. We observe bimodal hydration in this system as well. With polar interactions, the (de)wetting process happens at smaller domain separations. This is illustrated by the general shifting of charged curves in a and b by ∼2 Å toward smaller separation. In addition, the bimodal hydration is damped by partial charges since they assist the nucleation process in the interdomain region. The sharp transitions seen in a and b for tight initials indicate a strong dewetting occurring around 6–8 Å for uncharged case, capturing the hydrophobic collapse that has been also reported in MD simulations.[39] It is noticeable that at both small and large separations the loose and tight initial conditions yield slightly different energetics. This is because the small differences remain in the final optimized surfaces when drastically different initial surfaces are used for complex systems. The optimized surface is only one of multiple local minimizers of the VISM free-energy functional. The coexistence of such multiple local minima represents the fluctuating nature of the solvent–solute interface.[21]

Figure 9

Different parts of the PMF of BphC with respect to the domain separations, with or without partial charges, and with loose and tight initial surfaces.

Conclusions

In this work, we introduce the Coulomb-field approximation (CFA) of electrostatic free energy into the variational implicit-solvent model (VISM) and implement the level-set numerical method for relaxing the resulting VISM free-energy functional. We apply our theory and method to molecular systems of various complexities, including the protein BphC. These applications demonstrate that our level-set VISM with the CFA is an efficient approach to qualitatively capturing the stable equilibrium solute–solvent interfaces with the correct corresponding free energies.

We have found from our extensive numerical data and detailed analysis using potentials of mean force (PMF) the following:

The coupling of the geometrical, van der Waals, and electrostatic contributions in the VISM is essential in qualitatively estimating solvation free energy, and capturing multiple equilibrium states of wet and dry. Such multiple states exist in biomolecules in solution. They lead to the system hysteresis and fluctuations. But they are hard to describe using implicit-solvent models of fixed-surface type.

Charges impact strongly on the hydration/dehydration process from a simple molecular dimer to a complex protein. The detail of local electrostatic field is critical for wet and dry transitions, as shown in the example of two paraffin plates. In the more complex BphC case, polar interactions enhance wetting, which dampens the strong hydrophobic collapse propensity posed by hydrophobic residues at the center of the domain interface. While these have been observed and explained in molecular dynamics (MD) simulations, our theory and method seem to be the first with an implicit solvent to capture such properties.

Compared with MD simulations, our level-set VISM with the CFA is computationally much more efficient, and also qualitatively, often even quantitatively, accurate. It is therefore promising to apply our theory and method to study efficiently the hydration of large and realistic systems.

We notice that the CFA is known to be an efficient but not necessarily accurate approximation of the electrostatic free energy. We emphasize that the purpose of our work is mainly to couple the CFA with our level-set VISM to demonstrate that the coupling of electrostatics and nonpolar contributions in VISM can capture qualitatively some important features, such as the capillary evaporation in hydrophobic confinement, that are otherwise hard to capture by other fixed-surface type implicit-solvent models. More accurate continuum descriptions of electrostatics, such as the Yukawa-field approximation[42] and the Poisson–Boltzmann theory,[43] need to and can be included in our level-set VISM. We are currently working on these improvements.

Finally, we remark that finding an energy-minimizing solute–solvent interface using our level-set VISM approach is usually computationally more expensive and less efficient than generating a van der Waals surface, solvent-excluded surface, or solvent-accessible surface. To completely relax a system with our method, it can take minutes to hours depending very much on the choice of initial surface and the numerical resolution. We are currently working on speeding up our calculations so that our approach can be combined with MD simulations. We note, however, that the application of our theory and methods is not limited to the calculation of electrostatics for MD simulations. As we have shown in our previous and current work, we can efficiently capture various kinds of equilibrium conformations of biomolecules and estimate the solvation free energies. Therefore, our approach can be useful for efficient molecular recognition. We are also developing related theory and methods to include solvent fluctuations for the real dynamics of biomolecular solvation.