Volume 13 Supplement 4
Italian Society of Bioinformatics (BITS): Annual Meeting 2011
Bluues: a program for the analysis of the electrostatic properties of proteins based on generalized Born radii
 Federico Fogolari^{1, 2}Email author,
 Alessandra Corazza^{1, 2},
 Vijaylakshmi Yarra^{1},
 Anusha Jalaru^{1},
 Paolo Viglino^{1, 2} and
 Gennaro Esposito^{1, 2}
DOI: 10.1186/1471210513S4S18
© Fogolari et al.; licensee BioMed Central Ltd. 2012
Published: 28 March 2012
Abstract
Background
The PoissonBoltzmann (PB) equation and its linear approximation have been widely used to describe biomolecular electrostatics. Generalized Born (GB) models offer a convenient computational approximation for the more fundamental approach based on the PoissonBoltzmann equation, and allows estimation of pairwise contributions to electrostatic effects in the molecular context.
Results
We have implemented in a single program most common analyses of the electrostatic properties of proteins. The program first computes generalized Born radii, via a surface integral and then it uses generalized Born radii (using a finite radius test particle) to perform electrostic analyses. In particular the ouput of the program entails, depending on user's requirement:
1) the generalized Born radius of each atom;
2) the electrostatic solvation free energy;
3) the electrostatic forces on each atom (currently in a dvelopmental stage);
4) the pHdependent properties (total charge and pHdependent free energy of folding in the pH range 2 to 18;
5) the pKa of all ionizable groups;
6) the electrostatic potential at the surface of the molecule;
7) the electrostatic potential in a volume surrounding the molecule;
Conclusions
Although at the expense of limited flexibility the program provides most common analyses with requirement of a single input file in PQR format. The results obtained are comparable to those obtained using stateoftheart PoissonBoltzmann solvers. A Linux executable with example input and output files is provided as supplementary material.
Background
Generalized Born models
Electrostatic effects arising due to the interaction of solute charges among themselevs and with solvent and ion charges, are of utmost importance for biomolecular structure and function. Continuum methods based on the PoissonBoltzmann equation have been widely used for calculating electrostatic effects [1–5]. In the last decades much interest has been devoted to developing fast approximations to the solution of the PoissonBoltzmann (PB) equation.
Onufriev and coworkers have developed an analytical approximation to the exact potential inside and outside the low dielectric region of a sphere [6], that performs surprisingly well also for the complex shape of proteins [7–10] and is therefore more general than the generalized Born (GB) models.
When only self and interaction energies and forces are sought, the most widely used approach is based on generalized Born (GB) models. Recent reviews summarize the approach and highlight most interesting recent results [3, 11–14].
Central to these models is the estimation of polarization charge contributions to: i) the selfenergy of each charge (embedded in the solute); ii) the interaction energy of each pair of charges.
This expression was found to be more accurate than the exact expression for two charges embedded in a conducting sphere, although the coefficient of 8.0 instead of 4.0 in the exponential has been suggested [16] and other similar forms have been proposed [17, 18].
Computation of Born radii
Born radii could be in principle computed by solving a system of linear equations for the polarization charges at the boundaries of the solute volume [19–25]. Under the approximation that the ionic solution provides complete screening, amounting to the assumption that the surface behaves like a grounded conductor, polarization charges at the surface are such that the integral of their electric field at any outer surface point is exactly the opposite of the source charge field. This approximation amounts to setting ε_{ out } = ∞. In a different context the same approximation has been proposed many years ago as the conducting surface model (COSMO) [26, 27] and successfully used since then. Under this condition it is possible to obtain simple formulae for the generalized Born radius for a low dielectric region delimited by a sphere or a plane. The reaction field satisfies Laplace equation inside the surface and the solution may be obtained solving the equation using suitable basis functions with Dirichlet boundary conditions at discretized points on the surface [28].
where $\overrightarrow{n}\left({\overrightarrow{r}}_{S}\right)$ is the unit vector normal to the surface at point ${\overrightarrow{r}}_{S}$ and pointing outwards.
The above expression which will be referred hereafter as GBR6 following Grycuk [30] and Tjong and Zhou [31, 32], has been analysed in detail by Mongan et al. [33]. Equation 7 was found to perform extremely well also for "very" non spherical shapes and in the context of biomolecular models. Occasional large differences with respect to PoissonBoltzmann calculations were found for inner cavities and local concavities at the surface, i.e. in conditions where a continuum model is anyway questionable. That study concluded that with a correction for a small systematic error, the GBR6 model is a sufficiently accurate continuum electrostatic model [33].
Tjong and Zhou [31, 32] used an analytical implementation for the estimation of Born radii based on volume integrals (corresponding to the above formula (7) that uses a surface integral instead) and showed its superior accuracy compared to other existing methods for a set of 55 very different proteins. In their approach it is made clear that standard estimations of Born radii compute in fact only geometric properties, and they provide empirical formulae for the correct Born energy depending on the inner and outer dielectric constants, ionic strength, total charge and number of atoms.
In the absence of a theory which could be cast in a fast computational framework, fitting approaches have been successfully followed and tested on large sets of proteins showing excellent agreement with PoissonBoltzmann calculation results. Romanov et al. [34] proposed that a linear combination of integrals I_{3} to I_{6}, in the present notation, and a constant term could fit the selfpolarization energy and thus be used to compute generalized Born radii (see also the discussion by Mongan et al. [33]).
Applications of the generalized Born model
Based on the correct estimation of Born radii, the solvation contribution to the interaction between any two charges may be computed by using equation (2). Derivation of equation (2) with respect to atomic positions gives the electrostatic solvation forces acting on atoms. The implicit dependence of Born radii on atomic positions makes computation of forces far from trivial [35–40] for the generalized Born model as well as for the parent PoissonBoltzmann model where the boundaries depend on atomic positions [41]. The possibility of computing energies and forces faster with respect to the reference PoissonBoltzmann equation has made the Generalized Born model the choice of election for implicit solvent molecular dynamics simulations. Also, the computation of pairwise solvation energies allows for fast computation of pKa of multiple titrating groups as we discuss in the Methods section.
Applications of generalized Born model (and other implicit solvent methods) have been reviewed elsewhere [3, 11–14]. At variance with the reference PoissonBoltzmann model, the computation of electrostatic potential in space is outside the scope of the Generalized Born model where only interactions are considered thorugh equation 2. Here we use a finite radius test charge in order to define a potential within the frame of the generalized Born model.
Aim of this work
 1)
