Stochastic level-set variational implicit-solvent approach to solute-solvent interfacial fluctuations

Aqueous solvent plays a critical role in dynamical processes of biomolecules, such as protein folding, molecular assembly, and molecular recognition.^1–3 An accurate description of solvent effects is therefore crucial to understanding the underlying mechanisms of such biomolecular processes. Implicit-solvent models are efficient approaches to studying biomolecules solvated in an aqueous medium.^4–8 In a large class of implicit-solvent models, the solvent is described as a dielectric continuum with a uniform dielectric coefficient and the effect of solvent is described through the solute-solvent interface or dielectric boundary. However, thermal fluctuations from both solute atoms and solvent molecules are often ignored in such implicit-solvent models. Yet such fluctuations are known to be crucial in biomolecular dynamic processes.^9–13 

Within the framework of dielectric boundary implicit-solvent modeling, accurate predictions of solvation free energies and hydrophobic interactions of biomolecules rely sensitively on the definition of solute-solvent interface.^14–17 The commonly used Poisson–Boltzmann/Surface Area (PB/SA) implicit-solvent model is based on fixed solute-solvent interfaces, such as the van der Waals surface (vdWS), solvent-excluded surface (SES), and solvent-accessible surface (SAS).^18–22 While successful in many cases, these a priori defined surfaces may fail in capturing hydrophobic interactions, such as capillary evaporation^23–36 and hydrophobic-hydrophilic coupling.^36–38 

To address these issues, Dzubiella et al. developed a variational implicit-solvent model (VISM),^37,38 coupling contributions from surface energy, solute-solvent van der Waals (vdW) interaction energy, and electrostatic solvation free energy. Implemented with a robust level-set method, such a model is able to capture polymodal hydration states and provide qualitatively accurate estimates of solvation free energies.^{29,30,39–47} More precisely, in the VISM, one minimizes a macroscopic solvation free-energy functional of all possible solute-solvent interfaces to find stable, equilibrium biomolecular conformations that are often local minimizers of the functional as well as to determine the corresponding solvation free energies. The level-set method numerically implements this minimization process by moving an initially guessed surface with the “normal velocity” that is the VISM boundary force on the solute-solvent interface. With such an approach, for instance, we find bimodal hydration states of a host-guest system with different initial guesses of minimization.^44,45 In a receptor-ligand system, our level-set VISM calculations capture distinct metastable states that correspond to topologically different solute-solvent interfaces.^29,46 Each predicted hydration state, a local minimum of the VISM free-energy functional, contains much information of the solvated system in thermodynamic equilibrium, with contributions proportional to exp(−G[Γ_min]/k_BT).^29,30,44 Here G[Γ_min] is the VISM solvation free energy of the local minimum that has the solute-solvent interface Γ_min, and k_B is the Boltzmann constant and T temperature. Therefore, the local minima with lower free energies contribute more to the ensemble average of a quantity under consideration. The level-set method used here is to relax the solvation free energy with the steepest descent path and can only predict local minima of an underlying system. In order for an underlying system to jump out of a local minimum, one has to include thermal fluctuations to overcome energy barriers.

In this work, we develop a theory to include the solute-solvent interfacial fluctuations in the VISM framework. This theory is in the form of a Langevin geometrical flow: the motion of the solute-solvent interface is driven by the combination of both the VISM boundary force and the stochastic force arising from the solute-solvent interfacial fluctuations. We also develop a stochastic level-set method to implement our theory. Our numerical treatment of interfacial fluctuations with the level-set method originates from the recent work on computer visions by Juan et al.;⁴⁸ cf. also Ref. 49. We add the noise in the normal velocity in the level-set relaxation. We employ the Stratonovich stochastic calculus to formulate the stochastic level-set equation to ensure the mathematical well-posedness of the stochastic surface evolution.^50–52 To numerically solve the resulting stochastic partial differential equation (PDE), we discretize the time derivative by the Euler–Maruyama scheme. We discretize the parabolic-like spatial differential operators by the central differencing and hyperbolic-like spatial differential operators by a fifth-order weighted-essentially no-oscillations (WENO) scheme.⁵³ We combine such a numerical algorithm with a simulated annealing strategy for energy barrier crossing. In order to calculate geometrical quantities such as the curvature, we need to smooth out the fluctuating boundary. To do so, we design a numerical filter using the Fourier transform.

