Subject: Responses to my query and a question

From: Viswanadhan <[email-masked]>
Date: Mon, 4 Mar 1996 21:05:51 +48000

Dear Netters,

In response to my question, I received some responses summarized below.
The question was:

-------------------
We have a fairly rigid cyclic molecule (C) which can exist in two 
conformations (puckers), C1 and C2.

Its acyclic analog (A) can also exist in different conformations, but
they are not separated by high barriers. C1 and some conformers of A
appear to have good interactions with a protein, but not C2.

In terms of binding geometry/energies C1 is better than A but A is 
much better than C2. C1 and C2 are separated by a high barrier, but 
they are isoenergetic (internal energy from QM calculations). 

TCP (free energy perturbation) calculations for the  mutation 
C1 -> A (or C2 -> A) will not reveal much because of the high barrier 
between C1 and C2 (making interconversion almost impossible in a 
reasonable time). 

Do you think it is possible to qualitatively assess the binding of C vs. 
A ?

-----------------

The question I have for today is what is the best and most efficient method
to calculate the transion state barrier (and free energy profile) between 
two puckers for a cyclic molecule. I am aware of LST in Gaussian. 
I need a method that is accurate and works well for cyclic structures even if it 
takes a good deal of computer time. I will summarize the responses.

-----------------

This is what I think about my previous query and responses. 

While the answers are interesting, it seems to me that if I can correctly
compute the barrier and the torsion profile, I may guestimate the relative
frequency (or probability) of interconversion (C1 to C2) relative to the 
acyclic molecule with a smaller barrier. Computing interaction energies in the complex
and in the solvent for each molecule (one conformation in the complex as the start,
and both conformers in the solvent) would enable the estimation of relative
difference in binding enthalpy. If C is much better than A, and the torsional barrier
is not too high, I would expect that C should be better than A as an inhibitor,
because A would lose out on the entropic argument as well. While there is some subjectivity
in this process, it will lead to better design ideas such as how the binding conformation
can be further restricted, and what changes in the structure would lower the barrier etc.

Monte Carlo methods are promising however, they are not easily applied to endocyclic
torsions when dealing with protein-ligand complexes. The work of Dr. Still and coworkers
is useful as a start in this respect.

Vellarkad Viswanadhan

[email-masked]

------------------------
Responses and Comments.

Regarding your C1, C2, A problem: if you could
compute the difference in chemical potential (free energy)
between C1 and C2, you'd be all set.  It might
be possible to compute the work of forcing this
transition, via a free energy calculation.  

Mike Gilson

Computing the barrier would certainly help. However, it is not easy to accomplish
the free energy calculation via umbrella sampling etc, though in principle this 
should work.

From: [email-masked] (Michael K. Gilson)
To: [email-masked], [email-masked]
Message-Id: <[email-masked]>
Content-Transfer-Encoding: 7BIT
Status: RO
X-Status: 

Dear Viswanadhan,

If you know the equilibrium constants for binding of C1 and C2 to the
receptor, and also that for the interconversion of C1 and C2, then you
can use algebra to compute the "overall" or "apparent" equilibrium for
binding of C to the receptor.

Does this sound right?

Mike

To: Viswanadhan <[email-masked]>
Status: RO
X-Status: 

  It is theoretically possible to do what you want directly with Monte
Carlo calculations rather than FEP theory.  For example, you can sample
conformations of C1 and C2 in the bound form and determine the ratio
existing at equilibrium and do the same without the protein.  This can
be done with continuum solvation models such as Clark-Still's.  This
kind of work is going on here at Mt. Sinai for peptides in water but
could be extended to include protein.  For the other problem of the
binding of A versus C you could do a Monte Carlo calculation sampling
AP versus A + P and another sampling CP versus C + P.  You can use
internal coordinates in the sampling and restrict the range of values
sampled to include only those which exist in the free and bound forms
of these molecules.  This improves the sampling efficiency of Monte
Carlo greatly over blind approaches and local sampling methods such
as force bias.  I'm not aware of anyone having tried exactly what I've
described yet but it might work.
                               Sincerely, Dave Garmer
                               Physiology & Biophysics Dept.
                               Mt. Sinai Medical School
                               New York, NY 10029

