Comment on the Validation of Continuum

Electrostatics Models
MARCO SCARSI, AMEDEO CAFLISCH
Department of Biochemistry, University of Zürich, Winterthurerstrasse 190,
CH-8057 Zürich, Switzerland

Received 4 January 1999; accepted 1 June 1999

ABSTRACT: A validation based on solvation energies (vacuum to water

transfer) is not sufficient to justify the use of approximated models of
electrostatics to rank ligand/protein complexes. A full validation should be based
on energies in solution, i.e., solvation plus vacuum Coulomb energies, because of
the anticorrelation between solvation and vacuum energies. The energy in
solution is the relevant quantity in simulations of biological macromolecules
and complexes. c 1999 John Wiley & Sons, Inc. J Comput Chem 20:
1533–1536, 1999

Keywords: continuum electrostatics; solvation energy; Poisson equation;

thrombin; docking

F ast methodologies are needed to estimate sol-

vation effects for docking flexible ligands in re-
ceptor binding sites. Approximated continuum elec-
solute–solvent energy that is lost when part of the
surrounding solvent is replaced by molecules with
lower polarizability. Desolvation is the only energy
trostatic approaches are often used, but sometimes a term that disfavors the association between two op-
good agreement with the exact solution of the Pois- posite charges in solution.
son equation (finite-difference method1 ) can origi- In this comment, we focus on the screened inter-
nate uniquely from an inappropriate validation. action between two molecules (e.g., a receptor and
A comprehensive continuum treatment of the a ligand) in a solvent. It can be calculated from the
electrostatic effects of a high-dielectric solvent in potential solution of the Poisson equation, assuming
the calculation of ligand-receptor binding energies that one of the two molecules (in the present case the
should include the evaluation of the screened inter- ligand) is an uncharged low dielectric cavity:
molecular interaction, as well as the receptor and X
N


ligand desolvations. The former is the intermolec- solution =
Eint qi φ Eri (1)
ular energy of the complex in solution, whereas the i=1
desolvation of a solute molecule is defined as the where N is the number of ligand atoms, qi is the
partial charge of the ligand atom i, and φ(Eri ) is
Correspondence to: A. Caflisch; e-mail: caflisch@bioc.unizh.ch
the electrostatic potential generated by the charges
Contract/grant sponsor: Swiss National Science Foundation; of the receptor at the position of atom i. This ap-
contract/grant number: 31-53604.98 proach requires the solution of the Poisson equation

Journal of Computational Chemistry, Vol. 20, No. 14, 1533–1536 (1999)

c 1999 John Wiley & Sons, Inc. CCC 0192-8651 / 99 / 141533-04
SCARSI AND CAFLISCH

(a) (b)
FIGURE 1. (a) Screened interaction calculated according to eq. (1) (x-axis) and to eq. (2) (y-axis) for 1025
ligand–thrombin complexes. (b) Screening energy calculated according to eq. (6) (x-axis) and to eq. (7) (y-axis) for 1025
ligand–thrombin complexes. All values are in kcal/mol.

for every ligand-receptor conformation and, con- where Eint eint

solution and Esolution are the exact and the
sequently, is very time consuming. For efficiency approximate screened interaction, and Eint solvation and
reasons, it has been assumed that the receptor is the e
Eint
solvation are the exact and the approximate screen-
only low-dielectric medium that needs to be consid- ing. The screening is the change in intermolecular
ered and the additional contribution of a ligand to electrostatic energy upon solvation of the complex.
the polarization and the potential is negligible.2, 3 It is a solvation energy lacking the intramolecular
According to this assumption, the electrostatic po- contributions. With a solute dielectric constant of 1,
tential of the isolated receptor (φ̃) is calculated only one has:
once, and the screened interaction energy between X X qi qj X
N


the receptor and a bound ligand is obtained as: vacuum =
Eint = qi φvacuum Eri (5)
rij
i ∈ ligand j ∈ receptor i=1
X
N

X
N




e qi φ̃ Eri
solution =
Eint (2) solvation =
Eint qi φ Eri − φvacuum Eri (6)
i=1 i=1

This assumes that different positions of the lig- X

N