In the context of biomolecular solvation, dewetting transition corresponds to jumping out of a local minimum of a solvation free-energy functional. Water molecules in hydrophobic confinement are metastable because they cannot form a stable H-bond network. Induced by fluctuations, water molecules spontaneously evaporate to the bulk, generating void vapor adjacent to hydrophobic surfaces. Such a drying process, i.e., dewetting transition or capillary evaporation, has been found in experiment and computer simulations^{24–26,28,30,33–36,54–56} and is known to play a crucial role in many important biological processes, such as receptor-ligand binding.^3,32,57–59

We apply our stochastic level-set method to study the dewetting transition in two hydrophobic systems: two parallel hydrophobic plates, and the synthetic host molecule cucurbit[7]uril (CB[7]). By using the deterministic level-set method developed in our previous works,^{29,39,41,42,45,46} we are able to find multiple local minima with different initial guesses. We then start from one of the local minima with higher solvation energy and perform simulations many times to get different realizations of the stochastic level-set equation. Our extensive numerical results show that our method is indeed able to help the system jump out of an initial local minimum, capture the evolution of solute-solvent interface in dewetting transitions, and estimate well the corresponding activation energy barrier. In particular, for the two-plate system, we find that the estimated activation energy barrier of dewetting transitions increases quadratically with respect to inter-plate distances. This agrees very well with the results obtained by molecular dynamics (MD) simulations in the literature.^34,35,54 In addition, we also study how the fluctuations with different wavelengths can affect the dewetting process. We find that, as the wavelength of fluctuation decreases, the dewetting probability first increases and then decreases slightly and the dewetting energy barrier first decreases very little and then increases largely. Some of these studies are also carried out for the solvation of CB[7], a nonpolar system. Our numerical simulations successfully capture the dewetting transition in the toroidal shaped cavity of the molecule and provide the corresponding activation energy barrier from a wet state to a dry state.

We organize the rest of our paper as follows: In Section II, we describe our theory of VISM with solute-solvent interfacial fluctuations. In Section III, we present the stochastic level-set method that implements our theory. In Section IV, we apply our theory and method to the two-plate system and CB[7]. Finally, in Section V, we draw conclusions.

We consider the solvation of biomolecules in an aqueous solvent. The geometry of such a solvation system consists of three parts: the solute region Ω_m (m for solute molecule), the solvent region Ω_w (w for water), and the solute-solvent interface Γ, cf. Fig. 1. The solute molecules have a total of N atoms inside Ω_m and their positions are labeled by x₁, …, x_N.

We treat the solvent as a structureless continuum and its thermodynamic contributions are described by a few parameters, such as the bulk solvent density, the solute-solvent interfacial surface tension, and the solvent dielectric coefficient. The underlying biomolecular system is therefore determined well by the solute atomic positions x₁, …, x_N and the solute-solvent interface Γ which we shall identify as the dielectric boundary. We describe the fluctuating relaxation dynamics of the biomolecules with implicit solvent by the following system of generalized Langevin equations:

Here a dot denotes the time derivative and ∇ the spatial gradient. The effective Hamiltonian H = H[x₁, …, x_N; Γ], detailed below, combines the solute atomic molecular mechanical interactions and the VISM macroscopic solvation free energy through the interface Γ. Eq. (2.1) is the Langevin equation for the fluctuating relaxation dynamics of the jth solute particle x_j = x_j(t), with M_j the mobility coefficient and η_j the noise. Eq. (2.2) is the Langevin geometrical flow describing the fluctuating relaxation of the solute-solvent interface Γ = Γ(t), where M_Γ is a mobility constant, η_Γ is the noise term, and δ_Γ denotes the variation with respect to the location change of boundary Γ.

The effective Hamiltonian is given by

The first part V = V[x₁, …, x_N] is the potential of solute molecular mechanical interactions, consisting of bonding, bending, torsion, Lennard-Jone (LJ) potentials, and Coulomb potentials, as in a usual forcefield for solute atomic interactions. The second part is the VISM solvation free-energy functional, given by

Here, the first term describes the work it takes to create a solute cavity volume in a solvent medium at hydrostatic pressure P₀. The second term is the solute-solvent interfacial energy, where γ is the local interfacial surface tension that is allowed to have a curvature correction. Practically, we often set γ = γ₀(1 − 2τH), where γ₀ is the constant surface tension for a planar interface, τ the Tolman curvature correction coefficient,⁶⁰ and H the local mean curvature (i.e., the average of the two principal curvatures). We call the sum of the first two terms the geometrical part of the solvation free energy and denote it by G_geom[Γ]. These terms do not depend explicitly on x₁, …, x_N.

The third term describes the solute excluded-volume effect and the solute-solvent van der Waals (vdW) type dispersive interaction. The constant ρ_w is the bulk solvent density. For each i, the term U_i(|x − x_i|) describes the interaction between the solute atom at x_i and a solvent molecule or ion at x. We often use the LJ potentials

