Fast Newton Iterative Method for Local Steric Poisson–Boltzmann Theories in Biomolecular Solvation

11/01/2021
by   Minhong Chen, et al.
Shanghai Jiao Tong University
0

This work proposes a fast iterative method for local steric Poisson–Boltzmann (PB) theories, in which the electrostatic potential is governed by the Poisson's equation and ionic concentrations satisfy equilibrium conditions. To present the method, we focus on a local steric PB theory derived from a lattice-gas model, as an example. The advantages of the proposed method in efficiency lies on a key idea that ionic concentrations as scalar implicit functions of the electrostatic potential, i.e, generalized Boltzmann distributions, are numerically available. The existence, uniqueness, boundness, and smoothness of such distributions are rigorously established. A Newton iteration method with truncation is proposed to solve a nonlinear system discretized from the generalized PB equations. The existence and uniqueness of the solution to the discretized nonlinear system are established by showing that it is a unique minimizer of a constructed convex energy. Thanks to the boundness of ionic concentrations, truncation bounds for the potential are obtained by using the extremum principle. The truncation step in iterations is shown to be energy, residual, and error decreasing. To further speed-up computations, we propose a novel precomputing-interpolation strategy, which is applicable to other local steric PB theories and makes the proposed methods for solving steric PB theories as efficient as for solving the classical PB theory. Analysis on the Newton iteration method with truncation shows local quadratic convergence for the proposed numerical methods. Applications to realistic biomolecular solvation systems reveal that counterions with steric hindrance stratify in an order prescribed by the parameter of ionic valence-to-volume ratio. Finally, we remark that the proposed iterative methods for local steric PB theories can be readily incorporated in well-known classical PB solvers.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 18

page 19

page 20

page 21

page 22

09/11/2020

Analysis of a new implicit solver for a semiconductor model

We present and analyze a new iterative solver for implicit discretizatio...
06/10/2021

Existence of Strong Solution for the Complexified Non-linear Poisson Boltzmann Equation

We prove the existence and uniqueness of the complexified Nonlinear Pois...
08/31/2021

An Efficient Finite Element Iterative Method for Solving a Nonuniform Size Modified Poisson-Boltzmann Ion Channel Model

In this paper, a nonuniform size modified Poisson-Boltzmann ion channel ...
04/17/2022

On an interior-exterior nonoverlapping domain decomposition method for the Poisson–Boltzmann equation

A nonoverlapping domain decomposition method is studied for the lineariz...
09/28/2021

Solution decomposition for the nonlinear Poisson-Boltzmann equation using the range-separated tensor format

The Poisson-Boltzmann equation (PBE) is an implicit solvent continuum mo...
12/26/2019

Travelling wave mathematical analysis and efficient numerical resolution for a one-dimensional model of solid propellant combustion

We investigate a model of solid propellant combustion involving surface ...
02/27/2021

Reduced basis method for the nonlinear Poisson-Boltzmann equation regularized by the range-separated canonical tensor format

The Poisson-Boltzmann equation (PBE) is a fundamental implicit solvent c...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 long-range 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 mean-field 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 ion-ion 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 mean-field theory is known to neglect ionic steric effects and ion-ion 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 hard-sphere interactions between ions using the Lennard-Jones (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 hard-sphere 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 free-energy 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 memory-saving 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 lattice-gas 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 valence-to-volume ratio.

We remark that the proposed numerical method has salient features in efficiency and memory-saving. Several reasons account for such features: First, the iterations only involve the unknowns of electrostatic potential, rather than both potential and (possible multiple-species) ionic concentrations; Second, the derivative of concentrations with respect to the potential is used in iterations, making the super-linear convergence of the whole algorithm possible. Third, to further speed-up computations, we propose a precomputing-interpolation 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 well-established 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

Figure 1: A schematic view of a charged biomolecule solvated in an ionic solution with ionic species. The system occupies a region , which is divided into a solute region and a solvent region by a solute-solvent interface . The solute contains atoms, located at , carry partial charges (). The ions are shown individually to highlight the steric effects of ions in electric double layers next to charged biomolecules.

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 solute-solvent 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 boundary-value 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 free-energy 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 mean-field 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 hard-sphere interactions between ions can be described by the Lennard-Jones (LJ) potential energy, which gives rise to a nonlocal model [26]. To avoid computationally intractable integro-differential 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 hard-sphere 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 finite-difference 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 signed-distance level-set function whose zero level set represents the solute-solvent 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 over-relaxation) 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 un-updated 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 memory-saving 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 Speed-Up 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 Newton-type iteration methods are considered.

To speed-up 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 high-order 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 high-order interpolation schemes to interpolate and for each given with stored data and , when solving (3.10) iteratively. We remark that such a precomputing-interpolation 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 off-diagonal 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

  1. If satisfies , , then .

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

1:  Compute the upper and lower estimates and . Initialize so that . Let the iteration step . Choose a stopping tolerance .
2:  while  do
3:     Compute