Protein Engineering vol.8 no.7 pp.663-675, 1995
Free energy simulations of the HyHEL-10/HEL antibody-antigen
complex
R.Pomfes1-2-3, R.C.Willson4 and J.A.McCammon1-5
'Department of Chemistry and 4Department of Chemical Engineering and
Department of Biochemical and Biophysical Sciences, University of
Houston, Houston, TX 77204, USA
2Present address: Department of Chemistry, University of Montreal,
CP 6128, Succ. A, Montreal (Quebec) H3C 3J7, Canada
'Present address: Department of Chemistry and Biochemistry and
Department of Pharmacology, University of California at San Diego,
La Jolla, CA 92093, USA
'To whom correspondence should be addressed
Free energy simulations are reported for the NSl'-D muta-
tion, both in the HyHEL-10-HEL antibody-lysozyme
complex and in the unliganded antibody, using the thermo-
dynamic-cycle perturbation method. The present study
suggests that the mutation would change the free energy
of binding of the complex by -5.6 kcal/mol (unrestrained
free energy simulations), by -0.5 kcal/mol (free energy
simulations with a restrained backbone) and by 1.8 kcal/
mol (Poisson—Boltzmann calculations, which also use a
restrained geometry model). A detailed structural analysis
helps in estimating the contributions from various residues
and regions of the system. Enhanced recognition of HEL
by the mutant HyHEL-10 would arise from the combination
of thermodynamically more favorable conformational
changes of the CDR loops upon association and subsequent
charge pairing with Lys96 in the antigen.
Key words: antibody-antigen complex/free energy simulations/
single-point mutant/thermodynamic cycle
Introduction
The molecular recognition of large substrates is central to the
immunological response. The recognition and binding of
antibodies to antigenic molecular sites helps inactivate viruses
and bacteria and initiates a succession of events that results in
the destruction of foreign molecules. This process requires
both a very highly selective discrimination in order that the
integrity of the self be preserved and a dazzingly diverse
response, so that virtually any non-self molecule can be
identified. The elucidation of the 3-D structures of antibody-
antigen complexes revealed that these properties arise from
the association with a large number of antibody residues on
an extended surface made up of as many as six antigen-binding
loops whose primary sequence displays high variability (the
hypervariable loops) (Alzari et al., 1988).
The association of proteins is thought to result from thermo-
dynamic forces involving a large number of intermolecular
interactions that may be described as hydrogen bonds,
van der Waals contacts, and electrostatic interactions, purely
entropic contributions such as those arising from the restriction
of translational, rotational and internal conformational freedom
upon association, and solvent effects. Whereas the nature and
© Oxford University Press
origins of these various factors are generally understood, their
relative importance is a matter of current debate. For instance,
Dill (1990) argues that hydrophobic forces dominate protein
folding, but that ion pairing can also stabilize proteins, and he
notes that the magnitude of hydrogen bonds and van der
Waals interactions among polar amino acids remains poorly
understood. A detailed understanding, however, is paramount
to the effective development of protein engineering strategies,
in protein folding and stability and molecular recognition
problems alike (Novotny et al., 1989; Sauer and Lim, 1992).
Among the techniques that can be used to probe the stability
of protein-substrate interactions, site-directed mutagenesis
offers wide prospects towards a detailed assessment of the
factors at play in molecular recognition: the replacement of
individual side chains potentially allows for very specific
perturbations to the association. Thermodynamic measure-
ments such as calorimetric measurements can then be used to
determine the corresponding change in stability. In order to
extract unambiguous information from such studies (Sauer and
Lim, 1989), one needs detailed information on the nature of
the physical changes induced by the point mutation. This
is sometimes difficult to obtain in cases where structural
information on the mutant protein is lacking, because significant
conformational rearrangements may result from the mutation,
which in turn can complicate predictions made with reference
to the native state. In order to minimize this uncertainty, it
may be preferable to choose mutations that conserve the
general structure and activity of the native protein. In such
cases, one might hope that small and localized changes in the
physical properties of the mutant give rise to small and
localized changes in the structure and affinity of the protein
for its ligand.
Computer techniques appear well suited to address the
problem of structure and affinity changes upon small modifica-
tions of macromolecules, because they allow one to monitor
both effects on a microscopic scale. In particular, molecular
dynamics simulations used in conjunction with free energy
calculations are the technique of choice if one wants to monitor
the results of a point mutation (Subramaniam et al., 1989).
Such methods have been developed and used successfully in
a number of recent studies pertaining to the binding of ligands
(Lee, 1992 and references therein; Straatsma and McCammon,
1992), including the effect of mutations on the free energy of
dimerization of hemoglobin (Kuczera et al, 1990), on the
stability of T4 lysozyme (Prevost et al., 1991) and barnase
(Tidor and Karplus, 1991), and on the binding of tryptophan
by the trp repressor (Komeiji et al., 1992).
The work presented here attempts to extend the methodology
to an antibody-antigen complex as an example of a protein-
protein complex with a large interface. The engineered modi-
fication of antibodies is currently a very active field and
research groups have reported successful attempts at fabricating
antibodies with designed catalytic properties and human anti-
bodies have been designed with mouse hypervariable loops
(Presta, 1992 and references therein).
663
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

R.Pomfcs, R.CWIIlson and J-A.McCammon
We present the application of free energy calculations to
determine the relative antigen-binding affinity of a mutant
antibody with respect to the wild type HyHEL-10 (Smith-Gill
etai, 1984). The crystal structure of that antibody's antigen
binding fragment (Fab) complexed with its antigen protein,
hen egg white lysozyme (HEL), has been solved (Padlan et al,
1989). The purpose of the mutation considered, the replacement
of an asparagine side chain, Asn31L, by Asp, is to putatively
introduce a salt link between that residue and the positively
charged terminal of a lysine of lysozyme, Lys96W£L, which is
within contact distance in the crystallographic structure and
appears to be shielded from solvent. The latter distinction is
significant because electrostatic screening by water might
weaken the ion pair. In fact, recent mutational studies of salt
links on the surface of proteins converge to the conclusion
that surface salt links contribute marginally to the stability of
proteins (Horovitz et al., 1990; Dao-Pin etai, 1991). On the
other hand, the weakening of a solvent-inaccessible salt link in
HyHEL-5-HEL, another antibody-lysozyme complex, through
replacement of Arg48W£L by Lys snowed a 1000-fold loss in
affinity (Lavoie et al, 1989). Furthermore, the stability of the
association of the constant domains of immunoglobulins has
been found to correlate well with the number of buried salt
bridges between them (Schiffer etai, 1988).
Some care was taken in the solvation of our protein model.
This is of importance since there are a number of possibly
solvated cavities within the protein assembly and because
buried waters can be of critical importance to the structure
and activity of polypeptide molecules, as discussed in a
computational study (Wade et al., 1990, 1991) and a review
article (Meyer, 1992). However, the crystallographic structure
is of little help in the present case, since its resolution allows
for the detection of only one water molecule (Padlan etai,
1989). On the other hand, the association of the antibody Fab
HyHEL-5 with hen egg lysozyme has been believed to result
in large part from excellent complementarity that excludes
interface water molecules (Sheriff etai, 1988), an assumption
that will be tested in the present system.
One of the factors likely to affect the thermodynamics of
antibody-antigen binding is the extent of conformational
rearrangement of the protein(s) upon association. The conven-
tional view of antibody-antigen association regards it as a
'lock and key' process (Amit et al, 1986), whereby the proteins
are presented to each other in the conformation in which they
will bind. Whereas this was inferred from the highly conserved
sequence of most of the antigen-binding domains of antibodies
(the variability of the so-called 'variable' domain being in fact
essentially restricted to the hypervariable loops making up the
complementarity-determining regions, or CDRs, that bind the
antigen), evidence for some degree of adjustment of the
antigen-binding regions upon complexation or 'induced-fit'
has begun to emerge (Rini et al, 1992 and references therein).
One expects limited flexibility of the CDR loops, however,
because a high degree of flexibility could in principle allow
for non-specific recognition of various antigenic surfaces and
would augment the entropic cost of binding, thus limiting both
the specificity and affinity. Limited flexibility is, in turn,
compatible with the concept of canonical structures (Chothia
and Lesk, 1987). This assumption is essential to the current
tractability of computational studies, since in order to allow
for the convergence and accuracy of the free energy perturba-
tions, all accessible conformations must be sampled during
free energy calculations. Furthermore, we expect the specificity
of the antibody to be very sensitive to the precise relative
orientation of the CDR loops. The matter of CDR loop main
chain conformations will therefore be addressed in some detail
in the present study.
Finally, an outstanding problem concerns the treatment of
long-range electrostatic interactions. This aspect of molecular
dynamics calculations can be very demanding if the solvent
molecules are included explicitly in the model. For that reason,
alternatives to that approach have been sought. Among them,
electrostatic continuum methods, in which the solvent is treated
as a dielectric continuum, have been used recently in a number
of calculations (Sharp and Honig, 1990; Honig etai, 1993;
Roush et al, 1994). In particular, the application of continuum
electrostatic theory to the calculation of the complexation
energies of another antilysozyme antibody, HyHEL-5, with
HEL, for a large number of single-point mutants (Slagle et al,
1994), suggests that the method may be useful in the rational
design of antibodies.
Materials and methods
Thermodynamic-cycle perturbation method
The thermodynamic cycle used to calculate the relative
Helmholtz free energy of antigen-binding of a mutant Fab'
with respect to the wild type Fab is depicted below.
Fab +Ag
Fab' +Ag
Fab/Ag
Fab'/Ag
(1)
The theoretical basis used for the calculation of AAA =
A/42 — AAj is the thermodynamic-cycle perturbation method
(Tembe and McCammon, 1984), whereby AA4 is calculated
as A/t4 —AA3. Because the free energy is a state function, the
two paths are equivalent in theory. Computationally, however,
the latter path is preferred because it does not require the
simulation of bottle-necks such as the lengthy diffusional
encounter of the solutes and the desolvation which occurs
upon docking (equations 1 and 2). Instead, the unphysical
transformations represented by equations (3) and (4) correspond
to the numerical substitution of the wild type residue by its
mutant counterpart. On the assumption that both the wild type
and the mutant antibody variable domains have the same
folding, the calculations of &A3 and AA4 can be carried out
with configurational sampling restricted to local changes at
and around the site of the mutation in the unliganded and
liganded cases respectively.
Two methods have been used extensively that allow compu-
tation of the free energy differences between a molecular
system and a perturbed or mutated counterpart in the canonical
ensemble (Beveridge and DiCapua, 1989; van Gunsteren and
Berendsen, 1990; Straatsma and McCammon, 1991). In both
approaches, a perturbation is introduced into a molecular
dynamics simulation. The so-called 'perturbation method'
computes the free energy difference as
M = -/t71n[<«-AWAr>0] (2)
where k is the Boltzmann constant, AH is the perturbation to
the Hamiltonian of the system and the subscript indicates that
664
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