where the parameters ε_i of energy and σ_i of length can vary with i. We call the third term the van der Waals (vdW) part of the solvation free energy and denote it by G_vdW[x₁, …, x_N; Γ].

The last term G_elec[x₁, …, x_N; Γ] is the electrostatic part of the solvation free energy. It is determined by the Poisson–Boltzmann (PB) theory^61–65 with Γ the dielectric boundary and solute partial charges at x₁, …, x_N. More precisely, given the boundary Γ, we can solve the linear or nonlinear PB equation to get the electrostatic potential, with the dielectric coefficient being close to 80 in the solvent region and 1 in the solute region, respectively, and with the fixed charge density the superposition of all the partial point charges at x_i (1 ≤ i ≤ N). Then we can obtain the electrostatic free energy through a usual integral formula. See details in our previous work.^45,64

As the first step to tackle the problem of solute-solvent interfacial fluctuations with implicit solvent, in this work, we shall focus on the solute-solvent interfacial fluctuations for nonpolar, hydrophobic molecular systems. Therefore, we shall fix all the solute atomic positions x₁, …, x_N and neglect the electrostatic part in the Hamiltonian. We shall also set the Tolman correction parameter τ = 0 for our initial studies and set the mobility constant M_Γ = 1 as for the general case can be treated by choosing a suitable numerical time step. For notational simplicity, we shall write G[Γ] for G[x₁, …, x_N; Γ] and G_vdW[Γ] for G_vdW[x₁, …, x_N; Γ], respectively. Now our governing equation is the following simplified version of Eq. (2.2):

With the VISM free-energy functional G[Γ] now only including the geometrical part (without the curvature correction) and the vdW part, the negative variation of the VISM free energy G[Γ] with respect to the location change of Γ is given by^39,41

This is the effective VISM force on the boundary Γ.

In our previous works,^{29,39–42,45} we developed a robust and efficient level-set method to minimize the free-energy functional (2.4). We begin with an initial surface that encloses all the solute atoms and then move the surface in the direction of steepest descent of free energy until a steady state is reached. The surface at time t is represented by a zero level-set function ϕ = ϕ(x, t), i.e., the points x on the surface are exactly those such that ϕ(x, t) = 0. With such a level-set function, geometric quantities of the surface, such as the unit normal, mean curvature, and Gaussian curvature, can be explicitly obtained.⁴¹ For instance, the local mean curvature is given by

where Δϕ and ∇²ϕ are the Laplacian and Hessian of ϕ, respectively. In what follows, we shall write sometimes H = H(∇ϕ) to indicate the dependence of mean curvature H on the gradient of a level-set function ϕ. We shall also express the VISM solvation free-energy functional as a functional of the level-set function ϕ and denote it by G[ϕ], when no confusion arises.

We track the relaxation motion of surface by solving the level-set equation

for the level-set function ϕ(x, t). In this equation, F_n is the normal velocity (i.e., the normal component of velocity) with which the surface moves. The normal velocity is defined only on the surface but is extended to the entire computational box so that Eq. (3.2) can be solved. To minimize the solvation free-energy functional G[Γ], we choose the normal velocity to be the boundary force F_n = − δ_ΓG[Γ] as given by Eq. (2.7), i.e.,

Note that there is a pre-factor, which we choose to be the unity, with suitable units that converts the normal velocity to the normal component of force.

We model the Langevin geometrical flow (2.6) and (2.7) by adding the following noise term:⁴⁸ 

to the normal velocity F_n defined in (3.3) of the level-set equation (3.2). Here M is a numerical value, the functions W_i(t) (1 ≤ i ≤ M) are the components of a standard M-dimensional Wiener process, and the functions ψ_i(x) are smooth functions with small supports (i.e., each ψ_i vanishes outside a small region). We obtain therefore the stochastic level-set equation

This stochastic partial differential equation (PDE) is expressed in the sense of Stratonovich (indicated by the symbol ∘) so that the surface representation is independent of the choice of level-set function.⁴⁸ Direct numerical integration of this equation requires an implicit scheme, which is numerically inefficient due to the highly nonlinear terms in the drift coefficient. As in Refs. 48 and 49, we therefore convert Eq. (3.5) into an equivalent form of Itô that has a modified drift term,

For convenience of numerical discretization, we further write this equation as⁴¹ 

where

We use the forward Euler–Maruyama method to discretize the time derivative in (3.7),

