From: Jan Saam (saam_at_charite.de)
Date: Tue May 30 2006 - 12:38:05 CDT
Dear Jerome,
the discussion is getting quite interesting.
>>May I bug you with one more question?
>>I'm studying the diffusion of oxygen along a channel in a protein. The
>>channel is not straight but follows despite some kinks roughly along a
>>certain direction from the active center to the protein surface. There
>>are also some "bays" or dead end roads along the tunnel that can be
>>explored by oxygen.
>>
>>Thus, most of the time the gas molecule will not travel exactly along
>>the defined reaction coordinate.
>>Does that mean the following?:
>>
>>a) If distance active center -- O2 is the RC:
>>Oxygen could be temorarily trapped in a dead end road, its hammering
>>against the wall and the adapdive biasing force is cranked up until
>>oxygen breaks through that wall? And even when it's travelling back,
>>finding the right way, the time spent in the side channel will distort
>>my PMF?
>>
>>
>That is, to a certain extent, true. In such cases where the pathway is roughly
>described by the RC, but some details are not resolved, it usually boils down
>to a problem of timescales. If oxygen is able to diffuse fast enough between
>the main route and small dead ends surrounding it, then it will find its wasy
>and the PMF should not be distorted, because the time spent on the right
>pathway will dominate the average. The freedom to move back and forth along
>the way is crucial there. In an SMD simulation, it would be extremely hard
>for a system to return to the right track.
>If on the contrary, diffusion is a bit too slow, then yes, the system may
>become trapped in the dead end, and the adaptive bias will then grow to
>unrealistic values to try to pull oxygen through the wall.
>One parameter that may determine which of the two cases prevails is
>fullSamples. You want the system to find the right track *before* the bias is
>fully applied. That's where things become really system-dependent.
>
>
>
>
Within the wildtype channel oxygen diffuses more or less freely. But
since I try to investigate the effect of mutations on the transport
properties of the channel it could well be that the dead end is chosen
to be the better way and oxygen will be trapped. Probably something
better is needed.
>>b) If abscissa along main direction is used as RC:
>>When I'm defining a direction of the RC there will always be a certain
>>angle (say between 0 and 50 deg.) between the real channel V and the RC
>>due to the crookedness of the channel. I.e. the adaptive biasing force
>>is pushing with that angle against the channel wall, which will also
>>distort my PMF. Correct?
>>
>>
>I would say, no, at least not if the simulation is in near-equilibrium. Again,
>fullSamples will play a role on that. What matters is that there is a
>one-to-one mapping between your RC and the actual progression of the system
>along the pathway, that's enough to makes the RC relevant. And that's also
>why alternate pathways or dead-ends are very hard to take into account.
>
>Another problem with this, though, is to define an absolute direction. Is your
>protein able to tumble or rotate in some way?
>
>
>
Unfortunately yes, it will tumble. We should urge Jim to provide a
function to remove the COM motion constantly throughout the simulation,
not just once at the beginning. That would heö a lot defining RCs.
>>c) If I'm partitioning the channel into several shorter, almost straight
>>portions, I can stitch the PMF together and probably improve on these
>>problems. But I'll have to know the RC very well before and it is not
>>easy to define the directions in a system where everything moves.
>>
>>
>I would say the only useful coordinates would then be defined relative to the
>protein itself. The available RCs are probably not adapted for that, but we
>can talk about what would suit your needs. It would probably be helpful to
>many other people, too.
>
>
I think it would be cool to spiff up ABF (which I really have to say is
strikingly elegant!) so that it automatically find the reaction
coordinate and hence transition also states.
This might not be possible for all types of reactions but in my case I
could envision the following:
You probably know "conformational flooding" due to H. Grubmüller or the
similar method "metadynamics" by Parinello. I'll quickly outline
metadynamics since it's suited better to our problem: The basic idea is
that you obtain information about your free energy profile along a
certail RC by adding a a little gaussian function to the total force
field potential for the current state. Each sample adds one more little
gaussian contribution to a potential that builds up slowly filling the
energy minimum you are sitting in. Figuratively speaking you are
dropping litter into your well until you are carried out of the pit by
all that trash. The potential that you built up by then is something
like the inverse of the PMF along the RC (very elegant method, too!).
I think this gaussian accumulation technique could also be used to
actually *find* the most likely RC. At least in such a simple case as a
diffusiing particle in a channel this should be possible.
But may be there's even a simpler way without running a simulation
before to get the RC, since ABF has the same effect as the gaussian
accumulation.
Wouldn't it be possible to replace Xi by a 3-dimensional analog Xi_j ;
j=1,2,3 so we would get:
dA(xi_j)/dXi_j = <partial(V(x))/partial(Xi_j) -
1/beta*partial(ln|J_j|)/partialXi_j>_(Xi_j) = -<F_Xi_j>_Xi_j
(Underscores denote subscripts)
We would have j^3 bins but we would only use a small fraction of them,
namely the ones that have more than a certain number of samples. So we
would get a 3D potential map of the accessible conformations and the
longer the simulation is run the more the 3D subspace is explored.
Since the bins with no samples in them are even in the 1D case treated
as having the highest energy while all other bins have a lower energy
this is ideal to be visualize using a volume denstity isosurface in VMD
(you'd have to invert the potential first). Then, depending on the
chosen free energy isovalue you'd see blobs or tubes with an energy
lower than that.
The problem here is again, that one would have to make the RC grid
relative to a protein coordinate system or constrain the COM motions and
rotations of the system.
I also thought about case a) again. When the system has spent some time
ina dead-end even if that's in a potential well this minimum will be
leveled out after enough samples have been taken and the particle can
travel back onto the main road it will finally first break the lowest
potential barrier wether it's on the branch 1 or 2 of the way.
What do you think about this. Do you think it could work? Would it be
terribly hard to implement this?
Gosh, that was a long mail...
Cheers,
Jan
>Best,
>Jerome
>
>
>
>
-- --------------------------- Jan Saam Institute of Biochemistry Charite Berlin Monbijoustr. 2 10117 Berlin Germany +49 30 450-528-446 saam_at_charite.de
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:42:07 CST