Electrostatics interaction plays a very important role in almost all biomolecular systems.  The Poisson-Boltzmann equation is widely used to treat this electrostatic effect in an ionic solution.  In this work, a simple mixed discrete-continuum model is considered and boundary element method is used to solve for the solution.Singapore-MIT Alliance (SMA

Lim, Kian Meng

Tan, Lian Hing

White, Jacob K.

DSpace@MIT

     Abstract – Electrostatics interaction plays a very important role in almost all biomolecular systems.  The Poisson-Boltzmann equation is widely used to treat this electrostatic effect in an ionic solution.  In this work, a simple mixed discrete-continuum model is considered and boundary element method is used to solve for the solution.  Keywords – Boundary Element Method, Biomolecular electrostatics, Poisson-Boltzmann Equation.  1. INTRODUCTION Electrostatics interaction plays a very important role in almost all biomolecular systems.  Some of the areas where electrostatics has been very useful include protein structural stability, enzyme catalysis, biomecular recoginition and biomolecular encounter rates [1].  Thus, there is a need for an accurate modelling and simulation of the biomolecular electrostatiscs. There are few methods to simulate a biomolecule in an ionic solution.  In this paper, a mixed discrete-continuum approach based on combining a continuum description of the macromolecules and solvent with a discrete description of the atomic charges is adopted [2].  Instead of considering the non-linear Poisson-Boltzmann equation, the linearized Poisson-Boltzman equation is considered by assuming that the electrostatic energy of the ions is much lesser than their thermal energy [3]. The boundary element method is used to solve for solution of the equations generated from the above model.  It is the preferred method as it treats the infinite domain and point charges more naturally, and is able to reduce the dimensional of the problem by one. The next section briefly discusses the mixed discrete-continuum model and the boundary integral equations formulation.  Then, numerical solutions and some sample 2D computational result are presented.  Finally, some of the future investigations are discussed.    2. PROBLEM FORMULATION 2.1 Mixed Discrete-continuum Model Fig. 1 shows a simplified mixed discrete-continuum model for the problem. Region 1Ω  corresponds to the interior of the molecule and region 2Ω  corresponds to the surrounding solvent. In region 1Ω , the electrostatic potential is governed by a Poisson equation. 1Ω : 2 111( ) ( )cNkkkqr r rϕ δε=∇ = − −∑r r r     (1) where 1ϕ  is the electrostatic potential, rr is an evaluation position, krr is the location of the point charge, kq  is the strength of the point charge, cN  is the number of point charges, and 1ε  is the dielectric constant in 1Ω . In region 2Ω , Debye-Hückle theory suggests that the electrostatic potential should satisfy a nonlinear Poisson-Boltzmann equation.  This equation can be linearized under the assumption that the electrostatic energy of the ions is much lesser than their thermal energy.  The linearized Poisson-Boltzmann is also a Helmholtz equation. 2Ω : 2 22 2( ) ( ) 0Dr rϕ κ ϕ∇ − =r r      (2) where 2ϕ  is the electrostatic potential, Dκ  is the reciprocal of the Debye length and 2ε  is the dielectric constant in 2Ω .  Fig. 1: The mixed discrete-continuum model Numerical Study of the Poisson-Boltzmann Equation for Biomolecular Electrostatics  Lian Hing TAN 1, Kian Meng LIM1, Jacob WHITE2 1 Singapore-MIT Alliance, E4-04-10, 4 Engineering Drive 3, Singapore 117576 2 Singapore-MIT Alliance, Room 8-407, 77 Massachusetts Avenue, Cambridge, MA 02139 USA Ω1, ε1krr  kr r−r r  rrnr  orrΩ2, ε2S : Boundary surface Point charge    2.2 Boundary Integral Equations Formulation The fundamental solution to equations (1) and (2) are 2D: 11( , ') ln( ' )2G r r r rπ= − −r r r r         (3.1) (1)2 0( , ') ( ' )4iG r r H i r rκ= −r r r r        (3.2) 3D: 11( , ')4 'G r rr rπ= −r r r r           (4.1) '2 ( , ') 4 'r reG r rr rκπ− −== −r rr r r r          (4,2) By applying the direct formulation of the boundary integral equations to equations (1) and (2), we obtain The integral equation for region 1Ω  is 1 11 1 1( ') ( , ')( ) ( ) ( , ') ( ')Sr G r rc r r G r r r dSn nϕϕ ϕ∂ ∂⎡ ⎤= −⎢ ⎥∂ ∂⎣ ⎦∫r r rr r r r r   11( , )cNkkkqG r rε=+∑ r r  (5) The integral equation for region 2Ω  is 2 22 2 1( ') ( , ')( ) ( ) ( , ') ( ')Sr G r rc r r G r r r dSn nϕϕ ϕ∂ ∂⎡ ⎤= − +⎢ ⎥∂ ∂⎣ ⎦∫r r rr r r r r(6) where nr  is the outward pointing normal as shown in Fig. 1, S  is the boundary surface separating 1Ω  and 2Ω , and 0.5, if  is a point on a smooth boundary( ) 1 , if  is a point inside the domain          0 , if  is outside the domain                   or rc r rr=⎧⎪= Ω⎨⎪ Ω⎩r rr rr The potential 1ϕ  and 2ϕ  must satisfy the following boundary conditions. 1 2( ) ( )o or rϕ ϕ=r r         (7) 1 21 2( ) ( )o or rn nϕ ϕε ε∂ ∂=∂ ∂r r       (8) where or S∈r . Substitute (7) and (8) into (5) and (6), we obtain 1 11 1 1( , ') ( ')1 ( ) ( ') ( , ')2 SG r r rr r G r r dSn nϕϕ ϕ ∂ ∂⎡ ⎤+ −⎢ ⎥∂ ∂⎣ ⎦∫r r rr r r r  11( , )cNio iiqG r rε== ∑ r r  (9) 2 1 22 1 22( , ') ( ')1 ( ) ( ') ( , )2oo o o oSG r r rr r G r r dSn nε ϕϕ ϕ ε⎡ ⎤∂ ∂+ − +⎢ ⎥∂ ∂⎣ ⎦∫r r rr r r r0=  (10) Equations (9) and (10) can be used to compute 1ϕ  and 1nϕ∂∂  on the boundary surface S .   3. NUMERICAL SOLUTION 3.1 Discretization Method A standard piecewise-constant centroid collocation scheme is used to discretize (9) and (10).  The boundary surface S  is discretized into N  elements jS  and the midpoint of each element irr  is chosen as the collocation point.  This results in the following 2N system of equations. 111112 12212 ( , )102cj jj jNj j kS S i kkj jS SGI dS G dS qn G r rGI dS G dS nnϕεϕεε=∂⎡ ⎤ ⎡ ⎤+ − ⎡ ⎤⎢ ⎥∂ ⎢ ⎥⎢ ⎥⎢ ⎥ =∂ ⎢ ⎥⎢ ⎥⎢ ⎥∂ ⎢ ⎥− ⎢ ⎥∂⎣ ⎦⎢ ⎥ ⎣ ⎦∂⎣ ⎦∫ ∫ ∑∫ ∫r r(11) where jS∫ corresponds to an integration over jth element.   4. COMPUTATIONAL RESULTS In this paper, the problem being considered has a cylindrical geometry.  Therefore, the original 3D problem can be simplified into a 2D problem.  The simulation results for 2 test cases are presented below.  First, a hypothetical cylindrical molecule with a unit charge at the centre of the molecule is considered.  It is then followed by considering a hypothetical cylindrical molecule with a unit dipole around the centre of the molecule.  4.1 Potential of a hypothetical cylindrical molecule with a unit charge A hypothetical cylindrical molecule of radius a  with a unit charge located at the centre of the molecule in a solution is simulated. The case where 1 1ε = , 2 20ε =  and 1/D aκ =  is considered.  The potential at the radial distance 0.5Cr a= , 0.9Cr a= , Cr a= and 1.1Cr a=  is tabulated in Table 1.  4.2 Potential of a hypothetical cylindrical molecule with a unit dipole A hypothetical cylindrical molecule of radius a  with a unit dipole around the centre of the molecule in a solution is simulated. The case where 1 1ε = , 2 20ε =  and 1/D aκ =   is considered.  The locations of the unit dipole and the locations of the point where the potential is to be evaluated is shown in Fig. 2.  The values of these potentials are tabulated in Table 2.      N  (Number of elements) Potential 50 100 200 0.5Cr a=  0.1168367 0.1169774 0.1170018 0.9Cr a=  0.0232880 0.0234282 0.0234527 Cr a=  0.0096749 0.0096089 0.0095776 1.1Cr a=  0.0058682 0.0058344 0.0058170 Table 1: Potential of a hypothetical cylindrical molecule with a unit charge    Fig. 2: Location of the unit dipoles and the points where the potential is to be evaluated   N  (Number of elements) Potential 50 100 200 0.5Cr a=  0.0494663 0.0495021 0.0495109 0.9Cr a=  0.0083473 0.0084118 0.0084278 Cr a=  0.0031264 0.0031237 0.0031231 1.1Cr a=  0.0015498 0.0015483 0.0015480 Table 2: Potential of a hypothetical cylindrical molecule with a unit dipole   5. CONCLUSION In this paper, a 2D numerical simulation to the mixed discrete-continuum model of the biomolecular electrostatics based on boundary integral equation approach is presented.  Some computational results are presented to show the effectiveness of such model.  More work is required to extend the simulation to the 3D case which allows the simulation of a real biomolecular structure.  However, the tradeoff is that the problem size will become much larger.  Thus, it is necessary to develop a faster and more memory-efficient numerical algorithm.  6. REFERENCES [1] F.Fogolari, A.Brigo and H.Molinari. The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology. Journal of Molecular Biology, 15:377-392, 2002 [2] Shihhsien S.Kuo, Michael D. Altman, Jaydeep P.Bardhan, Bruce Tidor and Jacob K. White. Fast Methods for Simulation of Biomolecule Electrostatics, IEEE/ACM International Conference on 10-14 Nov. 2002, 466-473, 2002 [3] A. J. Bordner, G. A. Huber. Boundary Element Solution of the Linear Poisson-Boltzmann Equation and a Multipole Method for the Rapid Calculation of Forces on Macromolecules in Solution. Journal of Computational Chemistry, 24:353-367, 2003.  +1 -1 O rC 0.1 0.5 0.9 1.1 r=1 : Location of the dipole : Locations of the point where potential is to be evaluated 

Numerical Study of the Poisson-Boltzmann Equation for Biomolecular Electrostatics

http://dspace.mit.edu/bitstream/1721.1/30372/1/HPCES001.pdf

Numerical Study of the Poisson-Boltzmann Equation for Biomolecular Electrostatics

Abstract

Similar works

Full text

Available Versions

DSpace@MIT