where ϕ^m(x) is the approximation of ϕ(x, t_m) at time t_m = t_m−1 + Δt (m = 1, 2, …) and ΔW = (ΔW₁, …, ΔW_M) is a random vector with $ΔWi∼Δt N(0,1)(1≤i≤M)$⁠. Here and below Δt is the time step and $N(0,1)$ the standard normal distribution.

We denote by h the spatial step size. By the Courant–Friedrichs–Lewy (CFL)⁶⁶ condition, we set $Δt=C̄h2$ for some constant $C̄$ that is determined by trials. We have by (3.1) that A(ϕ) = D(∇ϕ) : ∇²ϕ, where the matrix

Here I is the 3 × 3 identity matrix and P = I − ∇ϕ ⊗ ∇ϕ/|∇ϕ|². Consequently, the eigenvalues of D(∇ϕ) are $(1/2)∑i=1Mψi2$⁠, γ₀, and γ₀, all positive. This means that the PDE ∂_tϕ = A(ϕ) is parabolic. Hence, we discretize A(ϕ) using the central differencing.⁴¹ For the hyperbolic-like term B(ϕ), we use the central differencing to discretize ∇ψ_i and a fifth-order WENO scheme⁵³ to discretize ∇ϕ. We also use the WENO scheme to discretize |∇ϕ| in C(dt, dW)|∇ϕ|. For the C(dt, dW) term, we first use the level-set function to interpolate the interface points and calculate the values of C on these points. We then use the fast marching method to extend these values in the normal direction to all the grid points.

To prevent the gradient of the level-set function ϕ being too large or too small, we reinitialize the level-set function ϕ every few time steps by solving

where ϕ₀ is the level-set function before reinitialization, sign (ϕ₀) is the sign of ϕ₀, and the time t is in general different from the one in the original level-set equation. We start from the initial value ϕ₀ at t = 0 and solve the equation with a few steps of time iteration.

In our implementation, we choose the functions ψ_i (1 ≤ i ≤ M) in the stochastic level-set equation (3.5) to be

where x_i is a finite-difference grid point, indicating the location of a source of noise. The parameter μ_i is the maximum strength of noise at x_i, and the parameter σ_i, related to the size of support of ψ_i(x), is usually chosen much less than the finite-difference grid size. We distribute the noise at grid points that are likely to have enhanced fluctuations. For instance, we put noise at regions that are inside the solute-solvent interface obtained with a loose initial surface (i.e., an initial, large surface) but outside the solute-solvent interface obtained with a tight initial surface (i.e., an initial surface that tightly wraps the solute atoms).

Since each function ψ_i is localized at the grid point x_i, the noise we add is spatially uncorrelated, leading to oscillatory solute-solvent interfaces. As a result, large numerical errors occur in direct finite-difference approximations of spatial derivatives near the interface and the solvation free energies. In addition, the dewetting transition described by a highly oscillatory interface may have much larger energy barrier. Therefore, in our implementation, we smooth out the interface by numerically filtering out the high frequency part of noise. More precisely, we first transform the noise in real space to the frequency space by using the Fast Fourier Transform (FFT); we then keep all the modes with wavenumber up to k_p, some critical wavenumber, and do Inverse Fourier Transform (iFFT) to get a smoothed noise.

In order for an underlying system to jump out of local minima, we apply a simulated annealing scheme to implement the stochastic geometrical flow (3.7). In such a scheme, we accept surface motions that decrease the free-energy functional and accept the other cases with probabilities depending on the amount of increased energy. This method also helps to prevent irregular surface motions that may lead to the overlap of the solvent and solute atoms, which can cause a dramatical increase in the vdW interaction energy. The simulated annealing scheme becomes more and more deterministic as the annealing temperature decreases. In the end, it almost only accepts interface motions that decrease the free energy. To accelerate the iterations, we remove the noise after certain iteration steps and switch to the deterministic level-set calculation.

• Input all the parameters P₀, γ₀, ρ_w, and atomic parameters x_i, ε_i, and σ_i for all i = 1, …, N. Initialize the annealing temperature T₀. Choose a computational box and discretize the box uniformly with a prescribed computational grid size h. Determine the location x_i, strength μ_i, and size of support σ_i of the noise function ψ_i (i = 1, …, M). Input the numerical time step Δt, numerical error tolerance, and the noise filtering parameter k_p. Generate an initial interface and initialize the corresponding level-set function ϕ⁰. Set the counter of time iteration m = 0.

• Generate the tight and loose VISM equilibrium surfaces by solving the deterministic level-set equation with a tight initial surface and a loose initial surface, respectively.

