Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC4922513

Formats

Article sections

- Abstract
- 1 Introduction
- 2 A Solvent Fluid Model
- 3 Numerical Methods
- 4 Numerical Tests
- 5 Applications
- 6 Conclusions
- References

Authors

Related links

J Sci Comput. Author manuscript; available in PMC 2017 May 1.

Published in final edited form as:

Published online 2015 September 12. doi: 10.1007/s10915-015-0099-z

PMCID: PMC4922513

NIHMSID: NIHMS722846

Hui Sun: ude.dscu@300suh; Shenggao Zhou: nc.ude.adus@uohzgs; David K. Moore: ude.dscu@eroomkd; Li-Tien Cheng: ude.dscu.htam@gnehcl; Bo Li: ude.dscu.htam@ilb

The publisher's final edited version of this article is available at J Sci Comput

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.

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, 19–21, 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:

- We discretize the Stokes equation using a ghost fluid finite difference scheme. We show that our scheme is second-order accurate;
- 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.

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 Γ = Ω_{−} ∩ Ω_{+}; cf. Figure 2.1. We assume all Ω, Ω_{−}, and Ω_{+} are open sets, and Ω_{−} is completely inside Ω, i.e., Ω_{−} ∩ Ω = . 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*) : Ω_{+} → ^{2} and *p* : Ω_{+} → the velocity and pressure of the solvent fluid, respectively. We also denote by *p*_{−} : Ω_{−} → *R* the pressure inside the solute region Ω_{−}.

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:$$\dot{\mathbf{x}}=\mathbf{u}(\mathbf{x}(t))\phantom{\rule{0.2em}{0ex}}\text{for}\phantom{\rule{0.2em}{0ex}}\mathbf{x}=\mathbf{x}(t)\in \mathrm{\Gamma},$$(2.1)where a dot donates the time derivative.
- The Stokes equation for incompressible flow:$$\{\begin{array}{c}\mu \mathrm{\Delta}\mathbf{u}-\nabla p+\mathbf{G}=\mathbf{0}\phantom{\rule{0.2em}{0ex}}\text{in}\phantom{\rule{0.2em}{0ex}}{\mathrm{\Omega}}_{+},\\ \nabla \xb7\mathbf{u}=0\phantom{\rule{0.2em}{0ex}}\text{in}\phantom{\rule{0.2em}{0ex}}{\mathrm{\Omega}}_{+}.\end{array}$$(2.2)
- The ideal-gas law:$${p}_{-}\text{Vol}({\mathrm{\Omega}}_{-})={C}_{m}.$$(2.3)
- Traction interface conditions for the fluid velocity and pressure:$$\{\begin{array}{c}\mu \mathbf{n}\xb7D(\mathbf{u})\mathbf{n}-p+{p}_{-}=-\mathbf{f}\xb7\mathbf{n}\phantom{\rule{0.2em}{0ex}}\text{for}\phantom{\rule{0.2em}{0ex}}\mathbf{x}\in \mathrm{\Gamma},\\ \mu \mathit{\tau}\xb7D(\mathbf{u})\mathbf{n}=-\mathbf{f}\xb7\mathit{\tau}\phantom{\rule{0.2em}{0ex}}\text{for}\phantom{\rule{0.2em}{0ex}}\mathbf{x}\in \mathrm{\Gamma}.\end{array}$$(2.4)
- Boundary conditions on Ω for the velocity and pressure:$$\mathbf{u}={\mathbf{u}}_{0}\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}p={p}_{\infty}\phantom{\rule{0.2em}{0ex}}\text{on}\phantom{\rule{0.2em}{0ex}}\partial \mathrm{\Omega}.$$(2.5)

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), *C _{m}* is a constant, proportional to the temperature and the number of solute atoms. If there are multiple components
${\mathrm{\Omega}}_{-}^{(i)}$ (1 ≤

$$\mathbf{f}(\mathbf{x})=-\gamma H(\mathbf{x})\mathbf{n}+{\mathbf{f}}_{\text{vdW}}(\mathbf{x})\phantom{\rule{0.2em}{0ex}}\text{for}\phantom{\rule{0.2em}{0ex}}\mathbf{x}\in \mathrm{\Gamma},$$

where *γ* is the constant surface tension and *H* is the mean curvature. The surface vdW force **f**_{vdw} is defined by [7, 8, 13, 29]

$${\mathbf{f}}_{\text{vdW}}(\mathbf{x})=(\sum _{i=1}^{I}{U}_{LJ}^{(i)}(\mathbf{x}))\mathbf{n},$$

with each
${U}_{LJ}^{(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 **u**_{0} will be specified later. The boundary pressure *p _{∞}* is a given function.

Due to the incompressibility of the solvent fluid, usual boundary conditions on Ω, 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

$${\int}_{\partial \mathrm{\Omega}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l={\int}_{\mathrm{\Gamma}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l.$$

(2.6)

This equates the fluid flux along Ω to that along Γ. It results from the incompressibility condition. The volume change of Ω_ is thus determined by the fluid flux along Ω. The calculation of such volume change depends on the way u_{0} is specified in (2.5). Here we design the boundary velocity **u**_{0} to achieve such volume change. Let Ω = (0, *l _{x}*) × (0,

- A parabolic profile:$$\{\begin{array}{cc}u(0,y)=C(t)({l}_{y}-y)y,& v(0,y)=0,\\ u({l}_{x},y)=-C(t)({l}_{y}-y)y,& v({l}_{x},y)=0,\\ u(x,0)=0,& v(x,0)=C(t)({l}_{x}-x)x,\\ u(x,{l}_{y})=0,& v(x,{l}_{y})=-C(t)({l}_{x}-x)x.\end{array}$$(2.7)
- A circular profile:$$\{\begin{array}{cc}u(0,y)=C(t)\frac{2{l}_{x}}{{l}_{x}^{2}+{(2y-{l}_{y})}^{2}},& v(0,y)=C(t)\frac{2{l}_{y}-4y}{{l}_{x}^{2}+{(2y-{l}_{y})}^{2}},\\ u({l}_{x},y)=-C(t)\frac{2{l}_{x}}{{l}_{x}^{2}+{(2y-{l}_{y})}^{2}},& v({l}_{x},y)=C(t)\frac{2{l}_{y}-4y}{{l}_{x}^{2}+{(2y-{l}_{y})}^{2}},\\ u(x,0)=C(t)\frac{4x-2{l}_{x}}{{(2x-{l}_{x})}^{2}+{l}_{y}^{2}},& v(x,0)=C(t)\frac{2{l}_{y}}{{(2x-{l}_{x})}^{2}+{l}_{y}^{2}},\\ u(x,{l}_{y})=C(t)\frac{4x-2{l}_{x}}{{(2x-{l}_{x})}^{2}+{l}_{y}^{2}},& v(x,{l}_{y})=-C(t)\frac{2{l}_{y}}{{(2x-{l}_{x})}^{2}+{l}_{y}^{2}}.\end{array}$$(2.8)

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** Ω : *ϕ*(**x**, *t*) = 0}. We also assume that Ω_{−} = {**x** Ω : *ϕ*(**x**, *t*) < 0} and Ω_{+} = {**x** Ω : *ϕ*(**x**, *t*) > 0}. The level-set equation is

$${\varphi}_{t}+{u}_{n}|\nabla \varphi |=0\phantom{\rule{0.2em}{0ex}}\text{in}\phantom{\rule{0.2em}{0ex}}\mathrm{\Omega},$$

(2.9)

where *u _{n}* =

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).

In this section, we introduce our numerical methods. We divide our computational domain Ω = (0, *l _{x}*) × (

We use a marker-and-cell (MAC) grid for the unknown fluid components *u, v*, and *p*. We approximate *p* at the center **x*** _{i, j}* of each cell,

- The regular velocity points are those points on the edges of the cells and are located inside the fluid region Ω
_{+}. - The regular pressure points are those points on the center of the cells, of which at least one edge contains a regular velocity point.
- The boundary velocity points are those points on Ω, or on an edge that intersects Ω.
- The boundary pressure points are the center points
**x**of those boundary cells, each of which has at least one edge entirely on Ω._{i, j} - 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.

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

$$\begin{array}{c}{\mathcal{C}}_{u}=\{u\u2010\text{points}\}=\{\text{regular points for}\phantom{\rule{0.2em}{0ex}}u\}\cup \{\text{boundary points for}\phantom{\rule{0.2em}{0ex}}u\}\cup \{\text{ghost points for}\phantom{\rule{0.2em}{0ex}}u\},\\ {\mathcal{C}}_{v}=\{v\u2010\text{points}\}=\{\text{regular points for}\phantom{\rule{0.2em}{0ex}}v\}\cup \{\text{boundary points for}\phantom{\rule{0.2em}{0ex}}v\}\cup \{\text{ghost points for}\phantom{\rule{0.2em}{0ex}}v\},\\ {\mathcal{C}}_{p}=\{p\u2010\text{points}\}=\{\text{regular points for}\phantom{\rule{0.2em}{0ex}}p\}\cup \{\text{boundary points for}\phantom{\rule{0.2em}{0ex}}p\}.\end{array}$$

For any non-boundary fluid point for velocity **x**_{i}_{−1/2, j}
_{u} or **x**_{i, j}_{−1/2}
* _{v}*, or any non-boundary regular point for pressure

$$\{\begin{array}{c}\mu ({D}_{+}^{x}{D}_{-}^{x}{u}_{i-1/2,j}+{D}_{+}^{y}{D}_{-}^{y}{u}_{i-1/2,j})+{D}_{-}^{x}{p}_{i,j}+{g}_{i-1/2,j}^{(1)}=0,\\ \mu ({D}_{+}^{x}{D}_{-}^{x}{v}_{i,j-1/2}+{D}_{+}^{y}{D}_{-}^{y}{v}_{i,j-1/2})+{D}_{-}^{y}{p}_{i,j}+{g}_{i,j-1/2}^{(2)}=0,\\ {D}_{+}^{x}{u}_{i-1/2,j}+{D}_{+}^{y}{v}_{i,j-1/2}=0,\end{array}$$

(3.1)

Where

$$\{\begin{array}{cc}{D}_{+}^{x}{(\xb7)}_{i,j}=\frac{{(\xb7)}_{i+1,j}-{(\xb7)}_{i,j}}{{h}_{x}},& {D}_{+}^{y}{(\xb7)}_{i,j}=\frac{{(\xb7)}_{i,j+1}-{(\xb7)}_{i,j}}{{h}_{y}},\\ {D}_{-}^{x}{(\xb7)}_{i,j}=\frac{{(\xb7)}_{i,j}-{(\xb7)}_{i-1,j}}{{h}_{x}},& {D}_{-}^{y}{(\xb7)}_{i,j}=\frac{{(\xb7)}_{i,j}-{(\xb7)}_{i,j-1}}{{h}_{y}}.\end{array}$$

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, *l _{x}* and

$${\int}_{\partial \mathrm{\Omega}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l\approx \sum _{j=1}^{{n}_{y}}({u}_{1/2,j}+{u}_{{n}_{x}+1/2,j}){h}_{x}+\sum _{i=1}^{{n}_{x}}({v}_{i,1/2}+{v}_{i,{n}_{y}+1/2}){h}_{y}.$$

(3.2)

To approximate the right-hand side of (2.6), we consider the region *ω* bounded by Γ and the lines *x* = *h _{x}, l_{x}* −

$$0={\int}_{\omega}\nabla \xb7\mathbf{u}\mathrm{d}A={\int}_{\mathrm{\Gamma}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l+{\int}_{\partial \omega \backslash \mathrm{\Gamma}}\mathbf{u}\xb7\mathbf{n}dl.$$

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

$${\int}_{\mathrm{\Gamma}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l=-{\int}_{\partial \omega \backslash \mathrm{\Gamma}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l\approx \sum _{j=2}^{{n}_{y}-1}({u}_{3/2,j}+{u}_{{n}_{x}-1/2,j}){h}_{x}+\sum _{i=2}^{{n}_{x}-1}({v}_{i,3/2}+{v}_{i,{n}_{y}-1/2}){h}_{y}.$$

(3.3)

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

$$C(t)=\frac{1}{\zeta}[\sum _{j=2}^{{n}_{y}-1}({u}_{3/2,j}+{u}_{{n}_{x}-1/2,j}){h}_{x}+\sum _{i=2}^{{n}_{x}-1}({v}_{i,3/2}+{v}_{i,{n}_{y}-1/2}){h}_{y}],$$

(3.4)

Where

$$\zeta =\{\begin{array}{cc}\frac{1}{6}[(2{n}_{x}^{3}+{n}_{x}){h}_{x}^{3}+(2{n}_{y}^{3}+{n}_{y}){h}_{y}^{3}]& \text{for parabolic profile}\phantom{\rule{0.2em}{0ex}}(2.7),\\ \sum _{i=1}^{{n}_{x}}\frac{4{l}_{y}{h}_{x}}{{l}_{y}^{2}+{l}_{x}^{2}{(2i-1-{n}_{x})}^{2}}+\sum _{j=1}^{{n}_{y}}\frac{4{l}_{x}{h}_{y}}{{l}_{x}^{2}+{l}_{y}^{2}{(2j-1-{n}_{y})}^{2}}& \text{for circular profile}\phantom{\rule{0.2em}{0ex}}(2.8).\end{array}$$

(3.5)

For any ghost velocity point **x**_{i}_{+1/2}_{, j}* _{u}* or

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
is asymmetric and sparse. It takes the following form

$$\mathcal{A}=\left[\begin{array}{cccc}{A}_{u}& {O}_{u}& {A}_{{p}_{x}}& {\mathbf{b}}_{u}^{t}\\ {O}_{v}& {A}_{v}& {A}_{{p}_{y}}& {\mathbf{b}}_{v}^{t}\\ {A}_{{u}_{x}}& {A}_{{v}_{y}}& {O}_{p}& {\mathbf{0}}^{t}\\ {\mathbf{c}}_{u}& {\mathbf{c}}_{v}& \mathbf{0}& 1\end{array}\right],$$

(3.6)

where the last row and last column correspond to (3.4) and (3.5), respectively. Rows of *A _{u}, A_{v}* correspond mainly to the five-point stencils for Laplacian. All entries of

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

$$\frac{{\varphi}^{(k+1)}(\mathbf{x})-{\varphi}^{(k)}(\mathbf{x})}{\mathrm{\Delta}t}+{u}_{n}^{(k)}(\mathbf{x})|\nabla {\varphi}^{(k)}(\mathbf{x})|=0,$$

(3.7)

where *ϕ*^{(k)}(**x**) and
${u}_{n}^{(k)}$ (**x**) are the approximations of *ϕ*(**x**, *t _{k}*) and

$$\mathrm{\Delta}t=\frac{0.25h}{{max}_{\mathbf{x}}|{u}_{n}(\mathbf{x},{t}_{k})|}.$$

(3.8)

If *u _{n}*(

To approximate the normal velocity *u _{n}*(

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

$${\varphi}_{t}+\text{sign}({\varphi}_{0})(|\nabla \varphi |-1)=0$$

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.

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*) Ω : *ϕ*(*x, y*) < 0}, where the level-set function is given by

$$\varphi (x,y)=\sqrt{{(x-0.5)}^{2}+{(y-0.5)}^{2}}-0.13.$$

(4.1)

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 **u**_{0}, 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:

$$\{\begin{array}{c}u(x,y)=sin(2\pi x)cos(2\pi y),\\ v(x,y)=-cos(2\pi x)sin(2\pi y),\\ p(x,y)=-cos2\pi xcos(2\pi y).\end{array}$$

(4.2)

We then solve the Stokes equation (2.2) on Ω_{+} = {(*x, y*) Ω : *ϕ*(*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 *L*^{2}-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 *L*^{2}-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*(*N ^{3}*). We believe that such spikes can be reduced by discretizing the interface condition using a least-squares approach.

Numerical solution to the Stokes equation (2.2). From left to right: *u, v* component of fluid velocity and the fluid pressure *p*.

This example is the same as that in Subsection 4.1, except that Ω_{−} = (1/2, 1/2) + Ω* _{c}*, where Ω

$${\mathrm{\Omega}}_{c}=\{(r,\theta ):r<\frac{1}{4}+\frac{1}{16}cos(3\theta )\}.$$

(4.3)

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 *L*^{2}-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*.

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*.

We now test on our numerical boundary conditions. Define

$$\varphi =\sqrt{{(x-0.5)}^{2}+{(y-0.5)}^{2}}-0.2.$$

(4.4)

Define also the circular region Ω_{−} = {(*x, y*) ^{2} : *ϕ*(*x, y*) < 0}. One can verify that

$$p(\mathbf{x})=0\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}\mathbf{u}(\mathbf{x})=-\frac{0.02(x-0.5,y-0.5)}{{(x-0.5)}^{2}+{(y-0.5)}^{2}}$$

(4.5)

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

$$-{\int}_{\mathrm{\Gamma}}\mathbf{u}\xb7\mathbf{n}\mathrm{d}l=0.04\pi .$$

(4.6)

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

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).

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

$$p(\mathbf{x})=0\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}\mathbf{u}(\mathbf{x})=-\frac{r(t)(2x-1,2y-1)}{{(2x-1)}^{2}+{(2y-1)}^{2}},$$

(4.7)

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
${\mathbb{R}}^{2}/\overline{{\mathrm{\Omega}}_{-}(t)}$, where

$${\mathrm{\Omega}}_{-}(t)=\{(x,y)\in {\mathbb{R}}^{2}:{(x-0.5)}^{2}+{(y-0.5)}^{2}<r{(t)}^{2}\}.$$

Moreover, *p* and **u** satisfy the interface condition (2.4) on Γ(*t*) = Ω_{−}(*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

$$\mathrm{\Gamma}(0)=\{(x,y)\in \mathrm{\Omega}:{(x-0.5)}^{2}+{(y-0.5)}^{2}={0.2}^{2}\}.$$

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** = −*H***n**, where *H* is the curvature, and the boundary conditions *p _{∞}* = 0 and boundary profile (2.7) and (2.8). Our numerical solution for Γ(

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.

We consider Ω = (0, *l _{x}*) × (0,

$$\varphi (x,y)=\{\begin{array}{cc}\sqrt{{(x-{l}_{x}/2+c)}^{2}+(y-{l}_{y}/2)}& \text{if}\phantom{\rule{0.2em}{0ex}}x\le {l}_{x}/2-c,\\ |y-{l}_{y}/2|& \text{if}\phantom{\rule{0.2em}{0ex}}{l}_{x}/2-c\le x\le {l}_{x}/2+c,\\ \sqrt{{(x-{l}_{x}/2-c)}^{2}+(y-{l}_{y}/2)}& \text{if}\phantom{\rule{0.2em}{0ex}}x\ge {l}_{x}/2+c.\end{array}$$

(5.1)

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

$$\mathbf{f}=(-\gamma \nabla \xb7(\frac{\nabla \varphi}{|\nabla \varphi |})+4\varepsilon \sum _{i=1}^{2}(\frac{{\sigma}^{12}}{{|\mathbf{x}-{\mathbf{x}}_{i}|}^{12}}-\frac{{\sigma}^{6}}{{|\mathbf{x}-{\mathbf{x}}_{i}|}^{6}}\left)\right)\mathbf{n},$$

(5.2)

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

We then solve the Stokes equation coupled with level-set equation (2.9), with *l _{y}* = 1 and various choices of

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.1*k*) and (0.5 + *c*, 0.25 + 0.1*k*) with *k* {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, *C _{m}* = 0.001 in (5.2). Moreover, for a loose initial condition of Γ, we define

$${\mathrm{\Omega}}_{-}=\{(x,y)\in \mathrm{\Omega}:0.15\le x\le 0.85,0.2\le y\le 0.8\}.$$

For a tight initial condition on Γ, we define

$${\mathrm{\Omega}}_{-}=\{(x,y)\in \mathrm{\Omega}:0.25\le x\le 0.45,0.2\le y\le 0.8\}\cup \{(x,y)\in \mathrm{\Omega}:0.55\le x\le 0.75,0.2\le y\le 0.8\}.$$

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.

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** = (*n*_{1}, *n*_{2}) and ** τ** = (–

$$\{\begin{array}{c}2{n}_{1}^{2}{u}_{x}+2{n}_{1}{n}_{2}({u}_{y}+{v}_{x})+2{n}_{2}^{2}{v}_{y}-p={f}_{\perp},\\ -2{n}_{1}{n}_{2}{u}_{x}+({n}_{1}^{2}-{n}_{2}^{2})({u}_{y}+{v}_{x})+2{n}_{1}{n}_{2}{v}_{y}={f}_{\parallel},\end{array}$$

where *f*_{} = **f** · **n** and *f*_{‖}= **f** · ** τ**. Some straightforward algebraic calculations together with the incompressibility condition (2.2) lead to

$$2{u}_{x}-({n}_{1}^{2}-{n}_{2}^{2})p=({n}_{1}^{2}-{n}_{2}^{2}){f}_{\perp}-2{n}_{1}{n}_{2}{f}_{\parallel},$$

(A.1)

$$2({n}_{2}^{2}-{n}_{1}^{2}){v}_{y}+2{n}_{1}{n}_{2}({u}_{y}+{v}_{x})-p={f}_{\perp}.$$

(A.2)

For any ghost velocity point **x**, we find a point **x*** Γ, such that |**x** – **x***| = *dist*(**x**, **Γ**). We call **x*** a projection point of **x** onto Γ. We then discretize equation (A.1) at each projection point
${\mathbf{x}}_{i-1/2,j}^{\ast}$ corresponding to each ghost velocity point **x**_{i}_{−1/2,}* _{j}* of

To obtain a second-order convergence scheme for *u, v*, and *p* up to Γ, we design a third-order discretization of *u _{x}, u_{y}, v_{x}, v_{y}*, and

$$|S(u,{\mathbf{x}}^{\ast},r)|=\frac{(r+1)(r+2)}{2},|S(v,{\mathbf{x}}^{\ast},r)|=\frac{(r+1)(r+2)}{2},\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}|S(p,{\mathbf{x}}^{\ast},r)|=\frac{r(r+1)}{2}.$$

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** *S*(*u*, **x***, *r*) *S*(*v*, **x***,*r*)*S*(*p*, **x***,*r*); (3) All *S*(*u*, **x***, *r*), *S*(*v*, **x***, *r*), and *S*(*p*, **x***, *r*) satisfy

$$\{\begin{array}{c}\mathit{max}\{i:{\mathbf{x}}_{i,j}\in S(u,{\mathbf{x}}^{\ast},r)\}-\mathit{min}\{i:{\mathbf{x}}_{i,j}\in S(u,{\mathbf{x}}^{\ast},r)\}=r,\\ \mathit{max}\{j:{\mathbf{x}}_{i,j}\in S(u,{\mathbf{x}}^{\ast},r)\}-\mathit{min}\{j:{\mathbf{x}}_{i,j}\in S(u,{\mathbf{x}}^{\ast},r)\}=r,\\ \mathit{max}\{i:{\mathbf{x}}_{i,j}\in S(v,{\mathbf{x}}^{\ast},r)\}-\mathit{min}\{i:{\mathbf{x}}_{i,j}\in S(v,{\mathbf{x}}^{\ast},r)\}=r,\\ \mathit{max}\{j:{\mathbf{x}}_{i,j}\in S(v,{\mathbf{x}}^{\ast},r)\}-\mathit{min}\{j:{\mathbf{x}}_{i,j}\in S(v,{\mathbf{x}}^{\ast},r)\}=r,\\ \mathit{max}\{i:{\mathbf{x}}_{i,j}\in S(p,{\mathbf{x}}^{\ast},r)\}-\mathit{min}\{i:{\mathbf{x}}_{i,j}\in S(p,{\mathbf{x}}^{\ast},r)\}=r-1,\\ \mathit{max}\{j:{\mathbf{x}}_{i,j}\in S(p,{\mathbf{x}}^{\ast},r)\}-\mathit{min}\{j:{\mathbf{x}}_{i,j}\in S(p,{\mathbf{x}}^{\ast},r)\}=r-1.\end{array}$$

(A.3)

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

$$({\lambda}_{1}{h}_{x},{\lambda}_{2}{h}_{y})={\mathbf{x}}^{\ast}-\mathbf{x},({s}_{1},{s}_{2})=(\text{sign}({\lambda}_{1}),\text{sign}({\lambda}_{2})).$$

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** = **x**_{i}_{−1/2,}* _{j}* and

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 **x**_{i−1/2j} with its projection point **x***. Let us denote

$$\mathbf{X}={({\mathbf{x}}_{0},{\mathbf{x}}_{1},{\mathbf{x}}_{2},{\mathbf{x}}_{3},{\mathbf{x}}_{4},{\mathbf{x}}_{5},{\mathbf{x}}_{6},{\mathbf{x}}_{7},{\mathbf{x}}_{8},{\mathbf{x}}_{9})}^{T}={({\mathbf{x}}_{i-\frac{1}{2},j},{\mathbf{x}}_{i-\frac{1}{2}+{s}_{1},j},{\mathbf{x}}_{i-\frac{1}{2}-{s}_{1},j+{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2},j+{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2}+{s}_{1},j+{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2}+2{s}_{1},j+{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2},j+2{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2}+{s}_{1},j+2{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2}+2{s}_{1},j+2{s}_{2}},{\mathbf{x}}_{i-\frac{1}{2},j+3{s}_{2}})}^{T},$$

and **F** = (*u, u _{x}, u_{y}, u_{xx}, u_{xy}, u_{yy}, u_{xxx}, u_{xyy}, u_{xxy}, u_{yyy}*)

$$A=\left[\begin{array}{cccccccccc}1& {h}_{x,1}& {h}_{y,1}& {h}_{x,1}^{2}/2& {h}_{x,1}{h}_{y,1}& {h}_{y,1}^{2}/2& {h}_{x,1}^{3}/6& {h}_{x,1}^{2}{h}_{y,1}/2& {h}_{x,1}{h}_{y,1}^{2}/2& {h}_{y,1}^{3}/6\\ 1& {h}_{x,2}& {h}_{y,2}& {h}_{x,2}^{2}/2& {h}_{x,2}{h}_{y,2}& {h}_{y,2}^{2}/2& {h}_{x,2}^{3}/6& {h}_{x,2}^{2}{h}_{y,2}/2& {h}_{x,2}{h}_{y,2}^{2}/2& {h}_{y,2}^{3}/6\\ 1& {h}_{x,3}& {h}_{y,3}& {h}_{x,3}^{2}/2& {h}_{x,3}{h}_{y,3}& {h}_{y,3}^{2}/2& {h}_{x,3}^{3}/6& {h}_{x,3}^{2}{h}_{y,3}/2& {h}_{x,3}{h}_{y,3}^{2}/2& {h}_{y,3}^{3}/6\\ 1& {h}_{x,4}& {h}_{y,4}& {h}_{x,4}^{2}/2& {h}_{x,4}{h}_{y,4}& {h}_{y,4}^{2}/2& {h}_{x,4}^{3}/6& {h}_{x,4}^{2}{h}_{y,4}/2& {h}_{x,4}{h}_{y,4}^{2}/2& {h}_{y,4}^{3}/6\\ 1& {h}_{x,5}& {h}_{y,5}& {h}_{x,5}^{2}/2& {h}_{x,5}{h}_{y,5}& {h}_{y,5}^{2}/2& {h}_{x,5}^{3}/6& {h}_{x,5}^{2}{h}_{y,5}/2& {h}_{x,5}{h}_{y,5}^{2}/2& {h}_{y,5}^{3}/6\\ 1& {h}_{x,6}& {h}_{y,6}& {h}_{x,6}^{2}/2& {h}_{x,6}{h}_{y,6}& {h}_{y,6}^{2}/2& {h}_{x,6}^{3}/6& {h}_{x,6}^{2}{h}_{y,6}/2& {h}_{x,6}{h}_{y,6}^{2}/2& {h}_{y,6}^{3}/6\\ 1& {h}_{x,7}& {h}_{y,7}& {h}_{x,7}^{2}/2& {h}_{x,7}{h}_{y,7}& {h}_{y,7}^{2}/2& {h}_{x,7}^{3}/6& {h}_{x,7}^{2}{h}_{y,7}/2& {h}_{x,7}{h}_{y,7}^{2}/2& {h}_{y,7}^{3}/6\\ 1& {h}_{x,8}& {h}_{y,8}& {h}_{x,8}^{2}/2& {h}_{x,8}{h}_{y,8}& {h}_{y,8}^{2}/2& {h}_{x,8}^{3}/6& {h}_{x,8}^{2}{h}_{y,8}/2& {h}_{x,8}{h}_{y,8}^{2}/2& {h}_{y,8}^{3}/6\\ 1& {h}_{x,9}& {h}_{y,9}& {h}_{x,9}^{2}/2& {h}_{x,9}{h}_{y,9}& {h}_{y,9}^{2}/2& {h}_{x,9}^{3}/6& {h}_{x,9}^{2}{h}_{y,9}/2& {h}_{x,9}{h}_{y,9}^{2}/2& {h}_{y,9}^{3}/6\end{array}\right]$$

(A.4)

with (*h _{x}*

$${u}_{x}({\mathbf{x}}^{\ast})=-\frac{({\lambda}_{2}-{s}_{2})({\lambda}_{2}-2{s}_{2})}{2{s}_{1}{h}_{x}}{u}_{0}+\frac{({\lambda}_{2}-{s}_{2})({\lambda}_{2}-2{s}_{2})}{2{s}_{1}{h}_{x}}{u}_{1}-\frac{3{\lambda}_{1}^{2}-6{\lambda}_{1}{s}_{1}+2}{6{s}_{1}{h}_{x}}{u}_{2}-\frac{{s}_{1}({\lambda}_{1}+{\lambda}_{1}{\lambda}_{2}{s}_{2})+3{\lambda}_{2}{s}_{2}/2-{\lambda}_{2}^{2}-3{\lambda}_{1}^{2}/2}{{s}_{1}{h}_{x}}{u}_{3}-\frac{{s}_{1}({\lambda}_{1}-2{\lambda}_{1}{\lambda}_{2}{s}_{2})+{\lambda}_{2}^{2}-{\lambda}_{2}{s}_{2}-1+3{\lambda}_{1}^{2}/2}{{s}_{1}{h}_{x}}{u}_{4}+\frac{{\lambda}_{2}/2-2{s}_{2}/3+{\lambda}_{1}^{2}{s}_{2}/2-{s}_{1}{\lambda}_{1}{\lambda}_{2}+{s}_{1}{\lambda}_{1}{s}_{2}}{{s}_{1}{s}_{2}{h}_{x}}{u}_{5}-\frac{({\lambda}_{2}-{s}_{2})({\lambda}_{2}{s}_{1}-2{\lambda}_{1}{s}_{2}+{s}_{1}{s}_{2})}{2{h}_{x}}{u}_{6}+\frac{({\lambda}_{2}-{s}_{2})({\lambda}_{2}{s}_{1}-4{\lambda}_{1}{s}_{2}+2{s}_{1}{s}_{2})}{2{h}_{x}}{u}_{7}+\frac{({\lambda}_{2}-{s}_{2})(2{\lambda}_{1}-{s}_{1})}{2{s}_{2}{h}_{x}}{u}_{8},$$

where the notation *u*(**X**) = (*u*_{0}, *u*_{1}, *u*_{2}, *u*_{3}, *u*_{4}, *u*_{5}, *u*_{6}, *u*_{7}, *u*_{8}, *u*_{9})_{T} is used. The third entry of **F**(**x***)^{t} gives a third-order discretization of *u _{y}*(

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

$$p({\mathbf{x}}^{\ast})=(\frac{3}{8}-{\lambda}_{1}{s}_{1}+\frac{{\lambda}_{1}^{2}}{2}){p}_{0}-\frac{(2{\lambda}_{1}-3{s}_{1})(2{\lambda}_{1}{s}_{2}-2{\lambda}_{2}{s}_{1}+3{s}_{1}{s}_{2})}{4{s}_{2}}{p}_{1}+\frac{((2{\lambda}_{1}{s}_{2}-2{\lambda}_{2}{s}_{1}+{s}_{1}{s}_{2})(2{\lambda}_{1}{s}_{2}-2{\lambda}_{2}{s}_{1}+3{s}_{1}{s}_{2}))}{8}{p}_{2}-\frac{({\lambda}_{2}-{s}_{2})(2{\lambda}_{1}-3{s}_{1})}{2{s}_{1}{s}_{2}}{p}_{3}+\frac{({\lambda}_{2}-{s}_{2})(2{\lambda}_{1}{s}_{2}-2{\lambda}_{2}{s}_{1}+3{s}_{1}{s}_{2})}{2\ast {s}_{1}}{p}_{4}+\frac{({\lambda}_{2}-{s}_{2})({\lambda}_{2}-2{s}_{2})}{2}{p}_{5},$$

Where ${p}_{0}=p({\mathbf{x}}_{i-\frac{{s}_{1}+1}{2},j+{s}_{2}})$, ${p}_{1}=p({\mathbf{x}}_{i-\frac{{s}_{1}-1}{2},j+{s}_{2}})$, ${p}_{2}=p({\mathbf{x}}_{i+\frac{3{s}_{1}-1}{2},j+{s}_{2}})$, ${p}_{3}=p({\mathbf{x}}_{i+\frac{{s}_{1}-1}{2},j+2{s}_{2}})$, ${p}_{4}=p({\mathbf{x}}_{i+\frac{3{s}_{1}-1}{2},j+2{s}_{2}})$, ${p}_{5}=p({\mathbf{x}}_{i+\frac{3{s}_{1}-1}{2},j+3{s}_{2}})$.

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*(**x**_{k}), *v*(**y**_{k}), and *p*(**z**_{k}) at **x***, with **x**_{k}*S*(*u*, **x***, *r*), **y**_{k}*S*(*v*, **x***, *r*), and **z**_{k}*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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |