c – The GROMACS development teams at the Royal Institute of Berendsen, Gromacs User Manual version beta1, (). Refer to the GROMACS 4 publication and the manual for details. As of version , GROMACS supports the use of GPU accelerators for running MD. This tutorial focuses specifically on issues related to dealing with the ligand, This tutorial assumes you are using a GROMACS version in the x or x.
|Published (Last):||17 May 2005|
|PDF File Size:||11.5 Mb|
|ePub File Size:||18.18 Mb|
|Price:||Free* [*Free Regsitration Required]|
In this tutorial we will try to obtain the free energy of binding of n-phenylglycinonitrile to T4 lysozyme using an alchemical pathway, in order to reproduce the result obtained in the published work of Boyce et al.
groacs The tutorial assumes knowledge of Gromacs 4. The thermodynamic cycle used to obtain the results will be reviewed first, and then we will go through the principal steps of the process. The input files needed to run the calculations will be provided, however you can also set everything up from scratch by yourself if you wish it can be a good exercise in fact.
For enquiries about the tutorial please feel free to email me. Here is a graphical representation of the cycle we will use in this tutorial, as I find it helpful to visualise the steps needed to carry out the calculations. It’s purpose is simply to summarise the steps we gromac need to consider – if this seems confusing, I would suggest visiting the following pages before proceeding with the tutorial: In the cycle above, the systems we need to simulate are indicated by having a black box around groamcs, restraints are indicated by a red circle, the transparent ligand means it is not interacting with the environment and the light blue background is reminding that water is present.
The right column has the simulations involving the complex, whereas the left column the simulations involving only the ligand. Starting from the top-right corner we have the complex, with the ligand n-phenylglycinonitrile and protein T4 lysozyme fully mabual as in a normal MD simulation. Now we want to decouple the ligand from the system in order to get to the bottom-right corner of the cycle.
The presence of restraints is indicated in the cycle scheme above by a red circle, which is trying to represent the fact that the ligand is being confined to a certain volume. The set of restraints described by Boresch is used for this work. In practice we will be doing the opposite, that is decoupling the ligand from the water box; however, note how this means running basically the same set of simulations.
This also results in the same free energy difference, just with opposite sign.
At this point we are finally at the top-left corner of the cycle, which means that summing up all the steps done so far we are going to obtain the quantity we are after: Or, in fact, its opposite since we went from bound to unbound state.
In order to have everything ready to run the windows at different lambda values we need to have the input files ready. Briefly, for this exercise the setup is as follows: More details about the parameters are in the. You can find the input files needed to complete the tutorial here: Note that in the. This is explained also later on in the text.
For consistence with the description of the cycle above, we will start taking care of the simulations involving the complex. This set of simulations involve a total of 30 windows, and for each of them energy minimization, NVT and NPT equilibration, and production run 1 ns have to be performed. This .46 we can use 30 identical.
As it is possible to see, we first apply the restraints, then we remove the coulombic and finally the Van der Waals interactions. We are now ready to run all the simulations and collect all the dhdl. Then we can analyse the results with the alchemical-gromacs.
GROMACS 4.6 example: n-phenylglycinonitrile binding to T4 lysozyme
To see the help for the script just run it with the ‘-h’ flag, so that it will show the meaning of all the options too.
Here a link to a summary of the results: Note that this might not be exactly the gromcas you will get even if you use exactly the same input files, parameters and GROMACS version I used but it gives an idea of what sort of value one should expect.
In fact, running two additional complex calculations, I obtained a standard deviation for the three runs of 0. Now the ligand is decoupled from the protein and the solvent. The next step in the cycle is to remove the restraints. As mentioned before, with this choice of restraints we can derive the free energy of restraints removal analytically.
The formula to use is the following,  but I would suggest to first read the relevant paper to understand what it means. If you are using exactly the set of restraints provided with the input gromads, this should give:. The same considerations discussed for the complex simulations apply for the ligand as well. State B has zero charges and dummy atom types and we need to run 20 simulations of the ligand in water, at the different states defined by the lambda vectors in the.
In this case we are only turning off coulombic and Lennard-Jones interactions, since the restraints have just been accounted for analytically.
The results can be obtained again in the same way as we did for the complex. Here a link to what I obtained after running alchemical-gromacs.
GROMACS Online Reference
At this point we gromacd carried out all the simulations we need to reproduce the result obtained by Boyce However, the authors also applied a long range dispersion correction to their final results.
Again, the best thing to do to have an idea of why and how this is done is probably to read the original paper. This can be done using the -rerun flag for grmacs and providing the trajectory file, for example:. Since we are reprocessing with a long cut-off only the infrequently stored configurations, and we would like to compare the same configurations with long and short cutoffs, we can rerun the two lambda states also with the original.
This value will be added to the free energy computed for the complex and will contribute to a higher affinity. Summing everything we obtain the free energy of dissociation, so if we want the free energy of binding we simply take the negative of it. This, in my case, results in:. The published result is However, it is possible you’ll get something slightly different. In fact, considering the final free energies I would obtain if I used the results from the other two complex simulations and their respective EXP-LR corrections Free Energy Fundamentals Theory.