• Generate noises on the finite-difference grid points bound by the tight and loose VISM equilibrium surfaces, and filter the high frequency part by using the FFT.

• Solve the stochastic level-set equation (3.7) with the initial ϕ^m to obtain a new level-set function $ϕ̄$⁠.

• Solve the Equation (3.8) to obtain a reinitialized level-set function ϕ^∗ from $ϕ̄$⁠.

• Evaluate the free energy G[ϕ^∗].

• If G[ϕ^∗] < G[ϕ^m], accept the level-set function and set ϕ^m+1 = ϕ^∗; else accept the ϕ^∗ with the probability exp(−(G[ϕ^∗] − G[ϕ^m])/k_BT(m)), where $T(m)=T0/m$⁠.

• Check the stopping criteria. If failed, go back to Step 2.

We use the TIP4P water model to determine the parameters for water and use the Lorentz–Berthelot mixing rules to determine the parameters of LJ potentials of interaction between water and individual solute atoms. In our simulations, we use k_BT and Ångström (Å) for the units of energy and length, respectively. Unless specified otherwise, we fix T = 300 K, P₀ = 0 bar, γ₀ = 0.1315 k_BT/Ų, and ρ_w = 0.0331 Å⁻³.

We study the solvation of a strongly hydrophobic system of two parallel paraffin-like model plates.^{36,39,42,45,67} Each plate consists of 6 × 6 CH₂ atoms and has a side about 30 Å in length. The LJ parameters for each of these solute atoms are ε = 0.265 k_BT and σ = 3.538 Å. We denote by d (Å) the center-to-center distance between the two parallel plates and choose it as the reaction coordinate of the two-plate system.

In our previous works,^{29,39,41,42,45} we used the deterministic level-set method to study the equilibrium conformations and potential of mean forces of this system. We observed two types of hydration states that are local minima of the VISM solvation free-energy functional: wet state and dry state; cf. Fig. 2. The wet state is obtained by using a tight initial surface, which is close to a vdW surface of the two plates. The dry state is obtained by using a loose initial surface, which is a loose and large surface enclosing all the atoms in two plates. The solvation free energy of the wet state is higher than that of the dry state. For instance, it is 520.5 k_BT vs. 458.2 k_BT for d = 12 Å. Clearly, water molecules between the two plates with certain inter-plate distances are metastable due to the enthalpy penalty that stems from less formed hydrogen bonds compared to the bulk phase. It is well known that the water molecules in such a hydrophobic confinement evaporate spontaneously,^{23–35,54,56,68} transitioning from the wet state to the thermodynamically more favorable dry state. Our deterministic level-set VISM calculations can capture these different states, but, due to the steepest descent nature of the method, are unable to help the system jump out of a local minimum, transitioning from one state to another. Here we apply our newly developed stochastic approach to investigate the transition pathway of the solute-solvent interface.

Fig. 3 displays the pathway of dewetting transition of the two-plate system with the inter-plate distance d = 12 Å, predicted by our stochastic level-set calculations, in which the noise filtering parameter is set to be k_p = 0.17 Å⁻¹. Our simulations begin with the wet state. With the noise, the solute-solvent interface fluctuates, generating a vapor bubble in between the two plates. With the simulated annealing, trials of the interface motion that are too close to the solute atoms are rejected due to the energy penalty caused by the overlap of water molecules and solute atoms. Subsequently, the vapor bubbles from two sides come closer and begin to connect, creating a thin vapor tube, spanning the two plates; cf. the fourth snapshot in Fig. 3. Such a state has the highest free energy along this transition pathway and is believed to be a transition state of the dewetting transition. Our finding is in a good agreement with a well-known result that the transition state involves the nucleation of a critical vapor tube spanning the region between the two plates.^34,35,54,55 We also note that the water fluctuation plays a critical role in forming such a vapor tube.^54,56 Clearly, the interfacial fluctuation is the main driving force behind the formation of the vapor tube. Since we used the simulated annealing strategy, the probability of accepting interface motions that increase the energy is very low in the late stage of interface evolution. After certain time steps, we turn off all fluctuations to accelerate the convergence of interfacial evolution to the dry state. Notice that the number of steps labeled on the horizontal axis in Fig. 3 is not uniform. Overall, the figure illustrates that our stochastic level-set approach can help the system jump out of a local minimum and capture the pathway of capillary evaporation.

