Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Sci Comput. Author manuscript; available in PMC 2017 May 1.
Published in final edited form as:
PMCID: PMC4922513

Numerical Treatment of Stokes Solvent Flow and Solute-Solvent Interfacial Dynamics for Nonpolar Molecules


We design and implement numerical methods for the incompressible Stokes solvent flow and solute-solvent interface motion for nonpolar molecules in aqueous solvent. The balance of viscous force, surface tension, and van der Waals type dispersive force leads to a traction boundary condition on the solute-solvent interface. To allow the change of solute volume, we design special numerical boundary conditions on the boundary of a computational domain through a consistency condition. We use a finite difference ghost fluid scheme to discretize the Stokes equation with such boundary conditions. The method is tested to have a second-order accuracy. We combine this ghost fluid method with the level-set method to simulate the motion of the solute-solvent interface that is governed by the solvent fluid velocity. Numerical examples show that our method can predict accurately the blow up time for a test example of curvature flow and reproduce the polymodal (e.g., dry and wet) states of hydration of some simple model molecular systems.

Keywords: Nonpolar molecules, solute-solvent interface, the Stokes equation, ghost fluid method, level-set method, interface motion, change of volume, traction boundary conditions

1 Introduction

Aqueous solvent plays a significant role in dynamical processes of biological molecules, such as conformational changes, molecular recognition, and molecular assembly, that control cellular functions of underlying biological systems [2, 12, 21, 22]. Implicit-solvent models are efficient descriptions of such dynamical processes. In such descriptions, the solvent is treated implicitly as a continuum and the effect of individual solvent molecules is coarse grained [3, 15, 23]. One of the successful dielectric boundary based implicit-solvent approaches is the variational implicit-solvent model (VISM) [7, 8]. In VISM, one minimizes a macroscopic solvation free-energy functional of all possible dielectric boundaries, coupling both nonpolar and polar contributions, and the solute-solvent van der Waals (vdW) interaction. Implemented by a robust level-set numerical method [4, 5], VISM can predict polymodal states of hydration, such as wet and dry states, subtle electrostatic effects, and solvation free energies of an underlying bimolecular system [4, 9, 18, 25, 28, 29].

While dielectric boundary based implicit-solvent models, including VISM, have been successful in many cases, they treat the solvent as a structureless dielectric medium, neglecting other solvent effects, such as the solvent hydrodynamic effect. Recent experimental and theoretical studies have indicated that the solvent shear motion can induce protein conformational changes and the solvent viscosity can affect the kinetics of such changes [1, 10, 11, 16, 17, 1921, 24].

In several recent works [13, 14, 26, 27], the authors have initiated the development of a fluid mechanics approach to treat the solvent fluid in molecular systems. The key features of such a new approach include: (1) the aqueous solvent (i.e., water or salted water) is treated as an incompressible fluid and its motion is by the Stokes or Navier-Stokes equation; (2) the solute pressure is simply described by the ideal-gas law; (3) the electrostatic interactions are modeled by the Poisson or Poisson–Boltzmann equation; and (4) all viscous force, electrostatic force, and vdW force are balanced on the solute-solvent interface that moves with solvent velocity. White [26] proposed to add the Landau–Lifshitz random stress tensor in the Stokes equation to model the solvent fluctuations. They also propose to describe the electrostatic interaction through the dielectric boundary force, without introducing ionic charge densities in the solvent. They further applied their model, termed dynamical implicit-solvent model (DISM), to a charged spherical molecule to derive a generalized Rayleigh–Plesset equation, a stochastic ordinary differential equation for the fluctuating radius. With the same deterministic model, Li et al. [13] study the linear stability of a cylindrical solute-solvent interface, and conclude that the viscosity can modify the critical wavelength of such stabilities. Luo et al. [14, 27] make a connection of solvent fluid mechanics model with statistical mechanics theory. They also develop numerical methods to solve the solvent fluid equations. In particular, they design boundary conditions on the boundary of a computational domain to allow solutes to change their volumes. They also implement an augmented immersed interface method for the Navier–Stokes flow with moving interface.

In this work, we develop numerical methods to solve the governing equations of the solvent fluid mechanics model. As the full model is quite complicated, we shall focus here on a two-dimensional setting for nonpolar molecular systems. Our main results are as follows:

  1. We discretize the Stokes equation using a ghost fluid finite difference scheme. We show that our scheme is second-order accurate;
  2. We refine the method proposed in [14, 27] to design artificial boundary conditions on the boundary of computational domain, allowing the change of solute volume;
  3. We couple the level-set method with our Stokes solver to track the motion of solute-solvent interface. Our method captures both dry and wet states for some simple model systems.

