[Prev][Next][Index][Thread]
WSN: RE: calculation of binding energies
From: Mill Lambert <mhl8097@glaxo.com>
In response to Bruce Bush's request, I can say something about our
experiences at Glaxo here in North Carolina. We have applied extensive
binding energy calculations in several structure-based drug discovery
projects. In all three cases, solvation was represented with a model
similar to that of Ooi and Scheraga (PNAS 84, 3086, 1987) where each
atom type is assigned a solvation parameter measuring the free energy
per square angstrom of solvent accessible surface area. The original
Ooi and Scheraga parameters were based on free energies of transfer from
vacuum to water for about 25 model compounds. I have used the original
parameters, and have also experimented with alternative
parameterizations incorporating more model compounds, as well as the
controversial volume entropy correction (Sharp, Nicholls, Friedman and
Honig, Biochemistry 30, 9686, 1991). Interactions with the protein are
represented with the Discover CFF91 force field excluding cross terms,
as implemented in my moelcular mechanics program MVP. The MVP program
uses a generalization of Scheraga's build-up procedure to generate
500-1000 alternative low energy conformations of the ligand within the
protein binding site, and then uses a filtering minimization procedure
to select subsets of 50-100 structures for minimization to rms gradients
of 0.1 kcal/angstrom. To include the solvation model in the
minimizations, I have coded an analytical algorithm within MVP for the
solvent accessible surface areas and their gradients. To get meaningful
binding energies, I carry out two build-up calculations, one in the
binding site and one free in solution. The estimated binding energy is
then the difference
dG = - ktlog[sum1 exp(-E/kt)] + ktlog[sum2 exp(-E/kt)]
where the first sum is over structures in the complex and the second sum
is over structures in the noninteracting protein-ligand system. If the
protein is held fixed, then the second sum reduces to a sum over
structures of the the ligand free in solution plus the solvation energy
of the protein by itself. However, we typically subject one or more
side chains in the protein to simultaneous conformational search, in
which case both sums must include the energies of the movable atoms in
the protein. The energies of movable atoms, in either the protein or
the ligand, include all interactions with other movable atoms, as well
as interactions with the fixed atoms of the protein. A calculation of
this sort requires 1-8 hours on an R4400, depending on the number of
structures requested and the restraints imposed. We have automated the
determination of ionization states, charges and atom types for the
functional groups we needed in our projects, so the setup time is fairly
quick. Normally, we build and submit molecules one at a time, possibly
running three or four molecules per night on idle Silicon Graphics
machines, depending on the interest and suggestions from medicinal
chemists. However, in one collaboration with Paul Charifson, we used 26
Silcon Graphics machines to carry out docking calculations on 1200
candidate compounds. The compounds were generated by selecting 1200
compounds from the fine chemicals database, generating structures, and
then substituting each into the same position on the lead compound. The
resulting binding energies were used to help select monomers for a
chemical library.
So what kind of results do we get? Of the three projects, our best
results, obtained very recently, gave a correlation coefficient of about
0.6 between calculated binding energy and measured Ki. This calculation
involved about 15 compounds within a particular series, with Ki values
ranging over 5 log units. More typically, calculations give correlation
coefficients of 0.2-0.5 for sets of 10-60 compounds within a series.
The correlation generally deteriorates if the range of Ki values is less
than 3 log units, as is often the case with our data. The correlation
breaks down completely if compounds from different series are included
in the same plot of energy against IC50. The correlation also generally
breaks down when different heteroatoms are substituted at functionally
important positions. I would presume that this is due to inaccuracies
in the solvation model. We have experimented with more sophisticated
solvation models, including the use of DelPhi reaction field energies,
but have had no success to date.
Mill Lambert mhl8097@glaxo.com (919)-941-3344