examples/ammo_basic/run_from_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "acc03043",
   "metadata": {},
   "source": [
    "# Allostery example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b19d3952",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from ammo import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4607630e",
   "metadata": {},
   "source": [
    "## System setup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64efd477",
   "metadata": {},
   "source": [
    "A common starting point in MD simulations is a PDB file. `ammo.setup.setup_system()` takes a (clean) input PDB file, and prepares it for further MD simulations. More information is in the description:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddaf8aae",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(setup.setup_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "352b9da7",
   "metadata": {},
   "source": [
    "If a `topology=` parameter is provided, the system will not be re-parameterised. Additionally, if `solvation=None`, no additional preparation will be done, and the function will go straight to minimisation. Otherwise, `solvation` can specify the size of box or shell to solvate the system in."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f8d834",
   "metadata": {},
   "source": [
    "An example structure of PTP1B is provided. We will run system setup for 7500 minimisation steps, 100 ps of heating and 250 ps of further equilibration. The file also contains an 11 residue peptide with a phosphorylated tyrosine, and so an additional parameter command is passed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93b0115e",
   "metadata": {},
   "outputs": [],
   "source": [
    "equilibrated_system = setup.setup_system('example_data/input_protein.pdb', \n",
    "                                         (7500, 100, 250),\n",
    "                                         parameters=['source leaprc.phosaa10'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2e76dac",
   "metadata": {},
   "source": [
    "## Run steered MD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eac7745",
   "metadata": {},
   "source": [
    "Once the system is prepared, the next step is to run steered MD simulations. This allows for better sampling of intermediate conformations which are unstable and therefore short-lived."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fecfac80",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(steering.run_smd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f989b5",
   "metadata": {},
   "source": [
    "The `topology` and `coordinate` parameters are self-explanatory. The `input` requires as pseudo PLUMED file, containing all the required information for steering, except the specific `ATOM` indices are replaced with [AMBER selection masks](https://amberhub.chpc.utah.edu/atom-mask-selection-syntax/), and in case of RMSD collective variables, an additional reference `FILE` parameter is added, which will be removed during PLUMED input preparation. An example pseudo PLUMED input file is given below (as well as in the specific use case examples), and more information can be found on the [PLUMED website](https://www.plumed.org/doc-v2.8/user-doc/html/_m_o_v_i_n_g_r_e_s_t_r_a_i_n_t.html). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c865758a",
   "metadata": {},
   "source": [
    "#### Single steering step"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb3ee551",
   "metadata": {},
   "source": [
    "Here is an example of multi-CV steering. The CVs will be the heavy atom RMSD of residues 178-184, the $\\chi$1 angle of Tyr152, the stacking of residues 185 and 179, and the distance between C$\\gamma$ atoms of residues 196 and 280. Note that in the masks below, the residue numbers are offset by 1. The system includes an ACE cap at the start, and the mask selection indices starting from 1. The steering will be carried out in 100 ns. The target values and forces used are based on knowledge of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28bd265a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plumed_input = ['rmsd: RMSD REFERENCE=:179-185&(!@/H) TYPE=OPTIMAL FILE=example_data/reference.pdb\\n',\n",
    "                'tyr: TORSION ATOMS=:153@N:153@CA:153@CB:153@CG\\n',\n",
    "                'pro1: DISTANCE ATOMS=:180@CE2:186@CD\\n',\n",
    "                'pro2: DISTANCE ATOMS=:180@CD1:185@CA\\n',\n",
    "                'stacking: CUSTOM ARG=pro1,pro2 FUNC=abs(x-y) PERIODIC=NO\\n',\n",
    "                'phe: DISTANCE ATOMS=:197@CG:281@CG\\n',\n",
    "                'MOVINGRESTRAINT ...\\n',\n",
    "                '  ARG=rmsd,tyr,stacking,phe\\n',\n",
    "                '  STEP0=0    AT0=initial,initial,initial,initial    KAPPA0=0.0,0.0,0.0,0.0\\n',\n",
    "                '  STEP1=2000    AT1=initial,initial,initial,initial    KAPPA1=3500.0,3500.0,3500.0,3500.0\\n',\n",
    "                '  STEP2=75000000    AT2=0.0,1.047,0.0,0.45    KAPPA2=3500.0,3500.0,3500.0,3500.0\\n',\n",
    "                '  STEP3=76000000    AT3=0.0,1.047,0.0,0.45    KAPPA3=0.0,0.0,0.0,0.0\\n',\n",
    "                '... MOVINGRESTRAINT\\n',\n",
    "                'PRINT STRIDE=2500 ARG=* FILE=steering.dat']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f0f2df5",
   "metadata": {},
   "source": [
    "The `\"initial\"` values for the CVs at steps 0 and 1 will be computed using PLUMED and filled in during final file setup. This, together with the use of AMBER atom masks, allows for easier steering preparation while still using the whole range of CVs in PLUMED."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b59b602",
   "metadata": {},
   "outputs": [],
   "source": [
    "steering = steering.run_smd('system.prm7', 'system_equilibrated.rst7', input=plumed_input)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3dbb0af",
   "metadata": {},
   "source": [
    "After the steering process is finished, the output files are copied over from the working directory, and additionally a dry copy of the trajectory (no waters or ions) is saved (in case of further analysis required)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "140069da",
   "metadata": {},
   "source": [
    "#### Custom steered MD protocols"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c279e377",
   "metadata": {},
   "source": [
    "`run_smd()` uses BioSimSpace to manage steered MD simulations. For more customization, sMD can be run [directly with BSS](https://github.com/michellab/BioSimSpaceTutorials/tree/main/03_steered_md)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aec72db4",
   "metadata": {},
   "source": [
    "## Analysing steered MD data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9242acd9",
   "metadata": {},
   "source": [
    "Once a steered MD trajectory is produced, it has to be checked to ensure steering has been successful, and snapshots need to be saved for seeded MD simulations. This simple trajectory analysis can be done however the user choses. There is an example notebook in `$ALLOSTERYHOME/data/sMD_analysis.ipynb` which can be a good starting point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbf21c6d",
   "metadata": {},
   "source": [
    "## Seeded MD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53686ecf",
   "metadata": {},
   "source": [
    "With snapshot saved from the sMD trajectory, they can be used as \"seeds\" to run equilibrium MD simulations. Since they are indeed just equilibrium MD simulations, the `ammo.equilibrium.run_eq_md()` function is most appropriate in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ace6cf78",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(equilibrium.run_eq_md)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "986ba041",
   "metadata": {},
   "source": [
    "The `coordinates=` parameter would simply be the saved snapshot coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b622250",
   "metadata": {},
   "source": [
    "Please note that this step is highly recommended to be done on a computing cluster with multiple GPU access and running multiple seeded MD simulations in parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08be1d10",
   "metadata": {},
   "source": [
    "## Trajectory featurization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4582e1d3",
   "metadata": {},
   "source": [
    "Once seeded MD simulations are finished, they can be used to build a Markov State Model. However, that requires dimensionality reduction, which starts by reducing trajectory data from all atom coordinates to select features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c25347e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(analysis.featurize)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59363560",
   "metadata": {},
   "source": [
    "`featurize()` computes a distance, dihedral or RMSD values for the trajectory specified, using an AMBER selection mask ([documentation](https://amberhub.chpc.utah.edu/atom-mask-selection-syntax/)). If RMSD is being calculated, `reference` has to be provided as well, and `shared` is the selection mask for atoms used for alignment. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a5c48d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "featurized_trajectory = analysis.featurize('example_data/steering.nc', 'example_data/system.prm7', \n",
    "                                           'distance', ':197@CG :281@CG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68dfb6c8",
   "metadata": {},
   "source": [
    "will calculate the distance between C$\\gamma$ atoms for residues 196 and 280. Since the trajectory path was given, topology is also provided. The results are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3acc48aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(10,5))\n",
    "\n",
    "ax.plot(featurized_trajectory, alpha=0.7)\n",
    "ax.set_xlabel('frames')\n",
    "ax.set_ylabel('distance/$\\AA$')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e3df3e",
   "metadata": {},
   "source": [
    "## MSM building"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c10721ca",
   "metadata": {},
   "source": [
    "An example of Markov State Modelling is available in a separate [notebook](msm.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32e4cc4b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "2e311702793db1e2cbdbca1b4af01943df9920d13b5c28949f5d0bcf0a016590"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}