Modification of the Generalized Born Model Suitable for Macromolecules
Alexey Onufriev, Donald Bashford, and David A. Case*
Department of Molecular Biology, The Scripps Research Institute, La Jolla, California 92037
ReceiVed: NoVember 16, 1999; In Final Form: January 27, 2000
The analytic generalized Born approximation is an efficient electrostatic model that describes molecules in
solution. Here it is modified to permit a more accurate description of large macromolecules, while its established
performance on small compounds is nearly unaffected. The modified model is also adapted to describe
molecules with an interior dielectric constant not equal to unity. The model is tested by computations of pK
shifts for a number of titratable residues in lysozyme, myoglobin, and bacteriorhodopsin. In general, except
for some deeply buried residues of bacteriorhodopsin, the results show reasonable agreement with both
experimental data and calculations based on numerical solution of the Poisson-Boltzmann equation. A very
close agreement between the two models is obtained in an application to the prediction of the pK shifts
associated with conformational change. The calculations based on this version of the generalized Born
approximation are much faster than finite difference solutions of the Poisson-Boltzmann equation, which
makes the present method useful for a variety of other applications where computational time is a critical
factor. The model may also be integrated into molecular dynamics programs to replace explicit solvent
simulations which are particularly time-consuming for large molecules.
1. Introduction
Over the past ten years classical electrostatic models based
upon numerical solution of the Poisson-Boltzmann (PB)
equation have been successfully applied to compute various
properties of macromolecules.
1-9
However, since solving a
system of partial differential equations is a computationally
costly procedure, the method may become quite time-consum-
ing, especially if applied to a large set of independent conforma-
tions of a macromolecule, or if it is incorporated into molecular
dynamics programs.
Recently, several fast analytic versions of the so-called
generalized Born (GB) approximation have been proposed as
an alternative to the computationally intensive Poisson-
Boltzmann approach.
10-15
For small molecules, the method has
been shown to reproduce solvation energies and individual
charge-charge interactions very well when compared to solu-
tions of the Poisson-Boltzmann equation.
11,12,14-16
However,
in similar comparisons on larger molecules, the agreement has
not been as close, the discrepancy generally being more pro-
nounced for molecules having large interior regions.
15,16
There,
the generalized Born approximation tends to overestimate the
solvation energy of deeply buried atoms and to underestimate
the interaction between them, as compared to numerical solution
of the Poisson-Boltzmann (PB) equation.
15,16
This appears to
be a general property of this type of GB approximation, in-
dependent of a particular parametrization. The resulting error
in the calculation of solvation energies is often acceptable, since
the atoms contributing the most to the solvent polarization
energy are the ones on the surface, and they are the most
accurately treated. Also, for an overall neutral molecule, the
individual errors in like-charge interactions largely cancel those
coming from opposite-charge ones.
15
However, when an ac-
curate estimate of the interactions between individual atoms
becomes important, such as in calculations of pK shifts, the
method is likely to show poorer performance.
12
For example,
the GB approximation was recently applied
17
to the estimation
of pK shifts of the active-site aspartic dyad in HIV-1 protease.
It was found that the set of atomic radii that worked well for
small molecules did not produce accurate estimates for the two
residues considered, and had to be adjusted. Proteins represent
a significant challenge to the GB approximation as they
generally contain biologically important functional groups which
may lie both on the surface and deep in the molecular interior.
Hence, a simple adjustment of atomic radii that may improve
the method’s performance for buried residues is likely to be
inadequate for the well-exposed ones, and vice versa. The theory
clearly needs to be improved to correctly describe interior
regions of large molecules while preserving its remarkable
accuracy already established for small compounds.
In this work we develop an analytic GB theory which is
suitable for proteins. We begin by re-deriving the GB ap-
proximation for the general case of a molecular interior dielectric
that is not necessarily equal to unity. We then formulate a simple
criterion which any GB theory must satisfy in order to correctly
describe interior atoms. Application of the criterion to the
pairwise GB method leads us to introduce a single new
parameter into the model, to account for the nonzero size of
solvent molecules. We optimize its value by comparing the
charge-charge interactions in myoglobin calculated by the GB
model with those obtained by numerical solution of the PB
equation. We then study the performance of the modified GB
model by computing pK
1/2
values for ionizable groups in
lysozyme, myoglobin, and bacteriorhodopsin. The predicted
values agree reasonably well with both experiment and calcula-
tions based upon solution of the PB equation. The agreement
between the two models becomes even better when the GB
approach is used to evaluate the difference in titration behavior
associated with conformational change. We also show that the
corrections have little effect on GB accuracy for small mol-
ecules, allowing one to use the same set of parameters for all
* Corresponding author. E-mail: case@scripps.edu.
3712 J. Phys. Chem. B 2000, 104, 3712-3720
10.1021/jp994072s CCC: $19.00 © 2000 American Chemical Society
Published on Web 03/22/2000
compounds, regardless of size. Our current implementation of
the GB method is significantly faster than its finite-difference
PB counterpart when applied to titration calculations on
medium-sized proteins such as lysozyme.
2. Theory
The continuum electrostatic model,
1
which is the starting point
for both the GB and PB models, subdivides the entire space
into two regions separated by a dielectric boundary: the solute,
characterized by a low dielectric value, and the solvent, which
has a high dielectric value. The electrostatic potential can then
be calculated directly, without any other approximations, by
solving the Poisson-Boltzmann equation numerically to a
desired degree of accuracy. Alternative methods, such as the
modified generalized Born model considered below, are based
on various further approximations that allow one to obtain an
analytical form of the potential or the electrostatic component
of the solvation free energy for a complex molecule. The
accuracy of the GB approximation can be assessed by comparing
its predictions to the corresponding PB results.
2.1. GB Model with Variable Interior Dielectric. We
represent each atom in the molecule as a sphere of radius F
i
with a charge q
i
at its center; the interior of the atom is assumed
to be filled uniformly with dielectric material of low dielectric
constant
p
. The molecule is surrounded by a solvent of a high
dielectric value
w
. The electrostatic component of the free
energy is the work W of creating a given charge distribution.
Its calculation is made nontrivial by the presence of surface
polarization that develops at the dielectric boundary between
the solute and the solvent if
p
*
w
. We single out the part of
W that is due to the polarization charges by considering a process
in which the molecule is transfered from a uniform medium of
dielectric value
p
into the solvent with dielectric constant
w
:
where W
0
is the energetic cost of creating the charge distribution
in a uniform dielectric
p
and ∆W is the energy cost of the
transfer, i.e., the solvation energy. The evaluation of W
0
is
straightforward since no surface polarization is present; it is
simply the Coulomb charge-charge interaction:
An efficient way to compute ∆W is given by the generalized
Born model which is conventionally applied
10,14-16
to the case
p
) 1. The theory is, however, not specific
19
to this particular
value of
p
, and one can easily modify it to describe the more
general case of
p
* 1:
where f
ij
gb
is a certain smooth function which is assumed to
depend only upon atomic radii F
i
and interatomic distances r
ij
.
A detailed analysis of the approximations on which the model
is based is presented in ref 19. Here we motivate eq 3 by
considering the exact analytical form of ∆W for a pair of atoms
i and j in the limiting cases of r
ij
f ∞ and r
ij
f 0, and requiring
that f
ij
gb
interpolate between the two extremes.
When the atoms are far apart, the work done on transferring
the pair from the medium of dielectric constant
p
into the one
of
w
must be
The first term is the Born solvation energy of transferring
the two independent spherical ions between the two media, and
the second term corresponds to the difference in the interaction
energy of two point charges in the two media. In the other
extreme r
ij
f 0, the two spheres merge into one, and if we also
assume that F
i
≈ F
j
, the system may be approximated by a sphere
of radius (F
i
F
j
)
1/2
with total charge (q
i
+ q
j
). Therefore, ∆W in
this case is just the Born solvation energy of transferring the
sphere between the two media:
The function f
ij
gb
in eq 3 should hence be such that f
ii
gb
)F
i
.As
r
ij
f ∞, we require f
ij
gb
f r
ij
, and as r
ij
f 0, f
ij
gb
f (F
i
F
j
)
1/2
. The
fact that these requirements on f
ij
gb
depend only on the geom-
etry of the system and not on the dielectric constant allows us
to use exactly the same f
ij
gb
here as in the standard formulation
of GB with
p
) 1, which is conventionally considered in
literature. This is important, since f
ij
gb
and the associated set of
atomic parameters have already been optimized by others. The
standard
11
form of f
ij
gb
(r
ij
) is particularly simple:
where R
i
g F
i
is the so-called effective Born radius which
replaces F
i
and accounts for the fact that atoms i and j may be
surrounded by neighbors that displace the solvent and therefore
decrease the polarization energy. A method to estimate R
i
, which
is the crux of our modification to the GB theory, is presented
in section 2.2. Larger effective radii result in smaller contribu-
tions to ∆W; their values reflect the extent of burial and depend
only on the mutual positions, radii, and types of the surrounding
atoms. In the limiting cases described above, and with other
neighbors far away, we require that R
i
f F
i
, and one can check
that eqs 4 and 5 are satisfied with f
ij
gb
from eq 6.
Given ∆W from eq 3, the total electrostatic work W of
creating the given charge distribution in the solvent is
2.2. Evaluation of Effective Radii. One way to estimate the
effective Born radius R
i
is to consider the change ∆W
i
in the
self-energy of atom i upon solvation. According to eqs 3
and 6:
On the other hand, this quantity can also be calculated directly
on the basis of classical electrostatics.
18,20
The work done on
creating a given charge distribution in an arbitrary dielectric
environment is
W ) W
0
+ ∆W (1)
W
0
)
1
2
∑
i*j
q
i
q
j
p
r
ij
(2)
∆W )-
1
2
(
1
p
-
1
w
)
∑
i,j
q
i
q
j
f
ij
gb
(r
ij
)
(3)
lim
r
ij
f∞
∆W )-
1
2
(
1
p
-
1
w
)(
q
i
2
F
i
+
q
j
2
F
j
)
-
(
1
p
-
1
w
)
q
i
q
j
r
ij
(4)
lim
r
ij
f0
∆W )-
1
2
(
1
p
-
1
w
)
(q
i
+ q
j
)
2
(F
i
F
j
)
-1/2
(5)
f
ij
gb
) [r
ij
2
+ R
i
R
j
exp(-r
ij
2
/4R
i
R
j
)]
1/2
(6)
W )-
1
2
(
1
p
-
1
w
)
∑
i,j
q
i
q
j
f
ij
gb
(r
ij
)
+
1
2
∑
i*j
1
p
q
i
q
j
r
ij
(7)
∆W
i
)-
1
2
(
1
p
-
1
w
)
q
i
2
R
i
(8)
W )
1
8π
∫
R
3
[DB(rb)]
2
(rb)
d
3
rb (9)
Modification of the Generalized Born Model J. Phys. Chem. B, Vol. 104, No. 15, 2000 3713
where DB(rb) is the dielectric displacement vector, and (rb)isthe
position-dependent dielectric constant. Therefore, the work ∆W
i
of transferring the atom i from a medium of uniform dielectric
constant,
p
, to the two-dielectric solute/solvent system is
where DB
i
(rb) is the total dielectric displacement due to charge i,
and
is the Coulomb field created by point charge q
i
in the uniform
dielectric environment. For convenience we have placed the
origin of the coordinate system at the center of atom i. So far,
eq 10 is an exact result. We now make the approximation:
The validity of the approximation above is considered in detail
in refs 19 and 20. On substituting DB
i
(rb) from eq 12 into eq 10,
the integrals over the solute volume cancel, while the solvent
volume integrals combine:
where R
i
is an effective radius defined by
The integration domain in the above equation can be changed
to the solute volume, which is computationally more convenient.
Note that
where θ(|rb| -F
i
) is the step function. Therefore, the effective
Born radius R
i
is also given by
As expected, R
i
)F
i
if the solute consists of only one atom i.
2.2.1. Integration Domain. Up until this point, the derivation
has not depended on the exact definition of the solute or solvent
volume. To proceed with the evaluation of the effective Born
radius, we need to specify the boundaries of the integration
domain in eq 16. In the PB calculations, it is common to define
the low dielectric solute as the region bounded by the molecular
surface.
21
With this definition, most small crevices between the
van der Waals (VDW) spheres of protein atoms fall in the solute
region because a solvent-sized probe cannot fit into them. This
implies that one could implement eq 16 by numerical integration
over the volume bounded by the molecular surface.
11,22
The
procedure would, however, undermine the main advantages of
the analytic GB model that were outlined in the Introduction,
as it would involve finding the molecular surface and performing
a costly numerical integration. To circumvent the problem, the
analytic GB method performs the integral in eq 16 over the
VDW spheres of atoms. This implies a definition of the solute
volume in terms of a set of spheres, rather than the molecular
surface, and allows one to obtain
14
an analytical form for R
i
.
The approach is known to work well in the case of small
molecules, but for large compounds such as proteins it has been
shown
15,16
to overestimate the solvation energy of deeply buried
atoms as well as to underestimate the interaction between them,
as compared to the PB calculations that use the molecular
surface-based definition of solute. One plausible reason for this
behavior is that, in the sphere-based GB method, only the inside
of each sphere has a low dielectric value and all crevices in
between are filled with high dielectric solvent. For small
molecules, this distinction is unimportant, but for large mac-
romolecules that have considerable interior regions from which
solvent is completely expelled, it results in the effective dielectric
constant of the molecular interior being too high. In other words,
the effective Born radii are underestimated for deeply buried
atoms. The theory apparently needs to be modified to correctly
describe, within the accuracy of the continuum approximation,
the electrostatics of the interior regions of large molecules, while
retaining its good performance for small molecules.
2.2.2. Packing Correction Factor. We begin by analyzing
the correct behavior of the effective Born radius as an atom
becomes buried deep inside the solute, and then modify the
analytic GB model to reproduce it. Consider a large macro-
molecule and assume that its interior is totally inaccessible to
the solvent, the situation corresponding to the molecular surface-
based definition of solute volume. The effective Born radius of
any buried atom i must then be no smaller than the shortest
distance L
i
between the atom and the molecule-surface
interface. For the roughly spherical molecular surface of Figure
1, the integral in eq 16 can be approximated analytically and
we obtain R
i
≈ L
i
. For a nonspherical molecule, the total volume
of the low dielectric region around atom i is always larger than
that corresponding to the inscribed sphere of radius L
i
, and
therefore the generalized Born radius of the atom must satisfy
the condition
∆W
i
)
1
8π
w
∫
solvent
[DB
i
(rb)]
2
d
3
rb+
1
8π
p
∫
solute
[DB
i
(rb)]
2
d
3
rb-
1
8π
p
∫
solvent
[DB
i
0
(rb)]
2
d
3
rb-
1
8π
p
∫
solute
[DB
i
0
(rb)]
2
d
3
rb (10)
DB
i
0
(rb) ≡
q
i
r
3
rb (11)
DB
i
(rb) ≈ DB
i
0
(rb) ≡
q
i
r
3
rb (12)
∆W
i
)
1
8π
(
1
w
-
1
p
)
∫
solvent
[DB
i
0
(rb)]
2
d
3
rb)
(
1
w
-
1
p
)
q
i
2
8π
∫
solvent
1
r
4
d
3
rb)-
1
2
(
1
p
-
1
w
)
q
i
2
R
i
(13)
R
i
-1
)
1
4π
∫
solvent
1
r
4
d
3
rb (14)
1
F
i
)
1
4π
∫
R
3
1
r
4
θ(| rb| -F
i
)d
3
rb
)
1
4π
∫
solute
1
r
4
θ(| rb| -F
i
)d
3
rb+
1
4π
∫
solvent
1
r
4
d
3
rb (15)
R
i
-1
)F
i
-1
-
1
4π
∫
solute
θ(| rb| -F
i
)
1
r
4
d
3
rb (16)
Figure 1. Schematic of a macromolecule, with the circles representing
the atomic VDW spheres. L
i
is the shortest distance from the atom i to
the molecular surface. For a roughly spherical molecule with the atom
i in its center, the generalized Born radius should be R
i
≈ L
i
.
3714 J. Phys. Chem. B, Vol. 104, No. 15, 2000 Onufriev et al.
Physically this condition implies that if L
i
f ∞ the atom does
not “see” the solvent and its Born solvation energy ∆W
i
in eq
13 must tend to zero. Any GB theory must satisfy this general
criterion in order to correctly describe large, densely packed
molecules.
A GB approximation does not satisfy eq 17 if the integral in
eq 16 is taken only over the solute volume based on VDW
spheres, instead of the molecular surface-based volume; we miss
the interatomic spaces inaccessible to solvent in Figure 1,
resulting in an underestimation of R
i
. To correct this, we
introduce a new parameter λ > 1 into eq 16:
where the integral is now taken over the VDW volume of the
spheres (shaded regions in Figure 1), and the λ factor compen-
sates for the missing volume. For an initial guess of λ, consider
the hypothetical case of a very large molecule made up of
identical atoms, the atom i being in its center, as in Figure 1.
Since even in the densest packing (without overlaps) of identical
spheres, only three-fourths of the total space is occupied,
choosing λ ∼ 4/3 should approximately compensate for the
missed volume, and bring eq 18 in better accord with condition
eq 17. The assumption of the dense packing is, in fact, a very
reasonable one for proteins.
23,24
Physically, the value of the
packing correction factor λ > 1 accounts for nonzero size of
the solvent molecule.
As for the effect of this correction on surface atoms and small
molecules, note that dR
i
/dλ ∼ R
i
(R
i
/F
i
- 1) for λ ∼ 1, so that
only large effective radii of deeply buried atoms change
significantly, while the small ones of the surface atoms change
only slightly, since for them R
i
≈ F
i
. Therefore, for our
calculations we can use nearly the same set of atomic parameters
as in earlier works, such as in ref 15, while making the GB
model more suitable for interior regions by optimizing the value
of λ; see below.
2.2.3. Parameterization of the Model. We employ the
approach of Hawkins et al.
14
and approximate the effective Born
radius of atom i by calculating analytically the contribution of
each atom j to the integral in eq 18 and adding the contributions
together. Since the procedure ignores overlaps among the atoms
surrounding atom i, we follow Hawkins et al.
15
and introduce
empirical scaling factors S
j
that partially account for this
behavior. Here we adopt the model that the S
j
’s depend only
on the identity of atom j, with values given elsewhere.
15
Following Still et al.
11
we begin the calculation of effective Born
radii with atomic radii reduced slightly from those used in the
corresponding numerical PB calculations; the offset is F
0
) 0.09
Å. We have already mentioned that setting λ ∼ 1.33 is expected
to have very little effect on individual charge-charge interac-
tions in a small molecule. However, since all of the effective
radii increase with λ, and since we also wish to retain the
remarkable accuracy of the GB in solvation energy calculations,
we shift all of the effective radii calculated via eq 19 downward
by a small term δ ) 0.15 Å in the end of the calculation: R
i
f
R
i
- δ. This term is purely empirical and allows us to avoid
complete reparametrization of the model, such as readjustment
of S
j
or F
0
. Note that setting δ ) 0.15 Å has little effect on the
large radii R
i
. 1 Å of buried atoms.
Introducing the sum over i and j and the above parameters
into eq 18, we write
where H(r
ij
,S
j
(F
j
-F
0
)) represents the result of integration (scaled
by r
-4
) over the VDW volume of atom j; its analytical form is
given in ref 14 and we use it here.
We model the effects of ionic strength by simple substitution
in eq 7, where κ is the Debye-Hu¨ckel screening parameter.
This approximation was introduced earlier
15
and was shown to
work reasonably well. Finally, we notice that since the free
energy W in eq 7 is a quadratic function of charges, we can
express it in terms of the GB analogue of the electrostatic
potential φ(rb
i
) defined at the position of each atom i:
These potentials enter into the titration calculations described
below. Unless otherwise specified, we use
w
) 80 and
p
) 4.
To model physiological conditions in our test calculations, we
set κ
-1
) 10 Å which corresponds to about a 0.1 M solution of
a monovalent salt. The present form of the GB approximation
permits straightforward introduction of a cutoff, where only the
atoms within a specified cutoff distance are taken into account
when computing both the effective Born radius and the atom-
atom interactions for a given atom.
2.2.4. Titration Calculations. The approach used here to
calculate the electrostatic contributions to titration in proteins
has been described in detail elsewhere.
2,29-31
It is assumed that
the difference in the titration behavior of an ionizable group in
a protein and in a model compound can be accounted for by
calculating the difference in the electrostatic work of altering
the charges from the unprotonated to the protonated state in
the protein and the work of making the same alteration in the
model compound. In its original formulation for the PB model,
this quantity is expressed through the values of the electrostatic
potential at the atoms’ positions. Therefore, the same formalism
can be applied verbatim for the GB method once the GB
potential analogue is defined via eq 21.
The protonation fraction of each site at any particular value
of the pH can be obtained by considering a Boltzmann-weighted
sum over all possible protonation states of the protein, or in
the case of a large number of sites, by a suitable approximation
method. In the present work, either the reduced site method
32
or a Monte Carlo method
33
is used. The pH at which the
protonation fraction of a site is 0.5 is then reported as the
calculated pK (or pK
1/2
) of the site.
We have incorporated the GB method into the MEAD suite
of programs,
26,31
where it optionally replaces the finite difference
solver for calculating electrostatic interactions. The codes and
input files necessary to reproduce the calculations presented in
this work will be made available for download through the Web
site, http://www.scripps.edu/bashford, with the next public
release of MEAD.
3. Results and Discussion
3.1. Optimization of λ. We begin by optimizing the value
of the correction factor λ introduced above by comparing the
charge-charge interaction matrix W
ij
calculated by the GB
R
i
g L
i
(17)
R
i
-1
)F
i
-1
-
λ
4π
∫
VDW
θ(| rb| -F
i
)
1
r
4
d
3
rb (18)
R
i
) [(F
i
-F
0
)
-1
- λ
∑
j
H(r
ij
,S
j
(F
j
-F
0
))]
-1
- δ (19)
(
1
p
-
1
w
)
f
(
1
p
-
e
-κf
ij
gb
w
)
(20)
φ(rb
i
) )-
(
1
p
-
e
-κf
ij
gb
w
)
∑
j
q
j
f
ij
gb
(r
ij
)
+
∑
j*i
1
p
q
j
r
ij
(21)
Modification of the Generalized Born Model J. Phys. Chem. B, Vol. 104, No. 15, 2000 3715
approach with the one obtained by solving the PB equation using
the same set of coordinates, charges, radii, and internal dielectric
constant. The interaction is calculated as
where φ
j
(rb
i
) is the electrostatic potential due to atom j at the
position of atom i; for the GB model it is given by eq 21 with
a fixed value of j.
We choose myoglobin as a test case, since it has a large
interior. The methodology used for obtaining the pairwise
interactions by the PB approach is identical to that of ref 15.
As in the previous work, we use the partial atomic charges
from the AMBER force field
27
and the standard Bondi radii.
28
Atomic coordinates are taken from the 1.5 Å X-ray structure.
34
The interior dielectric constant
p
is set to 1.0. The GB
calculations for different values of λ are performed, and the
root-mean-square deviation (rmsd) between the PB and GB
values for the W
ij
are calculated (Figure 2). The two models
are in the best agreement when λ ≈ 1.4, which is close to our
initial guess of λ ) 1.33 derived from simple packing
considerations
3.2. Pairwise Interaction Energies. We have assessed the
accuracy with which the modified GB model reproduces the
PB results for charge-charge interactions in small molecules.
Figure 3 illustrates the performance of the λ ) 1.4 GB model
relative to the PB calculations by plotting the charge-charge
interaction energies in aspartic acid calculated using the GB
theory vs those calculated using the PB model for the same
pairs of atoms. The numbers predicted by the two methods are
nearly identical. Similar results for the original (λ ) 1) GB
model were reported earlier.
15
Solvation energies of small
molecules calculated using the modified GB theory approximate
those obtained by the PB model reasonably well: for aspartic
acid, glutamic acid, and a GLU-GLU dipeptide in a standard
extended conformation the relative differences are 5%, 4%, and
2.5% respectively.
For larger compounds, we expect
15
larger deviations between
the GB and PB approaches, and this is indeed seen in Figure 4
(left panel); but as one would expect from Figure 2, the modified
GB model (Figure 4, right panel) agrees with the PB calculations
more closely than the original GB model. The distribution of
points around the x ) y line (perfect match) is narrower,
especially for larger energies, in the case λ ) 1.4 than for the
original model corresponding to λ ) 1.0.
It is also instructive to compare the values of the so-called
effective dielectric constant
ij
between two atoms i and j,
computed using different models. This quantity is defined so
that W
ij
) q
i
q
j
/(
ij
r
ij
). Smaller values of
ij
correspond to buried
atoms characterized by larger charge-charge interactions. We
compute this quantity using the same set of W
ij
as before, and
plot the resulting
ij
as a function of charge-charge separation
(Figure 5). As expected, the original GB model with λ ) 1.0
significantly overestimates
ij
(underestimates the charge-charge
interactions) for buried atoms relative to the PB calculations,
while the λ ) 1.4 theory shows considerable improvement. Note,
however, that setting λ ) 1.4 does not necessarily bring the
GB model into exact agreement with the general criterion
formulated earlier, eq 17. This value of λ merely represents the
optimum for the current model. Increasing λ even further to λ
) 1.7 (Figure 5, bottom left) brings the charge-charge
interactions of the most deeply buried atoms (low
ij
) closer to
the corresponding PB values, but worsens the overall agreement
between the two models, as shown in Figure 2. This is most
likely due to an overestimation of the effective Born radii of
surface atoms; the GB model originally was parametrized to
perform well only on small molecules.
3.3. pK Calculations. We have assessed the performance of
the modified GB model by calculating titration curves for a
number of individual residues in lysozyme, sperm whale
carbonmonoxymyoglobin, and bacteriorhodopsin. All of these
proteins have titratable groups both on the surface and in the
interior regions, and in the case of bacteriorhodopsin some of
the residues exhibit extremely large pK shifts. Theoretical pK
calculations based on the PB model for these proteins have been
reported previously,
2,30,31
and we follow the general methodol-
ogy presented in those references. Our primary objective is to
assess the performance of the modified GB model, mainly
relative to the established PB approach. The results are reported
in terms of a single number for each residue, the pK
1/2
, which
is the midpoint of its titration curve.
31
In all titration calculations we use the standard Bondi radii
set, and use
p
) 4.0,
w
) 80.0. Unless otherwise specified
we choose λ ) 1.4, δ ) 0.15 in the GB calculations, and a
probe radius of 1.4 Å to compute the dielectric boundary in the
PB method. Only single-conformer calculations are performed
for each structure.
3.3.1. Lysozyme. A first set of calculations was made using
the coordinates of the triclinic form of hen egg lysozyme as
determined by neutron scattering
35
(PDB Accession No. 0LZ5).
Since hydrogen atom positions are available in this case, we
perform no further manipulation on the structure and assign
atomic charges according to the standard AMBER classification
scheme. In contrast to the original lysozyme pK calculations,
which used a single-charge model of the titrating sites,
2
we use
a full set of partial charges for both the protonated and
deprotonated forms of the sites. We calculate pK
1/2
values for
21 titratable residues in lysozyme using the GB model, and
compare the results with experimental data and PB calcula-
tions in Figure 6a. The GB method gives a good overall
agreement with experiment. It correctly predicts the pK
1/2
values for most residues and reproduces the trends in the
experimentally observed pK
1/2
shifts, such as the pronounced
downward shift for ASP-48 and ASP-66 or an upward shift of
TYR-53. It is also important to notice that the GB pK
1/2
values
are highly correlated with the PB ones for all residues, even
Figure 2. rmsd between charge-charge interaction energies W
ij
in
myoglobin calculated within GB and PB models for different values
of λ. The rmsd is computed as ((1/N)∑
ij
(W
(ij)
GB
- W
(ij)
PB
)
2
)
1/2
over a
randomly selected set of N ≈ 30 000 pairs of atoms.
W
ij
) q
i
φ
j
(rb
i
) (22)
3716 J. Phys. Chem. B, Vol. 104, No. 15, 2000 Onufriev et al.