We remark that our ghost fluid scheme is quite different from the augmented immersed interface method used in [14, 27]. We hope that our method can be better coupled with some other compact schemes, such as the compact coupling interface method [29] for electrostatic interactions.

The rest of the paper is organized as follows: In Section 2, we present the solvent fluid model. In particular, we describe the boundary conditions for the Stokes equation. In Section 3, we describe our numerical schemes. In Section 4, four test examples are provided to show the accuracy of our numerical schemes. In Section 5, we compute the interface around two nonpolar spherical solute atoms and the interface around two nonpolar plates. In Section 6, we draw conclusions and discuss our future work.

2 A Solvent Fluid Model

We denote by Ω the region of an underlying salvation system. It is divided into the solute region Ω and the solvent region Ω+, separated by the solute-solvent interface Γ = [partial differential]Ω[partial differential]Ω+; cf. Figure 2.1. We assume all Ω, Ω, and Ω+ are open sets, and Ω is completely inside Ω, i.e., [partial differential]Ω[partial differential]Ω = [empty]. The interface Γ, the solute region Ω, and the solvent region Ω+ can all depend on t. At a given time t, we denote by u = (u, v) : Ω+R2 and p : Ω+R the velocity and pressure of the solvent fluid, respectively. We also denote by p : ΩR the pressure inside the solute region Ω.

Figure 2.1
The geometry of a salvation system. The solute-solvent interface Γ separates the solute region Ω from the solvent region Ω+. The unit normal and unit tangent vectors at Γ are denoted by n and τ, respectively. ...

Our governing equations and boundary conditions are as follows [13, 26]:

  • Interface motion:
    where a dot donates the time derivative.
  • The Stokes equation for incompressible flow:
  • The ideal-gas law:
  • Traction interface conditions for the fluid velocity and pressure:
  • Boundary conditions on [partial differential]Ω for the velocity and pressure:

Here, μ denotes the viscosity of solvent fluid and G = (g(1), g(2)) is the density of a body force exerted in Ω+. In the ideal-gas law (2.3), Cm is a constant, proportional to the temperature and the number of solute atoms. If there are multiple components Ω(i) (1 ≤ iI) of solutes then the ideal-gas law should be applied to each Ω(i) and the constants Cm(i) for Ω(i) can vary with i. As usual, we denote D(u) = [nabla]u + [nabla]uT to be the rate-of-strain tensor. The force f in the traction boundary condition is given by


where γ is the constant surface tension and H is the mean curvature. The surface vdW force fvdw is defined by [7, 8, 13, 29]


with each ULJ(i) the Lennard-Jones potential for atom i. We have implicitly assumed that there are I solute atoms inside the solute region Ω. The boundary velocity u0 will be specified later. The boundary pressure p is a given function.

Due to the incompressibility of the solvent fluid, usual boundary conditions on [partial differential]Ω, such as no-slip or periodic boundary conditions, will lead to a constant volume of the solute region Ω. In order to simulate molecular conformational changes and multiple hydration states, which often result the volume change of solutes, we need to design numerical boundary conditions. To do so, we use the consistency condition


This equates the fluid flux along [partial differential]Ω to that along Γ. It results from the incompressibility condition. The volume change of Ω_ is thus determined by the fluid flux along [partial differential]Ω. The calculation of such volume change depends on the way u0 is specified in (2.5). Here we design the boundary velocity u0 to achieve such volume change. Let Ω = (0, lx) × (0, ly) for some lx > 0 and ly > 0, we specify the boundary velocity profile in one of the following two ways:

  • A parabolic profile:
  • A circular profile:

In these equations, the function C(t) is to be determined by the consistency condition (2.6). The idea is that, under the assumption of the absence of external flow, the error of the boundary velocity prescribed by the profiles given in (2.7) and (2.8) would have relatively minor effect on the interface motion. This assumption is justified in Section 4 for a shrinking circle, whose analytical solution is well known.

Finally, to numerically track the motion of the solute-solvent interface Γ, we employ the level-set method. We denote by ϕ = ϕ(x, t) a level-set function of the interface Γ = Γ(t) at time t, i.e., Γ = {x [set membership] Ω : ϕ(x, t) = 0}. We also assume that Ω = {x [set membership] Ω : ϕ(x, t) < 0} and Ω+ = {x [set membership] Ω : ϕ(x, t) > 0}. The level-set equation is


where un = u · n is the normal component of the velocity field u at the interface Γ, but is suitably extended to Ω.

In summary, we couple the Stokes equation (2.2) in Ω+ with the traction interface conditions (2.4) and the ideal-gas law (2.3), the consistency condition (2.6), and the Dirichlet boundary condition (2.5) for p, and (2.7) or (2.8) for u. The fluid velocity dictates the motion of Γ, which is tracked by the level-set equation (2.9).