the generalized Born radius of each atom;
 2)
the solvation electrostatic free energy;
 3)
the electrostatic forces on each atom (the theory of electrostatic forces based on surface integrals is developed here and implementation is still at a developmental stage);
 4)
the pHdependent properties (total charge and pHdependent free energy of folding in the pH range 2 to 18);
 5)
the pKa of all ionizable groups;
 6)
the electrostatic potential at the surface of the molecule;
 7)
the electrostatic potential in a volume surrounding the molecule.
Results and discussion
Generalized Born radii
Self energies
model  RMSD  corr. coef. 

GBR6 SAS  10.0  0.995 
LC5 SAS  3.5  0.996 
LC9 SAS  3.5  0.996 
GBR6 MS  57.6  0.949 
LC5 MS  23.9  0.950 
LC9 MS  22.7  0.955 
Electrostatic solvation free energy
Electrostatic solvation forces
Since solvation energies depend, according to equation 2, on atomic positions explicitly and implicitly through Born radii, the computation of solvation electrostatic forces is quite complex, and strongly dependent on the interface model chosen [38–40] (see also for a general discussion [41]).
In order to estimate how important are effects due to the dependence of Born radii on atomic positions, solvation electrostatic forces have been computed for each atom of the 55 proteins as the derivative of the solvation energy under the approximation that Born radii are constant. Albeit approximate this way of computing electrostatic forces preserve by definition zero total electrostatic force as expected for isotropic media and as found by the correct expression for the force [41]. Due to the different way of computing forces we did not attempt comparing ionic boundary, dielectric boundary and charge times electrostatic field components of the electrostatic force, but we rather compared the total solvation force. The results are not very accurate with the average root mean square deviation, with respect to the forces computed using APBS using the same parameters, equal to 2.7 kJ/(Åmol) compared to an average square root value of the force of 5.0 kJ/(Åmol). The correlation coefficient is 0.84.
pHdependent properties (total charge, pHdependent free energy of folding, and pKa of ionizable groups
The test set provided by Gunner and coworkers has been used to test the prediction of pKa's in proteins [43]. We compared the results obtained with the results obtained by running the program Propka2.0 [44–47]. The latter program is very fast and provides predictions whose accuracy is comparable to that of more computationally intensive programs. We chose this program as a reference because it is widely used and it is seemingly the fastest available with a very good accuracy. A wide range of programs and methods have been compared recently [48] and the reader is referred to that work for more extensive comparisons of existing softwares.
Predicted vs. experimental pKa shifts
BLUUES  Propka v. 2.0  

RMSD  corr. coef.  RMSD  corr. coef.  
ASP  0.78  0.42  1.50  0.27 
GLU  0.90  0.54  1.75  0.24 
HIS  0.98  0.74  1.50  0.59 
LYS  0.79  0.45  0.79  0.25 
TYR  1.78  0.41  2.12  0.03 
NTR  0.54  0.93  0.96  1.00 
CTR  1.01  0.14  1.06  0.36 
All  1.00  0.51  1.59  0.36 
On the other hand some scaling of contributions is done and therefore large deviations from experimental results are not found. In this respect, Propka that uses ten adjustable parameters and other parameters chosen for best performance [46] apparently does not prevent very large shifts to be predicted. As a consequence the global RMSD from experimental data is large and could thus easily be reduced. The execution time of Propka is in the range of seconds, while our program runs in minutes. No optimization or approximation has been implemented as yet.
Surface electrostatic potential
The evaluation of the potential in regions outside the molecule is beyond the scope of the generalized Born approach which focuses instead on self and interaction energies. This is at variance with approaches that approximate the potential computed using the PoissonBoltzmann equation. In particular Onufriev and coworkers [7, 8] have found an approximation (not simply amounting to a truncation) based on Kirkwood's series expansion of the solution for a sphere [6]. Their analytical model performs surprisingly well for a large set of proteins and enables fast calculation of the potential at the surface and in the volume surrounding the molecule.
Electrostatic potential in a volume surrounding the molecule
Conclusions
A program for the analysis of the electrostatic properties of proteins based on the surface integral computation of generalized Born radii has been presented. Further work will be devoted to improve the efficiency of calculations in particular for what concerns electrostatic solvation forces. The program, together with examples is given as supplementary material (Additional file 1).
Methods
Reference APBS and UHBD calculations and parameters
Reference electrostatic calculations were performed with the programs APBS [53] and UHBD [54, 55]. Except where noted, standard parameters were used: inner dielectric constant 1.0, outer dielectric constant 78.54, temperature 298.15 K, ionic strength 0.150 M. In APBS the mesh of the grid was 0.25 Å which resulted in very large grid size (257^{3} points). In UHBD, used for surface potential calculations, a focusing procedure was used with the final mesh size of 0.8 Å.
The boundary surface was defined in the two programs in order to match as close as possible the definition used in our program, i.e. the molecular surface for atoms with radii inflated by the radius of the solvent. The surface point density was 10 Å^{2}.
Except where noted the same temperature, dielectric constant and ionic strength were used in our program.
Solvent accessible surface
Solvent accessible surface points and surface normal vectors have been generated by considering the van der Waals sphere of each atom inflated by the radius of the solvent (e.g., for an atom of radius 1.9Å and a solvent radius of 1.4Å the inflated radius was 3.3Å).
n_{ pt } was 200 corresponding for an atom with radius 1.9Å to a density of 4.4 pts Å^{2}.
Atoms are mapped on a grid and the position of each surface point is compared only with the position of the list of atoms associated with neighboring grid points. This procedure was found to be efficient and robust without any failure.
Each surface point within the inflated van der Waals volume of a different atom is considered buried. All surface integrals have been performed as the corresponding finite sums over exposed (non buried) surface points.
Generalized Born radii
Reference Born radii have been obtained from PoissonBoltzmann selfpolarization energies computed using APBS [53] according to eq. 1.
where I_{ n } is defined by equation 8, or by fitting the reference solvation self energy by linear combinations of the solvation energies corresponding to different α_{ n }'s (n ranging here from 3 to 6) as done by Romanov et al. [34]. Fitting was necessary because the published coefficients did not provide good results, as a result of the different procedures used to compute surfaces and solvation free energies.
An equation for the generalized Born radius
where the constant K depends on the coefficient used in the exponential and is 0.75 for the commonly adopted Still formula.
For the conducting grounded surface approximation the boundary condition is that the Born radius approaches 0 as the boundary surface is approached. It is easy to check that this equation is satisfied for the sphere and the plane where K is equal 1.
where J_{0}(x) is the Bessel function of order zero. We may consider only order 0 functions due to cylindrical symmetry. The coefficients are obtained by imposing that at the boundary the potential (U_{ react } + U_{ source }) is equal to zero and using the orthogonality relation: $\int x{J}_{0}\left(kx\right){J}_{0}\left({k}^{\prime}x\right)=\frac{1}{k}\delta \left(k{k}^{\prime}\right)$. The parameters k are chosen such that J_{0}(kR) = 0.
The solution of equation (10) is obtained by bisection, guessing first a value at the midpoint of the axis of symmetry and integrating the equation using an adaptive RungeKutta fourth order method [56] and requiring that the value of the Born radius at the boundary be 0.
Electrostatic solvation energies in ionic solutions
where the functions ϕ_{i,j}are the Green's functions at points i and j but do not include the self potential of each charge in the inner dielectric medium and δ_{i,j}is Kronecker's delta.
Electrostatic solvation forces
For a unit charge placed at 7.5 Å from the center of a globular protein of radius 10 Å, assuming an inner dielectric constant 4.0, the resulting force is 14 kJ/(mol Å). An opposite force of the same magnitude is applied on the center of the sphere defining the boundaries. This simple example shows that forces arising from variations of Born radii are not negligible.
The derivative of ${I}_{5}^{i}$ is not straightforward as the domain of integration depends on atomic coordinates.
Note that in the above equations the operator $\overrightarrow{\nabla}$ operates only on r_{ S } coordinates, i.e. at the point where the integrand function is evaluated.
where the apices i and j in I_{5} means that the integral involves the vector from the surface points to the coordinates of atoms i and j respectively. The computation of forces via surface integral according to the above equations is not practical, unless a cutoff on the distance between pairs of atoms and surface points is used. Indeed for each of the $\frac{N\left(N+1\right)}{2}$ pairs of atoms j and i, where N is the number of atoms, surface integrals must be computed (see equation (21)). This is the result of the brute force application of our surface integral definition of Born radius. Work is under way in our laboratory to make this computation efficient.
pHdependent properties (total charge and pHdependent free energy of folding in the pH range 2 to 18)
The method described by Antosiewicz et al. [60] is implemented here with some modifications. PDB entries corresponding to proteins for which data are available in the pKa database made available by Gunner and coworkers [43] were prepared with pdb2pqr [61].
The method implemented here is a single conformer method, at variance with other methods like the QM/MM/TI method by Simonson et al. [62], MCCE [43] and the constant pH molecular dynamics method by McCammon and coworkers [63]. These methods are expected to be more accurate, at the expense of large computational time.
Approaches similar to that implemented here have been used by Onufriev and coworkers [64, 65] where a significant speedup for larger proteins, with respect to our implementation, is due to clustering of interacting sites.
Significant improvements over continuum electrostatics (and more in general molecular mechanics) methods have been achieved by parametrization of interactions [49, 50, 66–69].
Other approaches which use a continuum method have been proposed recently which achieve very good agreement with experimental values, by using a pentapetide reference model for the titratable group, and by optimizing the hydrogen positions [70].
Improvements similar to those mentioned in the above paragraphs will be implemented in the future in our program.
We summarize here the theory that is discussed at length by Antosiewicz et al. [60].
Energy of ionization
where z is the charge of the group upon ionization, R is the gas constant, T is the temperature (298.15 K in the present calculations).
It is assumed that pKa of protein ionizable groups in an unfolded protein are the same as those of model compound. We assume here that the the pKa are the following: 4.5 for Glu, 3.8 for Asp, 6.5 for His, 12.5 for Arg, 10.5 for Lys, 9.0 for Tyr, 8.0 for Cys, 8.0 for the Nterminal ammine and 3.2 for the Cterminal carboxyl.
where each term describes a contribution to the free energy of ionization: the first term is the free energy of ionization in the model compound, the second term is the difference in ionization selfenergy, the third term describes the difference in interaction of ionization charges with partial charges in the unionized protein and model compound and the fourth term describes the interaction among ionization charges. Following Antosiewicz et al. the inner dielectric constant is set to 20.0. A wide range of inner dielectric constants has been used for pKa calculations. When a single dielectric constant has been used it has been chosen mostly larger than 1, e.g. 4 in MCCE [43], 8 in EGAD [68] and 11 in GB/IMC [70]. A large value of the dielectric constant accounts in an empirical way for the missing degrees of freedom implied by a single conformer method (see the discussion by Schutz and Warshel [71]).
Monte Carlo simulation of the ionization state
The protonation state of the ionizable sidechains is changed according to a Monte Carlo procedure. Protonation or deprotonation are simulated by adding or subtracting a unit charge to an atom representing the ionizable group. The following atoms have been considered as representative for each ionizable group: CG for Asp, CD for Glu, NZ for Lys, CE for Arg, SG for Cys, OH for Tyr, N for the Nterminal ammine and C for the Cterminal carboxyl. For histidine protonation may occur alternatively at the atom NE2 and ND1.
Monte Carlo simulations are performed at 0.5 pH intervals between 2.0 and 18.0 and at each pH value average ionization values and average components of ionization free energy are computed. At each pH value the number of equilibration steps is 100 times the number of ionizable sites and the number of Monte Carlo steps is 1000 times the number of ionizable sites.
Based on Monte Carlo simulation the output of the program includes: i) the titration curve for each ionizable site; ii) the list of pKa values for each ionizable site; iii) the charge state of the protein; iv) the pHdependent component of the free energy of folding. Compared to the scheme of Antosiewicz et al. we are using the solvent accessible surface (as detailed in the subsection "Solvent accessible surface") instead of the molecular surface because results are more stable.
Heuristic corrections (detailed hereafter) to the scheme of Antosewicz et al [60] have been done mainly to remove large contributions which are typically overestimated due to neglection of molecular flexibility. For instance for surface residues large unfavorable interactions may be accomodated by conformational relaxation. Also large chargecharge interactions are likely to lead to partial unfolding that relieves the strong unfavorable energy. These considerations are consistent with the idea of using different dielectric constants for buried and solvent exposed residues, as suggested many years ago by Demchuck and Wade [72].
Scaling self energies
Here we consider as solvent exposed all those sites that have generalized Born radius smaller than 4.0 Å and buried those that have generalized Born larger than 7.0 Å. Selfenergy interactions (second righthand term in Equation 32) are multiplied by a factor 2.0 for buried sites and by 0.25 for exposed residues. For all intermediate situations the multiplying factor is a linear combination of the two extreme values.
Scaling background interactions
For exposed residues unfavorable background interactions (third righthand term in Equation 32) are weighted by the empirical function $\mathrm{exp}\left(\frac{E}{2RT}\right)$. The factor 2 in the Boltzmannlike function is purely empirical to downweigh a contribution that is likely to be relaxed by protein flexibility.
For each ionizable site the pKa is obtained as the midpoint of the titration curve. The value is further corrected by separating the shift in pKa due to the desolvation selfenergy $\left(dp{K}_{a}^{self}\right)$, background interactions $\left(dp{K}_{a}^{bg}\right)$ and sitecharge  sitecharge interactions $\left(dp{K}_{a}^{ii}\right)$.
Final scaling of contributions
The sitecharge  sitecharge interactions shift is scaled by the function $dp{K}_{a}^{ii}=\frac{dp{K}_{a}^{ii}\times 3.0}{\leftdp{K}_{a}^{ii}\right+3.0}$ which sets 3 pKa units as the maximum contribution arising from interactions between titratable sites, consistent with the idea that large unfavorable interactions lead to partial or global unfolding. Similarly, background interaction shifts are scaled by the function $dp{K}_{a}^{bg}=\frac{dp{K}_{a}^{bg}\times 5.0}{\leftdp{K}_{a}^{bg}\right+5.0}$, which sets 5 pKa units as the maximum contribution due to background interactions, in order to avoid few very large shifts.
Test sets and conditions
The 55 proteins selected by Tjong and Zhou [31] for testing the GBR6 model have been used here. Atoms of the PDB structures have been assigned charges and radii taken from the CHARMM forcefield, using the program PDB2PQR [61, 73]. The PQR file format is described at [52]. The radius of hydrogen atoms was reset to 1.0 Å in order to avoid numerical inaccuracies in the analysis linked to the small radius of polar hydrogens in the CHARMM forcefield. The surface of the molecule was defined alternatively as the surface accessible to the center of a solvent probe sphere with radius 1.5 Å or as the surface accessible by the surface of a solvent probe sphere with radius 1.5 Å, except where noted. The first type of surface, here referred to as solvent accessible surface, is computed by default by the program while the second type of surface, here referred to as molecular surface, was computed using the program MSMS [42] and given as input to the program in the form of a list of vertices and normal vectors.
For the computation of self energies 1000 sites were randomly chosen in the 55 proteins providing an unbiased sample of different environments.
The inner dielectric constant was assumed to be 1.0, the outer dielectric constant was assumed to be 78.54, the temperature is 298.15 K, ionic strength 0.150 M, except where noted. For the computation of pKa of ionizable sites in proteins we used the database developed by Gunner and coworkers [43] available at [74]. The test set includes pKa for structures with PDB id: 1a2p, 1a6k, 1beg, 1bf4, 1bhc, 1bus, 1bvi, 1bvv, 1cdc, 1coa, 1cvo, 1de3, 1dg9, 1dwr, 1egf, 1gb1, 1goa, 1h4g, 1ig5, 1igc, 1kf3, 1lni and 1lse.
Exact generalized Born radii for a conducting sphere
where α_{ i } is the generalized Born radius at point i.
The work of Mongan et al. [33] reports a compact form for the Born radii corresponding to the integrals above expressed in term of GBR6 multiplied by a correction factor.
List of abbreviations used
 PB:

