Skip to content

Commit 1ed0713

Browse files
jmichel80jmichel80
authored
Merge pull request #14 from fjclark/feature_update_abfe_mdr
Update ABFE setup notebook to include multiple distance restraints
2 parents 6193bdc + 8a3f4f7 commit 1ed0713

2 files changed

Lines changed: 225 additions & 11 deletions

File tree

‎04_fep/03_ABFE/01_setup_abfe.ipynb‎

Lines changed: 115 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,9 @@
6262
"metadata": {},
6363
"outputs": [],
6464
"source": [
65-
"import BioSimSpace.Sandpit.Exscientia as BSS"
65+
"import BioSimSpace.Sandpit.Exscientia as BSS\n",
66+
"import numpy as np\n",
67+
"import pandas as pd"
6668
]
6769
},
6870
{
@@ -317,14 +319,24 @@
317319
"\n",
318320
"# Search for the optimal Boresch restraints\n",
319321
"restraint = BSS.FreeEnergy.RestraintSearch.analyse(\n",
320-
" \"output/restraint_search\",\n",
321-
" system,\n",
322-
" traj,\n",
323-
" 298 * BSS.Units.Temperature.kelvin,\n",
322+
" work_dir = \"output/restraint_search\",\n",
323+
" system = system,\n",
324+
" traj = traj,\n",
325+
" temperature = 298 * BSS.Units.Temperature.kelvin,\n",
324326
" method=\"BSS\",\n",
327+
" restraint_type = \"Boresch\",\n",
325328
")"
326329
]
327330
},
331+
{
332+
"cell_type": "code",
333+
"execution_count": null,
334+
"metadata": {},
335+
"outputs": [],
336+
"source": [
337+
"restraint"
338+
]
339+
},
328340
{
329341
"cell_type": "markdown",
330342
"metadata": {},
@@ -333,7 +345,9 @@
333345
"\n",
334346
"In this algorithm candidate sets of Boresch restraints are generated and the fluctuations of the associated Boresch degrees of freedom are tracked over the unrestrained simulation. The optimum Boresch degrees of freedom are then selected based on minimum configurational volume accessible to the restrained non-interacting ligand, and the force constants are calculated to mimic the native protein-ligand interactions. Likely stability of the restraints is assessed by ensuring that the restraints impose an energy penalty of at least 10 $k_\\mathrm{B}T$ for unstable anchor point geometries. See S6 and S7 [here](https://ndownloader.figstatic.com/files/41107959) for a discussion of some of these points.\n",
335347
"\n",
336-
"It's always a good idea to check the distributions of the Boresch degrees of freedom over the course of the simulations and their values against simulation time. Check the output directory to see the plots.\n",
348+
"As well as Boresch restraints, BioSimSpace supports using [multiple distance restraints](https://pubs.acs.org/doi/full/10.1021/acs.jctc.3c00139) for ABFE calculations with both SOMD and GROMACS - simply set `restraint_type = \"multiple_distance\"` in `RestraintSearch.analyse`. Note that you'll have to specify `method = \"numerical` when runnning `restraint.getCorrection` and you'll need to run an extra calculation stage to release all but one restraint after decoupling the ligand. Examples of how to set up ABFE calculations for both Boresch and multiple distance restraints are given at the end of this notebook for both GROMACS and SOMD.\n",
349+
"\n",
350+
"It's always a good idea to check the distributions of the restrained degrees of freedom over the course of the simulations and their values against simulation time. Check the output directory to see the plots.\n",
337351
" \n",
338352
"<div class=\"alert alert-success\">\n",
339353
"<b>Exercise 1: Issues with Restraint Selection</b>\n",
@@ -712,6 +726,101 @@
712726
" gromacs_protocol, \n",
713727
" engine='gromacs', restraint=restraint, \n",
714728
" work_dir='output/gromacs')\n",
729+
"```\n",
730+
"\n",
731+
"<div class=\"alert alert-success\">\n",
732+
"<b>Example: Setting up the bound leg with GROMACS using multiple distance restraints</b>\n",
733+
"</div>\n",
734+
"\n",
735+
"The two main differences when using multiple distance restraints compared to Boresch restraints with GROMACS are a) the restraint strenght is scaled with `restraint` lambda, rather than `bonded` lambda, and b) we need to run a seperate \"release_restraint\" stage, where we calculate the free energy change for turning off all but one of the distance restraints. This is because we cannot calculate the entropic correction for releasing multiple distance restraints without simulation, but we can easily calculate this for a single distance restraint. It should also be noted that the \"release_restraint\" stage is effectively run in the opposite direction to the other states (lambda 0 - > 1 turns on the restraints, but we want the free energy of turning off the restraints), so the sign of the free energy change need to be reversed before being added to the \"full\" stage free energy change. An example of how to set up these stages is given below:\n",
736+
"\n",
737+
"```Python\n",
738+
"\n",
739+
"# First, set up the \"full\" stage, as above. During this stage,\n",
740+
"# we turn on the restraints with the ligand interacting, discharge\n",
741+
"# the ligand, and vanish the ligand.\n",
742+
"initial_lam_full=pd.Series(data={'restraint': 0.0, 'coul': 0.0, 'vdw': 0.0})\n",
743+
"full_stage_lam_vals=pd.DataFrame(\n",
744+
" data={'restraint': [0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.2, \n",
745+
" 0.35, 0.5, 0.75, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
746+
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, \n",
747+
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],\n",
748+
" 'coul': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
749+
" 0.0, 0.0, 0.16, 0.33, 0.5, 0.67, 0.83, 1.0,\n",
750+
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
751+
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],\n",
752+
" 'vdw': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
753+
" 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, \n",
754+
" 0.2, 0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75,\n",
755+
" 0.8, 0.85, 0.9, 0.95, 1.0]}\n",
756+
" )\n",
757+
"\n",
758+
"gromacs_full_protocol = BSS.Protocol.FreeEnergy(runtime=6*BSS.Units.Time.nanosecond, \n",
759+
" lam=initial_lam_full,\n",
760+
" lam_vals=full_stage_lam_vals, \n",
761+
" perturbation_type=\"full\")\n",
762+
"\n",
763+
"# Make sure that you have created a multiple distance restraint object.\n",
764+
"gromacs_full_stage = BSS.FreeEnergy.AlchemicalFreeEnergy(restraint.system, \n",
765+
" gromacs_full_protocol,\n",
766+
" engine='gromacs', restraint=multiple_distance_restraint,\n",
767+
" work_dir='output/gromacs_mdr_full')\n",
768+
"\n",
769+
"# Now, we set up the \"release_restraint\" stage. This involves releasing all but one of the\n",
770+
"# distance restraints while the ligand is decoupled.\n",
771+
"initial_lam_release_restraint=pd.Series(data={'restraint': 0.0, 'coul': 1.0, 'vdw': 1.0})\n",
772+
"\n",
773+
"# It is recommended scale lambda^5 to smoothen the removal of the harmonic bonds. This is\n",
774+
"# done automatically by SOMD, but not by GROMACS. The coul and wdw lambdas are set to 1.0\n",
775+
"# at all times to keep the ligand decoupled.\n",
776+
"release_restraint_stage_lam_vals=pd.DataFrame(\n",
777+
" data={\n",
778+
" \"restraint\": [x**5 for x in np.linspace(0.0, 1.0, 10)],\n",
779+
" \"coul\": [1.0 for _ in range(10)],\n",
780+
" \"vdw\": [1.0 for _ in range(10)],\n",
781+
" }\n",
782+
" )\n",
783+
"\n",
784+
"gromacs_release_restraint_protocol = BSS.Protocol.FreeEnergy(runtime=6*BSS.Units.Time.nanosecond, \n",
785+
" lam=initial_lam_release_restraint,\n",
786+
" lam_vals=release_restraint_stage_lam_vals,\n",
787+
" perturbation_type=\"release_restraint\")\n",
788+
"\n",
789+
"# Make sure that you have created a multiple distance restraint object.\n",
790+
"gromacs_full_stage = BSS.FreeEnergy.AlchemicalFreeEnergy(restraint.system, \n",
791+
" gromacs_release_restraint_protocol,\n",
792+
" engine='gromacs', restraint=multiple_distance_restraint,\n",
793+
" work_dir='output/gromacs_mdr_release_restraint')\n",
794+
"\n",
795+
"# The free energy change of releasing the final remaining distance restraint can then be calculated\n",
796+
"# without simulation.\n",
797+
"print(f\"Correction for releasing the final distance restraint: {multiple_distance_restraint.getCorrection(method='numerical')}\")\n",
798+
"\n",
799+
"# The free energy change of releasing the final remaining distance restraint can then be calculated\n",
800+
"# without simulation.\n",
801+
"print(f\"Correction for releasing the final distance restraint: {multiple_distance_restraint.getCorrection(method='numerical')}\")\n",
802+
"\n",
803+
"```\n",
804+
"<div class=\"alert alert-success\">\n",
805+
"<b>Example: Setting up the bound leg with SOMD using multiple distance restraints</b>\n",
806+
"</div>\n",
807+
"\n",
808+
"An ABFE calculation with SOMD using multiple distance restraints can be set up in the same way as for Boresch restraints, except we need to add an extra stage where we release all but one of the restraints, as discussed above. The additional \"restraint_restraint\" stage can be set up as shown below:\n",
809+
"\n",
810+
"```Python\n",
811+
"\n",
812+
"lam_vals_release_restraint = [0.000, 0.100, 0.200, 0.300, 0.400, 0.500, 0.600, 0.700, 0.800, 0.900, 1.000]\n",
813+
"\n",
814+
"release_restraints_protocol = BSS.Protocol.FreeEnergy(runtime=6*BSS.Units.Time.nanosecond, \n",
815+
" timestep=1*BSS.Units.Time.femtosecond, \n",
816+
" lam_vals=lam_vals_release_restraint, \n",
817+
" perturbation_type=\"release_restraint\")\n",
818+
"\n",
819+
"release_restraints_fe_calc = BSS.FreeEnergy.AlchemicalFreeEnergy(restraint.system, \n",
820+
" release_restraints_protocol, \n",
821+
" engine='somd', \n",
822+
" restraint=multiple_distance_restraint,\n",
823+
" work_dir='output/release_restraint')\n",
715824
"```"
716825
]
717826
},

0 commit comments

Comments
 (0)