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

[chris.neale at utoronto.ca]

> I am trying to calculate the PMF of a drug inside a membrane protein
> which forms a channel. The reaction coordinate is the z-axis. In  
> order to do so, I would like to use Umbrella sampling of the drug in
> many positions along the z-axis, followed by WHAM. That should give me
> the PMF for each position along the pore, inside-out, right?

yup. Assuming you want to see the free energy surface for extraction  
along the z-axis.

> 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.

> Then Umbrella sampling should be performed on each of those structures,
> so I should also have 30 files of pull.ppa. The difference between them
> is the position to which the pulled group (drug) is restrained to. Does
> that sound reasonable so far?

yup

> This is how my pull.ppa looks like: --
> verbose                  = no
> runtype                  = umbrella
> ngroups                  = 1
> group_1                  = AMAN
> ;reference_group          = weights_1                =  
> reference_weights        = reftype                  = com
> reflag                   = 0
> pulldim                  = Y Y N
> ; DYNAMIC REFERENCE GROUP OPTIONS
> r                        = 0
> rc                       = 0
> update                   = 1
> ; UMBRELLA SAMPLING OPTIONS       K1                       = 1000
> Pos1                     = 33.482  31.225  33.173
> --

looks fine, assuming that you are pulling out in XY. At the top you  
said that the rx coord is Z, so in that case your pulldim here is  
incorrect. I didn't inspect it thouroughly though (that's your job) so  
don't take my word that it is error free.

> K1 is the force constant acting on the pulled group in order to restrain
> it, right?

yup.

> 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.

> For Pos1 - is that OK to input the COM (center of mass) of the ligand,
> as it is in each starting structure? Since this is the position I would
> like to restrain it to.

Yup, assuming that you have made your starting structures appropriatly.

> 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.

> 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.

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.