A new finite element and finite difference hybrid method for computing electrostatics of ionic solvated biomoleculeJournal of Computational Physics


Jinyong Ying, Dexuan Xie
Physics and Astronomy (miscellaneous) / Computer Science Applications


Journal of Computational Physics 298 (2015) 636–651

Contents lists available at ScienceDirect

A fo


De a













So 1. ac co ha an

AP st us le so ˜ tio po * ht 00Journal of Computational Physics www.elsevier.com/locate/jcp new finite element and finite difference hybrid method r computing electrostatics of ionic solvated biomolecule nyong Ying, Dexuan Xie ∗ partment of Mathematical Sciences, University of Wisconsin-Milwaukee, Milwaukee, WI, 53201, USA r t i c l e i n f o a b s t r a c t ticle history: ceived 24 November 2014 ceived in revised form 10 June 2015 cepted 24 June 2015 ailable online 29 June 2015 ywords: isson–Boltzmann equation ite element method ite difference method molecular electrostatics main decomposition ultigrid preconditioning lvation free energy

The Poisson–Boltzmann equation (PBE) is one widely-used implicit solvent continuum model for calculating electrostatics of ionic solvated biomolecule. In this paper, a new finite element and finite difference hybrid method is presented to solve PBE efficiently based on a special seven-overlapped box partition with one central box containing the solute region and surrounded by six neighboring boxes. In particular, an efficient finite element solver is applied to the central box while a fast preconditioned conjugate gradient method using a multigrid V-cycle preconditioning is constructed for solving a system of finite difference equations defined on a uniform mesh of each neighboring box. Moreover, the PBE domain, the box partition, and an interface fitted tetrahedral mesh of the central box can be generated adaptively for a given PQR file of a biomolecule. This new hybrid

PBE solver is programmed in C, Fortran, and Python as a software tool for predicting electrostatics of a biomolecule in a symmetric 1:1 ionic solvent. Numerical results on two test models with analytical solutions and 12 proteins validate this new software tool, and demonstrate its high performance in terms of CPU time and memory usage. © 2015 Elsevier Inc. All rights reserved.


Electrostatic interactions are important in understanding the structure and biological function of biomolecules, catalytic tivity, ligand association, and pKa values [1–4]. The Poisson–Boltzmann equation (PBE) is one widely-used implicit solvent ntinuum model for calculating electrostatics of biomolecules (protein or nucleic acids) in ionic solvent [5–7]. A lot of work s been done in the development of PBE numerical solvers and program packages based on finite difference, finite element, d boundary integral equation approaches [8–13]. The popular PBE program packages and web-based resources, such as

BS [10], DelPhi [14], PBEQ [12,15], and UHBD [16], have become powerful simulation aides in the study of biomolecular ructure and function, biomolecule-ligand association, ion channel modeling, and computer-aided drug design.

To improve the accuracy of PBE numerical solutions, we recently developed a solution decomposition PBE solver (SDPB) ing finite element and minimization techniques [13]. Instead of solving PBE directly, in SDPB, two other interface probms — one linear interface problem for  and one nonlinear interface problem for ˜ — are solved to yield a numerical lution u of PBE as a sum of a given function G with  and ˜. Since G collects all the singularity points of u, both  and become twice continuously differentiable within both the solute and solvent regions. Thus, their numerical approximans can be calculated much more easily and less expensively than directly calculating u, which is singular at each atomic sition r j , to reach a required numerical accuracy. We note that in DelPhi, APBS and PBEQ, a finite difference method

Corresponding author.

E-mail address: dxie@uwm.edu (D. Xie). tp://dx.doi.org/10.1016/j.jcp.2015.06.016 21-9991/© 2015 Elsevier Inc. All rights reserved.

J. Ying, D. Xie / Journal of Computational Physics 298 (2015) 636–651 637 de

It fin in tw to m of el hy p va im of di of co ne lin en (P el it to bo nu ra m al on of gi pa se si ge th w on po

Th of ca th ge hy ic pr th w 11 sa ne m thfined on a uniform mesh accelerated by geometric multigrid techniques is the major algorithm for solving PBE directly. s mesh size h may become very small to catch the solution behavior of PBE around each singularity point r j , making the ite difference system too large to be solved due to the memory limitation of computers. Hence, it is critical to solve PBE directly through computing  and ˜.

Another feature of SDPB is to use an unstructured interface fitted tetrahedral mesh to approximate the interface  beeen the solute and solvent regions. Such a treatment can yield a numerical solution of PBE in high accuracy in comparison the case of a uniform finite difference mesh [17]. However, using an unstructured mesh may cause efficient geometric ultigrid techniques no longer to work, and require extra memory locations to store the mesh data and the nonzero entries the coefficient matrix of each involved linear system.

In order to confine the disadvantages while retaining the advantages of SDPB, in this paper, we present a new finite ement and finite difference hybrid method based on the Schwartz domain decomposition approach [18,19]. This new brid PBE solver was motivated from the fact that the linear variational problem for determining the search direction k of the modified Newton minimization algorithm from SDPB can be reformulated as a linear elliptic interface boundary lue problem (see Theorem 3.1). From this fact it implies that the modified Newton minimization algorithm can also be plemented based on the finite difference approach. The only change to be made is to solve the linear interface problems  (see (5)) and pk (see (11)) by a finite difference method. To combine the advantages of finite element and finite fference approaches together, we propose an overlapped box iterative method for solving the linear interface problems  and pk based on a special partition of the whole domain  into seven overlapped boxes, in which one central box ntains the solute region Dp and is surrounded by the six neighboring boxes (see Fig. 1 for an illustration). Since each ighboring box is a part of the solvent region Ds only, the corresponding equations of  and pk are reduced to the regular ear elliptic boundary value problems. Thus, they can be solved by a finite difference method defined on a uniform mesh to able the usage of efficient geometric multigrid techniques. To do so, we construct the preconditioned conjugate gradient