e qi φ̃ Eri − φvacuum Eri
solvation =
Eint (7)
ands will not significantly perturb the potential.2, 3 i=1
Because the solvent effects are used as a correc-
tion to a vacuum energy function, one can rewrite where φvacuum (Eri ) is the electrostatic potential of the
isolated receptor in vacuo. It can be evaluated ana-
eqs. (1) and (2) as:
lytically exactly on every position, Eri , or on the same
grid used for φ(Eri ). Using a grid spacing of 0.5 Å,
solution = Evacuum + Esolvation
Eint int int
(3)
the two approaches practically coincide (correlation
e
Eint = Eint +e
Eint (4) 0.9999, slope 1.0022) and we adopted the first one.
solution vacuum solvation

1534 VOL. 20, NO. 14

 COMMENT ON VALIDATION OF CONTINUUM ELECTROSTATICS MODELS

(a) (b)
FIGURE 2. (a) Screened interaction calculated according to eq. (1) (x-axis) and to the “scrambled” model (y-axis) for
1025 ligand–thrombin complexes. (b) Screening energy calculated according to eq. (6) (x-axis) and to the “scrambled”
model (y-axis) for 1025 ligand–thrombin complexes. All values are in kcal/mol.

One is tempted to validate the solvation model by thrombin active site from P3 to P20 . A solute dielec-
comparing the approximate with the exact screen- tric constant of 1 and a solvent dielectric constant
ing (e Eint int
solvation with Esolvation ). Such validation is not of 78.5 were employed in all calculations together
correct because the energy function used for dock- with a grid spacing of 0.5 Å. The comparisons be-
ing ultimately contains only e Eint
solution , the screened tween the screened interactions and between the
interaction. Solvation energies of molecular com- screenings are shown in Figure 1a and b, respec-
plexes often show a strong anticorrelation with the tively. As is evident from Figure 1a, the approximate
vacuum energies P of the complexes themselves.4 – 6 In method does not show good agreement with the ex-
eqs. (6) and (7), i qi φvacuum (Eri ) can, in general, be act solution of the Poisson equation. Nevertheless,
P P
much higher in module than i qi φ(Eri ) or i qi φ̃(Eri ). a validation based on comparison between screen-
Consequently, the correlation between e Eint
solvation and ings (Fig. 1b) gives the opposite impression (see also
int Fig. 2 of ref. 3).
E
Psolvation can result mainly from the correlation of
i qi φ vacuum (E
ri ) with itself. As a paradox, we generated another set of
The UHBD program7 was used to calculate screened interactions for the same complexes by
the approximate and the exact screened interac- randomly mixing the list of 1025 exact screened
tion (and, consequently, also the corresponding interactions. As expected, the “scrambled” model
screening) for 1025 complexes of the rigid struc- showed no agreement with the exact solution of the
ture of thrombin with small organic compounds Poisson equation (Fig. 2a). Yet, the screening result-
(acetate ion, benzoate ion, methylsulfonate ion, ing from the “scrambled” model still showed a very
methylammonium ion, methylguanidinium ion, good correlation with the exact screening (Fig. 2b)
2,5-diketopiperazine, and benzene). These com- because of the aforementioned correlation of Eintvacuum
pounds were distributed on the surface of the with itself.

JOURNAL OF COMPUTATIONAL CHEMISTRY 1535

SCARSI AND CAFLISCH

This indicates that approximate continuum elec-

trostatic models should be validated on the energies References
in solution rather than on the solvation energies.4
1. Warwicker, J.; Watson, H. C. J Mol Biol 1982, 157, 671.
The latter quantities often tend to show a much bet-
2. Shoichet, B. K.; Stroud, R. M.; Santi, D. V.; Kuntz, I. D.; Perry,
ter agreement with the exact solution of the Poisson
K. M. Science 1993, 259, 1445.
equation than the former ones, but only energies in
3. Vieth, M; Hirst, J. D.; Kolinski, A.; Brooks, C. L. III J Comput
solution are useful for discriminating between fa- Chem 1998, 19, 1612.
vorable and unfavorable binding modes. 4. Scarsi, M.; Apostolakis, J.; Caflisch, A. J Phys Chem A 1997,
101, 8098.
5. Fujinaga, M.; Huang, K.; Bateman, K. S.; James, M. N. G. J Mol
Acknowledgment Biol 1998, 284, 1683.
6. Vorobjev, Y. N.; Almagro, J. C.; Hermans, J. Prot Struct Func-
We thank Prof. J. A. McCammon for providing tion Genet 1998, 32, 399.
the UHBD program, which was used for all finite- 7. Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A.
difference calculations. Comp Phys Commun 1991, 62, 187.

1536 VOL. 20, NO. 14