1 Introduction
Electrostatic interactions between biomolecules, water molecules, and mobile ions are fundamental to the function and stability of solvated biomolecules, such as membranes and proteins [41]. Ionic steric effects in electrostatic interactions have profound impact on the dynamics of an underlying biomolecular system. For example, ionic steric effect plays a key role in the selection of ions that permeate through transmembrane channels [22]. Due to the longrange nature of the Coulomb force, electrostatic interactions are difficult to deal with, not only in continuum level but also in molecular simulations. As one of the most popular theories, the classical Poisson–Boltzmann (PB) theory has achieved great success in describing electrostatic interactions in biomolecular systems, electrochemical systems, etc. [14, 2]. Based on the meanfield approximation, the classical PB theory treats ions as point charges, takes the screening effect of water solvent into account as a continuum with a certain dielectric coefficient, and ignores direct ionion interaction details, resulting in ions interacting only with the background electrostatic potential. With such approximations, the PB equation, which couples the Poisson’s equation for the electrostatic potential and Boltzmann distributions for ionic concentrations, becomes a mathematically simple model that can be efficiently solved by various numerical methods [40, 32].
The PB theory has been very successful in many applications. Nonetheless, the meanfield theory is known to neglect ionic steric effects and ionion correlations [21]. Without steric hindrance, the PB theory could predict unphysically high counterion concentration next to a charged surface [7, 53]. Recent years have seen a growing interest in incorporating the ionic steric effects within the framework of the PB theory [7, 8, 30, 29, 39, 53]. The steric effect can be effectively incorporated in a variational approach by adding an excess chemical potential energy to the classical PB free energy [5]. One of the most popular local steric PB models is based on the statistical mechanics of ions and solvent molecules on lattices. In such a model, the excess chemical potential is described by the entropy of solvent molecules [6, 7, 29, 39, 24]. The excess chemical potential can also represent hardsphere interactions between ions using the LennardJones (LJ) potential energy [26]. This leads to a nonlocal model, which can be further approximated using the Fourier analysis to obtain a computationally more tractable local model [25, 35, 16]. Alternatively, the excess chemical potential can be determined by the Carnahan–Starling (CS) equation of state for hardsphere liquids of a uniform size or the Boublik–Mansoori–Carnahan–Starling–Leland (BMCSL) equation of state for unequal sizes [5]. The excess chemical potential determined by the equation of state of CS or BMCSL can also derive local steric PB theories.
Extensive efforts have been devoted to the development of numerical methods for the classical PB equation, ranging from finite difference methods [15, 9, 27, 56, 44], finite element methods [3, 23, 13, 48, 46] to boundary element methods [37, 20, 19, 38]. Such advances have led to the development of several successful software packages, such as APBS [4, 18], DelPhi [43, 34], MIBPB [56, 12], and AFMPB [38]. Nonetheless, not many numerical methods have been developed for steric PB models, especially the models with nonuniform ionic sizes. One of the main obstacles to the development of an efficient solver for steric PB models with nonuniform ionic sizes is that the Boltzmann distributions in the classical PB theory are no longer available. An augmented Lagrange multiplier method has been proposed to minimize a freeenergy functional subject to constraints for a steric PB model in the work [53]. Newton iterative relaxation finite element methods have been developed to solve steric PB models with nonuniform ionic sizes, which, after discretization, become a large nonlinear system coupling unknowns of both ionic concentrations and the electrostatic potential on computational grid points [47, 33, 49]. The unknowns of each ionic concentrations on grid points were updated in each relaxation step in a coupled way [47]. Based on an observation made in the work [28] that the unknowns of concentrations are spatially decoupled for each grid point with a given electrostatic potential, efficient and memorysaving iteration methods could be developed to solve steric PB models with nonuniform ionic sizes.
In this work, we first briefly review several versions of local steric PB theories. We then focus on a widely used steric PB model based on the latticegas theory, and introduce the corresponding numerical discretization method for the Poisson’s equation and equilibrium conditions for ionic concentrations. Since ionic concentrations are spatially decoupled for each grid point with a given electrostatic potential [28]
, we treat ionic concentrations as implicit scalar functions of the electrostatic potential, i.e., generalized Boltzmann distributions. The existence, uniqueness, boundness, and smoothness of such distributions are rigorously established for a given electrostatic potential. We develop a Newton iterative method with truncation for the nonlinear system resulting from discretization of the generalized PB equation, which couples the Poisson’s equation with the generalized Boltzmann distributions. The existence and uniqueness of the solution to the discretized nonlinear system are proved by showing that the solution is a unique minimizer of a constructed convex energy functional. By the boundness of ionic concentrations, we are able to establish upper and lower estimates on the electrostatic potential using the extremum principle. Such estimates are employed to truncate the solution in Newton iterations. The truncation step is proved to be error, energy, and residual decreasing. Further analysis on the iteration method establishes that the method has local quadratic convergence rate. Numerical simulations are performed to demonstrate the efficiency and robustness of the proposed numerical methods. Applications to biomolecular solvation systems reveal that ions stratify next to the surface of biomolecules when the electrostatic potential is strong. Also, it is confirmed that the order of stratification is prescribed by the parameter of ionic valencetovolume ratio.
We remark that the proposed numerical method has salient features in efficiency and memorysaving. Several reasons account for such features: First, the iterations only involve the unknowns of electrostatic potential, rather than both potential and (possible multiplespecies) ionic concentrations; Second, the derivative of concentrations with respect to the potential is used in iterations, making the superlinear convergence of the whole algorithm possible. Third, to further speedup computations, we propose a precomputinginterpolation strategy, in which the generalized Boltzmann distributions and their derivatives are interpolated from precomputed and stored data. Such a strategy is directly applicable to other local steric PB theories, and makes the proposed algorithm for solving steric PB theories as efficient as for solving the classical PB theory. Finally, we remark that our efficient iterative methods for the local steric PB theory can be readily incorporated to wellestablished PB solvers, such as the APBS [4, 18], DelPhi [43, 34], MIBPB [56, 12], and AFMPB [38].
The rest of this paper is organized as follows: In Section 2, we introduce local steric PB theories. In Section 3, we propose numerical methods, and present some analysis on the model and numerical methods. In Section 4, we report some numerical simulations to demonstrate the effectiveness of the numerical methods. Finally, we draw conclusions in Section 5.
2 Local Steric Poisson–Boltzmann Theories
Consider the solvation of charged biomolecules in an ionic solution, in which the aqueous solvent is treated implicitly as a continuum and the distribution of ions is described by ionic concentrations. As shown in the schematic view of Fig. 1, the solvation system occupies a region , which is an open and connected subset of with a smooth boundary . A solutesolvent interface separates the solute region, , from the solvent region, . Assume that the solute consists of atoms located at , carrying fixed partial charges (). Suppose that there are ionic species in the ionic solution. Deonte by and the valence and local ionic concentration of the th () ionic species at position , respectively. Let The volume associated to an ion of the th species is denoted by .
The electrostatic potential of the charged system is governed by a boundaryvalue problem (BVP) of the Poisson’s equation
(2.1) 
where is the vacuum permittivity, the dielectric coefficient is if and is if , and is a given bounded and smooth function, representing the electrostatic potential on the boundary. The total charge density is given by
where
is a characteristic function of the solvent region
, is the elementary charge, and the fixed charge distribution represents the point charges carried by solute atoms.To incorporate ionic steric effects in the framework of the Poisson–Boltzmann (PB) theory, we consider the electrostatic solvation freeenergy functional of ionic concentrations:
(2.2) 
Here the electrostatic solvation energy is given by [51, 36]
where is determined by the BVP (2.1), denotes the unit exterior normal on , and is the reaction potential, with being the electrostatic potential of the same charged biomolecules before solvation:
(2.3) 
The entropy of ions is given by
where with the Boltzmann constant and the temperature. The excess potential energy accounts for ionic steric effects. First variation of with respect to leads to chemical potentials
Here, the sum of first two terms is the chemical potential of ideal gas, and the excess chemical potential accounts for the ionic steric effects. See the references [30, 29, 36] for formal calculations on the first variation. In the equilibrium, the total chemical potential can be further determined by the bulk ionic concentrations, , which are achieved when the electrostatic potential vanishes. Thus, we have equilibrium conditions for ionic concentrations:
(2.4) 
Such nonlinear equilibrium conditions determine the generalized Boltzmann distributions . As such, we obtain generalized PB equations with steric effects:
(2.5) 
Ions are treated as point charges in the meanfield approximations, which break down when the steric repulsion and correlation are no longer ignorable. To account for steric effects, one simple approach is using local density approximations, in which the excess chemical potential is approximated as a local function of ionic concentrations [5]. Here, we briefly review several commonly used local models in literature. One popular model is based on the statistical mechanics of ions and solvent molecules on lattices [6, 7, 29, 39]. The entropy of solvent molecules accounts for the excess chemical potential
(2.6) 
where is the volume of a solvent molecule and is the solvent concentration defined by
(2.7) 
The hardsphere interactions between ions can be described by the LennardJones (LJ) potential energy, which gives rise to a nonlocal model [26]. To avoid computationally intractable integrodifferential equations, local approximations of nonlocal integrals can be employed to obtain local models [25, 35, 16]. In such a local model, the excess chemical potential is given by
(2.8) 
where is a symmetric matrix to represent the cross interactions between different species and self interactions between ions of the same species. Alternatively, the excess chemical potential can be determined by the Carnahan–Starling (CS) equation of state for hardsphere liquids of a uniform size [5]:
(2.9) 
where is volume fraction of ions. To extend the model to mixtures of unequal size, one can choose the Boublik–Mansoori–Carnahan–Starling–Leland (BMCSL) equation of state to derive the excess potential [5].
For simplicity of presentation, we shall focus on the most popular local steric PB theory, which is the generalized PB equations (2.5) with the excess chemical potential given by (2.6), i.e.,
(2.10) 
We remark that the numerical methods presented below are applicable to other local steric PB theories with the excess chemical potential derived from the approximation of LJ potential (2.8), the CS equation of state (2.9), or the BMCSL equation of state.
3 Numerical Methods
3.1 Discretization
We choose a computational domain , and cover the domain with a uniform mesh
where is a positive number and is the uniform grid spacing. Define
and . For a given known function , we denote . For an unknown function , we denote as the numerical approximation of on the grid point .
The dielectric coefficients is determined by a given molecular surface . The numerical instability in the calculations of solvation energies and forces, due to relative location and orientation of biomolecules with respect to the finitedifference grid, has been realized in literature [44, 27]. Mesh refinement could be a solution, but too small mesh resolution is still not affordable in 3D calculations. One approach to reduce the instability due to grid dependence is to smooth the sharp transition of dielectric coefficient across the solute and solvent regions [15, 9, 44, 18, 10, 11]. We here introduce a smoothed characteristic function of the solvent region:
where is a signeddistance levelset function whose zero level set represents the solutesolvent interface , and is a smeared Heaviside function
Here is a parameter to control the width of transition layer. With such a smooth transition, the dielectric coefficient function becomes
The electrostatic potential governed by the Poisson’s equation has singularities due to fixed charges carried by solute atoms. To get accurate approximations of , we consider the reaction potential instead. It follows from (2.3) that
Therefore, satisfies the BVP
(3.1) 
With the standard central differencing, the BVP (3.1) is approximated by
(3.2) 
where is the numerical approximation of the concentration of the th species, and the discrete Laplacian operator is defined by
(3.3) 
The term is defined analogously. To improve the convergence of reaction field energies [15, 9], dielectric coefficients on half grid points are approximated by the harmonic average
Given concentrations , the electrostatic potential can be obtained by , where the reaction potential is obtained by solving the linear system (3.2). The nonlinear equilibrium conditions in (2.10) can be discretized as
(3.4) 
In summary, we have the following nonlinear system after discretization:
(3.5) 
Numerical iterative methods for solving the coupled system (3.5) have been proposed in literature to obtain the unknowns . For instance, a nonlinear SOR (successive overrelaxation) scheme has been developed in [47], in which the linear system (3.2) was solved with concentrations given in a previous iteration step, and the nonlinear system (3.4) was solved with a Newton iteration method for each species in a Gauss–Seidel manner, with the electrostatic potential and unupdated concentrations of other species given in a previous step. The unknowns for each ionic concentration on grid points were coupled together in the Newton iteration method. Thus, a linear system of was needed to solve for each species in each Newton iteration for (3.4). As observed in the work [28], the unknowns of concentrations in (3.4) are spatially decoupled for each grid point with a given electrostatic potential. Based on such an observation, much more efficient and memorysaving iteration methods for (3.4) on every single grid point can be designed. Interested readers are referred to the work [28] for more details. In this work, we shall propose a novel and super efficient Newton iteration method with truncation for the coupled system (3.5), by treating the concentrations as implicit scalar functions of the electrostatic potential, i.e., generalized Boltzmann distributions.
3.2 Generalized Boltzmann Distribution
This section focuses on solving the nonlinear equilibrium conditions in (2.10) with a given electrostatic potential , to obtain generalized Boltzmann distributions . The existence, uniqueness, boundness, and smoothness of such distributions are established in the following theorem.
Theorem 3.1.
For each , there exists a unique solution to the equilibrium conditions in (2.10). Moreover, each is a smooth function.
Proof.
The nonlinear equilibrium conditions in (2.10) can be simplified into
(3.6) 
by using the volume fraction of solvent [31], which is defined by
(3.7) 
Analogously, is defined with given bulk concentrations . A combination of (3.6) and (3.7) leads to an equation for :
(3.8) 
There is a unique root for in the interval . Indeed, it is easy to verify that and
Simple calculations show that
Therefore, is a continuous, increasing function with a unique root in . It follows that for each , there exists a unique corresponding . By (3.6), one can see that there exists a unique for each . This establishes the function and the generalized Boltzmann distributions . It follows from (3.7) that . From the equation (3.8), one can show by the implicit function theorem that is a smooth function. Therefore, is a smooth function by (3.6). ∎
For each , the corresponding can be calculated efficiently with a Newton iteration scheme
(3.9) 
With the obtained , the numerical concentrations for each can be calculated via (3.6).
3.3 Precomputing SpeedUp Strategy
With the implicit generalized Boltzmann distributions numerically available for each given , the unknowns on governed by the coupled system (3.5) can be obtained more efficiently by solving
(3.10) 
Here, instead of , is served as the iterative variable. To solve the nonlinear system (3.10) iteratively, we need to evaluate for many times the generalized Boltzmann distributions and , if Newtontype iteration methods are considered.
To speedup computations, we propose to precompute the functions and , and store the functions on a mesh of in a certain interval. For instance, we consider and assume that the interval is large enough to cover the range we are interested in. The two end points and can be estimated in specific applications; See (3.17) for their estimation. We cover the interval with a mesh for , where the mesh spacing . For each grid point , we first compute with the iteration scheme (3.9), then compute with (3.6). Note that can be an initial guess for the iterations (3.9) when computing . Such a continuation method can provide very good initial guesses for the Newton iterations.
For each , it follows from (3.6) that
(3.11) 
where
(3.12) 
By Theorem 3.1, we know that and are smooth functions. Thus, we can compute the derivatives alternatively with highorder difference schemes, e.g.,
with the stored data
. Therefore, we can precompute and store the vectors
and for . To save memory, we can alternatively store the vectors and instead, especially when the number of ionic species under consideration is large. With and stored, we then can compute and directly with (3.6) and (3.11), respectively.By the smoothness of and , we can use highorder interpolation schemes to interpolate and for each given with stored data and , when solving (3.10) iteratively. We remark that such a precomputinginterpolation strategy is applicable to other local steric PB theories with implicit generalized Boltzmann distributions determined by equilibrium conditions (2.4), in which the excess chemical potential can be derived from the approximation of LJ potential (2.8), the CS equation of state (2.9), or the BMCSL equation of state. Furthermore, such a precomputing and interpolation strategy makes our algorithm for solving steric PB theories as efficient as for solving the classical PB theory.
3.4 Newton Iteration Method with Truncation
The nonlinear system (3.10) can be rewritten in a matrix form
(3.13) 
where is the unknown vector, is the coefficient matrix corresponding to the discrete Laplacian, and is a column vector. Here and the column vector results from boundary conditions and the known term . It is standard to show that is a symmetric positive definite matrix with positive diagonal elements and negative offdiagonal elements. Define the residual function by , and its Jacobian matrix by , where . The existence and uniqueness of the solution to the nonlinear system (3.13) can be established in the following theorem.
Theorem 3.2.
There exists a unique solution to the nonlinear system (3.13).
Proof.
Define an energy by
(3.14) 
By direction calculations, we have
It follows from (3.11) and (3.12) that for any ,
By the Cauchy–Schwarz inequality, we have
(3.15) 
Therefore, for any ,
where the positive definiteness of , the definition of , and the conclusion of (3.15) have been used in the last step. Hence, the function is a convex function. Therefore, there exists a unique minimizer of satisfying , i.e., the nonlinear system (3.13). ∎
3.4.1 Upper and Lower Estimates
It follows from Theorem 3.1 that and for . Such bounds can help establish upper and lower estimates for the unknown . Define grid functions and by the following problems
(3.16) 
where are components of the vector . From the discretization of the Laplacian, we have the following standard discrete extremum principle without giving its proof.
Lemma 3.1.
(Extremum Principle) Let be a grid function defined on . Let be the discrete Laplacian operator defined by (3.3). Then

If satisfies , , then .

If satisfies , , then .
Clearly, we have
where the fact that the volume fraction has been used. By the extremum principle, we have the solution , where the set is given by
(3.17) 
with . We now propose the following Newton iteration method with truncation for , with computed upper and lower estimates and .
Comments
There are no comments yet.