We have had positive experience with GB/SA model in water for small 
molecules. However, protein-ligand calculations have strong limitations 
with continuum models, especially if the ligand is at the solvent 
interface and it is charged. If we let the side chains be flexible (some 
of which are charged), you will see big movements and large rms 
deviations. Our own work on a database of twenty proteins convinced us 
that improvements are necessary in the methods to apply MC 
methods+continuum models successfully (for the bound protein).

However, your approach should work with good params and models.

--------
Also,

Well, you expect the polar and charged groups in solution to be 
surrounded by ions in solution. Since we do not put those, the sidechains 
tend to make internal salt bridges which are unrealistic. This will also 
be at the expense of right geometry for ligand and protein sidechains 
(interactions) since these are at the surface. With regard to the 
internal sites, if they are far away from solvent, why do you think 
solvent interactions would play an important role in binding? It is a 
quite different matter if you are modeling the whole flap opening or 
something like that.

From: Dave Garmer <[email-masked]>
To: Viswanadham <[email-masked]>
Status: RO
X-Status: 

Dear Viswanadham,
  To answer the easy question first, there is evidence of ordered water
molecules around the retinal in rhodopsin and the modeling work also
suggests that there are enough unsatisfied hydrogen bonds to generate
a water network around GPCR ligands.  I would guess that an exposed
site would have less ordering and so better described as a polarizable
continuum.
  Regarding crystallographic protein surface details, the requirements
for crystallization usually mean that the local ion concentration is
high compared to an organism.  Also, crystallization is converting
molecules which exist in a liquid environment into a solid, thus
probably eliminating some inherent flexibility.  Also, crystallographers
try to assign what they think is the best single structure and very
seldom are able to identify what might be minority structures for small
parts of the protein.  Theoretical studies have indicated that there
probably are minority structures often existing even for internal
packed sidechains.  Therefore, I would trust crystallography only to
suggest particular accessible structures for flexible units.  I would
also think that a decent continuum model like Clark-Stills probably
can more often than not give you a sample of some of the thermally
accessible structures.  They may not be the same as the ones being
sampled in a crystal structure refinement.
                                Dave Garmer

From: "E. Lewars" <[email-masked]>
Message-Id: <[email-masked]>
To: [email-masked]
Subject: SUBSTRATE CONFORMERS/PROTEIN
Status: RO
X-Status: 

It seems to me you should first try to find the global minimum. If this is 
also the best-binding conformer, then the barriers to conformer interconversion
are not relevant to the problem (I think). If the global minimum is _not_ the
the conformer that binds strongest, then you might have to write out
a scheme involving all the species and analyze it as a problem in (mainly)
kinetics.  Sorry, but I haven't looked at this in detail.
Errol Lewars
====

This is rigid molecule (C) with only two minima which are equivalent 
(internal energy-wise). Only one is the binding conf. The question is C 
vs. A. I am not sure what you precisely mean by 'problem in kinetics'. Is 
there any similar problem that is analysed? 

From: "E. Lewars" <[email-masked]>
Message-Id: <[email-masked]>
To: [email-masked]
Subject: CONFORMERS
Status: RO
X-Status: 

Hello,  I thought that the acyclic molecule A was in equilibrium with C;
evidently this is not so, i.e. A is a quite separate species.  If C1 and C2
are in equilibrium and only C1 binds to the protein, then if enough protein
is present all the C1 will get tied up and then C2 will (slowly, if the 
C2--> C1 barrier is high) drain away to C1.protein. To see if this 
happens on a useful timescale, you have to calc the C1/C2 barrier and estimate 
the half-life of C2. Ab initio calcs give more reliable results for
barriers than do semiempirical calcs.
Good luck
Errol Lewars
===