Hi,
I was hoping to perform RBFE calculations where I change Phenylalanine (PHE) to an unnatural amino acid 2-methylPhenylalanine(MPH2) using CharmmFF. The mutations supported by pmx are all natural amino acids. I searched for materials on creating dual topologies for non-standard amino AA mutatios but could not find one. That’s why I hoped I could generate dual topologies for PHE–>MPH2 using the ligand approach and use that in my calculations. While using the ligand approach I am running into problems at the ligandHybrid step.
What I did is, I generated the topology files for ACE-PHE-NME and ACE-MPH2-NME using pdb2gmx and added the parameters explicitly in the itp file. Next I followed the steps for ligand mutation and got stuck at the ligandHybrid step.
The last few lines of my hybrid.log is as follows:
ligandHybridTop__log_> Dummies in stateA:
ligandHybridTop__log_> Dummy…: 7 DUM_HGA3 | 0.00 | 1.01 → HGA3 | 0.09 | 1.01
ligandHybridTop__log_> Dummy…: 8 DUM_HGA3 | 0.00 | 1.01 → HGA3 | 0.09 | 1.01
ligandHybridTop__log_> Dummy…: 9 DUM_HGA3 | 0.00 | 1.01 → HGA3 | 0.09 | 1.01
ligandHybridTop__log_> Dummies in stateB:
ligandHybridTop__log_> Construct bonds…
ligandHybridTop__log_> Error: Something went wrong while assigning bonds!
ligandHybridTop__log_> A-> Atom1: 11-HN Atom2: 12-CA
ligandHybridTop__log_> B-> Atom1: 8-HN Atom2: 9-CA
ligandHybridTop__log_> Exiting…
I believe it is trying to make a bond between HN and CA but cannot figure out why. The command I used
pmx ligandHybrid -i1 phe.pdb -i2 mph2.pdb -itp1 phe_updated_cleaned.itp -itp2 mph2_updated_cleaned.itp -pairs pair1.dat -oA mergedA.pdb -oitp merged.itp --fit
I am attaching the pdb and itp files for reference. Any help will be much appreciated. I could not attach the topology files with .top or .itp extension and hence changed the extension to .txt.
phe.pdb (2.6 KB)
mph2.pdb (2.8 KB)
phe.txt (13.4 KB)
mph2.txt (14.5 KB)
Also if anyone can share some resources on how to perform non-standard amino acid mutations, it will be extremely helpful.
Hi,
The pdb and itp files have different sets of resname, ACE-PHE-NME and ACE-MPH-NME. Could you just rename them to the same resname (say MOL for both tripeptides) and try again?
Also, if the error persists, please provide the same itp and pdb files used for us to reproduce the error.
Another way of dealing with non-standard amino acid would be to parameterize an amino acid first, add it to .rtp, adjust genlib script (pmx genlib) and create an entry in the mutation library.
All the best.
Hi Sudarshan,
Thank you very much. Since the previous method did not work, what I did was to change the mutres.rtp and mutres.mtp file. Then I make the mutant pdb manually do gmx pdb2gmx to generate topology and then pmx gentop seems to work fine. It generates the necessary hybrid topologies required to do the simulations. However there is one more concern that I have.
In these simulations, I sample several lambda windows which generate the xvg files. Then I use gmx bar to calculate the deltaG of transformation. In the output I get a warning which says
WARNING: Some of these results violate the Second Law of Thermodynamics:
This is can be the result of severe undersampling, or (more likely)
there is something wrong with the simulations.
At first I thought it may be because I made some mistake while changing the mutres files but then I noticed that this error is in some of the mutations of the natural amino acids. For example I am working with a 15 residue peptide binding to a protein. When I change Leu to Val I have no error but Leu to Ala mutation has this error. Similar is the error for changing Phe to 2-methyl Phe. Do you have any idea why this error may arise? Is it simply an undersampling issue? I am uploading the mdp file that I am using. Let me know if you think there is something wrong with the parameters I am using.
And also if there is some way of understanding which window is undersampled if undersampling is the issue.
md_11.mdp (5.0 KB)
There are issues with your protocol. For RBFE you cannot use coul-lambdas and vdw-lambdas with a single topology file, because charges from one state will get swapped into the charges of another state, but the vdw will not be appeared yet. To make it work you can either create three steps with “pmx gentop --split” or use a concerted perturbation where coulomb and lambda are changed simultaneously (for this approach you might need to use sc-function=gapsys in mdp)
Hi Vytas,
Thank you very much for pointing out the flaw in my protocol. I appreciate it very much. I do have a few questions regarding this if you don’t mind. I have seen the tutorial for the Trp cage where we do end state simulations and then do rapid non-equilibrium transitions to get the dG value. I believe this protocol should also be correct for my type of systems where I aim to look into Protein peptide interactions using mutations at certain amino acids?
But I was wondering in case I want to do RBFE using equilibrium simulations where I sample each window separately, how do I do that? I used pmx gentop --split as you suggested but after that I am a bit lost. I see it generates separate TOPOLOGY files for qon, qoff and vdw_on. But these topology files include the old topology file for the system (the mutant topology files generated using pdb2gmx) and I get another new topology file which is just the dual topology. I am confused as to how to use these files to do the simulations. I could not find any tutorial which does equilibrium simulations for RBFE.
As for your second suggestion, I believe what you mean is, if I am infact using the dual topology approach, I should change the coulomb-lambda and vdw-lambda simultaneously at the same time and sample each window. But for this method, I need to use use sc-coupl=yes and sc-function=gapsys? If I am using this approach, can I just use fep-lambdas instead of vdw and coulomb-lambdas separately?
Thanks,
Shounak
Yes, you can also use non-equilibrium calculation protocol for your case.
For the split topology files, you get three steps: 1) qoff turns off the relevant charges, 2) vdw makes the vdw transition, i.e. dummies appearing as real atoms and real atoms turned into dummies, 3) qon turns the charges of the B-state on. You should treat each step as an independent free energy calculations and afterwards add the results from the three steps to get the overall dG. To control the alchemical transition use fep-lambdas in mdp.
This is not a dual topology approach; the topology is a hybrid between single and double. For a concerted alchemical transition of coulomb and vdw you need to switch on softcore interactions for both coulomb and vdw. You can use fep-lambdas,
Vytas
Hi Vytas,
Thank you very much for replying. I would like to clear a few confusions about the split topology approach.
I understand the approach you described but where I am confused is how does gromacs understand which parameter to use for each of the three independent free energy calculation. The reason I ask is using pmx gentop --split the three topology files I get include the old topolgy files not the hybrid. In that case how would gromacs understand which charges or vdw parameters it should use for a particular lambda window or an atom. For example, if an atom changes from type CT2(State A) to CT3(State B), the itp files included in the split topology files only have the parameters when the atom type is CT2(StateA). In that case do I need to change some parameters in the mdp file?
The only way I can think of doing that is using the following parameters in mdp. Please let me know if I am on the right track
Step 1: Turning off charges ==> couple-lambda0=q couple-lambda1=none and use soft core interactions and coulomb-lambdas change from 0 to 1
Step2: vdw transition ==> couple-lambda0=vdw, couple-lambda1=vdw and vdw-lambdas change from 0 to 1.
Step3: turning on charges ==> couple-lambda0=none and couple-lambda1=q and changing couple-lambdas from 0 to 1 and using soft-core potentials.
Then add the free energy results I get for each state to get the overall free energy change. Do it for PPI and the peptide in water box and take the difference to get ddG values.
But how do I make sure that these changes occur only for the particular residue I want. From what I understand is that I cannot use couple-moltype if I am just changing one amino acid in the peptide.
I am new to these calculations and apologize if some of the questions or concerns I have are trivial.
Best,
Shounak