[gmx-users] Umbrella sampling of a ligand inside a pore

[Hadas Leonov]

First of all - Thank you very much for your thorough reply!

> > I am having trouble defining the parameters for the simulations. As  
> > I understand, I should generate many starting structures, in which
> > the position of the drug differ by, say, 0.5 Angstrom apart on the z
> > axis. So let's say I have 30 starting structures (protein + drug + lipid
> > bilayer + water).
> 
> Be sure that you generate these carefully. Don't yank them all from  
> one initial starting structure or your huge forces bay severly change  
> the structure in some instances. A stepwise method of generation would  
> be better.
What I did was to position the ligand in the center of the protein pore,
and consecutively move it by 0.5 steps along Z (both up and down until
it reaches the bulk on both sides). In each step I re-center the ligands
x,y inside the pore to avoid clashes (using nearby residues). Then, each
of the generated structure is energy-minimized. Hope it's careful
enough.


> > how do I determine the best value for it?
> 
> trial and error. You must have enough overlap to reduce WHAM errors.  
> However, too much overlap and you will end up prolonging simulation  
> time. But that's not too much of an issue if you have lots of cpus  
> since it won't be wasted sampling. If you understand WHAM and can code  
> yourself then you won't need to waste any data if you change force  
> constants or nominal positions half-way through. Basically if you have  
> tons of cpus then just overkill the overlap. If not, spread as many as  
> you can and then make histograms of P(rx) vs. rx and measure overlap  
> (where rx is your reaction coord. e.g. z) Then in places where there  
> is not enough overlap you should do new simulations with umbrellas  
> that are in intermediate positions.
> 
> Don't like trial and error? I am sure there are better ways. Look into  
> adaptive umbrella sampling.
> 
I'll probably perform trial and error before I move to more advanced
methods. Do I have to optimize the force constant for each simulation or
is it OK to use the same for all?

> > Pulldim = Y Y N, since for each starting structure, the point is to see
> > if it moves along x,y due to the forces acting on it by the protein.  
> > As I understand, since the ligand's z position is different in each
> > starting structure, we can analyze it with WHAM and calculate the PMF.
> > Does that sound right, or am I completely wrong?
> 
> Completely wrong. Pick a reaction coordinate and then have different  
> starting structures along that reaction coordinate and also  
> harmonically restrain to that reaction coordinate. If the reaction  
> coordinate that you want is xy then your starting structures should  
> have the xy close to their Force=0 value. In this case either give all  
> z the same value or whatever you want, but be sure to make sure that  
> this has converged. If the reaction coordinate is z then "as it moves  
> along x y" makes no sense since you won't get any information about  
> this.
> 
That's the part I'm confused about.
The reaction coordinate is definitely the z axis. That's why I have 30
different starting structures of the ligand along Z inside the pore.
But when you configure pulldim=N N Y, you enable it to move only in the
z-axis. So you won't see movements of the ligand in the xy plane inside
the pore (surely the harmonic restraint of x and y wasn't the optimal
x,y). So isn't it better to use the default pulldim=Y Y Y in that case?

The confusing issue is that it seems we're already performing the
pulling by generating overlapping starting structures, and the remaining
question is what's the free energy when the ligand is around that
specific Z (for each Z), yet still let it move freely within XY.

> > The result of a simulation with the aforementioned parameters is that
> > the ligand escapes the protein *into* the lipid environment, which
> > doesn't make sense. When I use K1=10, the drug does not escape the pore.
> > So I'm confused here. I thought K1=1000 was supposed to put more force
> > on the drug to stay in its original position.
> 
> Yes, larger Fc puts more force. But what happens depends on where  
> 33.482  31.225  33.173 is. Your simulation is definately in error if  
> the protein pore is at X=33.482, y=31.225 and a K1=100kJ/mol/nm^2  
> pulls it far away from that position. If so, are you sure that you  
> don't have a negative sign in the force constant?
> > Any help and ideas regarding why this is wrong?
> 
> You are definately confused about z vs. xy and what is your reaction  
> coordinate. I suspect that although you have ensured that z=33.173 is  
> in your pore, you did not ensure that X=33.482, y=31.225 is in your  
> pore. If this is true and you can't understand why that's a problem  
> then you need to focus on the definition of a reaction coordinate.
> 
Now this feels stupid. I Forgot that gromacs uses nm units, and gave it
coordinates in Angstroms.

Thanks again,
Hadas. 

> If this is not the case and you can't figure it out, then reduce your  
> system to something more managable for a test. A single water molecule  
> in vaccuo using the sd integrator should be a very fast test system  
> for spatial umbrella testing. If you don't like that idea then put  
> your ligand in a small box of water or some other idea that simulates  
> faster than a bilayer system.
> 
> > Thanks in advance,
> > Hadas.
> 
> You're welcome,
> Chris.