3 Numerical Methods

In this section, we introduce our numerical methods. We divide our computational domain Ω = (0, lx) × (0, ly) into nx × ny grid cells, with nx and ny two positive integers. We denote hx = lx/nx and hy = ly/ny, and define xi, j = ((i + 1/2)hx, (j + 1/2)hy), xi±1/2, j = ((i + 1/2 ± 1/2)hx, (j + 1/2)hy), and xi, j±1/2 = ((i + 1/2)hx, (j + 1/2 ± 1/2)hy).

3.1 Discretization of the Stokes Equation

We use a marker-and-cell (MAC) grid for the unknown fluid components u, v, and p. We approximate p at the center xi, j of each cell, u at the midpoints of vertical cell edges xi−1/2, j, and v at the midpoints of horizontal cell edges xi, j−1/2; cf. Figure 3.1. For convenience, we define the regular points, boundary points, and ghost points as follows:

  1. The regular velocity points are those points on the edges of the cells and are located inside the fluid region Ω+.
  2. The regular pressure points are those points on the center of the cells, of which at least one edge contains a regular velocity point.
  3. The boundary velocity points are those points on [partial differential]Ω, or on an edge that intersects [partial differential]Ω.
  4. The boundary pressure points are the center points xi, j of those boundary cells, each of which has at least one edge entirely on [partial differential]Ω.
  5. A ghost velocity point is a point located either inside Ω or on the interface Γ, and is a neighbor to a regular velocity point. Two velocity points are neighbors of each other if the edges they are on share a same vertex or if they are on the edges of the same cell.
Figure 3.1
A schematic MAC grid.

With the assumption that Γ is far away from [partial differential]Ω, we can simply assume that any ghost point is not a boundary point, and any boundary point is a fluid point. We then denote

Cu={upoints}={regular points foru}{boundary points foru}{ghost points foru},Cv={vpoints}={regular points forv}{boundary points forv}{ghost points forv},Cp={ppoints}={regular points forp}{boundary points forp}.

