Abstract
Calculations of electrostatic potential and solvation free energy of macromolecules are essential for understanding the mechanism of many biological processes. In the classical implicit solvent Poisson–Boltzmann (PB) model, the macromolecule and water are modeled as two-dielectric media with a sharp border. However, the dielectric property of interior cavities and ion-channels is difficult to model realistically in a two-dielectric setting. In fact, the detection of water molecules in a protein cavity remains to be an experimental challenge. This introduces an uncertainty, which affects the subsequent solvation free energy calculation. In order to compensate this uncertainty, a novel super-Gaussian dielectric PB model is introduced in this work, which devices an inhomogeneous dielectric distribution to represent the compactness of atoms and characterizes empty cavities via a gap dielectric value. Moreover, the minimal molecular surface level set function is adopted so that the dielectric profile remains to be smooth when the protein is transferred from water phase to vacuum. An important feature of this new model is that as the order of super-Gaussian function approaches the infinity, the dielectric distribution reduces to a piecewise constant of the two-dielectric model. Mathematically, an effective dielectric constant analysis is introduced in this work to benchmark the dielectric model and select optimal parameter values. Computationally, a pseudo-time alternative direction implicit (ADI) algorithm is utilized for solving the super-Gaussian PB equation, which is found to be unconditionally stable in a smooth dielectric setting. Solvation free energy calculation of a Kirkwood sphere and various proteins is carried out to validate the super-Gaussian model and ADI algorithm. One macromolecule with both water filled and empty cavities is employed to demonstrate how the cavity uncertainty in protein structure can be bypassed through dielectric modeling in biomolecular electrostatic analysis.
Similar content being viewed by others
References
Abrashkin A, Andelman D, Orland H (2007) Dipoloar Poisson–Boltzmann equation: ions and dipoles close to charge interface. Phys Rev Lett 99:077801
Alexov EG, Gunner MR (1997) Incorporating protein conformational flexibility into the calculation of pH-dependent protein properties. Biophys J 72:2075–2093
Alexov EG, Gunner MR (1999) Calculated protein and proton motions coupled to electron transfer: electron transfer from QA- to QB in bacterial photosynthetic reaction centers. Biochemistry 38:8253–8270
Baker NA, Sept D, Joseph S, Holst MJ, McCammon JA (2001) Electrostatics of nanosystems: application to microtubules and the ribosome. Proc Natl Acad Sci 98:10037–10041
Bates P, Wei GW, Zhao S (2008) Minimal molecular surfaces and their applications. J Comput Chem 29:380–391
Bates PW, Chen Z, Sun YH, Wei GW, Zhao S (2009) Geometric and potential driving formation and evolution of biomolecular surfaces. J Math Biol 59:193–231
Blinn JF (1982) A generalization of algegraic surface drawing. ACM Trans Graph 1:235–256
Bohinc K, Bossa GV, May S (2017) Incorporation of ion and solvent structure into mean-field modeling of the electric double layer. Adv Colloid Interface Sci 249:220–233
Chakravorty A, Jia Z, Li L, Zhao S, Alexov E (2018a) Reproducing the ensemble average polar solvation energy of a protein from a single structure: Gaussian-based smooth dielectric function for macromolecular modeling. J Chem Theory Comput 14:1020–1032
Chakravorty A, Jia Z, Peng Y, Tajielyato N, Wang L, Alexov E (2018b) Gaussian-based smooth dielectric function: a surface-free approach for modeling macromolecular binding in solvents. Front Mol Biosci 5:25
Che J, Dzubiella J, Li B, McCammon JA (2008) Electrostatic free energy and its variations in implicit solvent models. J Phys Chem B 112:3058–3069
Chen M, Lu B (2011) TMSmesh: a robust method for molecular surface mesh generation using a trace technique. J Chem Theory Comput 7:203–212
Chen DA, Chen Z, Chen CJ, Geng WH, Wei GW (2011) Software news and update MIBPB: a software package for electrostatic analysis. J Comput Chem 32:756–770
Cheng L-T, Dzubiella J, McCammon JA, Li B (2007) Application of the level-set method to the solvation of nonpolar molecules. J Chem Phys 127:084503
Connolly ML (1983) Analytical molecular surface calculation. J Appl Crystallogr 16:548–558
Dai S, Li B, Liu J (2018) Convergence of phase-field free energy and boundary force for molecular solvation. Arch Ration Mech Anal 227:105–147
Deng W, Xu J, Zhao S (2018) On developing stable finite element methods for pseudo-time simulation of biomolecular electrostatics. J Comput Appl Math 330:456–474
Duncan BS, Olson AJ (1993) Shape analysis of molecular surfaces. Biopolymers 33:231–238
Geng WH, Zhao S (2013) Fully implicit ADI schemes for solving the nonlinear Poisson–Boltzmann equation. Mol Math Biophys 1:109–123
Geng W, Zhao S (2017) A two-component matched interface and boundary (MIB) regularization for charge singularity in implicit solvation. J Comput Phys 351:25–39
Giard J, Macq B (2010) Molecular surface mesh generation by filtering electron density map. Int J Biomed Imaging 2010:923780
Grant JA, Pickup B (1995) A Gaussian description of molecular shape. J Phys Chem 99:3503–3510
Grant JA, Pickup BT, Nicholls A (2001) A smooth permittivity function for Poisson–Boltzmann solvation methods. J Comput Chem 22:608–640
Hage KE, Hedin F, Gupta PK, Meuwly M, Karplus M (2018) Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size. eLife 7:e35560
Hammel M (2012) Validation of macromolecular flexibility in solution by small-angle X-ray scattering (SAXS). Eur Biophys J 41:789–799
Hu L, Wei GW (2012) Nonlinear Poisson equation for heterogeneous media. Biophys J 103:758–766
Huggins DJ (2015) Quantifying the entropy of binding for water molecules in protein cavities by computing correlations. Biophys J 108:928–936
Im W, Beglov D, Roux B (1998) Continuum solvation model: computation of electrostatic forces from numerical solutions to the Poisson–Boltzmann equation. Comput Phys Commun 111:59–75
Jia Z, Li L, Chakravorty A, Alexov E (2017) Treating ion distribution with Gaussian-based smooth dielectric function in DelPhi. J Comput Chem 38:1974–1979
Koehl P, Orland H, Delarue M (2009) Beyond the Poisson–Boltzmann model: modeling biomolecular-water and water-water interactions. Phys Rev Lett 102:087801
Kokkinidis M, Glykos NM, Fadouloqlou VE (2012) Protein flexibility and enzymatic catalysis. Adv Protein Chem Struct Biol 87:181–218
Lee B, Richards FM (1973) Interpretation of protein structure: estimation of static accessibility. J Mol Biol 55:379–400
Li C, Li L, Zhang J, Alexov E (2012) Highly efficient and exact method for parallelization of gridbased algorithms and its implementation in DelPhi. J Comput Chem 33:1960–1966
Li C, Li L, Petukh M, Alexov E (2013a) Progress in developing Poisson–Boltzmann equation solvers. Mol Based Math Biol 1:42–62
Li L, Li C, Zhang Z, Alexov E (2013b) On the dielectric “constant” of proteins: smooth dielectric function for macromolecular modeling and its implementation in DelPhi. J Chem Theory Comput 9:2126–2136
Li L, Li C, Alexov E (2014) On the modeling of polar component of solvation energy using smooth Gaussian-based dielectric function. J Theory Comput Chem 13:1440002
Li L, Wang L, Alexov E (2015) On the energy components governing molecular recognition in the framework of continuum approaches. Front Mol Biosci 2:5
Lu BZ, Zhou YC, Holst MJ, McCammon JA (2008) Recent progress in numerical methods for the Poisson–Boltzmann equation in biophysical applications. Commun Comput Phys 3:973–1009
Mengistu DH, Bohing K, May S (2009) Poisson–Boltzmann model in a solvent of interacting Langevin dipoles. EPL (Europhys Lett) 88:14003
Ng J, Vora T, Krishnamurthy V, Chung S-H (2008) Estimating the dielectric constant of the channel protein and pore. Eur Biophys J 37:213–222
Nymeyer H, Zhou HX (2008) A method to determine dielectric constants in nonhomogeneous systems, application to biological membranes. Biophys J 94:1185–1193
Pang X, Zhou HX (2013) Poisson–Boltzmann calculations: van der Waals or molecular surface? Commun Comput Phys 13:1–12
Qiao ZH, Li ZL, Tang T (2006) A finite difference scheme for solving the nonlinear Poisson–Boltzmann equation modeling charged spheres. J Comput Math 24:252–264
Quillin ML, Wingfield PT, Matthews BW (2006) Determination of solvent content in cavities in IL-1\(\beta \) using experimentally phased electron density. Proc Natl Acad Sci 103:19749–19753
Richards FM (1977) Areas, volumes, packing and protein structure. Annu Rev Biophys Bioeng 6:151–176
Sanner M, Olson A, Spehner J (1996) Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 38:305–320
Simonson T, Perahia D (1995) Internal interfacial dielectric properties of cytochrome c from molecular dynamics in aqueous solution. Proc Natl Acad Sci 92:1082–1086
Song X (2002) An inhomogeneous model of protein dielectric properties: intrinsic polarizabilities of amino acids. J Chem Phys 116:9359
Takano K, Yamagata Y, Yutani K (2003) Buried water molecules contribute to the conformational stability of a protein. Protein Eng 16:5–9
Tian W, Zhao S (2014) A fast ADI algorithm for geometric flow equations in biomolecular surface generation. Int J Numer Method Biomed Eng 30:490–516
Voges D, Karshikoff A (1998) A model of a local dielectric constant in proteins. J Chem Phys 108:2219
Wang L, Li L, Alexov E (2015a) pKa predictions for proteins RNAs and DNAs with the Gaussian dielectric function using DelPhiPKa. Proteins Struct Funct Bioinform 83:2186–2197
Wang L, Zhang M, Alexov E (2015b) DelPhiPKa Web Server: predicting pKa of proteins RNAs and DNAs. Bioinformatics 32:614–615
Warshel A, Russell ST (1984) Calculations of electrostatic interactions in biological systems and in solutions. Q Rev Biophys 17:283–422
Warshel A, Sharma PK, Kato M, Parson WW (2006) Modeling electrostatic effects in proteins. Biochim Biophys Acta 1764:1647–1676
Wilson L, Zhao S (2016) Unconditionally stable time splitting methods for the electrostatic analysis of solvated biomolecules. Int J Numer Anal Modell 13:852–878
Yu Z, Holst MJ, Cheng Y, McCammon JA (2008) Feature-preserving adaptive mesh generation for molecular shape modeling and simulation. J Mol Graph Modell 26:1370–1380
Zhang Y, Xu G, Bajaj C (2006) Quality meshing of implicit solvation models of biomolecular structures. Comput Aided Geom Des 23:510–530
Zhao S (2011) Pseudo-time-coupled nonlinear models for biomolecular surface representation and solvation analysis. Int J Numer Method Biomed Eng 27:1964–1981
Zhao S (2014) Operator splitting ADI schemes for pseudo-time coupled nonlinear solvation simulations. J Comput Phys 257:1000–1021
Zhao Y, Kwan YY, Che J, Li B, McCammon JA (2013) Phase-field approach to implicit solvation of biomolecules with Coulomb-field approximation. J Chem Phys 139:024111
Zhou YC, Zhao S, Feig M, Wei GW (2006) High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources. J Comput Phys 213:1–30
Acknowledgements
The research of Alexov was supported in part by the National Institutes of Health (NIH) Grant R01GM093937 and the National Science Foundation (NSF) Grant DMS-1812597. The research of Zhao was supported in part by the National Science Foundation (NSF) Grant DMS-1812903 and the Simons Foundation award 524151.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
Theorem
The density function for the \(i^{th}\) atom is defined by
where \(r_i\) and \(R_i\) are the center and radius of the \(i\mathrm{th}\) atom, respectively. Also, here \(\vec {r}\) is the position vector, \(\sigma \) is the relative variance and m is the power of super-Gaussian function. Suppose \(\sigma =1\) for simplicity. Next, the total density function of a biomolecular system is defined as \(g_0^s=1-\prod (1-g_i^s(\vec {r}))\) and the dielectric function of that system is modeled as
Here \(\epsilon _m\) and \(\epsilon _s\) are the dielectric constants of the solute and solvent respectively. Then we have that \(\displaystyle {\lim _{m\rightarrow \infty }\epsilon _G^s=\epsilon _2}\) at the solute and solvent regions where \(\epsilon _2\) is the dielectric function of the classical two-dielectric model.
Proof
Let us consider three cases where the position vector is either inside or outside the solute, or on the Van der Walls (VDW) molecular surface.
-
Case I: There exists an atom (say ith atom) such that \(|\vec {r}-r_i|< R_i\) or, \(\displaystyle {\frac{|\vec {r}-r_i|}{ R_i}}<1\). In this case \(\displaystyle {\lim _{m\rightarrow \infty }\Big (\frac{|\vec {r}-r_i|}{ R_i}\Big )^{2m}}=0\). Hence \(\displaystyle {\lim _{m\rightarrow \infty }\exp \Big [-\Big (\frac{|\vec {r}-r_i|}{ R_i}\Big )^{2m}\Big ]}=1\), which means \(g^s_i(\vec {r}) =1\) and \(g^s_0(\vec {r}) =1\). Therefore, if \(|\vec {r}-r_i|< R_i\) for some i (inside the VDW surface), \(\epsilon _G^s=\epsilon _m\).
-
Case II: For all atoms, we have \(|\vec {r}-r_i|> R_i\) or \(\displaystyle {\frac{|\vec {r}-r_i|}{ R_i}}>1\) for any i. In this case \(\displaystyle {\lim _{m\rightarrow \infty }\Big (\frac{|\vec {r}-r_i|}{ R_i}\Big )^{2m}}=\infty \). So, \(\displaystyle {\lim _{m\rightarrow \infty }\exp \Big [-\Big (\frac{|\vec {r}-r_i|}{ R_i}\Big )^{2m}\Big ]}=0\), which means that \(g^s_i(\vec {r}) =0\) for all i. Hence \(g^s_0(\vec {r}) =0\). Therefore, if \(|\vec {r}-r_i|> R_i\) for all i (outside the VDW surface), \(\epsilon _G^s=\epsilon _s\).
-
Case III: In the last case, the position vector \(\vec {r}\) has to be located on the VDW surface. Without the loss of generality, we assume that \(\vec {r}\) is on the sphere boundary of the \(i^{th}\) atom and does not locate inside any other atoms. So, we have \(|\vec {r}-r_i|= R_i\) or \(\displaystyle {\frac{|\vec {r}-r_i|}{ R_i}}=1\). And, for any \(j \ne i\), \(|\vec {r}-r_j|> R_j\). In this case \(\displaystyle {\lim _{m\rightarrow \infty }\exp \Big [-\Big (\frac{|\vec {r}-r_i|}{ R_i}\Big )^{2m}\Big ]}=\frac{1}{e}\), which means \(\displaystyle g^s_i(\vec {r}) =\frac{1}{e}\). For any \(j \ne i\), \( g^s_j(\vec {r}) =0\). Hence, \(\displaystyle g^s_0(\vec {r}) =\frac{1}{e}\). Therefore, on the VDW surface, we have \(\displaystyle {\epsilon _G^s=\epsilon _m g_0^s+\epsilon _s (1-g_0^s)=\epsilon _m \frac{1}{e}+\epsilon _s (1-\frac{1}{e})}=50.9375\) for \(\epsilon _m=1\) and \(\epsilon _s=80\). In all cases, the new dielectric model converges to a two-dielectric model based on the VDW surface
$$\begin{aligned} \lim _{m\rightarrow \infty }\epsilon _G^s(\vec {r})=\epsilon _2(\vec {r}) = \left\{ \begin{array}{ll} \epsilon _m, &{}\quad \vec {r}~~\text {is inside the VDW surface}\\ \epsilon _m/e + \epsilon _s(e-1)/e, &{}\quad \vec {r}~~\text {is on the VDW surface}\\ \epsilon _s, &{}\quad \vec {r}~~\text {is outside the VDW surface}. \\ \end{array} \right. \end{aligned}$$(35)
Rights and permissions
About this article
Cite this article
Hazra, T., Ahmed Ullah, S., Wang, S. et al. A super-Gaussian Poisson–Boltzmann model for electrostatic free energy calculation: smooth dielectric distribution for protein cavities and in both water and vacuum states. J. Math. Biol. 79, 631–672 (2019). https://doi.org/10.1007/s00285-019-01372-1
Received:
Revised:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00285-019-01372-1
Keywords
- Poisson–Boltzmann equation
- Gaussian dielectric model
- Minimal molecular surface
- Alternating direction implicit (ADI)
- Protein cavity
- Electrostatic free energy