Repository navigation
Septop tutorial inquiry #315
Description
Hello,
I was working through the SepTop tutorial (openmm_septop/septop_tutorial.ipynb) with the example Tyk2 structures (ejm 31 --> 42), followed by simulation of the ERa system (with ligands 3b and 2e and aux paper setup). I ran triplicates of each system, but free energy results I received for were quite noisy for both:
Tyk2:
Run 1: 0.769 kcal/mol
Run 2: 1.540 kcal/mol
Run 3: 0.3927 kcal/mol
ERa:
Run 1 = -3.37 kcal/mol
Run 2 = -0.72 kcal/mol
Run 3 = -4.29 kcal/mol
I ran through the tutorial as is, which I think is meant to replicate the publication (https://pubs.acs.org/doi/10.1021/acs.jctc.3c00282) parameters. I was hoping you might be able to provide some insight as to whether there is something I am missing that might have caused the precision noise, as I was expecting higher precision for this specific transformation (referencing the paper). I checked the overlap matrices and they look fine, but the simulations were clearly not converged.
Happy to provide further details if needed. Thank you!
Activity
hannahbaumann commentedon Jun 23, 2026
Hi @priscilladuggan4 ,
thank you for raising this!
It's always tricky to figure out why some legs fail to converge as tightly as expected, but here are some thoughts I have:
Regarding replicating the publication: a key difference is that in the paper, restraint atom selection and the equilibration used to inform it were performed once per system and reused across all three repeats; in openfe, each
protocol_repeatsrun reselects restraint atoms independently (documented in the SepTop protocol docs), likely contributing to higher inter-repeat spread than reported. Also, the implementation from the paper combines HREX (Hamiltonian replica exchange) with LJ-epsilon softening (REST2-like enhanced sampling), while openfe only does HREX. In the paper the plain-HREX variant sometimes had higher stdev than the REST2-like one, so this could be another explanation for the slightly higher inter-repeat variability here.Different things you could check:
openfe gather-septop --report rawto see how the individual legs differ.Let me know what you find or if you have any more questions!
priscilladuggan4 commentedon Jun 29, 2026
Thank you so much for getting back to me and for all the helpful information! I've looked into each point you raised...
Inter-repeat noise analysis:
I've examined the results and it seems like the inter-repeat differences are mainly coming from the complex leg:
Structural changes:
I did structural analysis for the ERa runs and it seems like there is a lot of movement for Ligand B (2e) in the complex leg during Run 2:
As well as Ligand A (3b) in the complex leg during Run 3:
On Boresch restraints:
I went into the .json files and saw that there are different restraints for each replicate. The guest atoms are the same for each replicate but the host atoms are different each time. I am wondering if there is a way to rerun a simulation and force the same restraint selection across multiple runs? That way I would be able to rerun multiple replicates with the same restraints to see if the noise was coming from that or not?
Reverse transformations:
I haven't had a chance to run the reverse transformations yet, but I plan to do this in the near future.
Interestingly, I reran the ERa simulation in triplicate with a simulation time of 50ns (instead of the 10ns the tutorial settings defaults to) and results with much higher precision (note experimental value is –2.44 kcal/mol):
This would suggest that if the difference in Boresch restraints was adding to the noise in the 10ns runs, that effect becomes negligible when the production time is increased by a factor of five.
Thank you again,
Priscilla
hannahbaumann commentedon Jul 3, 2026
@priscilladuggan4 ,
thank you for looking into this so thoroughly, this is great to see that the precision is much higher when running the simulations for longer!
Regarding the structural changes: What I see in these plots makes sense, it's always the non-interacting lambda windows (state 18 for ligand A and state 0 for ligand B) that have the large RMSD/COM displacement. This means that the Boresch restraints are not sufficiently keeping the non-interacting ligand close to the bound conformation in the interacting state.
Regarding Boresch restraints across repeats: you're right, these are being re-calculated based on the repeated equilibration trajectories separately for each repeat. We're working on enabling user-supplied Boresch restraint atoms (first for ABFEs here: OpenFreeEnergy/openfe#2019) which would enable what you're looking for. But as of right now this is not easily possible. If it's important enough to be worth the effort, one option is to modify this part of the code (https://github.com/OpenFreeEnergy/openfe/blob/97bef39d3ccfabb259636728b7782490640f42ed/src/openfe/protocols/openmm_septop/equil_septop_method.py#L550-L564) to move the setup unit outside the repeat loop, such that the Boresch selection is only run once. Then you'd have to set the
protocol_repeats = 3in the settings and run all three repeats sequentially using a single quickrun command, different from how we're doing it in the example notebook SepTop tutorial where we're setting theprotocol_repeats = 1and submitting each transformation json 3x. Happy to help debug if you give it a go and hit issues.In the meantime, I had a quick go at a Boresch restraint analysis notebook. This is still a rough draft, but maybe it could help you nail down the sampling problem you observed in those runs: OpenFreeEnergy/openfe_analysis#120
Please let me know if you have any questions!
priscilladuggan4 commentedon Jul 31, 2026
Hi Hannah,
I apologize for the delayed response! Thank you so much for the information regarding Boresch restraints across repeats - this has been very helpful. My colleague helped me with our subsequent simulations by implementing your suggestions into a custom openfe software environment, allowing us to force the same restraints across 3 replicates. I then ran three sets of separate calculations for 3b --> 2e with
protocol repeats = 3. Please see the results below:These were interesting results because, as shown in the table, there is significant variance in inter-replicate precision (ranging from 0.199 to 0.949 kcal/mol), suggesting that restraint identity may strongly affect repeat-to-repeat precision. I am working on some analysis now but hopefully I get some time to look more into this in the future! Thank you for your help with all of this.
All the best,
Priscilla