search

Gallicchio-Lab / AToM-OpenMM

OpenMM-based framework for absolute and relative binding free energy calculations with the Alchemical Transfer Method
Other
118 stars 31 forks source link

Simulation crashes when reaching intermediate state #2

Closed qsabanes closed 1 year ago

qsabanes commented 2 years ago

Hi! In the past few weeks I started to use your plugin as we discussed. The tutorials run smoothly and without any issues but I started to encounter some troubles while trying to run other systems. As it happens in the other issue #1 my systems usually crash upon reaching the intermediate state. The thing is that in my case I'm using two small molecule ligands so I really don't understand what is happening at this point. For example when testing a system, 3 out of 24 ligand pairs crash while only 3 of them keep running without issue past replica 12. I tried to follow the comments from the previous issue, with specific care to the receptor residues selection, it solved some previous issed I had but no luck further. If you require further information and/or files I can provide them with no issue.

Many thanks

egallicc commented 2 years ago

Yes, absolutely, I would be happy to help debug the issues with your systems. Please share the input files here if possible, otherwise e-mail me at egallicchio@brooklyn.cuny.edu.

qsabanes commented 2 years ago

Hi Emilio, attached you can find the input files I started with (Inputs for FEP) as well as the ones I've been using in the end for CDK2 and TYK2. In the case of CDK2 the ligand pairs 1oiy - 1oi9, 20 - 1h1q and 1oiy - 1h1q seemed to run fine, the rest crash. You'll see that I put a longer displacement vector since the ligand that should be in corner was being placed close to the protein surface. In the case of TYK2 the few ligand pairs I tried crash. I also did tests with MCL-1 also without success so far.

Inputs_for_FEP.zip Ligand pairs data: https://docs.google.com/spreadsheets/d/1T1C2SsjOfhz3RJE3WLAaSkBDUoRTxh6UoLLx3AxoV6Q/edit?usp=sharing CDK2.zip TYK2.zip

Many thanks

egallicc commented 2 years ago

The CDK2 ligands do not seem to be docked into the protein receptor. image Maybe the protein got displaced during preparation? AToM requires a set of aligned ligands docked to the receptor.

Also, in setup-settings.sh there are some stray commas. For example:

"3,20,22 3,20,22", "3,20,22 3,20,22"

should be

"3,20,22 3,20,22"  "3,20,22 3,20,22"

(no commas between the two pairs of atom triplets; these are simply bash arrays whose elements are separated by white space.)

A displacement vector of (44 44 44) is probably overkill. It would result in a very large solvent box and slow MD performance. It seems to me that a shorter displacement in the xy-plane, (40 40 0) for example, would be sufficient to clear the protein with several layers of waters in between. It is not necessary to place the ligand in the corner of the box. The only requirement is that the ligand, when displaced, is separated by the protein by at least 3 hydration layers.

Let me know if this is sufficient to get you going for now. We can work on the other systems as well.

qsabanes commented 2 years ago

Hmmm that's weird. There should be one ligand inside the pocket and the other outside right? This is what I see after equilbration: image Which I thought it was the ligand placed inside the binding site but after a closer look it seems to be placed in the same region but in another pocket. But you're absolutely right that the starting structures dont have the ligand docked into the receptor. This is probably due to the protein preparation protocol, since the inputs provided were generating some errors in AnteChamber. I'll align the prepared structure to the starting one and see if it works!

About the bash thing, you're right. I wasn't doing that mistake for the other systems.

Many thanks again!

egallicc commented 2 years ago

Right, the workflow assumes that all ligands are placed in the binding pocket. During tleap, the second ligand in the pair is automatically displaced into the solvent.

qsabanes commented 2 years ago

Yeah, it works fine now! I was having some issues with the protein hydrogens, which tleap was not happy with some inputs and tampering with the protein led to this mess. Now it works fine for most proteins in my benchmark but not in CDK2 as showed in the example before. The presence of a TPO confuses the setup scripts as it ignores that residue as it is initially read as an HETATM.

Furthermore, for convergence analysis, do you use alchemlyb or another tool?

Many thanks!

egallicc commented 2 years ago

Thank you for the update. Right, the workflows we prepared for the tutorials are far from foolproof. They work for the specific systems considered there but they can easily get confused with other systems. The setup workflow (setup-atm.sh) is just a bash shell script meant to be tinkered with depending on the situation. In particular, as you found out, the logic for identifying the ligand, or ligands, is very narrow and it breaks easily.

Companies and collaborators we work with usually develop their own setup workflows. We are not experts in automated system setup workflows. We welcome any help or suggestion you can give us. Were you able to fix the CDK2 case?

Convergence analysis is possible but it is also not automated. (Sorry, I do not know alchemlyb. What is it?) For free energy analysis, we use UWHAM as implemented in R. The analysis.sh script takes optional arguments to control the minimum and maximum timestamps to analyze. By varying the trajectory window one can get an idea of the rate of convergence. Also, as you probably noticed, the analysis process produces a Rplots.pdf file with rough quality control graphs. Two of them display the probability densities of the perturbation energy as a function of lambda. If any of those displays a "gap" in the sequence of distributions--that is if neighboring probability densities do not overlap--the resulting free energy estimate cannot be trusted even if it appears converged. This often occurs when the transformation is too extensive. (We developed tricks to correct some of these occurrences and a bit of our work in this respect is published.) Conversely, if there is a continuous unbroken sequence of distributions, the free energy estimate is more trustworthy.

Oulfin commented 2 years ago

Hi,

I have a follow up question regarding this issue. I just tried to run PTP1B and the simulation crashed when trying replica 11. I have a general question, what are the the things I should pay attention to solve this issue? Is it the CM of the residues from the receptor? What other parameters I should play with in the *.cntl file that can allow me to avoid the simulation to crash?

I would like to test first to see if I can fix it and if not I will open a new issue with the files.

I run Eralpha without any issues.

Really nice tool, hype to explore what we can do with this.

Nice work Emilio and co-workers something new in the FEP topic :D

egallicc commented 2 years ago

Thank you. Yes, we also often incur similar problems when we work with a new system for the first time. In all cases I can recall, setup issues were to blame.

Assuming an RBFE project, things I recommend checking in the *_asyncre.cntl file are below. We tend to do some of these by loading the *_0.pdb file in VMD and use selections/representations. But of course, other tools work just as well.

We recently posted an improved workflow for RBFE setups. Among other things, it adds an alchemical annealing step (mdlambda) that should do a better job of preparing the system in the alchemical intermediate (lambda = 1/2). Try it out.

Oulfin commented 2 years ago

Hi,

I had some time and I could make it work using your new improved workflow.

The main changes I made was r0 and CM_TOL to 25, tol to 2.5 in restrain selected atoms.

Non-colinear selection of ref_atoms of ligands.

I run for 8 hours in one pair of ligand for PTP1B.

PTP1B 23467 23476

exp ddG = -2.07 calc ddG = -1.32 +/- 0.86

Sorry for taking over this issue, It did not crash in lambda=1/2, thanks for the information you provided me.

egallicc commented 2 years ago

Setting CM_TOL to a large value essentially removes the binding site restraint. Your experiment I think indicates that there is something wrong with the definition of the binding site. In any case, it is okay in most cases to ignore the binding site definition in RBFE calculations.

I feel we need to develop a tutorial on how to pick the 3 reference atoms ...