PoissonBoltzmann
 GB:

Generalized Born
 COSMO:

Conducting Surface Model
 RMSD:

Root mean square deviation.
Declarations
Acknowledgements
AJ and VY have been supported by Italian Ministry of Education and University  Borse giovani ricercatori indiani 2008. Dr. I.M. Lait is gratefully aknowledged.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 4, 2012: Italian Bioinformatics Society (BITS): Annual Meeting 2011. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/13/S4.
Authors’ Affiliations
References
 Fogolari F, Brigo A, Molinari H: The PoissonBoltzmann equation for biomolecular electrostatics: a tool for structural biology. J Mol Recogn 2002, 15: 377–392. 10.1002/jmr.577View ArticleGoogle Scholar
 NevesPetersen MT, Petersen SB: Protein electrostatics: a review of the equations and methods used to model electrostatic equations in biomoleculesapplications in biotechnology. Biotechnol Annu Rev 2003, 9: 315–395.View ArticlePubMedGoogle Scholar
 Baker NA: Improving implicit solvent simulations: a Poissoncentric view. Curr Opin Struct Biol 2005, 15: 137–143. 10.1016/j.sbi.2005.02.001View ArticlePubMedGoogle Scholar
 Lu BZ, Zhou YC, Holst MJ, McCammon JA: Recent progress in numerical methods for the PoissonBoltzmann equation in biophysical applications. Commun Comput Phys 2008, 3: 973–1009.Google Scholar
 Wang J, Tan C, Tan YH, Lu Q, Luo R: PoissonBoltzmann solvents in molecular dynamics simulations. Commun Comput Phys 2008, 3: 1010–1031.Google Scholar
 Kirkwood JG: Theory of Solutions of Molecules Containing Widely Separated Charges with Special Application to Zwitterions. J Chem Phys 1934, 2: 351–361. 10.1063/1.1749489View ArticleGoogle Scholar
 Fenley AT, Gordon JC, Onufriev A: An analytical approach to computing biomolecular electrostatic potential. II. Derivation and analysis. J Chem Phys 2008, 129: 075101. 10.1063/1.2956497PubMed CentralView ArticlePubMedGoogle Scholar
 Gordon JC, Fenley AT, Onufriev A: An analytical approach to computing biomolecular electrostatic potential. II. Validation and applications. J Chem Phys 2008, 129: 075102. 10.1063/1.2956499PubMed CentralView ArticlePubMedGoogle Scholar
 Sigalov G, Fenley A, Onufriev A: Analytical electrostatics for biomolecules: Beyond the generalized Born approximation. J Chem Phys 2006, 124: 124902. 10.1063/1.2177251View ArticlePubMedGoogle Scholar
 Sigalov G, Scheffel P, Onufriev A: Incorporating variable dielectric environments into the generalized Born model. J Chem Phys 2005, 122: 094511. 10.1063/1.1857811View ArticlePubMedGoogle Scholar
 Bashford D, Case DA: Generalized born models of macromolecular solvation effects. Annu Rev Phys Chem 2000, 51: 129–152. 10.1146/annurev.physchem.51.1.129View ArticlePubMedGoogle Scholar
 Feig M, Brooks CL: Recent advances in the development and application of implicit solvent models in biomolecule simulations. Curr Opin Struct Biol 2004, 14: 217–224. 10.1016/j.sbi.2004.03.009View ArticlePubMedGoogle Scholar
 Koehl P: Electrostatics calculations: latest methodological advances. Curr Opin Struct Biol 2006, 16: 142–151. 10.1016/j.sbi.2006.03.001View ArticlePubMedGoogle Scholar
 Chen J, Brooks CL, Khandogin J: Recent advances in implicit solventbased methods for biomolecular simulations. Curr Opin Struct Biol 2008, 18: 140–148. 10.1016/j.sbi.2008.01.003PubMed CentralView ArticlePubMedGoogle Scholar
 Still WC, Tempczyk A, Hawley RC, Hendrickson T: Semianalytical treatment of solvation for molecular mechanics and dynamics. J Am Chem Soc 1990, 112: 6127–6129. 10.1021/ja00172a038View ArticleGoogle Scholar
 Lee MS, Feig M, Salsbury FR, Brooks CL: New analytic approximation to the standard molecular volume and its application to generalzied Born calculations. J Comp Chem 2003, 24: 1348–1356. 10.1002/jcc.10272View ArticleGoogle Scholar
 Onufriev A, Case DA, Bashford D: Effective Born radii in the generalized Born approximation: the importance of being perfect. J Comput Chem 2002, 23: 1297–1304. 10.1002/jcc.10126View ArticlePubMedGoogle Scholar
 Lee MS, Salsbury FR, A OM: New analytic approximation to the standard molecular volume and its application to generalzied Born calculations. J Comp Chem 2004, 25: 1967–1978. 10.1002/jcc.20119View ArticleGoogle Scholar
 Miertus S, Scrocco E, Tomasi J: Electrostatic interaction of a solute with a continuum. A direct utilizaion of AB initio molecular potentials for the prevision of solvent effects. Chem Phys 1981, 55: 117–129. 10.1016/03010104(81)850902View ArticleGoogle Scholar
 Zauhar RJ, Morgan RS: A new method for computing the macromolecular electric potential. J Mol Biol 1985, 186: 815–820. 10.1016/00222836(85)903997View ArticlePubMedGoogle Scholar
 Rashin AA, Namboodiri K: A simple method for the calculation of hydration enthalpies of polar molecules with arbitrary shapes. J Phys Chem 1987, 91(23):6003–6012. 10.1021/j100307a038View ArticleGoogle Scholar
 Yoon BJ, Lenhoff AM: A boundary element method for molecular electrostatics with electrolyte effects. J Comp Chem 1990, 11: 1080–1086. 10.1002/jcc.540110911View ArticleGoogle Scholar
 Lu B, Cheng X, Huang J, McCammon JA: An Adaptive Fast Multipole Boundary Element Method for PoissonBoltzmann Electrostatics. J Chem Theory Comp 2009, 5: 1692–1699. 10.1021/ct900083kView ArticleGoogle Scholar
 Bardhan JP: Interpreting the Coulombfield approximation for generalizedBorn electrostatics using boundaryintegral equation theory. J Chem Phys 2008, 129: 144105. 10.1063/1.2987409View ArticlePubMedGoogle Scholar
 Bardhan JP: Numerical solution of boundaryintegral equations for molecular electrostatics. J Chem Phys 2009, 130: 094102. 10.1063/1.3080769View ArticlePubMedGoogle Scholar
 Klamt A, Schüürmann G: COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J Chem Soc Perkin Trans 1993, 2: 799–805.View ArticleGoogle Scholar
 Klamt A, Eckert F, Arlt W: COSMORS: an alternative to simulation for calculating thermodynamic properties of liquid mixtures. Ann Rev Chem Biomol Eng 2010, 1: 101–120. 10.1146/annurevchembioeng073009100903View ArticleGoogle Scholar
 Jackson JD: Classical Electrodynamics Third Edition. New York, NY, USA: Wiley and sons; 1998.Google Scholar
 Ghosh A, Rapp CS, Friesner RA: Generalized Born Model Based on a Surface Integral Formulation. J Phys Chem B 1998, 102: 10983–10990. 10.1021/jp982533oView ArticleGoogle Scholar
 Grycuk T: Deficiency of the Coulombfield approximation in the generalized Born model: An improved formula for Born radii evaluation. J Chem Phys 2003, 119: 4817–4826. 10.1063/1.1595641View ArticleGoogle Scholar
 Tjong H, Zhou HX: GBr^{ 6 }: a parametrization free, accurate, analytical generalized Born method. J Phys Chem 2007, 111: 3055–3061. 10.1021/jp066284cView ArticleGoogle Scholar
 Tjong H, Zhou HX: GBr^{ 6 }NL: a generalized Born method for accurately reproducing solvation energy of the nonlinear PoissonBoltzmann equation. J Chem Phys 2007, 126: 195102. 10.1063/1.2735322View ArticlePubMedGoogle Scholar
 Mongan J, SvrcekSeiler WA, Onufriev A: Analysis of integral expressions for effective Born radii. J Chem Phys 2007, 127: 185101. 10.1063/1.2783847View ArticlePubMedGoogle Scholar
 Romanov AN, Jabin SN, Martynov YB, Sulimov AV, Grigoriev FV, Sulimov VB: Surface generalized Born method: a simple, fast and precise implicit solvent model beyond the Coulomb approximation. J Phys Chem 2004, 108: 9323–9327. 10.1021/jp046721sView ArticleGoogle Scholar
 Hawkins GD, Cramer CJ, Truhlar DG: Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J Phys Chem 1996, 100: 19824–19839. 10.1021/jp961710nView ArticleGoogle Scholar
 Onufriev A, Bashford D, Case DA: Modification of the generalised Born model suitable for macromolecules. J Phys Chem 2000, 104: 3712–3720. 10.1021/jp994072sView ArticleGoogle Scholar
 Onufriev A, Bashford D, Case DA: Exploring protein native states and largescale conformational changes with a modified generalized Born model. Proteins: Struct Func Gen 2004, 55: 383–394. 10.1002/prot.20033View ArticleGoogle Scholar
 Dominy BN, Brooks CL: Development of a Generalized Born Model Parametrization for Proteins and Nucleic Acids. J Phys Chem B 1999, 103: 3765–3773. 10.1021/jp984440cView ArticleGoogle Scholar
 Im W, Lee MS, Brooks CL: Generalized born model with a simple smoothing function. J Comp Chem 2003, 24: 1691–1702. 10.1002/jcc.10321View ArticleGoogle Scholar
 Yu Z, Jacobson MP, Friesner RA: What role do surfaces play in GB models? A newgeneration of surfacegeneralized born model based on a novel gaussian surface for biomolecules. J Comp Chem 2006, 27: 72–89. 10.1002/jcc.20307View ArticleGoogle Scholar
 Gilson MK, Davis ME, Luty BA, McCammon JA: Computation of electrostatic forces on solvated molecules using the PoissonBoltzmann equation. J Phys Chem 1993, 97: 3591–3600. 10.1021/j100116a025View ArticleGoogle Scholar
 Sanner M, Spehner JC, Olson A: Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 1996, 38: 305–320. 10.1002/(SICI)10970282(199603)38:3<305::AIDBIP4>3.0.CO;2YView ArticlePubMedGoogle Scholar
 Song Y, Mao J, Gunner MR: MCCE2: Improving protein pKa calculations with extensive side chain rotamer sampling. J Comp Chem 2009, 30: 2231–2247.Google Scholar
 Li H, Robertson AD, Jensen JH: Very fast empirical prediction and rationalization of protein pKa values. Proteins 2005, 61: 704–721. 10.1002/prot.20660View ArticlePubMedGoogle Scholar
 Bas DC, Rogers DM, Jensen JH: Very fast prediction and rationalization of pKa values for proteinligand complexes. Proteins 2008, 73: 765–783. 10.1002/prot.22102View ArticlePubMedGoogle Scholar
 Olsson MHM, Søndergaard CR, Rostkowski M, Jensen JH: PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions. J Chem Theory Comput 2011, 7: 525–537. 10.1021/ct100578zView ArticleGoogle Scholar
 Søndergaard CR, Olsson MHM, Rostkowski M, Jensen JH: Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values. J Chem Theory Comput 2011, 7: 2284–2295. 10.1021/ct200133yView ArticleGoogle Scholar
 Stanton CL, Houk KN: benchmarking pKa prediction methods for residues in proteins. J Chem Theory Comput 2008, 4: 951–966. 10.1021/ct8000014View ArticleGoogle Scholar
 Huang RB, Du QS, Wang CH, Liao SM, Chou KC: A fast and accurate method for predicting pKa of residues in proteins. Protein Eng Des Sel 2010, 23: 35–42. 10.1093/protein/gzp067View ArticlePubMedGoogle Scholar
 Burger SK, Ayers PW: A parameterized, continuum electrostatic model for predicting protein pKa values. Proteins 2011, 79: 2044–2052. 10.1002/prot.23019View ArticlePubMedGoogle Scholar
 Humphrey W, Dalke A, Schulten K: VMD Visual Molecular Dynamics. J Mol Graph 1996, 14: 33–38. 10.1016/02637855(96)000185View ArticlePubMedGoogle Scholar
 APBS  File Formats[Http://www.poissonboltzmann.org/fileformats/]
 Baker NA, Sept D, Joseph S, Holst MJ, McCammon JA: Electrostatics of nanosystems: application to microtubules and the ribosome. Proc Natl Acad Sci USA 2001, 98: 10037–10041. 10.1073/pnas.181342398PubMed CentralView ArticlePubMedGoogle Scholar
 Madura JD, Davis ME, Gilson MK, Wade R, Luty BA, McCammon JA: Biological applications of electrostatics calculations and Brownian dynamics simulations. Rev Comp Chem 1994, 5: 229–267.Google Scholar
 Madura JD, Briggs JM, Wade R, Davis ME, Luty BA, Ilin A, Antosiewicz JA, Gilson MK, Bagheri B, Ridgway Scott L, McCammon JA: Electrostatics and diffusion of molecules in solution: simulations with the University of Houston Brownian Dynamics program. Comput Commun Phys 1995, 91: 57–95. 10.1016/00104655(95)00043FView ArticleGoogle Scholar
 Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical recipes in C (2nd ed.) the art of scientific computing. Cambridge, UK: Cambridge University Press; 1992.Google Scholar
 Srinivasan J, Trevathan MW, Beroza P, Case DA: Application of a pairwise generalized Born model to proteins and nucleic acids: inclusion of salt effects. Theor Chem Acc 1999, 101: 426–434. 10.1007/s002140050460View ArticleGoogle Scholar
 Flanders H: Differentiation Under the Integral Sign. Am Math Month 1973, 80: 615–627. 10.2307/2319163View ArticleGoogle Scholar
 Aris R: Vectors, tensors, and the basic equation of fluid mechanics. New York, NY, USA: Dover Publications; 1962.Google Scholar
 Antosiewicz J, McCammon JA, Gilson MK: Prediction of pHdependent properties of proteins. J Mol Biol 1994, 238: 415–436. 10.1006/jmbi.1994.1301View ArticlePubMedGoogle Scholar
 Dolinsky TJ, Czodrowski P, Li H, Nielsen JE, Jensen JH, Klebe G, Baker NA: PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res 2007, 35: W522–525. 10.1093/nar/gkm276PubMed CentralView ArticlePubMedGoogle Scholar
 Simonson T, Carlsson J, Case DA: Proton binding to proteins: pK(a) calculations with explicit and implicit solvent models. J Am Chem Soc 2004, 126: 4167–4180. 10.1021/ja039788mView ArticlePubMedGoogle Scholar
 Mongan J, Case DA, McCammon JA: Constant pH molecular dynamics in generalized Born implicit solvent. J Comput Chem 2004, 25: 2038–2048. 10.1002/jcc.20139View ArticlePubMedGoogle Scholar
 Gordon JC, Myers JB, Folta T, Shoja V, Heath LS, Onufriev A: H++: a server for estimating pKas and adding missing hydrogens to macromolecules. Nucleic Acids Res 2005, 33: W368–371. 10.1093/nar/gki464PubMed CentralView ArticlePubMedGoogle Scholar
 Myers J, Grothaus G, Narayanan S, Onufriev A: A simple clustering algorithm can be accurate enough for use in calculations of pKs in macromolecules. Proteins 2006, 63: 928–938. 10.1002/prot.20922View ArticlePubMedGoogle Scholar
 Mehler EL, Guarnieri F: A selfconsistent, microenvironment modulated screened coulomb potential approximation to calculate pHdependent electrostatic effects in proteins. Biophys J 1999, 77: 3–22. 10.1016/S00063495(99)768682PubMed CentralView ArticlePubMedGoogle Scholar
 Wisz MS, Hellinga HW: An empirical model for electrostatic interactions in proteins incorporating multiple geometrydependent dielectric constants. Proteins 2003, 51: 360–377. 10.1002/prot.10332View ArticlePubMedGoogle Scholar
 Pokala N, Handel TM: Energy functions for protein design I: efficient and accurate continuum electrostatics and solvation. Protein Sci 2004, 13: 925–936. 10.1110/ps.03486104PubMed CentralView ArticlePubMedGoogle Scholar
 He Y, Xu J, Pan XM: A statistical approach to the prediction of pK(a) values in proteins. Proteins 2007, 69: 75–82. 10.1002/prot.21478View ArticlePubMedGoogle Scholar
 Spassov VZ, Yan L: A fast and accurate computational approach to protein ionization. Protein Sci 2008, 17: 1955–1970. 10.1110/ps.036335.108PubMed CentralView ArticlePubMedGoogle Scholar
 Schutz CN, Warshel A: What are the dielectric "constants" of proteins and how to validate electrostatic models? Proteins 2001, 44: 400–417. 10.1002/prot.1106View ArticlePubMedGoogle Scholar
 Demchuk E, Wade RC: Improving the Continuum Dielectric Approach to Calculating pKas of Ionizable Groups in Proteins. J Phys Chem 1996, 100(43):17373–17387. 10.1021/jp960111dView ArticleGoogle Scholar
 Dolinsky TJ, Nielsen JE, McCammon JA, Baker NA: PDB2PQR: an automated pipeline for the setup of PoissonBoltzmann electrostatics calculations. Nucleic Acids Res 2004, 32: W665–667. 10.1093/nar/gkh381PubMed CentralView ArticlePubMedGoogle Scholar
 Protein pKa database[Http://pka.engr.ccny.cuny.edu]
 Wolfram Mathematica online integrator[Http://integrals.wolfram.com]
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.