For any non-boundary fluid point for velocity xi−1/2, j [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms722846ig1.jpgu or xi, j−1/2 [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms722846ig1.jpgv, or any non-boundary regular point for pressure xi, j [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms722846ig1.jpgp, we apply the following second-order discretization scheme to the incompressible Stokes equation:




Note that this scheme enforces the divergence free condition on each cell with center a non-boundary regular pressure point. For any boundary velocity point or pressure point, we apply the Dirichlet boundary condition (2.7) or (2.8), and (2.5), respectively.

The boundary velocity profile on Γ(t) depends on a new variable C(t), which is determined by discretizing (2.6). We approximate the left-hand side of (2.6) by numerical integrating along the lines x = 0, lx and y = 0, ly:


To approximate the right-hand side of (2.6), we consider the region ω bounded by Γ and the lines x = hx, lxhx, and y = hy, lyhy, and use the incompressibility condition to get


Here the right-hand side of (2.6) can be approximated by


By substituting (2.7) or (2.8) into (3.3), and combining (2.6) and (3.2), we arrive at



ζ={16[(2nx3+nx)hx3+(2ny3+ny)hy3]for parabolic profile(2.7),i=1nx4lyhxly2+lx2(2i1nx)2+j=1ny4lxhylx2+ly2(2j1ny)2for circular profile(2.8).

For any ghost velocity point xi+1/2, j [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms722846ig1.jpgu or xi, j+1/2 [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms722846ig1.jpgv, we first project it onto Γ and find the projection point yi+1/2, j or yi, j+1/2, respectively. Then we use the nearby fluid velocity points, ghost velocity points, and pressure points to form a locally third-order approximation to the derivative components ux, uy, vx, vy, and pressure component p at yi+1/2, j and yi, j+1/2, respectively. Thus, we obtain a third-order approximation to the traction interface conditions at the projection points. This third-order approximation scheme on the interface covers various cases of the local geometry of interface. In some other cases, where the geometry of Γ does not allow enough points for a third order scheme, we switch to a second order or even first order scheme, which causes the solution to be reverted to a lower order approximation. We include a detailed discussion and formulation of the discretization scheme for ghost velocity points in the appendix.

By discretizing the system as illustrated above, we obtain a linear system, with an array of unknown composed of u, v at the fluid points, p at the pressure points, and C(t). The discretization matrix An external file that holds a picture, illustration, etc.
Object name is nihms722846ig2.jpg is asymmetric and sparse. It takes the following form


where the last row and last column correspond to (3.4) and (3.5), respectively. Rows of Au, Av correspond mainly to the five-point stencils for Laplacian. All entries of Ou, Ov, and Op are zero, except a few that correspond to ghost points. All Apx, Apy, Aux, and Avy result mainly from the discretization of px, py, ux, and vy, respectively. In each of these submatrices, there are rows corresponding to the ghost points, where the discretization of the (2.4) couples the u, v, and p points and introduces many nonzero off-diagonal entries. This linear problem is similar in structure to a saddle point problem. It may suggests a solution method involving the Schur complement reduction. However, the highly coupled interface conditions together with the geometry dependent discretization make the conditional number of An external file that holds a picture, illustration, etc.
Object name is nihms722846ig2.jpg very high. Furthermore, submatrices such as Au, Av are far away from diagonal dominant. At a few ghost points, the diagonal entries have their magnitude smaller than off-diagonals. As a consequence, it is not efficient to apply any Krylov subspace solver to even a subproblem. Leaving the development of numerical algorithms to future work, as a first step, we use UMFPACK, an implementation of a direct multifrontal sparse LU factorization method, to solve this system [6]. The time cost for UMFPACK to solve a sparse linear system of size ~ 30,000-by-30,000 is O(0.1) seconds.

3.2 Solving the Level-Set Equation

We now consider the discretization of the level-set equation (2.9). For the time derivative, we use the explicit forward Euler scheme


where ϕ(k)(x) and un(k) (x) are the approximations of ϕ(x, tk) and un(x, tk), respectively, at time tk = kΔt (k = 1, 2,…) and Δt is the time step satisfying the CFL condition


If un(x, t) = 0 for all x on Γ at a certain time t, then Γ no longer moves, and a steady state interface is reached.

To approximate the normal velocity un(x), we use the level-set function, to interpolate the points on the interface, calculate the normal fluid velocity on these points, and use a fast marching algorithm to extend the values in the normal direction to all spatial grid points. Note that we here use the same grid as that for the discretization of the pressure p. For spatial derivatives, we use a fifth-order WENO method to approximate |[nabla]ϕ(k)(x)|. We use the homogeneous Neumann boundary conditions at the outer boundary of the computational box.

To keep the level-set function as a signed distance function, we reinitialize ϕ by solving the equation


for a few steps with the initial value ϕ = ϕ0 at t = 0. Here ϕ0 is the level-set function before reinitialization, and the time t is different from that in the original level-set equation.

4 Numerical Tests

4.1 Flow Outside a Circular Region

In this example, we test the convergence of our Stokes solver on flow outside a circular region. We set Ω = (0, 1) × (0, 1) and Ω = {(x, y) [set membership] Ω : ϕ(x, y) < 0}, where the level-set function is given by


Note that Ω is a circular region. Note also that the region Ω and the interface Γ are fixed all the time. So there is no moving interface in this test. Furthermore, we set μ = 1 and fix p =0. We select the boundary velocity u0, the boundary pressure p, the fluid body force G, and the interface force f to yield the following velocity and pressure fields for an incompressible flow:


We then solve the Stokes equation (2.2) on Ω+ = {(x, y) [set membership] Ω : ϕ(x, y) > 0} with such boundary conditions. The numerical solutions of u, v, and p to this test example are plotted in Figure 4.1 on an N × N spatial grid with N = 400. We analyze the error between the numerical solutions and the analytical solutions (4.2), by generating in Figure 4.2 six log-log plots of the L2-norm and L-norm of the error for u, v, and p versus N, the number of grid points in both x and y directions. The log-log plots show that the error for u, v, and p all decay in an order of O(N−2) in both L2-norm and L-norm, with spikes intermittently. These spikes arise due to the insufficient grid resolution for the curved interfacial geometry resulting unpredicted sudden increase of the conditional number of discretization matrix. In average, the conditional number increases in the order of O(N3). We believe that such spikes can be reduced by discretizing the interface condition using a least-squares approach.

Figure 4.1
Numerical solution to the Stokes equation (2.2). From left to right: u, v component of fluid velocity and the fluid pressure p.
Figure 4.2
Log-log plots of errors vs. the number of discretization (grid points). Top row from left to right: The L2-norm of the error for u, v and p, with the red straight line having slope −2. Bottom row from left to right: The L-norm of the ...

4.2 Flow Outside a Clover Shaped Region

This example is the same as that in Subsection 4.1, except that Ω = (1/2, 1/2) + Ωc, where Ωc is a three-fold clover region defined in polar coordinates as


Our numerical solutions of u, v, and p are plotted in Figure 4.3, with the spatial discretization parameter N = 400 in both x and y directions. An analysis of the error between the numerical and analytical solutions shows similar behavior than the previous example in terms of L2-norm and L-norm of the error on the three fluid components u, v, and p, as shown in Figure 4.4. We get an average second-order convergence in u, v, and p.

Figure 4.3
Numerical solution to the Stokes equation (2.2) described in Subsection 4.2. From left to right: u, v component of fluid velocity, and the fluid pressure p.
Figure 4.4
Top row from left to right: The L2-norm of the error for u, v and p, with the red straight line being of slope −2. Bottom row from left to right: The L-norm of the error for u, v and p, with the red straight line being of slope −2. ...

4.3 Flow around a Disk with Designed Numerical Boundary Conditions

We now test on our numerical boundary conditions. Define


Define also the circular region Ω = {(x, y) [set membership] R2 : ϕ(x, y) < 0}. One can verify that


solve the Stokes equation (2.2) with μ = 1, G = 0, together with the interface condition p = 0, and f = − n, and also that its flux along Γ is


We solve the Stokes equation numerically for u and p on Ω = (0, 1) × (0, 1) with p = 0, and the boundary velocity profile given by (2.7) or (2.8). In Figure 4.5, we plot the solution with profile (2.7) and nx = ny = N = 101. The tendency of shrinking of the volume of Ω can be observed from the inward velocity field.

Figure 4.5
Numerical solution to the Stokes equation (2.2) and the corresponding flux. Left: The quiver velocity field u in the whole domain Ω. Right: A zoom-in velocity field u along Γ.

As we increase N, the flux ζC(t) approaches a constant. In Figure 4.6, we plot the flux ζC(t) as a function of N. We see that for both (2.7) and (2.8), ζC(t) converges to a certain value. For (2.8), ζC(t) approaches the analytical value 0.04π, because the circular velocity profile is exact in this case. For (2.7), ζC(t) approaches to a value which is slightly less than, but reasonably close to 0.04π. This shows that even though the assumptions of boundary velocity profiles are artificial, they can be expected to work relatively accurately. To test the convergence of ζC(t) with respect to increasing N, we take the value of ζC(t) at N = 600 as a reference value for both (2.7) and (2.8). We then compute the absolute values of the differences between the fluxes at the other values of N and the reference values. This process gives us the convergence plot in Figure 4.6, from which we observe a convergence rate of roughly first order for (2.7) and second order for (2.8).

Figure 4.6
Left: The total fluid flux ζC(t) versus N. The blue horizontal line is 0.04π. The analytic value of ζC(t) : the black and red lines and dots correspond to the solutions with profile (2.7) and profile (2.8), respectively. Right: ...

4.4 Moving Interface Driven by Solvent Fluid Flow with Curvature Force

We now test our level-set method coupled with our Stokes solver on a moving circular interface. Let us define


where r(t) = 0.2 − t/2 with 0 ≤ t < 0.4. It is easy to verify that p and u solve the Stokes equation (2.2) with μ = 1 and G = 0 in the region R2/Ω(t)¯, where


Moreover, p and u satisfy the interface condition (2.4) on Γ(t) = [partial differential]Ω(t) with p = 0 and f = −n/r(t) (i.e., the curvature flow). The time needed for the circle Γ(t) to shrink to the center (0.5, 0.5) is t = 0.4.

We now set Ω = (0, 1) × (0, 1), and the initial interface to be


We use our level-set method and ghost fluid Stokes solver to solve numerically (2.1) and (2.2) with G = 0, (2.4) with p = 0 and f = −Hn, where H is the curvature, and the boundary conditions p = 0 and boundary profile (2.7) and (2.8). Our numerical solution for Γ(t) as plotted in Figure 4.7 exhibits circular shrinking with constant speed. The center is slightly shifted due to small numerical error. The numerical critical time tc for both boundary profiles (2.7) and (2.8) agree well with the analytic value.

Figure 4.7
Left: The solution to the curvature flow problem with parabolic profile (2.7), at time t = 0, 0.1, 0.2, 0.3, 0.4, from outside to inside, respectively. Right: The numerical radius r(t) versus t with parabolic profile (2.7) (black curve) and with circular ...

5 Applications

In this section, we present two examples of application of our method to two model nonpolar molecular systems. One is a two-particle system and the other is a two plate system.

5.1 A Two-Particle System

We consider Ω = (0, lx) × (0, ly) with lx and ly two positive numbers. We assume two particles located at x1 = (lx/2 – c, ly/2) and x2 = (lx/2 + c, ly/2), respectively, where c > 0 is an adjustable parameter. We define the initial regions Ω and Ω+ using this level-set function


We set the body force G = 0. We also combine the curvature force and Lenard-Jones force into the surface force


where γ, ε, and σ are all constants, and take the following values for this example: γ = 1, ε = 25, σ = 0.2. Furthermore, we take Cm = 0.001 in (2.3) and (2.4).

We then solve the Stokes equation coupled with level-set equation (2.9), with ly = 1 and various choices of lx and c : (1) lx = 1, c = 0.05; (2) lx = 1, c = 0.2; (3) lx = 2, c = 0.25. The final steady state for Γ is plotted in Figure 5.1. We see that the interface breaks into two parts as the distance between the two particles is large enough; cf. case (3).

Figure 5.1
The initial interface (red curve) and the numerically computed steady-state interface (black curve). The two blue dots are the positions of two particles. Left: Ω = (0, 1) × (0, 1), x1 = (0.45, 0.5), and x2 = (0.55, 0.5). Middle: Ω ...

5.2 A Two-Plate System

In this example, we set Ω = (0, 1) × (0, 1) and two plates composed of six atoms each. These atoms are located at coordinates (0.5 – c, 0.25 + 0.1k) and (0.5 + c, 0.25 + 0.1k) with k [set membership] {0, 1, (...), 5} and c an adjustable constant. With certain parameter choices, one observe dry/wet polymodal states depending on the initial solute region Ω. We choose c = 0.15, γ = 5, ε = 25, σ = 0.1, Cm = 0.001 in (5.2). Moreover, for a loose initial condition of Γ, we define


For a tight initial condition on Γ, we define


We solve the Stokes equation (2.2) coupled with level-set equation (2.9). The interface Γ goes to a dry final state from the loose initial condition, and a wet final state from the tight initial condition; cf. Figure 5.2. The evolution of Γ in between the two states captures the dynamics of the transitioning process, which cannot be obtained by an energy variational approach such as VISM. In our future work, we would like to develop a 3D fluid solver, and compare our numerical results with the results obtained from molecular dynamics simulations.

Figure 5.2
The numerical solutions of the interface Γ: an initial (red curve), an intermediate (blue curve), and a final, steady state (black curve) interfaces. The blue dots are the atom positions. Both plots are with γ = 5 and σ = 0.1. ...

6 Conclusions

In this paper, we model the aqueous solvent by the incompressible Stokes flow and treat the solute with the ideal-gas law. The solute-solvent interface moves with the solvent fluid flow. All the viscous, pressure, surface tension, and the solute-solvent vdW forces are balanced on the interface, leading to a traction boundary condition. To allow the change of solute volume, we design special numerical boundary conditions on the boundary of our computation domain. We design a second-order ghost fluid method for solving the Stokes equation. We also couple the level-set equation for the moving solute-solvent interface. Our methods accurately predict the blowup time of a shrinking bubble under surface tension. Moreover, our methods capture the dry and wet polymodal interfaces for the two-plate system.

Some existing issues of this model include: (1) The boundary conditions assume a velocity profile, which may not be realistic; (2) The numerical error is sensitive to the location of the interface; (3) We solve the linear system using a direct LU factorization, which can be problematic in extension to three dimensions.

As a first step in the future, we need to propose more realistic boundary conditions. We also plan to decrease the sensitivity of error dependence on the interface by considering a least square problem. Moreover, an extension of our two dimensional fluid solver to three dimensions is necessary. After we develop such a fluid solver, we can add noise to the solvent and observe the fluctuations of the interface. Through a combination of such a three dimensional fluid solver with fluctuations and a robust Poisson–Boltzmann solver, one can better observe and describe the dynamics of a solvent–solute system.


This work was supported by the US National Science Foundation (NSF) through grant DMS-1319731 and the US National Institutes of Health (NIH) through grant R01GM096188. Work in McCammon's group is supported in part by NSF, NIH, HHMI, and NBCR. The authors thank Dr. Robert Krasny, Dr. Ray Luo, and Mr. Li Xiao for helpful discussions.


In the appendix, we provide details of the ghost fluid discretization on the interface. First of all, with the notations n = (n1, n2) and τ = (–n1,n2), the traction boundary conditions (2.4) read


where f[perpendicular] = f · n and f= f · τ. Some straightforward algebraic calculations together with the incompressibility condition (2.2) lead to


For any ghost velocity point x, we find a point x* [set membership] Γ, such that |xx*| = dist(x, Γ). We call x* a projection point of x onto Γ. We then discretize equation (A.1) at each projection point xi1/2,j corresponding to each ghost velocity point xi−1/2,j of u. Similarly, we discretize equation (A.2) at each projection point xi,j1/2 corresponding to each ghost velocity point xi, j−1/2 of v.

To obtain a second-order convergence scheme for u, v, and p up to Γ, we design a third-order discretization of ux, uy, vx, vy, and p in the fluid region. This requires 10 stencil points for u, v, and 6 stencil points for p. We denote by S(u, x*, r), S(v, x*, r), and S(p, x*, r) respectively, the sets of u, v and p stencil points for discretizing [nabla]u, [nabla]v and p with an order r at the projection point xi,j1/2. Then


In choosing these stencil points, we follow three criteria: (1) Each of these stencil points needs to be either a ghost point or a fluid point; (2) The stencil points need to include x, that is x [set membership] S(u, x*, r) [union or logical sum] S(v, x*,r)[union or logical sum]S(p, x*,r); (3) All S(u, x*, r), S(v, x*, r), and S(p, x*, r) satisfy


The criterion (3) is important for the invertibility of matrix A in equation (A.4), as shall be discussed later.

We now describe the process of constructing S(u, x*, r), S(v, x*, r), and S(p, x*, r), and the corresponding schemes. In this process, we use the following notations


It is easy to see that |λ1| < 1 and |λ2| < 1. Since a ghost point is a neighbor to a fluid point, at least one of its neighbor needs to be in Ω+. We name a neighbor of x a check point, if that neighbor point is in Ω+. There are totally six different cases for the combination of ghost point and check point: (a) x = xi−1/2,j and ϕ(xi−1/2+s1,j) > 0; (b) x = xi−1/2j and ϕ(xi−1/2j+s2) > 0; (c) x = xi−1/2,j and ϕ(xi−1/2+s1/2, j+s2/2) > 0; (d) x = xi, j−1/2 and ϕ(xi+s1,j−1/2) > 0; (e) x = xi, j−1/2 and ϕ(xi, j−1/2+s2) > 0; (f) x = xi, j−1/2 and ϕ(xi+s1/2,j−1/2+s2/2) > 0. For each of these cases, we obtain the corresponding S(u, x*, r), S(v, x*, r), and S(p, x*, r) with r [set membership] {1, 2, 3}. In Figure A.1, we schematically plot all six cases and the corresponding S(u, x*, r), S(v, x*, r), and S(p, x*, r) for r = 3. Notice that for cases (a)–(c), equation (A.1) is discretized, and S(v, x*, r) is not needed.

Figure A.1
Different cases of combination of a ghost point and a check point.

We now describe the steps of constructing a third-order discretization scheme for case (a), whereas all other cases just follow tediously. We consider a u ghost point xi−1/2j with its projection point x*. Let us denote


and F = (u, ux, uy, uxx, uxy, uyy, uxxx, uxyy, uxxy, uyyy) T, where (·)T denotes transpose. Then a third-order Taylor expansion leads to a linear system X = AF(x*), where


with (hx,k, hy,k) := xkx*. A sufficient condition to have A invertible is (A.3). Inverting A symbolically, we get F(x*) = A−1X, and the second entry of F(x*) corresponds to a third-order discretization of ux(x*):


where the notation u(X) = (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9)T is used. The third entry of F(x*)t gives a third-order discretization of uy(x*). However, it is not necessary in discretizing equation (A.1).

Similarly, one can derive a third-order discretization scheme for p(x*):


Where p0=p(xis1+12,j+s2), p1=p(xis112,j+s2), p2=p(xi+3s112,j+s2), p3=p(xi+s112,j+2s2), p4=p(xi+3s112,j+2s2), p5=p(xi+3s112,j+3s2).

We use this symbolic inverse matrix method to find a third-order scheme for discretizing equation (A.1) around x*. To find a second-order or first-order scheme, we form a smaller matrix A by applying a second or first-order Taylor expansion to u(xk), v(yk), and p(zk) at x*, with xk [set membership] S(u, x*, r), yk [set membership] S(v, x*, r), and zk [set membership] S(p, x*, r), respectively. Using similar computations, we are able to construct the first, second, and third-order schemes for all six cases (a)–(f), but we do not need to record the ten-page formulas here.


1. Alexander-Katz A, Schneider MF, Schneider SW, Wixforth A, Netz RR. Shear–flow–induced unfolding of polymeric globules. Phys Rev Lett. 2006;97:138101. [PubMed]
2. Baron R, McCammon JA. Molecular recognition and ligand association. Annu Rev Phys Chem. 2013;64:151–175. [PubMed]
3. Chen J, Brooks CL, III, Khandogin J. Recent advances in implicit solvent based methods for biomolecular simulations. Curr Opin Struct Biol. 2008;18:140–148. [PMC free article] [PubMed]
4. Cheng L, Dzubiella J, McCammon JA, Li B. Application of the level–set method to the implicit solvation of nonpolar molecules. J Chem Phys. 2007;127:084503. [PubMed]
5. Cheng L, Li B, Wang Z. Level–set minimization of potential controlled hadwiger valuations for molecular solvation. J Comput Phys. 2010;229:8497–8510. [PMC free article] [PubMed]
6. Davis TA. Algorithm 832: UMFPACK v4.3–an unsymmetric-pattern multifrontal method. ACM Trans Math Softw. 2004;30:196–199.
7. Dzubiella J, Swanson J, McCammon J. Coupling hydrophobicity, dispersion, and electrostatics in continuum solvent models. Phys Rev Lett. 2006;96:087802. [PubMed]
8. Dzubiella J, Swanson J, McCammon J. Coupling nonpolar and polar solvation free energies in implicit solvent models. J Chem Phys. 2006;124:084905. [PubMed]
9. Guo Z, Li B, Dzubiella J, Cheng LT, McCammon JA, Che J. Heterogeneous hydration of p53/MDM2 complex. J Chem Theory Comput. 2014;10:1302–1313. [PMC free article] [PubMed]
10. Hagen SJ. Solvent viscosity and friction in protein folding dynamics. Curr Protein Pept Sci. 2010;11:385–395. [PubMed]
11. Klimov DK, Thirumalai D. Viscosity dependence of folding rates of protein. Phys Rev Lett. 1997;79:317–320.
12. Levy Y, Onuchic JN. Water mediation in protein folding and molecular recognition. Annu Rev Biophys Biomol Struct. 2006;35:389–415. [PubMed]
13. Li B, Sun H, Zhou S. Stability of a cylindrical solute–solvent interface: Effect of geometry, electrostatics, and hydrodynamics. SIAM J Appl Math. 2015;75:907–928. [PMC free article] [PubMed]
14. Li Z, Cai Q, Zhao H, Luo R. A semi–implicit augmented IIM for Navier–Stokes equations with open and traction boundary conditions. J Comput Phys. 2015;297:182–193. [PMC free article] [PubMed]
15. Roux B, Simonson T. Implicit solvent models. Biophys Chem. 1999;78:1–20. [PubMed]
16. Schneider SW, Nuschele S, Wixforth A, Gorzelanny C, Alexander-Katz A, Netz RR, Schneider MF. Shear–induced unfolding triggers adhesion of von Willebrand factor fibers. Proc Natl Acad Sci USA. 2007;104:7899–7903. [PubMed]
17. Sekhar A, Latham MP, Vallurupalli P, Kay LE. Viscosity-dependent kinetics of protein conformational exchange: Microviscosity effects and the need for a small viscogen. J Phys Chem B. 2014;118:4546–4551. [PubMed]
18. Setny P, Wang Z, Cheng LT, Li B, McCammon JA, Dzubiella J. Dewetting–controlled binding of ligands to hydrophobic pockets. Phys Rev Lett. 2009;103:187801. [PMC free article] [PubMed]
19. Siedlecki CA, Lestini BJ, Kottke-Marchant KK, Eppell SJ, Wilson DL, Marchant RE. Shear–dependent changes in the three–dimensional structure of human von Willebrand factor. Blood. 1996;88:2939–2950. [PubMed]
20. Singh I, Themistou E, Porcar L, Neelamegham S. Fluid shear induces conformation change in human blood protein von Willebrand factor in solution. Biophys J. 2009;96:2313–2320. [PubMed]
21. Szymczak P, Cieplak M. Hydrodynamic effects in proteins. J Phys: Condens Matter. 2011;23:033102. [PubMed]
22. Tanford C. The Hydrophobic Effect: Formation of Micelles and Biological Membranes. John Wiley; New York: 1973.
23. Tomasi J, Persico M. Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent. Chem Rev. 1994;94:2027–2094.
24. Vergauwe RMA, Uji–i H, De Ceunynck K, Vermant J, Vanhoorelbeke K, Hofkens J. Shear–stress–induced conformational changes of von Willebrand factor in water-glycerol mixture observed with single molecule microscopy. J Phys Chem B. 2014;118:5660–5669. [PubMed]
25. Wang Z, Che J, Cheng L, Dzubiella J, Li B, McCammon JA. Level–set variational implicit-solvent modeling of biomolecules with the coulomb–field approximation. J Chem Theory Comput. 2012;8:386–397. [PMC free article] [PubMed]
26. White M. PhD thesis. University of California; San Diego: 2015. Mathematical Theory and Numerical Methods for Biomolecular Modeling.
27. Xiao L, Cai Q, Li Z, Zhao H, Luo R. A multi–scale method for dynamics simulation in continuum solvents I: Finite-difference algorithm for Navier–Stokes equation. Chem Phys Lett. 2014;616:67–74. [PMC free article] [PubMed]
28. Zhou S, Cheng L, Sun H, Che J, Dzubiella J, Li B, McCammon JA. LS–VISM: A software package for analysis of biomolecular solvation. J Comput Chem. 2015;36:1047–1059. [PMC free article] [PubMed]
29. Zhou S, Cheng LT, Dzubiella J, Li B, McCammon JA. Variational implicit solvation with Poisson–Boltzmann theory. J Chem Theory Comput. 2014;10:1454–1467. [PMC free article] [PubMed]