Next, for each of the several inter-plate distances varying from d = 8 Å to d = 16 Å, we perform 200 realizations and evaluate the solvation free energy (2.4) (with no curvature correction and electrostatic energy). Fig. 4 (left) displays the VISM solvation free energy of the system against the time steps with the inter-plate distance d = 12 Å. The blue lines represent different realizations and the red line is the mean energy of these realizations. Note that the number of steps labeled on the horizontal axis is not linear. Clearly, we can see that it requires climbing an energy barrier to go to the state of lower solvation free energy. For each realization, there is an energy barrier that is defined as the difference between the maximum energy value along the pathway and the energy value of the (initial) wet state. For ease of reading, we have shifted all the profiles so that the initial state has zero energy value. We use the minimum energy barrier among all these 200 realizations to estimate the activation energy barrier of the dewetting transition. We find that this barrier is 8.6 k_BT.

Since the dewetting transition involves the formation of a vapor tube, the volume of such a vapor region is a quantity of interest.^35,54 In our level-set setting, this volume is given as the difference of volumes of the solute region (i.e., the solvent excluded region defined by our VISM surface) at certain time step and that of the wet (initial) state.⁴⁶ We use the volume of solute region as a parameter to probe the transition state in the dewetting process. Fig. 4 (right) shows the relation between the free energy and volume of solute region. The two insets are zoom-in details near the wet state and the dry state, respectively. Again, the energy values are shifted so that the values on the vertical axis represent the energy differences compared to the wet (initial) state. As shown in the first set, the transition state occurs when the volume of vapor tube is roughly about 238 ų.

Assuming a cylindrical tube spanning two plates, one can use the macroscopic theory of thermodynamics to explain the dependence of energy barrier on the inter-plate distances.^{34,35,54,56,68} However, MD simulations³⁴ and lattice-model simulations⁵⁴ reveal that the nucleated vapor tube spanning two plates is not necessarily a cylindrical shape right in the center of the two plates. Our stochastic level-set method captures irregularly shaped tubes spanning two plates.

Fig. 5 shows our stochastic level-set VISM prediction of the dependence of the energy barrier on the inter-plate distance d. We can see that the data points, results of our stochastic level-set calculations marked “StoLSM,” fit perfectly well with a quadratic polynomial

where β is the inverse of k_BT. This agrees well with the existing MD simulations.^34,35,56

The interfacial fluctuations of long and short wavelength (or, equivalently, small or large wavenumber) play different roles in the dewetting transition. In an early stage of such transition, long-wavelength fluctuations, with wavelength larger than the average diameter of vdW spheres of solute atoms, reduce the strong attraction between the solute-solvent interface and the solute atoms in each of the two plates. Such fluctuations therefore make it easy for the solute-solvent interface to get away from the solute atoms. Once small and scattered parts of the surfaces touch, in the next stage, short-wavelength fluctuations help the fast merging of the two branches of solute-solvent interface.^54,55 Too oscillatory surfaces, however, greatly increase the surface energy largely, giving rise to difficulties in accurate estimation of the solvation free energy and the energy barrier. In order to better understand the role of fluctuations in terms of their wavelengths, we perform 200 simulations with various values of k_p in the dewetting transition of the two-plate system with the inter-plate distance d = 12 Å. We recall that the parameter k_p is the largest wavenumber kept in the frequency space after the Fourier transform of spatial noises on finite-difference grid points (cf. Section III C).

Table I shows the estimated energy barrier and the probability of dewetting transitions. This probability is defined as the ratio of the number of simulations dewetted to the number of total simulations. In the second column, we can see that the dewetting probability first increases with a larger k_p up to 0.22 Å⁻¹. This result can be explained by the findings obtained using Monte Carlo simulations on lattice gas systems.⁵⁴ The long-wavelength fluctuations promote the interface to be away from the solute atoms. Later, the short-wavelength fluctuations come into effect by nucleating vapor bubbles, which further develop into vapor tubes connecting the two plates. As such, short-wavelength fluctuations are necessary in the dewetting transition. This explains why simulations with relatively small k_p values have a lower dewetting probability. As k_p gets larger (greater than 0.22 Å⁻¹), the dewetting probability slightly decreases. In our simulations, we find that large wavenumber fluctuations give rise to a very oscillatory interface which generates many very small nucleation bubbles. Some bubbles over a critical size may grow and develop into vapor tubes connecting two plates; however, small ones may shrink and disappear later. If the interface is too oscillatory, generated bubbles may be too small to grow. This fact is responsible for the slight decrease in the dewetting probability as k_p becomes large.

