02_funnel_metad/02_bss-fun-metad-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "875aaaed",
   "metadata": {},
   "source": [
    "# Part 2 - Analysis\n",
    "\n",
    "After running the simulations, let's see what information we can get out of it. The primary thing we are interested is the absolute binding free energy (ABFE) of ligand binding. What funnel metadynamics (fun-metaD) fundamentally does is deposit bias along a set of reaction coordinates (aka collective variables or CVs), forcing the ligand to unbind, explore the solvent phase, and eventually allow it to rebind back. \n",
    "\n",
    "Since we can track the amount of bias deposited along these CVs, we can construct a free energy surface by inverting the accumulated bias and projecting it along the CVs. To recover the ABFE, we can look at the projection CV and find the difference in free energy between the bound phase and the unbound/solvent phase.\n",
    "\n",
    "It is important to remember that we aren't truly samply the unbound phase, instead we limit the exploration of the ligand to a small volume at the 'tip' of the funnel. Correction for the loss of rotational and translational freedom of the ligand in the unbound phase due to the funnel restraints are\n",
    "calculated as described in [Rhys *et al* 2020](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7467642/).\n",
    "\n",
    "To determine whether the simulation has converged or not, we will have a look at the time evolution of the CVs and the free energy estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdce3ab1",
   "metadata": {},
   "source": [
    "Let's import the inputs for this tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01fe50ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "from get_tutorial import download\n",
    "download(\"02\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a26ad2d",
   "metadata": {},
   "source": [
    "Now we import the relevant libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "celtic-jurisdiction",
   "metadata": {},
   "outputs": [],
   "source": [
    "import BioSimSpace as BSS\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b760ab6",
   "metadata": {},
   "source": [
    "We will need the CV definitions that we used in running the fun-metaD simulation to reconstruct the free energy surface using PLUMED. We begin with loading the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "threatened-mistress",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = BSS.IO.readMolecules(\"inputs_02/0_run/fun-metaD/openmm.*7\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95a98c38",
   "metadata": {},
   "source": [
    "Now we redefine p0 and p1 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fatal-vault",
   "metadata": {},
   "outputs": [],
   "source": [
    "p0, p1 = BSS.Metadynamics.CollectiveVariable.makeFunnel(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bafb79b",
   "metadata": {},
   "source": [
    "Write the CV definitions as you put them into the BSS script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vital-imaging",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_upper_bound = BSS.Metadynamics.Bound(3.5 * BSS.Units.Length.nanometer)\n",
    "cv = BSS.Metadynamics.CollectiveVariable.Funnel(p0, p1, upper_bound=new_upper_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6689ac72",
   "metadata": {},
   "source": [
    "Use the CV definitions to reconstruct the protocol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pediatric-jewel",
   "metadata": {},
   "outputs": [],
   "source": [
    "protocol = BSS.Protocol.Metadynamics(cv)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6b25ca4",
   "metadata": {},
   "source": [
    "Create a PLUMED process that will run the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "infinite-python",
   "metadata": {},
   "outputs": [],
   "source": [
    "plumed = BSS.Process.Plumed(work_dir=\"inputs_02/0_run/fun-metaD/\")\n",
    "config = plumed.createConfig(system, protocol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ef72fe1",
   "metadata": {},
   "source": [
    "### Part 2.1 - CVs vs time\n",
    "\n",
    "Now we can plot the CVs across time. The first column, index 0, is the projection CV and the second column, index 1, is the extent CV.\n",
    "\n",
    "We can see that the extent CV is more or less diffusive - the ligand is rapidly exploring across the entire range of allowed CV values. \n",
    "\n",
    "The projection variable is much less diffusive, where the ligand is freely exploring the unbound phase, from about 2.0 nm to 3.5 nm, but is far slower in the bound phase. This is understandable, the ligand interacts with protein and is slowed down in its exploration.\n",
    "\n",
    "One of the most important things in determining whether the fun-metaD simulation has converged is the return of the ligand from the solvent back into the binding site. Looking at the projection CV, we see the ligand at start is around 1.0-1.2 nm from the origin. It unbinds and returns to that value multiple times, especially clearly at around 35-40 ns. At least one rebinding event is needed to construct an accurate FES."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "raising-closer",
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.Notebook.plot(\n",
    "    x=plumed.getTime(time_series=True),\n",
    "    y=plumed.getCollectiveVariable(0, time_series=True),\n",
    ")\n",
    "BSS.Notebook.plot(\n",
    "    x=plumed.getTime(time_series=True),\n",
    "    y=plumed.getCollectiveVariable(1, time_series=True),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aefbc72e",
   "metadata": {},
   "source": [
    "### Part 2.2 - 2D free energy surface\n",
    "\n",
    "Now let's have a look at the free energy surfaces. Here we plot the 2D surface.\n",
    "\n",
    "The 0th index in the `free_nrg` matrix contains the projection values, 1st index - the extent values and the 2nd index has the free energy values, in kJ/mol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "three-checklist",
   "metadata": {},
   "outputs": [],
   "source": [
    "free_nrg = plumed.getFreeEnergy()\n",
    "BSS.Notebook.plotContour(\n",
    "    free_nrg[0], free_nrg[1], free_nrg[2], xlabel=\"proj (nm)\", ylabel=\"ext (nm)\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bbf98b4",
   "metadata": {},
   "source": [
    "### Part 2.3 - 1D FES\n",
    "\n",
    "We can integrate the FE values for each CV. The smoothness of the FES in the solvent phase (last 0.5 nm of projection CV) is another indicator of convergence. \n",
    "\n",
    "Note - the extent CV shouldn't have any free energy below 0.0 nm, but because the bias potential is Gaussian-shaped, the Gaussian 'tails' overspill to non-real values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "assured-length",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the projection CV vs FE\n",
    "free_nrg0 = plumed.getFreeEnergy(0, kt=2.49 * BSS.Units.Energy.kt)\n",
    "BSS.Notebook.plot(free_nrg0[0], free_nrg0[1])\n",
    "\n",
    "# the extent CV vs FE\n",
    "free_nrg1 = plumed.getFreeEnergy(1, kt=2.49 * BSS.Units.Energy.kt)\n",
    "BSS.Notebook.plot(free_nrg1[0], free_nrg1[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d323494c",
   "metadata": {},
   "source": [
    "### Part 2.4 - Absolute binding free energy\n",
    "\n",
    "As I mentioned in the beginning of this tutorial, to recover the ABFE, we need to find the difference in free energy between the bound and the unbound/solvent phases. We will use a 1D FES along the projection CV and will add the funnel correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "794a11ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "funnel_correction = cv.getCorrection().kcal_per_mol()\n",
    "\n",
    "# the last 0.5 nm of the projection CV are between indices -75 to -25\n",
    "free_nrg_floats = [i.kj_per_mol().value() for i in free_nrg0[1][-75:-25]]\n",
    "\n",
    "estimate = -(np.mean(free_nrg_floats) / 4.184 + funnel_correction.value())\n",
    "\n",
    "print(f\"the funnel correction = {funnel_correction}\")\n",
    "print(f\"the ABFE estimate = {estimate:.2f} kcal/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb0efa89",
   "metadata": {},
   "source": [
    "### Part 2.5 - $\\Delta$G vs time\n",
    "\n",
    "Finally, we need to look at the time evolution of the free energy estimate. We can get the 1D FES every 1 ns, and since we recorded the CV every 2 ps, we select a 500 step stride. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5affd31",
   "metadata": {},
   "outputs": [],
   "source": [
    "free_nrg_series = plumed.getFreeEnergy(0, stride=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a71be213",
   "metadata": {},
   "source": [
    "Now let's find the estimate at every 1 ns of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bedb5b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "funnel_correction = cv.getCorrection().kcal_per_mol()\n",
    "\n",
    "estimates_list = []\n",
    "for item in range(np.shape(free_nrg_series)[0]):\n",
    "    free_nrg_floats = [i.value() for i in free_nrg_series[item][1][-75:-25]]\n",
    "    estimate = -(np.mean(free_nrg_floats) / 4.184 + funnel_correction.value())\n",
    "    estimates_list.append(estimate)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acb17949",
   "metadata": {},
   "source": [
    "One important aspect of fun-metaD, and all metadynamics simulations, is achieving convergence. If the fun-metaD simulation has fully converged i.e. the free energy surface has been fully explored, then the ABFE estimate should plateau at a certain value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2e10354",
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.Notebook.plot(estimates_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcfdea9a",
   "metadata": {},
   "source": [
    "### Part 2.6 - Multiple repeats\n",
    "\n",
    "We can see from the plot above that the free energy estimate still fluctuates quite a lot. The estimate never completely plateaus, just oscillates within a 1 kcal/mol range.\n",
    "\n",
    "In order to get a good error estimate, you should do multiple independent simulations. I've provided the data for 5 repeats. Let's plot them all together, while also estimating the mean of the 5 repeats at each timepoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dcd8adc",
   "metadata": {},
   "outputs": [],
   "source": [
    "funnel_correction = cv.getCorrection().kcal_per_mol()\n",
    "\n",
    "plt.title(\"$\\Delta$G vs time\")\n",
    "plt.xlabel(\"time(ns)\")\n",
    "plt.ylabel(\"$\\Delta$G (kcal/mol)\")\n",
    "all_estimates = []\n",
    "\n",
    "for run_index in range(0, 5):\n",
    "    plumed = BSS.Process.Plumed(work_dir=f\"inputs_02/{run_index}_run/fun-metaD/\")\n",
    "    config = plumed.createConfig(system, protocol)\n",
    "\n",
    "    free_nrg_series = plumed.getFreeEnergy(0, stride=500)\n",
    "\n",
    "    estimates_list = []\n",
    "    for item in range(np.shape(free_nrg_series)[0]):\n",
    "        free_nrg_floats = [i.value() for i in free_nrg_series[item][1][-75:-25]]\n",
    "        estimate = -(np.mean(free_nrg_floats) / 4.184 + funnel_correction.value())\n",
    "        estimates_list.append(estimate)\n",
    "\n",
    "    plt.plot(estimates_list)\n",
    "    all_estimates.append(estimates_list)\n",
    "\n",
    "all_estimates = np.array(all_estimates)\n",
    "mean_estimates = []\n",
    "for i in range(len(all_estimates[0])):\n",
    "    mean_estimate = np.mean(all_estimates[:, i])\n",
    "    mean_estimates.append(mean_estimate)\n",
    "plt.plot(mean_estimates, linewidth=3, color=\"black\", label=\"The mean\")\n",
    "exp_affinity = -4.7\n",
    "plt.axhline(y=exp_affinity, color=\"red\", linestyle=\"--\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "038d4ed9",
   "metadata": {},
   "source": [
    "## Closing thoughts\n",
    "\n",
    "You can see that each fun-metaD simulation explores the free energy surface very differently. The ABFE estimates diverge immediately and fluctuate wildly. However, the mean of the estimates is quite stable. This indicates that, while each fun-metaD simulation is exploring a different part of the FES at any given point, the full free energy surface is explored quickly. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}