Analysis of FEP-REST2 simulations

[galmas]

Hi All,
I'm experimenting with FEP-REST2 simulations using gromacs-2020.6 and plumed-2.7.2 for retrieving the relative free energy of binding of two ligands to a protein.

Does anyone has some experience with analyzing the results of such simulations?

I guess I cannot simply use gmx bar with the dhdl.xvg files, like I would if these were standard FEP simulations, since the topologies for each replica are scaled.

Any help would be much appreciated!

Cheers,
Gal.

[Giovanni Bussi]

Hi,

If you are using separate topologies and hrex flag (and not the lambda perturbation in GROMACS)  then you should be able to use the WHAM scripts as explained in PLUMED tutorials,

 though you will have to produce the equivalent of the bias files using a gmx mdrun -rerun. It's not really documented anywhere I think, but we used that approach several times and it works.

Giovanni

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/fba6221b-a995-49d4-90cc-9585a4fe19ban%40googlegroups.com.

[galmas]

Hi Giovanni,
Thanks for the quick response, it is very much appreciated!
I’m not sure that I understand though how to obtain the bias potential with gmx mdrun -rerun.
Should I rerun each lambda leg using the unscaled topology (i.e. the one used for lambda_0) and compute the potential energy or dH/dl?

Also, just to make sure that what I want to do is legit, this is what I had in mind:

1. Produce hybrid topology for ligands A and B using pmx.
2. Producing processed topologies using grompp -pp and modified plumed-partial_tempering that can process the hybrid topology file. In the processed topologies the effective temperature is increased half-way from A to B, and then decreased back, so that both A and B are simulated at the same effective temperature.
3. Doing FEP-REST2 with gromacs+plumed (using -hrex flag), with different topologies and mdp files for each lambda leg. Topologies should control the effective temperature and the mdp files the transition from ligand A to ligand B.
4. Post processing with WHAM – this part I’m trying to figure out.

Cheers,

Gal.

[Giovanni Bussi]

Hi,

so you plan a single ladder right? I mean: going from state A without REST2 to state B without REST2, passing through a series of intermediate states that are both alchemically perturbed (using lambda) and with modified parameters (using REST2).

I think the acceptance with -hrex would be correct but I am not 100% sure. What I did in the past was to run only -hrex (without lambda), e.g. to remove a single proton from the system (pKa calculations). In this case, you use rerun to compute the difference between the energy in each of the simulated topologies and the reference one, and you feed that to WHAM. If two replicas have both a different lambda AND different force field parameters you have to double check if the acceptance is correct. Once you checked this, I think you can use the same procedure, computing the energy associated to each topologies and its difference wrt to reference replica.

Notice that the WHAM scripts included with plumed might suffer numerical problems in cases where there are very large difference in the energies (that's the case when you do alchemical changes, with DG of the order of hundreds kcal/mol and much smaller DDGs). If you have that problem, I can point you to an alternative WHAM script that we use now in my group and that solves this problem.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/b1e03a33-2ba5-4d21-9115-2fde58ea43ffn%40googlegroups.com.

[galmas]

Hi Giovanni,
Thanks again for the reply, it is extremely helpful!

Yes, I’m attempting a single ladder.
I’ll attempt the approach you suggested – rerun with the reference topology, compute the difference in energy between the simulated topology and reference topology and feed it to WHAM.

So far, I did a test run for a pair of ligands with 15 lambda windows and effective temperature that ranges from 300K to 400K (increasing halfway and then decreasing as I wrote above). Simulations lasted for 5 ns attempting exchange every 1000 steps.
Attached are the exchange statistics that I got. Looks better for the protein-ligand complex system. Maybe I should add more lambdas at the edges…

Cheers,
Gal.

[Pallav Sengupta]

Hi Galmas and all,

Have you been able to do this FEP/REST2 with GROMACS,

if yes can you provide us with some insight in doing so.

Like in normal alchemical transformations:

1. I would go about making the hybrid topologies and structure with pmx tool

2. Then do the transformations in respective lambda windows.

will this be like replica exchange in corresponding lambdas in parallel simulations or this would just sample in that particular window?

It would be helpful if someone has performed this sort of calculations could comment on the matter.