Skip to content

Commit

Permalink
Update ABFE to include multiple distance restraints
Browse files Browse the repository at this point in the history 
This requires the changes to BioSimSpace implemented
here: OpenBioSim/biosimspace#178.
  • Loading branch information
@fjclark
fjclark committed Nov 18, 2023
1 parent 0f9196a commit 8a3f4f7
Show file tree
Hide file tree
Showing 2 changed files with 225 additions and 11 deletions.
121 changes: 115 additions & 6 deletions 04_fep/03_ABFE/01_setup_abfe.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,9 @@
"metadata": {},
"outputs": [],
"source": [
"import BioSimSpace.Sandpit.Exscientia as BSS"
"import BioSimSpace.Sandpit.Exscientia as BSS\n",
"import numpy as np\n",
"import pandas as pd"
]
},
{
Expand Down Expand Up @@ -317,14 +319,24 @@
"\n",
"# Search for the optimal Boresch restraints\n",
"restraint = BSS.FreeEnergy.RestraintSearch.analyse(\n",
" \"output/restraint_search\",\n",
" system,\n",
" traj,\n",
" 298 * BSS.Units.Temperature.kelvin,\n",
" work_dir = \"output/restraint_search\",\n",
" system = system,\n",
" traj = traj,\n",
" temperature = 298 * BSS.Units.Temperature.kelvin,\n",
" method=\"BSS\",\n",
" restraint_type = \"Boresch\",\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"restraint"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -333,7 +345,9 @@
"\n",
"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",
"\n",
"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",
"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",
"\n",
"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",
" \n",
"<div class=\"alert alert-success\">\n",
"<b>Exercise 1: Issues with Restraint Selection</b>\n",
Expand Down Expand Up @@ -712,6 +726,101 @@
" gromacs_protocol, \n",
" engine='gromacs', restraint=restraint, \n",
" work_dir='output/gromacs')\n",
"```\n",
"\n",
"<div class=\"alert alert-success\">\n",
"<b>Example: Setting up the bound leg with GROMACS using multiple distance restraints</b>\n",
"</div>\n",
"\n",
"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",
"\n",
"```Python\n",
"\n",
"# First, set up the \"full\" stage, as above. During this stage,\n",
"# we turn on the restraints with the ligand interacting, discharge\n",
"# the ligand, and vanish the ligand.\n",
"initial_lam_full=pd.Series(data={'restraint': 0.0, 'coul': 0.0, 'vdw': 0.0})\n",
"full_stage_lam_vals=pd.DataFrame(\n",
" data={'restraint': [0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.2, \n",
" 0.35, 0.5, 0.75, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, \n",
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],\n",
" 'coul': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
" 0.0, 0.0, 0.16, 0.33, 0.5, 0.67, 0.83, 1.0,\n",
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
" 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],\n",
" 'vdw': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
" 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, \n",
" 0.2, 0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75,\n",
" 0.8, 0.85, 0.9, 0.95, 1.0]}\n",
" )\n",
"\n",
"gromacs_full_protocol = BSS.Protocol.FreeEnergy(runtime=6*BSS.Units.Time.nanosecond, \n",
" lam=initial_lam_full,\n",
" lam_vals=full_stage_lam_vals, \n",
" perturbation_type=\"full\")\n",
"\n",
"# Make sure that you have created a multiple distance restraint object.\n",
"gromacs_full_stage = BSS.FreeEnergy.AlchemicalFreeEnergy(restraint.system, \n",
" gromacs_full_protocol,\n",
" engine='gromacs', restraint=multiple_distance_restraint,\n",
" work_dir='output/gromacs_mdr_full')\n",
"\n",
"# Now, we set up the \"release_restraint\" stage. This involves releasing all but one of the\n",
"# distance restraints while the ligand is decoupled.\n",
"initial_lam_release_restraint=pd.Series(data={'restraint': 0.0, 'coul': 1.0, 'vdw': 1.0})\n",
"\n",
"# It is recommended scale lambda^5 to smoothen the removal of the harmonic bonds. This is\n",
"# done automatically by SOMD, but not by GROMACS. The coul and wdw lambdas are set to 1.0\n",
"# at all times to keep the ligand decoupled.\n",
"release_restraint_stage_lam_vals=pd.DataFrame(\n",
" data={\n",
" \"restraint\": [x**5 for x in np.linspace(0.0, 1.0, 10)],\n",
" \"coul\": [1.0 for _ in range(10)],\n",
" \"vdw\": [1.0 for _ in range(10)],\n",
" }\n",
" )\n",
"\n",
"gromacs_release_restraint_protocol = BSS.Protocol.FreeEnergy(runtime=6*BSS.Units.Time.nanosecond, \n",
" lam=initial_lam_release_restraint,\n",
" lam_vals=release_restraint_stage_lam_vals,\n",
" perturbation_type=\"release_restraint\")\n",
"\n",
"# Make sure that you have created a multiple distance restraint object.\n",
"gromacs_full_stage = BSS.FreeEnergy.AlchemicalFreeEnergy(restraint.system, \n",
" gromacs_release_restraint_protocol,\n",
" engine='gromacs', restraint=multiple_distance_restraint,\n",
" work_dir='output/gromacs_mdr_release_restraint')\n",
"\n",
"# The free energy change of releasing the final remaining distance restraint can then be calculated\n",
"# without simulation.\n",
"print(f\"Correction for releasing the final distance restraint: {multiple_distance_restraint.getCorrection(method='numerical')}\")\n",
"\n",
"# The free energy change of releasing the final remaining distance restraint can then be calculated\n",
"# without simulation.\n",
"print(f\"Correction for releasing the final distance restraint: {multiple_distance_restraint.getCorrection(method='numerical')}\")\n",
"\n",
"```\n",
"<div class=\"alert alert-success\">\n",
"<b>Example: Setting up the bound leg with SOMD using multiple distance restraints</b>\n",
"</div>\n",
"\n",
"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",
"\n",
"```Python\n",
"\n",
"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",
"\n",
"release_restraints_protocol = BSS.Protocol.FreeEnergy(runtime=6*BSS.Units.Time.nanosecond, \n",
" timestep=1*BSS.Units.Time.femtosecond, \n",
" lam_vals=lam_vals_release_restraint, \n",
" perturbation_type=\"release_restraint\")\n",
"\n",
"release_restraints_fe_calc = BSS.FreeEnergy.AlchemicalFreeEnergy(restraint.system, \n",
" release_restraints_protocol, \n",
" engine='somd', \n",
" restraint=multiple_distance_restraint,\n",
" work_dir='output/release_restraint')\n",
"```"
]
},
Expand Down
Loading

0 comments on commit 8a3f4f7

Please sign in to comment.