Free energy simulations
averages are computed in the reference (unperturbed) ensemble.
Whenever the perturbation is large, it is desirable to break up
the simulation into smaller successive perturbation calculations
and subsequently sum up over such 'windows', in order to
improve the accuracy of the averaging.
In the second method, often referred to as 'thermodynamic
integration', the perturbed Hamiltonian is written as
H(k) = (1 - X)"H(0) + k"H(\) (3)
where n is an integer and X, a perturbation parameter that
varies from 0 to 1, gradually introduces the perturbation.
Free energy changes are computed as summations over the
ensemble-averaged derivative of the Hamiltonian with respect
to X. Thus, in the case of linear scaling (i.e. n = 1),
AA
-
<H(\) - i
(4)
X - 0
Here, the right-hand side of the equation is a computable
approximation for a finite number of successive integration
steps separated by an increment 8\.
Decomposition of the free energy change
Whereas the two methods described above are equivalent in
the limit of infinite sampling, the thermodynamic integration
method with linear scaling offers the possibility of approxi-
mately separating individual contributions to the free energy
change. That is, if AH = //(I) - H(0) is partitioned into a
sum of components Z(-AW = Z,{//*(1) - ff(0)}, then equation
(4) can be rewritten as
AA = A//1
A//1 = £ AA' (5)
/ \
where AA' represents the contribution to the free energy change
arising from the perturbation of//'. It should be noted, however,
that although A is a state function, the values of the AA' terms
are path dependent (Straatsma and McCammon, 1992; Mark
and van Gunsteren, 1994). Therefore, contributions to AA
should not be considered as physically exact values, but rather
as an indication of the relative magnitude of contributions
arising from different regions of the molecular system or,
alternatively, from different terms used to describe the force
field, in the particular path chosen for the perturbation.
Molecular dynamics simulations
This and the following sections describe successively the force-
field used in the computations, the solvation and equilibration
procedures and the perturbation calculations for Asn to Asp
mutation in the HyHEL-10-HEL complex and in the unli-
ganded antibody. Similar calculations were initially performed
on a reference solvated dipeptide system. The dipeptide model
with an Asn side chain is depicted in Figure 1 in the crystallo-
graphic conformation of residue 31^.
The starting point for the simulations is provided by the
crystallographic structure of the antigen binding fragment
(Fab) of the antibody HyHEL-10 with its antigen, hen egg
white lysozyme (HEL), solved at 3.0 A resolution (Padlan
etai, 1989). The coordinates of the complex are available
from the Brookhaven Protein Data Bank (Bernstein etal.,
1977). The methodology described in detail below was used
Fig. 1. Dipeptide model used for reference free energy calculations, shown
with an Asn side chain.
for simulations of the antibody-antigen complex, but similar
procedures were used for the unliganded antibody and for
reference simulations of dipeptide molecules, except where
otherwise noted.
The force-field used in all the simulations is that of the
CHARMM program (Brooks etal., 1983), version 20. An
extended-atom approach was taken for most carbon atoms
bearing non-polar hydrogen atoms. The all-atom representation
was used only for the side chain of Asn (Asp)31L. In that side
chain, y} is effectively confined by steric repulsion with the
neighboring peptide groups in their main chain conformation,
so that the principle degree of freedom is x2 (see Figure 1).
In the free energy perturbation calculations, the torsional
energy potential for x2 was replaced by a harmonic potential
which was fitted to the gauche torsional energy well corres-
ponding to the crystallographic conformation and restrained
the sampling to that well. In free energy perturbations of the
side chain, a restrained conformation avoids the need for
longer sampling that may arise from the conformational
isomerization of the mutating groups.
The mass of all explicit H atoms was weighted by a factor
of 10 so as to restrict the high-frequency vibration of covalent
bonds involving hydrogens, which allowed the use of a time
step increased to 4 fs (Pomes and McCammon, 1990). The
Verlet algorithm was used to propagate the equations of motion.
All non-bonded energy calculations were truncated, at 8 A
with a switching function in the case of van der Waals
interactions and at 8.5 A with a shifting function in the case
of Coulombic terms.
Because of the extended size of the system, a 'solvated
droplet' model (Shen et al., 1989) was set up. A sphere 26 A
in radius was selected, centered on C, of AsnSl'- of the
antibody. This sphere includes all four loops of the L chain,
all three CDR loops of the H chain of the antibody and most
of the antigen. Residues lying entirely beyond the 26 A
boundary from the center were deleted. As described below,
a spherical droplet of TIP3P water molecules (Jorgensen et al.,
1983) was used to solvate the protein system. During all
subsequent dynamical stages, only those residues lying within
15 A of the droplet center were allowed to move.
No positional restraints were applied to atoms inside the
dynamical system in the simulations of the complex, but a
665
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

R.Pomis, R.C.WIIIson and Jj\.McCammon
weak harmonic restraint was used on the C^ of Asn31 (Asp in
the mutant) in the dipeptide simulations, so as to keep the
small solute at the center of the droplet. The simulations of
the unliganded antibody were run twice, once with no restraints
and once with harmonic restraints on the backbone dihedral
angles of the loops of the L chain, in an attempt to gauge
the influence of backbone conformational flexibility and to
facilitate the comparison of molecular dynamics results with
the outcome of Poisson-Boltzmann calculations, performed in
the assumption of fixed geometry. The latter antibody system
is henceforth referred to as 'conformationally restrained Fab'.
Some care was taken in order to ensure realistic protonation
states of histidine side chains. The orientation and protonation
site of each of the buried His rings were assessed using
molecular mechanics energy calculations so as to identify
the combination leading to the best hydrogen bonding with
neighboring residues.
Solvation procedure
A large sphere of TIP3P waters (Jorgensen etai, 1983)
equilibrated at 300 K was superimposed on the protein model.
Solvent molecules whose oxygen centers lay within 2.4 A of
any protein heavy atom were deleted. This procedure left a
number of waters either buried in the protein interior or trapped
at the interface between the H, L, and HEL chains, molecules
whose fate was decided according to the following criteria.
Molecular mechanics calculations were performed for a water
probe placed successively on the points of a grid overlapping
the unsolvated protein complex, using the program GRID
(Goodford, 1985; Boobbyere/a/., 1989). Isoenergetic contours
were then calculated and displayed at hydrogen bonding energy
levels (typically, below -8.0 kcal/mol) in order to identify
potentially solvated cavities. T1P3P molecules were then placed
into these contours, oriented so as to minimize the hydrogen
bonding energy, and subjected to short energy minimizations
while the protein heavy atoms were held fixed. CHARMM
molecular mechanics energies were then calculated individually
for these buried solvent molecules and their values compared
to the previously obtained GRID values. Only those water
molecules participating in strong hydrogen bonds with sur-
rounding residues were retained. This procedure was made
necessary by the discrepancies between GRID and CHARMM
energies and its difficulty was compounded by the concurrent
need to reassess the protonation states of nearby histidine side
chains and by the possibility of water-water interactions within
cavities or channels. The ultimate test of assessing cavity
solvation was provided during thermalization of the solvent,
when a buried water molecule diffused through the protein
interior and was subsequently deleted.
Equilibration
Before proceeding with molecular dynamics, a number of
steps were required to ensure the removal of bad contacts
between atoms and the orientational relaxation of water molec-
ules at the interface with the protein(s). This is achieved
through successive energy minimizations and dynamical ther-
malization of the solvent and the protein residues of the inner
sphere. The algorithm chosen for these minimizations uses the
steepest descent procedure so as to leave the system in a local
minimum close to the crystallographic conformation. For the
same reason, the solvent is thermalized first, before the heavy
atoms of the protein are allowed to move.
First, the hydrogen atoms of water molecules inside the 19
A boundary were subjected to 100 steps of minimization,
during which the protein atoms were fixed and the positions
of water oxygen atoms were restrained via a steep harmonic
potential. The energy of water molecules inside the 15 A
boundary was then refined with 200 steps of minimization,
after which the solvent within that inner dynamical sphere was
thermalized in the following fashion: three stages of 2 ps
simulations with temperature gradients followed by 3 ps
equilibrations, at 100, 200 and 300 K successively. The solvent
was then taken through 5 ps of equilibration with velocity
reassignments at 300 K every 0.2 ps. At that stage, an ultimate
comparison of surface and buried water locations with the low-
energy contours constructed using GRID proved satisfactory,
which provided confidence in the solvation of the crystallo-
graphic structure.
The positions of solvent molecules were then fixed and the
protein residues inside the dynamical sphere were taken
through 200 steps of steepest-descent energy minimization.
The procedure followed to thermalize the protein ran over 18
ps and mirrored that used for the solvent.
All residues inside the dynamical sphere were then allowed
to move and the system was subjected to a 10 ps equilibration
with velocity reassignments every 0.2 ps, followed by 10 ps
with velocity reassignments every picosecond. Equilibration
was then completed with a 20 ps run with temperature
relaxation. Small r.m.s. deviations of the positions of inner
atoms confirmed the stability of the structure and the conver-
gence of the equilibration.
Free energy calculations
Free energy differences were calculated using the perturbation
algorithm of CHARMM. The perturbations were made on a
hybrid side chain with a unique topology up to and including
Cy and both Asn and Asp terminal atoms present (dual
topology). Consistent with the topology and parameters of the
CHARMM force-field, the types and parameters of all atoms
besides the terminal NH2 and O atoms of the Asn and Asp
side chains were identical during the simulation, except for
the partial charges attributed to Cy and to Cp and its H
substituents, which varied with X during the perturbation.
An initial perturbation run using the double sampling tech-
nique (Beveridge and Dicapua, 1989) with five reference points
at respectively X. = 0.1, 0.3, . . ., 0.9 was performed on the
dipeptide system and yielded a free energy profile characteristic
of a charge creation in aqueous solution. In order to limit the
free energy difference of each integration step to the order of
thermal fluctuations (kT) and thus reduce sampling problems
(McCammon and Harvey, 1987), all subsequent mutations
were broken down into 77 integration steps with X spacings
of 0.01 between 0.01 and 0.05 and between 0.49 and 0.99 and
spacing of 0.02 between 0.05 and 0.49. At each integration
step, a minimum of 2 ps of equilibration with velocity
reassignments every 0.2 ps and a minimum of 2 ps of data
collection were run. During the production stage, coordinate
sets and the energy values of the reference and perturbed states
were stored every 0.02 ps for post-processing. In the case of
the complex, an additional 2 ps of production was run at each
integration step so as to help assess convergence. The total
length of each perturbation thus ranged between 336 and 452 ps.
Data were then analyzed using the CHARMM algorithm
with both the integration and the perturbation methods. A
minimum of 42 500 configurations was used to compute each
overall free energy change, but 20 times fewer configurations
were used for the decomposition.
666
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

Poisson-Boltzmann calculations
The electrostatic energy contributions to the difference in the
solvation energy of the Asn and Asp dipeptides and to the
relative free energy of binding of wild type and mutant
HyHEL-10 with HEL were estimated by means of continuum
electrostatic calculations using the finite difference Poisson-
Boltzmann (PB) algorithm of the University of Houston
Brownian Dynamics (UHBD) program (Davis etal, 1991).
The theory and methodology of such calculations is described
elsewhere (Gilson etal., 1987).
The conformations used for the PB calculations were
obtained as averages from 20 ps molecular dynamics trajector-
A»n93»BL
Trp97B
TyrSO"-
Flg. 2. Cross-section of a hydrated channel at the core of the Fab-Ag
complex obtained from a 20 ps averaged trajectory showing two water
molecules labeled 1 and 2. The hydrogen bonds involving water 1 are
shown as dashed lines. They bridge all three polypeptide chains H, L and
HEL together.
Free energy simulations
ies of the Asn and Asp dipeptides and of the wild type and
mutant Fab-Ag complexes, after energy minimization of the
positions of hydrogen atoms. The solvation energy of the
dipeptides was calculated as the difference between electro-
static energies of the solute in a continuum with dielectric
constants of 1 and 78, successively. The electrostatic binding
energy of the complex was calculated as the difference between
the electrostatic energy of the Fab-Ag complex and its analogs
for the isolated Fab and Ag moieties in a solvent continuum
with dielectric constant 78. In these calculations, the conforma-
tion of the Fab and Ag moieties was assumed to be retained
upon decomplexation. The entire antigen-binding fragment of
the antibody, Fab, was included. The charges used were
identical to the CHARMM set used in the MD simulations.
The N- and C-termini were assumed to be charged. No buried
water molecule was included. All protonation states and partial
charges were consistent with those used in the MD simulations.
The ionic strength was set to zero; the solute dielectric was
taken to be 4. A grid size of 653 with a spacing between grid
points of 0.3 A was found to be adequate in the finite-
difference calculations of the dipeptide system. In the finite-
difference calculations of the complex system, a grid of
90X90X110 with a spacing of 1.0 A was used as a starting
point so as to include the entire Fab—Ag complex. Focusing
was then used; the focusing grid was made large enough to
include the dynamical MD system and contained 1O53 points
with a spacing of 0.5 A. Calculations were performed on the
same focusing grid so that the self-energy terms would cancel
for each binding energy moiety of the cycle (see equation 1).
Tests of convergence of the PB calculations included varying
the size of the grid by ± 1 unit and varying the grid spacing
by ±10% of its value.
Results
Solvation
The procedure used for solvation led to the hydration of six
cavities and, importantly, of a channel extending from the
A s n 3 1 L
Fig. 3. Cross-section of the Asn31 t-Lys96 t f£L region of the Fab-Ag interface obtained from a 20 ps averaged trajectory showing the proximity of a (buried)
water molecule. The hydrogen bonds involving the positively charged lysine terminal are shown in dashed lines, while the double arrow indicates the salt link
to be created with Asn to Asp mutation of residue 31L.
667
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

R.Pomes, R.CWlllson and J.A.McCammon
edge of the antibody-antigen interface into the core of the
assembly between the L and H chains on the one hand and
HEL on the other. This channel was found to accommodate
three water molecules, two of which are depicted in Figure 2
with a few surrounding residues.
One of the water molecules solvating this channel forms
hydrogen bonds between all three of the polypeptide macro-
molecules. Such hydration sites appear of structural importance
to the quaternary assembly since (i) they help conserve the
crystallographic conformation of channel residues around them
and (ii) their hydrogen bond network is well conserved over
molecular dynamics trajectories extending into the hundreds
of picoseconds.
Another hydrated pocket, in contact with the terminal of
Lys96W£L, was found to accommodate one water molecule that
makes a very stable hydrogen bond with the charged NH3
donor, as depicted in Figure 3. The presence of this water
molecule is important because it could compete with the
formation of the proposed salt link between Lys96//£i and
AspSl'-. In addition to the procedure described previously,
short simulations run with and without that water molecule
showed that it helps preserve the local crystallographic con-
formation around the Lys96W£Z- terminal (see Figure 3).
Antibody conformational flexibility
The backbone torsional angles of the loops of the L chain of
the antibody were calculated as averages from 20 ps molecular
dynamics trajectories computed at all four points of the
thermodynamic cycle (see equation 1) and are listed in Table
I. In this table, the names of residues making up the loop turns
proper appear in bold face and significantly large values or
changes are stressed with bold face type. The organization of
Table I. Main chain
coordinates and the
and Fab
Loop
CDR,
CDRj
Non-CDR
CDR3
dihedral angles ($,
dynamical averages
Residue
Cys23
Arg24
Ala25
Ser26
Gln27
Ser28
Ile29
Gly30
Asn31
Asn32
Uu33
His34
Trp35
Tyr36
Leu46
Leu47
Ilc48
Lys49
iyi-50
AlaJl
Ser52
Gln53
Ser54
Ile55
Ser56
Phe62
Ser63
Gly64
Ser65
Gly66
Ser<>7
Gly68
Thr«9
Asp70
Phe7l
Thr72
Leu73
Ser74
Phe87
Cys88
Gln89
Gln90
Ser91
Asn92
Ser93
V) of residues making up
for the wild type (wt) and
X-ray
wt
-155.
-104.
-85,
-88.
172.
-64.
9,
-63.
85.
-105.
-129,
-126,
-104.
-125,
-81,
-104,
-121.
-92.
30,
-158,
-134,
-86,
-73.
-70,
-97.
-80,
-149,
-137,
-161.
108,
163.
113,
-105,
-75,
-129,
-140,
-127,
-113,
-115,
-124,
-152,
-92,
-145,
-72,
-126,
107
109
135
64
160
-38
63
107
-15
72
119
145
138
172
138
-66
133
179
-72
-A\
7
102
134
167
136
135
141
155
175
177
141
-122
-48
112
134
151
123
135
166
161
126
166
16
-48
153
the loops of the iantibody L chain calculated respectively from the crystallographic
mutant (mu) proteins, successively for
Fab-Ag
wt
-159,
-91,
-81,
-88.
-157,
-99,
-129,
-63,
43,
-129,
-111,
-123,
-121.
-128,
-76,
-113.
-119,
-119,
73,
-177,
-142,
-77,
-71.
-71,
-110,
-80,
-152,
-136,
-144,
156,
-153,
62.
-118,
-111,
-130,
-140,
-134,
-134,
-126,
-114,
-144,
-107.
-110,
-85,
-88,
97
135
177
64
175
93
62
146
49
6
148
155
143
152
145
-73
145
148
-52
-4\
1
104
115
-177
136
135
150
165
169
137
150
-100
-18
132
149
157
149
139
156
156
156
145
10
-74
148
mu
-160,
-90.
-80.
-88,
-149,
-80,
-132,
-146,
36,
-158,
-123,
-118,
-116,
-127
-77,
-116,
-110,
-116,
70.
-165,
-135,
-73.
-79,
-97,
-53,
-92,
-156,
-156,
-151,
168,
-156,
80,
-144,
-103,
-121,
-140,
-125,
-129,
-120,
-105,
-143,
-120,
-111,
-88,
-95,
the complexed
%
137
176
64
158
99
160
156
59
38
145
148
140
151
148
-81
142
150
-61
-43
-14
122
151
-176
135
144
162
170
164
124
164
-85
-39
136
150
148
146
140
141
155
150
143
3
-67
153
and uncomplexed systems Fab-Ag
Fab
wt
-161,
-99,
-80.
-99,
-150,
-60,
34,
-37,
47,
-132,
-120.
-122,
-127,
-129,
-71,
-112,
-119,
-116,
47,
60,
-139,
-84.
-94,
-77,
-97,
-80,
-152,
-104,
-147,
99,
-144,
-87.
-79,
-104,
-127
-140,
-135,
-125,
-118,
-112,
-137,
-107.
-112,
-93,
-113,
107
127
164
15
167
-46
47
106
48
51
146
160
144
153
140
-73
139
155
36
-11
-22
147
148
173
136
135
151
155
180
147
-33
-144
-9
137
159
160
142
141
158
155
145
136
20
-58
147
mu
-133.
-142.
-96.
-59.
-104.
-75.
47,
60,
54,
-109,
-127,
-130,
-126,
-128.
-72,
-116,
-122.
-110,
47,
60.
-141,
-76.
-84,
-75,
-97,
-83,
-150,
-125,
-137,
146.
-142,
-149,
-73,
-101,
-132.
-135.
-116,
-128,
-100.
-82,
-142,
-105,
-130,
-97,
-121,
155
134
177
-68
170
-7
29
23
41
73
143
164
145
154
143
-70
137
155
43
-10
-35
139
145
173
136
134
161
149
164
169
-12
-136
-17
142
151
141
134
133
131
150
144
146
29
-46
150
The names of residues making up loop turns appear in bold face type, as do the values of dihedral angles that deviate largely from the crystallographic
reference or from each other. All values are in degrees.
668
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

R-Pomfes, R.C.Wlllson and J.A.McCammon
displays further changes in the conformation of its CDR, (see
residues 23-24, 26-27 and 30). More detailed observations
are needed in order to assess whether these further changes
are directly induced by the mutation or whether they simply
result from additional sampling, some 400 ps down the
dynamics trajectory.
The Ca deviations from the crystallographic conformation,
and the r.m.s. fluctuations of individual Cas are listed in Table
II. These data illustrate the nature of the limited structural
changes taking place upon mutation and/or upon desolvation
of the antigen binding surface and are consistent with the
peptide reorientations evidenced in Table I. In fact, all the Ca
displacements in excess of 1.0 A and most of the large r.m.s.
fluctuations of the Ca atoms (upwards of 0.40 A) occur at or
next to residues whose backbone underwent conformational
reorientation with respect to the crystallographic fold. In
particular, the CDR! residues of mutant Fab-Ag and Fab
systems affected by backbone transitions display large displace-
ments and record displacements of up to 4 A are reached in
the Fab cases by non-CDR loop residues that have unfolded
into the solvent. Larger r.m.s. fluctuations generally occur in
both of the wild type Fab-Ag and Fab systems compared to
their mutant counterparts, but this is more obvious in the
unliganded Fab case, suggesting that in spite of the screening
of Coulombic interactions by the high dielectric solvent, the
aspartate charge introduced upon the Asn to Asp mutation
might assist in restricting the mobility of nearby residues'
backbone dipoles, particularly Gly30, Asn31, Ala51, Ser52,
Ser67 and Gly68. To some extent, however, the enhanced
mobility of the wild type Fab may also result from incomplete
equilibration. In the following analysis, we try to address these
questions in more detail.
Binding free energy
The overall free energy changes for the Asn to Asp mutations
in the dipeptide, complex and unliganded Fab systems are
given in Table III. The relative solvation free energy of the
dipeptides, before electrostatic cut-off correction, is
— 54.9 kcal/mol (second column in the table). The relative
binding free energy of the mutant antibody is -5.6 kcal/mol
(columns 3 and 4). The calculation performed with harmonic
restraints on the unliganded antibody's main chain conforma-
tion (columns 3-5), is also included in Table III for later
comparison with continuum electrostatic calculations. There
is no straightforward way to estimate the statistical error
attached to these values, but a number of alternative calculations
using fewer perturbation steps, a different number of configura-
tions, extended sampling collection and different algorithms
for the calculation may help in assessing the convergence of
the results. The values obtained for the dipeptide mutation
using the perturbation method with five steps and the integration
method with 77 steps are respectively -54.1 and -54.9 kcal/
mol, within 2% of each other. One expects the calculation
with fewer integration steps to be subject to somewhat larger
statistical errors because the changes in free energy between
successive steps may be large compared to thermal fluctuations,
thus making the reference ensemble unrepresentative of the
perturbed one, which in turn leads to larger statistical uncer-
tainty. On the other hand, the values obtained using the
perturbation method should be the same as those obtained
from the same trajectory with the thermodynamic integration
method in the limit of infinite sampling. In our simulations,
they are separated by 0.2 and 0.1% respectively in the Fab-
Ag complex and the free, unliganded Fab cases. These errors
are small and concentrated at the end-points. Calculating the
free energy difference for the mutation in the Fab-Ag complex
using 20-fold fewer configurations and 50% more (from
additional sampling) yielded free energy differences of -52.8
and —52.2 kcal/mol respectively. These values lay respectively
within 0.6 and 0.1 kcal/mol from the value obtained with the
standard procedure described in the Materials and methods
section. Thus, the convergence of the free energy value seems
within reasonable reach of our simulations. The values retained
for the calculation of the relative free energy of binding of the
mutant antibody were chosen as the average of the values
obtained respectively with the perturbation and the thermodyn-
amic integration (equations 2 and 4). This yields a relative
binding free energy change of — 5.6 kcal/mol for the flexible
chain antibody and -0.5 kcal/mol in the case where the
unliganded antibody's main chain conformation is restricted.
Other sources of statistical error due to conformational
sampling include possibly insufficiently sampled conforma-
tions of the flexible antigen-binding antibody loops and the
sampling of the mutating side chains' conformations. The
former constitutes a potentially formidable problem, because
the time scale corresponding to loop conformational rearrange-
ments may extend into the thousands of picoseconds, a time
range not easily accessible with current molecular dynamics
of biomacromolecules. Among the other degrees of freedom
of the antibody molecule that may be restricted upon antigen
binding are those defining the torsional isomers of the mutating
side chain itself. Unless the native and mutant side chains
have identical rotamers and identical equilibrium distributions
of rotameric states (which is unlikely), a contribution to the
free energy will arise from substituting one side chain for the
other (Straatsma and McCammon, 1989). This contribution
was estimated elsewhere using potential of mean force calcula-
tions on solvated Asn and Asp dipeptide systems (Pomes,
1993) and was found to be small (-0.2 kcal/mol) compared
with the overall free energy change for the Asn to Asp mutation
reported in the present study. This justifies the neglect of
secondary rotamers of the mutating side chain imposed by the
Table III. Helmholtz free energy changes for Asitfl1 to Asp
antibody, the conformationally restricted HyHEL-10 antibody
respectively
Total
Van der Waals
Electrostatic
Covalent
Dipeptide
-54.9
-2.2
^»7.7
-4.5
Fab-Ag
-52.2
-3.4
-40.7
-8.6
mutation in a
(see text) and
Fab
-46.5
-2.3
-36.8
-7 1
dipeptide model, the HyHEL-10-HEL
the thermodynamic-cycle perturbation
Fab*
-51.7
-2.1
-41.8
-7.0
Columns :
-5.6
-1.1
-3 8
-1 5
complex, the unliganded HyHEL-10
values for columns 3 and 4 and for 3-5
i and 4 cycle Columns 3 and 5 cycle"
- 0 5
-1.3
1 2
-1.5
All values are in kcal/mol.
"Conformationally restrained antibody backbone.
670
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

Free energy simulations
harmonic potential on %2 in the free energy calculations
reported here.
The systematic error due to insufficient equilibration of the
system can be estimated by calculating free energy differences
using two different sets of energy values: (i) those collected
immediately after equilibration (0-2 ps) and (ii) those obtained
from additional sampling (2-4 ps) at every integration step.
These values were calculated for the Fab—Ag perturbation
calculation as -52.22 and -52.29 kcal/mol, respectively. The
relative error between these two amounts to less than 0.2%,
well within the statistical error. We therefore consider the extent
of equilibration sufficient in all the calculations reported here.
Other systematic errors include those arising from the force-
field and the truncation of non-bonded interactions. The
comparison of results obtained from our reference calculations
on dipeptide molecules with experimental data provides an
approximate measure of the model's performance in those
respects. The experimental value of the difference in the free
energy of hydration of acetate with respect to acetamide is
directly deduced from reported data (Pearson, 1986) as —73
± 2 kcal/mol. This compares fairly well to our calculated
value of - 7 5 kcal/mol for the relative hydration free energy
of an aspartate versus an asparagine dipeptide. The latter value
was obtained from the 77-step perturbation run, after correcting
for the truncation of long-range electrostatic interactions with
the Asp charge. This correction was approximated by subtrac-
tion of the Bom free energy for solvating a singly-charged 8.5
A radius spherical cavity into a continuum with a dielectric
constant of 78.5 (Bom, 1920; Straatsma and Berendsen, 1988).
The Bom correction amounts to —20.5 kcal/mol, which brings
our calculated Asn to Asp dipeptide relative solvation free
energy value from -54.9 to —75.4 kcal/mol. Consistently, the
relative solvation free energy of the Asp dipeptide with respect
to Asn obtained from continuum calculations with the finite-
difference Poisson-Boltzmann method was found to be
- 7 6 kcal/mol. These results seem to indicate that the force-
field is reasonably adequate for the calculation and that the
Bom correction provides a crude but fair approximation of the
truncated solvent-solute interactions.
In cases where protein residues extend beyond the cut-off
distance from the new charge, however, the uncertainty due
to the cut-off distance of long-range interactions constitutes
an outstanding problem because the Bom approximation cannot
be applied. This is the case of all four stages of the thermodyn-
amic cycle for Fab and Fab-Ag, where both antigen and
antibody proteins extend well beyond the cut-off radius. In
our protein free energy simulations, we note that the only full
charges within the cut-off from the center of charge of the
mutating residue are those of Lys96 in the antigen and
Asp31L itself. The free energy perturbations, therefore, involve
consistent changes in the net charge within the cut-off, respect-
ively from 0 to — 1 (unliganded Fab) and from +1 to 0 (Fab-
Ag complex). Although in principle the Poisson-Boltzmann
calculations might be preferred over the free energy simulations
because they take into account the long-range component of
the electrostatic energy, the large error that would result from
the imbalance of charges within the cut-off was avoided in
the present simulations.
In the future, the long-range component of the electrostatic
contribution could be estimated separately using continuum
techniques more sophisticated than the Bom model for solva-
tion, such as those utilizing the Poisson-Boltzmann equation
to calculate the electrostatic field on a grid superimposed on
the solute(s) via finite-difference approximations (Warwicker
and Watson, 1982; Davis and McCammon, 1989). This is a
challenging problem, especially because of the rapidly chan-
ging state of the PB art (Antosiewicz et al., 1994).
Decomposition of the free energy change
Overall physical contributions. The free energy results were
broken down approximately into different overall physical
contributions: internal (bonded), non-bonded electrostatic and
non-bonded van der Waals. These results are shown in Table III.
The electrostatic term dominates and contributes ~80% to
the free energy change for the mutation in all three systems.
Whereas both protein mutations yield an electrostatic term less
negative than their 'fully solvated' dipeptide counterpart, the
value for Fab—Ag is smaller than that of Fab alone, which
suggests that the penalty for burying a charge is more than
made up for by the salt link introduced by pairing the charges
of Asp31L and Lys96H£t. However, this does not hold in the
case of the conformationally restrained unliganded antibody
(column 4 in Table III), for which the electrostatic term is
more comparable to that of the complex. Accordingly, the
relative electrostatic binding free energy of the mutant Fab
calculated with the Poisson-Boltzmann method was found to
be 1.8 kcal/mol, in fair agreement with the molecular dynamics
free energy calculation result of 1.2 kcal/mol for the restrained-
Fab cycle. In the Poisson-Boltzmann calculations, the Fab and
Ag moieties retain their crystallographic conformations and
the agreement suggests that the electrostatic cut-off correction
to the molecular dynamics simulation results is likely to be
small for this system. The PB calculations also confirm that
with the assumption of a rigid Fab backbone, the single-point
mutation is unfavorable in electrostatic terms. However, the
favorable value of —3.8 kcal/mol obtained for the cycle in
which no restriction was imposed on the backbone of the
unliganded Fab indicates that assuming conservation of the
Fab loops conformation upon solvation or desolvation, as the
lock and key mechanism would imply, may not be valid in
this case.
The van der Waals contribution is small in all cases, which
reflects the isosteric nature of the mutation. Thus, the analysis
10
0
- 1 0
<* - 3 0
5. -30
- 4 0
- 6 0 -
- 6 0
-
-
r
1
t i 1 1
V
4
0.0 0.1 0.2 0.3 0.4 o.e o.a 0.7 o.a 0.9 1.0
Fig. 4. Helmholtz free energy change for Asn (K « 0) to Asp (X = 1)
mutation of residue 31^ in the Fab-Ag complex, as a function of the
perturbation parameter X, is indicated by a solid line. Its components are
also shown as follows: dashed for electrostatic, dash-dotted for covalent and
dotted for Lennard-Jones.
671
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

Free energy simulations
the AAA value from the total non-bonded energy thus appear
as the result of a balance between solvent and protein terms.
In the complex, the overall antibody term is small, but is
itself the balance of large residue contributions that nearly
cancel, mainly a large positive term from the mutating residue
itself (see Asn31 and CDR, values in Table IV) compensated
by favorable mid-range electrostatic contributions with nearby
CDR2. Interestingly, the small overall effect of the mutation
on the antibody's energy in Fab-Ag correlates with a well-
conserved secondary structure and conformation. The antigen
contribution is largely favorable to the mutation free energy,
arising essentially from the creation of a salt link between
AspSl'- and Lys96H£z-, which was the proposed goal of the
point mutation. The strand residues, however, contribute to
halving this favorable term. The solvent contribution compares
in size with that of Ag, although it amounts merely to one-
quarter of its counterpart in the Fab system. This simply
reflects the limited access of water molecules to the rim of the
Fab—Ag interface.
In the Fab system, by contrast, the solvent is responsible
for a very largely favorable contribution to AAA, whereas the
antibody segment reduces this favorable term by approximately
one half. This suggests that whatever protein conformational
rearrangements are taking place upon solvation of the antibody
surface are driven by the solvent and occur at the expense of
protein-protein interactions. Thus, in the Fab system, all three
L chain loops near the mutation site appear to oppose the path
of the mutation.
We now discuss the contributions of a number of individual
amino acid residues and solvent molecules near the site of
mutation in relation to the relative energy changes that they
induce during the mutation.
The non-bonded contribution from the mutating residue
itself, Asn31(L)Asp, is large, positive and nearly identical for
both Fab-Ag and Fab systems. Because this side chain was
conformationally restrained by construction, because the muta-
tion is practically isosteric, and because the backbone itself
retained a consistent conformation at y3 , all along the thermo-
dynamic cycle (see Table I), this large contribution results
mainly from the electrostatic repulsion between the new charge
and the backbone carbonyl.
Among the other amino acids from CDR[ of the L chain,
Ser28 contributes to the non-bonded energy only in the Fab
system. This is because during the mutation, the main chain
carbonyl oxygen gets closer to the new Asp charge center by
3 A.
Similarly, in the unliganded Fab, the y backbone transition
of GlySO*" results in closer proximity of the carbonyl group of
GlySO*- to the mutant group.
Adjacent Asn32Z- provides the only other significant term
from CDR| in the complex. The side chain reorients during
the course of the mutation, from (g~,t) to (g~,g~). There is no
hydrogen bond formed with the terminal of AsnSl^ at any
time (it would compete with Lys96W£i"), but the terminal NH2
orients in an electrostatically favorable way and ends up within
van der Waals contact to the mutant side chain. In the
unliganded case, the interaction is purely electrostatic and not
as strong because the side chain is kept away by strong
hydrogen bonds with the solvent.
CDR2 of the L chain provides many strong interactions in
the complex, but these tend to decrease in magnitude in
the unliganded case because the CDR( conformation shifts
Asn31(Asp) into the solvent somewhat, away from the protein
core, decreasing its proximity to CDR2. This is the case with
Lys49, whose carbonyl group points toward the new Asp
charge some 5-6 A away in the complex but is weakened in Fab.
The ring of Tyr50 lies within contact distance of the terminal
of Lys96H£/\ The backbone carbonyl points away from the
mutant residue, which makes for a favorable charge-dipole
interaction. Inversely, the backbone NH group of Ala51L,
whose methyl is in contact with the mutating side chain, points
toward the carboxylate of Asp31L, thus providing the other
half of that charge-peptide bond stabilization. On the contrary,
the peptide bond linking residues 50 and 51 reorients in the
wild type solvated Fab so that the carbonyl oxygen of TytfO*-
ends up only 4 A from the new charge.
In the complex, the backbone NH group of Ser52i- is within
4 A from Og2 of Asp3lL. It forms a stable and lasting hydrogen
bond with a water molecule which also hydrogen bonds
the mutating terminal. This water molecule reorients from
hydrogen acceptor to hydrogen donor during the mutation
without disrupting its hydrogen bond with S e ^ ^ .
The NE atom of Gln53L lies 5 A away from the new
carboxylate group and donates its hydrogens to two water
molecules: the buried water that bridges with Lys96W£L and a
surface water which bridges with Asn31(L)Asp.
The backbone carbonyl group of Gly66L provides the only
significant contribution from the non-CDR turn in the Fab-
Ag complex. The unfavorable orientation of that group with
respect to the charge introduced by the mutation is made
possible by the presence of solvent.
In the unliganded Fab case, loop turn reorientations tend to
bring CDR| and the non-CDR loop closer together. Several
large unfavorable energy contributions arise as the non-CDR
loop orients its backbone dipoles perpendicularly to the plane
of the turn, with peptide carbonyl moieties pointing at the new
charge. That is the case for residues 66, 68 and 70.
Argl4W£t is responsible for one of the two large unfavorable
terms arising from the 'strand' part of the antigen's epitope.
The backbone of Argl4H£/-, like that of neighboring residues,
does not reorient during the mutation. Its carbonyl points
straight at the side chain of Asn31L, some 3.5-4.5 A away.
As was the case with Gln53L, this carbonyl forms strong
hydrogen bonds with two water molecules that also bind
successively to the carboxylate of Asp31z", thus compensating
for this large positive energy term.
The carbonyl oxygen of His 15HEL forms a strong and stable
hydrogen bond with the ammonium group of Lys96W£L < 3
A away. This also results in its close proximity to the charge
created upon mutation.
The side chain carbonyl of Asn93H£i forms stable hydrogen
bonds with both the buried water molecule that also binds
Lys96W£Z- and Gln53L and with the terminal of Gln53L, which
locks it into a mildly unfavorable orientation with respect to
the mutant side chain.
Finally, the ammonium terminal of Lys96H£/' has its charge
paired with Asp31L (Ne ends up 2.6 A away from Og,), but it
also maintains its strong hydrogen bonds throughout the
mutation respectively with the backbone oxygens of AsnSl*"
and Hisl5W£L and with the neighboring buried water (see
Figure 3).
This buried water itself contributes little to the mutation
free energy, presumably owing to its sideways orientation with
respect to the carboxylate. This orientation is well maintained
during dynamics because of strong hydrogen bonds with
Lys96W£i., Q\n52L and Asn93H£L, as discussed above.
673
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

R.Pomis, R.C.WIllson and J.A.McCammon
For illustration, the energetic contributions of two other
water molecules that engage in hydrogen bonds with the
carboxylate during the course of dynamics in the Fab-Ag
system were also computed. The first one also binds Se^*-
and contributes -8.9 kcal/mol to the relative binding free
energy of the mutant, while the second one bridges the mutant
residue with Gln53L and contributes -10.2 kcal/mol. The
magnitudes of these contributions are comparable to the largest
protein residue terms, with the exception of the salt link itself.
Summary. In the antibody-antigen complex, the overall free
energy gain arising from the point mutation divides evenly
between solvent and antigen contributions, with a marginal
contribution from the antibody itself. The latter arises primarily
from the near cancellation of a positive mutant self-term with
negative contributions from the CDR2 loop. As to the antigen
term, it is dominated by a large salt link term, but tempered
by the adverse orientations of backbone carbonyl groups close
to the site of mutation.
Out of the 12 protein residues selected as giving significant
non-bonded contributions to the free energy change for mutat-
ing Asn31L into Asp in the Fab-Ag complex, one involves a
charge-charge term and the other 11 are dominated by charge-
dipole interactions: six unfavorable terms with carbonyl groups
(five of them on the main chain), one favorable term with a
peptide bond and three with uncharged amide or peptide
nitrogen groups (two of them in side chains).
As many as five of these residues are involved in hydrogen
bonds with waters that also bind the mutant side chain. Among
these are residues that contribute large positive energy terms,
which the bridging water terms presumably compensate for.
Finally, it seems important to emphasize that all protein
contributions to the Fab-Ag free energy change develop in a
regular and consistent way throughout the perturbation. This
correlates with the consistency of the model throughout the
dynamics and with the initial crystallographic data. Thus,
there are few backbone reorientations and limited side chain
conformational transitions in the Fab-Ag complex, despite the
absence of positional or internal constraints inside the dynam-
ical sphere of 30 A diameter. This solid-like behavior helps
in assigning structural origins to the free energy changes and
in understanding them from a simple knowledge of the
end-points.
In the unliganded Fab mutation, however, more than half a
dozen backbone carbonyl groups provide all the large, adverse
contributions to the mutation, emphasizing the importance of
the CDR and non-CDR loop turns' conformation to the free
energy outcome.
Discussion
The present study suggests that the mutation of Asn31L into
Asp in the antibody would lead to a 5.6 kcal/mol stabilization
of the HyHEL-10-HEL complex. If confirmed experimentally,
this would translate into a 10 000-fold increase in the binding
constant. This calculated value seems high, but is off only by
an order of magnitude from values observed in another
antibody-lysozyme complex (Lavoie etal, 1989) and in an
enzyme (Fersht et al., 1985, 1986) upon removal of salt links.
In the latter case, Fersht etal. (1985, 1986) observed a 1000-
fold decrease in enzyme specificity upon unpairing charges in
tyrosyl-tRNA synthetase, which corresponds to a destabiliza-
tion of 4 kcal/mol. This value is, in absolute terms, within the
estimated range of error for our calculations.
As far as the structure of antibody-antigen complexes is
concerned, the observations made in this study question two
assumptions regarding the mechanism of antibody-antigen
association. Firstly, the assumption of water exclusion from
the interface, since we found that a few water molecules are
instrumental in maintaining the crystallographic conformation
at the core of the complex, in agreement with recent experi-
mental observations on another HEL antibody complex (Braden
and Poljak, 1995). Otherwise, the observations made on
the basis of the crystallographic structure are confirmed by
simulations of the complex, as far as the nature and detail of
contacts involving the antibody L chain CDR loops and the
main chain conformation of these regions are concerned.
Secondly, the case for a lock and key association is weakened
by the observation of the consistent, if partial, reorientation of
a number of peptide bonds of the CDR2 and the non-CDR
loop turns in our model of the unliganded Fab. These findings
agree qualitatively with the small but significant CDR con-
formational changes observed in an important recent study of
the association of the antibody Fab 17/9 to a nonapeptide
antigen from influenza virus hemagglutinin (Rini etal., 1992
and references therein) and support the possibility of a limited
induced-fit in the mechanism of association of antibodies with
their antigen (Wilson and Stanfield, 1993).
The approximate decomposition of the free energy changes
provides some insight into the free energy change associated
with the point mutation. Mark and van Gunsteren (1994)
have warned against the overinterpretation of free energy
decompositions, since these depend on the path adopted for
the perturbation and should therefore be considered in relation
to the "k dependence of the Hamiltonian. For their part, Boresch
etal. (1994) have pointed out that as long as the integration
path is taken into account, such decompositions are useful in
the interpretation of results, because they shed light on the
underlying behavior of the system. In the present study,
a cautious approach was adopted. On the one hand, the
decomposition by energy type was given for the sake of
comparison with the Poisson-Boltzmann electrostatic calcula-
tions and in relation to the X dependence of the Hamiltonian,
so as to help evaluate possible sources of error. On the other
hand, the analysis by residue was presented in relation to the
conformation of the end-points.
Based on such correlation with structural information pro-
vided by simulations at the end-points of the thermodynamic
cycle, we can attribute the mutation free energy changes to
different factors: in the complex, the stabilization of the mutant
complex arises from the formation of a salt bridge with
Lys96W£L, with small or negligible changes due to limited
conformational flexibility. Charge formation does not lead to
a new, direct hydrogen bond between the charged groups.
Instead, Lys96W£i maintains three hydrogen bonds with main
chain carbonyl groups and a water molecule. In the unliganded
Fab, the solvent Coulombic term was seen to dominate
over its protein counterpart because the protein conformation
adopted upon full solvation of the Fab was largely unfavorable
to the point mutation. The competition of solvent for electro-
static interactions with protein main chain dipoles is well
illustrated in this system. The proximity of the solvent stabilizes
the charge, but it also forces nearby residues at the surface of
the antibody into adopting conformations that destabilize it
somewhat.
In short, the enhanced recognition of HEL by the mutant
HyHEL-10 appears to arise from the combination of thermo-
674
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from

Free energy simulations
dynamically more favorable conformational changes of the
CDR loops upon association and subsequent charge pairing in
the complex.
Reliable quantitative predictions of the molecular recognition
of large substrates using computational tools are seen to depend
critically upon a detailed knowledge of the structures at the
end-points of the thermodynamic cycle. Obtaining a full cut-
off correction to the electrostatic binding energy of flexible
systems using continuum electrostatics calculations remains a
computationally challenging problem. This may become more
feasible with the inclusion of recent developments of the
method (Luty et ai, 1992). More experimental data are needed
to verify the large increase in antibody-antigen affinity one
might expect from the calculations reported here upon forma-
tion of a salt link between an anti-lysozyme antibody and its
antigen. In the future, these could include structural data on
unliganded antibodies that may probe the limited extent of
conformational rearrangements suggested by our model and
thermodynamic studies of site-directed mutants of HyHEL-10.
A gene for a single-chain Fv of HyHEL-10 has been synthesized
(Tsumoto etal., 1994) and the HyHEL-10 antibody has been
expressed in Escherichia coli both as a chimeric Fab fragment
(K.Patel, S.Smith-Gill and R.C.Willson, unpublished results)
and as a single-chain antibody (Neri et ai, 1995), making the
results of the present study subject to experimental verification.
The energetics of HyHEL-10 association with HEL has also
been characterized by isothermal titration calorimetry (K.Shick
and R.C.Willson, manuscript in preparation).
Acknowledgements
We thank Drs Tjerk Straatsma, Shankar Subramaniam, Rebecca Wade and
Martin Zachanas for helpful suggestions. This work was supported in part by
grants from NIH, NSF, the Robert A.Welch Foundation, the Human Frontier
Science Program and the NSF Supercomputer Centers Metacenter Program.
References
Alzari.P.M., Lascombe,M.-B. and Poljak.RJ. (1988) Annu. Rev. Immunol, 6,
555-580.
Amil.A.C, Mariuzza.R.A., Phillips.S.E.V. and Poljak.R.J (1986) Science,
233, 747-753.
AntosiewiczJ., McCammonJ.A. and Gilson.M.K. (1994) J. Mol. Biol., 238,
415-436.
Bernstein.F.C, Koetzle.T.F., Williams.GJ.B., Meyer.E.F., Bnce.M.D.,
RodgersJ.R , Kennard.O., Shimanouchi.T. and Tasumi.M. (1977) J. Mol.
Biol., 112, 535-542.
Beveridge.D.L. and DiCapua.F.M. (1989) Annu. Rev. Biophys. Biophys. Chem.,
18,431^(92.
Boobbyer.D.N A., Goodford.P.J., McWhinnie.P.M. and Wade.R.C. (1989)
J. Med Chem., 32, 1083-1094.
Boresch.S., Archontis.G. and Karplus.M. (1994) Proteins: Struct. Fund.
Genet., 20, 25-33.
Bom.M. (1920) Z Phys., 1,45.
Braden.B.C. and Poljak.R.J. (1995) FASEB J., 9, 9-16.
Brooks.B., Bruccoleri.R.E., Olafson.B., States.D.J., Swaminathan.S. and
Karplus.M. (1983)7. Comput. Chem., 2, 187-217.
Chothia.C. and Lesk.A.M. (1987) J. Mol. Biol., 196, 901-917.
Dao-Pin.S , Sauer.U., Nicholson,H. and Matthews.B W. (1991) Biochemistry;
30,7142-7153.
Davis.M.E. and McCammonJ.A. (1989) / Comput. Chem., 10, 386-391.
Davis.M.E., MaduraJ.D., Luty.B.A. and McCammonJ.A. (1991) Comp. Phys.
Commun., 62, 187-197.
Dill.K.A. (1990) Biochemistry, 29, 7133-7155.
Fersht.A.R., ShiJ.-P, Knill-JonesJ., Lowe.D.M., Wilkinson.AJ., Blow.D M.,
Brick.P., Carter.P, Waye.M.M.Y. and Winter.G. (1985) Nature, 314, 235-
238.
Fersht.A.R., Leatherbarrow.R.J. and Wells.T.N.C. (1986) Phil. Trans. R. Soc.
Land., A317, 305-320.
Gilson.M.K., Sharp.K.A. and Honig.B.H. (1987) J. Comput. Chem., 9,
327-335.
Goodford.P.J. (1985) J. Med. Chem., 28, 849-857.
Honig.B.H., Sharp.K.A. and Yang,A.-S. (1993) J. Phys. Chem., 97, 1101-1109.
Horovitz,A., Serrano.L., Avron.B.. Bycroft.M. and Fersht.A.R. (1990) J. Mol.
Biol., 216, 1031-1044.
Jorgensen.W.L., Chandrasekhar.J., MaduraJ.D.. Impey.R.W. and Klein.M.L.
(1983) J. Chem. Phys., 79, 926-935.
Komeiji.Y., Uebayasi.M., SomeyaJ. and Yamato.I. (1992) Protein Bigng, 5,
759-767.
Kuczera.K., GaoJ., Tidor.B. and Karplus.M. (1990) Proc. Natl Acad. Sci.
USA.S1, 8481-8485.
Lavoie.T.B. etal. (1989) In Smith-Gill.S.J. and Sercarz.E. (eds), The Immune
Response to Structurally Defined Proteins: The Lysozyme Model, Adenine
Press, New York, pp. 151-168.
Lee.C. (1992) Curr. Opin. Struct. Biol., 2, 217-222.
Luty.B.A., Davis.M.E. and McCammonJ.A. (1992) J. Comput. Chem., 13,
768-771.
McCammonJ.A. and Harvey.S.C. (1987) Dynamics of Proteins and Nucleic
Acids. Cambridge University Press, Cambridge.
Mark.A.E. and van Gunsteren.W.F. (1994) J. Mol. Biol., 240, 167-176.
Meyer.E. (1992) Protein Sci., 1, 1543-1562.
Neri.D., Momo.M., ProsperoJ. and Winter.G. (1995) J. Mol. Biol., 246,
367-373.
NovotnyJ., Bruccoleri.R.E. and Saul.F.A. (1989) Biochemistry, 28,4735-4749.
Padlan.E.A., Silverton.E.W., Sheriff.S., Cohen.G.H., Smith-Gill.SJ. and
Davies.D.R. (1989) Proc. Natl Acad. Sci. USA, 86, 5938-5942.
Pearlman.D.A. and Kollman.P.A. (1991) J. Chem. Phys., 94, 4532^*545.
Pearson.R.G. (1986) J. Am. Chem. Soc., 108, 6109-6114.
Pomes.R. (1993) PhD thesis. University of Houston.
Pomes.R. and McCammonJ.A. (1990) Chem. Phys. Lett., 166. 425-428.
Presta.L.G. (1992) Curr. Opin. Struct. Biol., 2, 593.
Prevost.M., Wodak.SJ., Tidor.B. and Karplus.M. (1991) Proc. Natl Acad. Sci.
USA, 88, 10880-10884.
RiniJ.A., Schulze-Gahmen.U. and Wilson.l.A. (1992) Science, 255, 959-965.
Roush.D.J., Gill.D.S. and Wlllson.R.C. (1994) Biophys. J., 66, 1290-1300.
Sauer.R.T. and Lim.W.A. (1992) Curr. Opin. Struct. Biol., 2, 46.
Schiffer.M., Chang.C, Naik.V. and Stevens.FJ. (1988) J. Mol. Biol., 203,
799-802.
Sharp.K.A. and Honig.B.H. (1990) Annu. Rev. Biophys. Chem., 19, 301-332.
ShenJ., Subramaniam.S., Wong.C.F. and McCammonJ.A. (1989)
Biopolvmers, 28, 2085-2096.
Shenff.S., Silverton.E.W., Padlan.E.A., Cohen.G.H., Smith-Gill.S J., Finzel,B.
and Davies.D.R. (1988) In Sarma,R.H. and Sarma.M.H. (eds), Structure
and Expression. Volume I. From Proteins to Ribosomes. Adenine Press,
New York, pp. 49-53.
Slagle.S.R, Kozack.R.E. and Subramaniam.S. (1994) J. Biomol. Struct. Dvn.,
12, 439-456.
Smith-Gill.SJ., Lavoie.T.B. and Mainhart,C.R. (1984) J. Immunol.. 133,
384-393.
Straatsma.T.P. and Berendsen.HJ.C. (1988) J. Chem. Phys., 89, 5876-5886.
Straatsma.T.P. and McCammonJ.A. (1989) J. Chem. Phys., 90, 3300-3304.
Straatsma.T.P. and McCammonJ.A. (1991) Methods Enzymol., 202. 497-511.
Straatsma.T.P. and McCammonJ.A. (1992) Annu. Rev. Phys. Chem., 43,
407-435.
Straatsma.T.P., Zachanas.M. and McCammonJ.A. (1992) Chem. Phys. Lett.,
196, 207-302.
Subramaniam.S., McCammonJ.A. and Bacquet.RJ. (1989) In Smith-Gill.SJ.
and Sercarz.E. (eds), Tlie Immune Response to Structurally Defined Proteins:
The Lysozyme Model. Adenine Press, New York, pp. 169-176.
Tembe.B.L. and McCammonJ.A. (1984) Comput. Chem., 8, 281-283.
Tidor.B. and Karplusjvi. (1991) Biochemistry', 30, 3217-3228.
Tsumoto.K., Ueda.Y, Maenaka,K., Watanabe.K., Ogasahara.K., Yutani.K. and
Kumagai.I. (1994) J. Biol. Chem., 269, 28777-28782.
van Gunsteren.W.F. and Berendsen.HJ.C. (1990) Angew. Chem, Intl Ed. Engl.,
29, 992-1023.
Wade.R.C., Mazor.M.H., McCammonJ.A. and Quiocho.F.A. (1990) J. Am.
Chem. Soc., 112, 7057-7059.
Wade.R.C., Mazor.M.H., McCammonJ.A. and Quiocho.F.A. (1991)
Biopolymers, 31, 919-931.
WarwickerJ. and Watson,H.C. (1982) J. Mol. Biol., 157, 671-679.
Wilson.l.A. and Stanfield.R.L. (1993) Curr. Opin. Struct. Biol., 3, 113-118.
Received October 23, 1995: revised April 20, 1995: accepted May 8, 1995
675
at U
niversity of W
aterloo on D
ecem
ber 9, 2011
http://peds.oxfordjournals.org/
D
ow
nloaded from