Table I also displays the energy barriers estimated by the simulations with different k_p values. There are two factors that can contribute to the estimated energy barriers: the number of dewetted cases among a total of 200 simulations (or dewetting probability) and the wavenumber k_p of fluctuations. In general, more dewetted cases can give a better estimate, i.e., a lower energy barrier. Fluctuations with short wavelength lead to oscillatory interfaces, resulting in high surface-area penalty and energy barrier. The results shown in the table reflect the competition between these two factors. As k_p grows from 0.06 Å⁻¹, the estimated energy barrier first gently decreases to the lowest value when k_p = 0.17 Å⁻¹, due to the increasing dewetting probability. After that, the estimated energy barrier increases quickly to large values, on account of short-wavelength fluctuations. Our results verify that fluctuations of various wavelengths play different, important roles at different stages of deweting transitions.⁵⁴ 

We now apply our stochastic level-set method to study the solvation of Cucurbit[7]uril (CB[7]), a synthetic host molecule. This is a small molecule with relatively simple chemical complexity; yet has ultra high binding affinity of aromatic and metal guests,^{44,45,69–72} e.g., the bicyclo[2.2.2]octane (B2) derivative. Since water molecules cannot form a stable H-bond network due to the limited hydrophobic confinement, one of the main contributions to the binding affinity comes from the replacement of water molecules in the hydrophobic cavity of CB[7].^69,72 MD simulations show that the number of water molecules in the cavity exhibits enhanced fluctuations compared to the water molecules in the bulk phase. Induced by solvent fluctuations, dewetting transitions have remarkable impact on the binding kinetics of guests to the CB[7].^57,59 Here we use our stochastic level-set VISM approach to investigate the dewetting transition of the hydrophobic cavity of CB[7].

The solvation of nonpolar CB[7] has been systematically studied by the deterministic level-set method in our previous works,^44,45,47 where we found a wet state and a dry state, cf. Fig. 6. The color indicates the mean curvature of the surface: blue for concave and red for convex curvature. With the deterministic level-set method, if we start from a tight initial guess, then we get the wet state, where the water molecules can penetrate into central cavity; if we start from a loose initial guess, then we get the dry state, where the water molecules stay out of the central cavity. The solvation free energy of the wet state is higher than that of the dry state, with an energy difference of 4.98 k_BT. MD simulations confirm that there are enhanced fluctuations in the hydrophobic cavity and the molecule without electrostatics is mostly found in the dry state.^44,57

We first run the deterministic level-set method to find the equilibrium solute-solvent interface for the dry state and the wet state. Then we distribute noise sources on the finite-difference grid points that are inside the equilibrium interface of the dry state but are outside that of the wet state. To illustrate the pathway of dewetting transition, we present one realization in Fig. 7 as an example. We begin with the wet state which is one local minimum of the free-energy functional (2.4) (without the Tolman curvature correction and electrostatic part). To highlight the interface evolution along the pathway of dewetting transition, we display snapshots of cross sections of the solute-solvent interface in insets. We can see that the solute-solvent interface around the cavity fluctuates. As the energy increases, such an interface breaks into two parts, indicating that water molecules inside the cavity begin to evaporate. As in the two-plate system, the formation of vapor cavity seems to be a critical state of the dewetting transition. Again, we turn off all fluctuations after certain time steps to accelerate the convergence of interfacial evolution to the dry state. Also, we notice that the number of steps labeled on the horizontal axis in Fig. 7 is not uniform.

We perform 200 simulations on CB[7] in the same way as for the two-plate system; cf. Section IV A. Fig. 8 (left) displays the change of free energy against the time step. The blue lines are individual realizations and the red line represents the mean free energy of these realizations. We observe that the system needs to overcome an energy barrier to reach the dry state. We can estimate the activation energy barrier of the dewetting transition by finding the realization that has the lowest energy barrier. It is about 1.04 k_BT. Fig. 8 (right) shows the relation between the solvation free energy and the volume of solute region. Clearly, we can see that there is an energy barrier of 1.04 k_BT when the volume of solute region is about 2064 ų.

We have constructed a stochastic variational implicit-solvent model for molecular solvation. Our governing equation is the Langevin geometrical flow that describes the fluctuating relaxation of a solute-solvent interface. We have also developed a stochastic level-set method that numerically implements such a model. In our method, we add the noise into the normal velocity of the level-set equation. We have designed finite-difference schemes and combined them with numerical filters with the FFT and a simulated annealing strategy, to solve the stochastic level-set equation for the fluctuating relaxation of solute-solvent interface. With our theory and methods, we have studied the two-plate system and the synthetic host molecule CB[7].

Our computational results demonstrate that the interfacial fluctuation can help an underlying system jump out of a local minimum and reach a different local minimum with a lower free energy. Our theory and methods can capture the relaxing motion of a fluctuating solute-solvent interface during the dewetting transition of hydrophobic interactions. The evolution of corresponding solvation free energy reveals that an energy barrier is crossed in such interfacial motion. For the two-plate system, such an energy barrier occurs at the onset of a vapor tube between the two plates. Moreover, we predict successfully the quadratic dependence of the solvation free energy on the plate separation distance, agreeing well with existing MD simulations. We have also studied the role of interfacial fluctuations with different wavelengths in the dewetting transition and found that, as the wavelength of fluctuation decreases, the dewetting probability first increases and then decreases slightly and the dewetting energy barrier first decreases very little and then increases largely. All of these are achieved within our VISM framework and are computationally very efficient.

We now discuss several issues and possible further improvements of our approach. First, solute-solvent interfacial fluctuations in biomolecular systems are in general rather complicated.^14,15,73,74 For instance, an equilibrium interfacial structure under the influence of fluctuations should be described with the fluctuation-dissipation theorem. In our current setting, this is related to the detailed balance in the sense that the probability of the interfacial system getting into and out of a local minimum should correspond to the change in free energy. Moreover, the correct statistics of interfacial fluctuation is directly related to the kinetics of interface, which can be practically very important in biomolecular processes. To get the correct statistics of an equilibrium system, we need to sample all different states of interface motions with probability of transitioning to others satisfying the detailed balance. Here, we have taken a simple approach to treating the interfacial fluctuation by adding the noise in the normal velocity, or equivalently the boundary force. Our results show that this is a good approximation. However, we have used rather empirical descriptions of the location and strength of the noise. Moreover, the coupling or decoupling of the temporal and spatial noises is not well understood and characterized. More theoretical work is therefore needed. In particular, correct statistics can help to sample different hydration states and predict accurately the solvation free energies. We are currently working on the development of interface motions that satisfy the detailed balance in each motion step. With more reasonable treatment of the noises, quantitative characterization of the physical picture between the wavenumber k_p and the inter-plate distance d for the two-plate system will be another future work of importance.

Second, it is known that solvent density fluctuations contribute significantly to many dynamical processes in biomolecular systems, such as protein folding and protein-ligand binding.^{1,3,13,59,75–77} Assuming a Gaussian statistics, Lum–Chandler–Weeks⁷⁸ include such fluctuations in a continuum solvent theory for hydrophobic interactions.²³ In this theory, the solvent density is decomposed into a macroscopic part, which is governed by a van der Waals functional of gradient-squared term and a double-well potential for the solute-solvent interfacial structure, and a microscopic part, which describes the solvent fluctuations; cf. also Refs. 79–81. Within our dielectric-boundary based VISM, which is efficient in the description of electrostatic effects, we can possibly include the solvent fluctuation through a fluid mechanics approach using the Landau–Lifshitz random stress tensor,^82,83 which, though, can be more computationally expensive. Another possible approach is to design those spatial “basis” functions ψ_i(x) in the noise (3.4) to describe the solvent fluctuation, similar to a gas-lattice model,^76,77 somewhat ad hoc but can be straightforward and efficient.

Third, we have only considered nonpolar systems in our current work. It is natural to try to extend our approach and numerical methods to include the description of electrostatics using the Poisson–Boltzmann (PB) continuum theory, as electrostatic interactions in biomolecular systems are known to be crucial and yet can be computationally costly. Water molecules are often attracted to charged solute atoms, and such an attraction substantially suppresses the dewetting transition. A simple way to include the electrostatic interaction is to use the Coulomb-field approximation.⁴² Dielectric-boundary PB equation, linear or nonlinear, can also be solved with an additional step of smoothing out the fluctuating dielectric boundary for solving the PB equation.

Finally, we note that we have not been able to capture the drying transition in our simulations with a reasonable amount of computational cost. This is in fact a general problem of modeling and simulation of barrier crossing. There are several possible ways that we can pursue. One is to make our computation fast by using parallel computing. But a more promising way can be to use a better modeling strategy, such as the string method,^84–86 to efficiently connect different hydration states in the fluctuating relaxing motion of solute-solvent interface.

We thank the anonymous reviewers for helpful comments. This work was supported by Soochow University through a start-up Grant No. Q410700415 (S.Z.), NSF through Grant No. DMS-1319731 (L.-T.C. and B.L.), and NIH through Grant No. R01GM096188 (L.-T.C., B.L., and J.A.M.). J.D. acknowledges funding from the ERC (European Research Council) within the Consolidator Grant with Project No. 646659–NANOREACTOR. The work in the McCammon group is supported in part by NSF, NIH, HHMI, NBCR, and SDSC.