rna_md_sim_setup.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BioSimSpace\n",
    "Author: Rory Portman <br>\n",
    "Email:&nbsp;&nbsp; portmanrk@gmail.com\n",
    "\n",
    "_This notebook forms part of my BSc research project on Binding Free Energy Calculation protocols for Small molecule-RNA Complexes._\n",
    "\n",
    "# Molecular Dynamics Simulations Setup and Analysis for Small Molecule-RNA Complexes\n",
    "\n",
    "## Introduction\n",
    "\n",
    "__This notebook covers the necessary steps to setup, run and analyse Molecular Dynamics Simulations for Small Molecule-RNA complexes. Included in this notebook:__\n",
    "\n",
    "* Loading, parameterising, solvating and merging the Ligands, the RNA and the Waters.\n",
    "* Running Equilibration and Minimisation for the Small Molecule:RNA system.\n",
    "* Writing scripts to initiate 100 ns molecular dynamics simulations using various simulation engines.\n",
    "* Analysing MD simulations and producing plots for reports.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__For this project I will be using small molecule ligands and an RNA PreQ1 riboswitch reported by Connely et al.<sup>1</sup> and performing MD Simulations using AMBER, <sup>2</sup> GROMACS, <sup>3</sup> and OpenMM.<sup>4</sup>__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `[1]` _Connelly, C.M., Numata, T., Boer, R.E. et al. Synthetic ligands for PreQ1 riboswitches provide structural and mechanistic insights into targeting RNA tertiary structure. Nat Commun 10, 1501 (2019). https://doi.org/10.1038/s41467-019-09493-3_\n",
    "\n",
    "- `[2]` _Case DA, Cheatham TE 3rd, Darden T, Gohlke H, Luo R, Merz KM Jr, Onufriev A, Simmerling C, Wang B, Woods RJ. The Amber biomolecular simulation programs. J Comput Chem. 2005 Dec;26(16):1668-88. https://doi.org/10.1002/jcc.20290 PMID: 16200636; PMCID: PMC1989667._\n",
    "\n",
    "- `[3]` _M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX,2015, 1-2, 19–25.https://doi.org/10.1016/j.softx.2015.06.001_\n",
    "\n",
    "- `[4]` _Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, R. P. Wiewiora, B. R. Brooks and V. S. Pande,PLOS Computational Biology, 2017, 13, e1005659. https://doi.org/10.1371/journal.pcbi.1005659_\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1: Imports and Inputs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, all necessary directories were imported:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import BioSimSpace as BSS\n",
    "from BioSimSpace import _Exceptions\n",
    "import os\n",
    "import sys\n",
    "import csv\n",
    "import struct\n",
    "import glob\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import MDAnalysis as mda\n",
    "from MDAnalysis.analysis import rms, align\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the following small molecule ligand (see below) reported in Connely et al. was obtained from the Protein Data Bank (https://www.ebi.ac.uk/pdbe/):\n",
    "\n",
    "- __Ligand 1:__ 2-[(dibenzo[b,d]furan-2-yl)oxy]ethan-1-amine [PDB ID: 6E1S]\n",
    "\n",
    "This ligand, along with the PreQ1 riboswitch [PDB ID: 3Q50] and waters, is loaded below. \n",
    "\n",
    "__Note:__\n",
    "\n",
    "The 6E1S structure contained the PreQ1 riboswitch backbone with two omitted residues (13 and 14). To ensure a representative simulation of the ligand bound to the PreQ1 riboswitch, a structure with the cognate metabolite bound to the PreQ1 riboswitch (PDB ID: 3Q50), which included all residues, was also obtained from Connelly et al. This complete structure was then aligned with the binding pose of Ligand 1 (6E1S) using PyMOL, resulting in an RMSD of 0.4 Å. This signified that there was considerable overlap between the metabolite and the ligand, with residue 13 positioned slightly closer to the ligand than to the metabolite. The metabolite and the incomplete backbone were subsequently removed in PyMOL leaving just Ligand 1 (PDB ID: 6E1S) and the complete PreQ1 riboswitch backbone (PDB ID: 3Q50). $SO_4^{2-}$, which was present in the structure due to the crystallization and purification processes, was deleted using PyMOL as it is not typically observed under physiological conditions. The PreQ1 backbone and Ligand 1 were then protonated using PDBFixer and Open Babel (version 3.0.0) respectively. Since no parameters for the PreQ1 riboswitch were available in the AMBER GAFF2 force field, the PreQ1 backbone and the crystal structure waters were parameterized following the methodology of Hashem et al.<sup>5</sup> using tleap prior to loading.\n",
    "\n",
    "- `[5]` _Y. Hashem and P. Auffinger, Methods, 2009, 47, 187–197. https://doi.org/10.1016/j.ymeth.2008.09.020_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LOAD LIGAND, FREE_LIGAND, RNA AND WATER\n",
    "\n",
    "ligand = BSS.IO.readMolecules(\"input_rory/lig_1_protonated.pdb\")[0]\n",
    "free_ligand = ligand\n",
    "water=BSS.IO.readMolecules(\"input_rory/waters.*\")\n",
    "rna_bases = BSS.IO.readMolecules(\"input_rory/input_rory_3q50/rna.*\")\n",
    "\n",
    "print (\"ligand =\", ligand)\n",
    "print(\"free_ligand =\",free_ligand)\n",
    "print(\"water =\",water)\n",
    "print (\"rna_bases =\", rna_bases)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2: Ligand 1-PreQ1 riboswitch System Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ligand 1 was parametrised using the GAFF2 AMBER forcefield. As discussed above, the RNA and Waters were parameterised prior to loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PARAMATERISE LIGAND, FREE_LIGAND AND WATER\n",
    "\n",
    "gaff_free_ligand = BSS.Parameters.gaff2(free_ligand).getMolecule()\n",
    "gaff_ligand = BSS.Parameters.gaff2(ligand).getMolecule()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameterised inputs of Ligand 1, PreQ1 and the Waters were then combined into one system called `molecules`. It is sometimes useful to view a loaded system using `BSS.Notebook.View()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#COMBINE RNA, WATER AND LIGAND\n",
    "\n",
    "molecules= rna_bases + gaff_ligand + water\n",
    "view=BSS.Notebook.View(molecules)\n",
    "view.system()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the creation of the Ligand 1-PreQ1 it was solvated using `BSS.Solvent`. In this case, the `tip3p` solvating model was used and an automatically calculated solvent box sizewas determined using `.getAxisAlignedBoundingBox()`. It is important to specify the counter ion concentration. In this case 0.15 M of NaCl was specified to match the physiological concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DETERMINE BOX SIZE AND SOLVATE COMBINED SYSTEM (CHECK COUNTER IONS)\n",
    "\n",
    "box_min, box_max = rna.getAxisAlignedBoundingBox()\n",
    "box_size = [y - x for x, y in zip(box_min, box_max)]\n",
    "padding = 15 * BSS.Units.Length.angstrom\n",
    "box_length = max(box_size) + 2 * padding\n",
    "\n",
    "molecules_sol = BSS.Solvent.tip3p(molecule= molecules, box=3 * [box_length], ion_conc = 0.15)\n",
    "search = molecules_sol.search(\"not mols with atomidx 2\")\n",
    "for ion in search.atoms():\n",
    "    print(f\"element = {ion.element()}, charge = {ion.charge()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This process was also carried out for the unbound (free) Ligand for subsequent use in Relative Binding Free Energy calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DETERMINE BOX SIZE AND SOLVATE FREE LIGAND\n",
    "\n",
    "free_box_min, free_box_max = free_ligand.getAxisAlignedBoundingBox()\n",
    "free_box_size = [y - x for x, y in zip(free_box_min, free_box_max)]\n",
    "free_padding = 15 * BSS.Units.Length.angstrom\n",
    "free_box_length = max(free_box_size) + 2 * free_padding\n",
    "print(free_box_length)\n",
    "\n",
    "free_molecules_sol = BSS.Solvent.tip3p(molecule= gaff_free_ligand, box=3 * [free_box_length], ion_conc = 0.15)\n",
    "search = free_molecules_sol.search(\"not mols with atomidx 2\")\n",
    "for ion in search.atoms():\n",
    "    print(f\"element = {ion.element()}, charge = {ion.charge()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3: Ligand 1-PreQ1 riboswitch System Energy Minimisation and Equilibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ensure the stability of the Ligand 1-PreQ1 riboswitch system prior to molecular dynamics simulations, energy equilibration and minimisation steps were conducted using `PMEMD.CUDA`. The following approach for minimisation and equilibration is adapted from: https://github.com/michellab/BioSimSpaceTutorials/blob/69e2a9f97b2d2f5ae21c5c5b45a2992b2879621c/03_fep/execution_model/scripts/BSSligprep.py#L225 \n",
    "\n",
    "The Ligand 1-PreQ1 riboswitch and free ligand systems were each minimised in 250 steps using `PMEMD.CUDA` at 1 atm and gradually heated from 0 K to 300 K. The systems were then equilibrated in the NVT ensemble for 5 ps with all non-solvent heavy atoms restrained, followed by 50 ps with restraints, and 50 ps without restraints. Lastly, equilibration in the NPT ensemble was performed for 200 ps with non-solvent heavy atom restraints and for another 200 ps without restraints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, the settings, temporary directories and inputs are defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Settings.\n",
    "minim_steps = 250\n",
    "runtime_short_nvt = 5 # ps\n",
    "runtime_nvt = 50 # ps RESET to 50 after TESTING ! \n",
    "runtime_npt = 200 # ps RESET to 200 after TESTING !\n",
    "\n",
    "### preamble. tmp_dir should at some point be derived using os.environ.get(\"\")\n",
    "\n",
    "tmp_dir = \"temp\"\n",
    "\n",
    "amber_home = \"/usr/local/amber22\"\n",
    "pmemd_path = amber_home + \"/bin/pmemd.cuda\" \n",
    "\n",
    "#################\n",
    "lig_p_solvated = free_molecules_sol\n",
    "system_solvated = molecules_sol\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/free_lig_s\", lig_p_solvated, [\"PRM7\", \"RST7\"])\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/mol_sys_s\", system_solvated, [\"PRM7\", \"RST7\"])\n",
    "### Minimise and equilibrate our systems to make the ligand relax inside the pocket.\n",
    "\n",
    "lig_p_solvated = BSS.IO.readMolecules([f\"{tmp_dir}/free_lig_s.prm7\", f\"{tmp_dir}/free_lig_s.rst7\"])\n",
    "system_solvated = BSS.IO.readMolecules([f\"{tmp_dir}/mol_sys_s.prm7\", f\"{tmp_dir}/mol_sys_s.rst7\"])\n",
    "\n",
    "print (\"lig_p_solvated is =\", lig_p_solvated)\n",
    "print (\"original is =\", free_molecules_sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy equilibration and minimisation is then run for the unbound (free) ligand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def runProcess(system, protocol, pmemd=False):\n",
    "        \"\"\"\n",
    "        Given a solvated system (BSS object) and BSS protocol, run a process workflow with either \n",
    "        Sander (CPU) or pmemd.cuda (GPU). NPT is typically done with GPU to save computing time.\n",
    "        Returns the processed system.\n",
    "        \"\"\"\n",
    "        # Create the process passing a working directory.\n",
    "        if not pmemd:\n",
    "            process = BSS.Process.Amber(system, protocol)\n",
    "        elif pmemd:\n",
    "            process = BSS.Process.Amber(system, protocol, exe=pmemd_path)\n",
    "\n",
    "        # Start the process.\n",
    "        process.start()\n",
    "\n",
    "        # Wait for the process to exit.\n",
    "        process.wait()\n",
    "\n",
    "        # Check for errors.\n",
    "        if process.isError():\n",
    "            print(process.stdout())\n",
    "            print(process.stderr())\n",
    "            raise _Exceptions.ThirdPartyError(\"The process exited with an error!\")\n",
    "\n",
    "        # If it worked, try to get the system. No need to block, since it's already finished.\n",
    "        system = process.getSystem()\n",
    "\n",
    "        return system\n",
    "\n",
    "\n",
    "############# first minimise/equilibrate the solvated ligand.\n",
    "print(\"\\n#### Working on solvated ligand.\")\n",
    "\n",
    "print(f\"Minimising in {minim_steps} steps..\")\n",
    "protocol = BSS.Protocol.Minimisation(steps=minim_steps)\n",
    "minimised = runProcess(lig_p_solvated, protocol)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/min_free_lig_s\", lig_p_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "\n",
    "print(f\"PMEMD NVT equilibration for {runtime_short_nvt} ps while restraining all non-solvent atoms..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_short_nvt*BSS.Units.Time.picosecond, \n",
    "                                temperature_start=0*BSS.Units.Temperature.kelvin, \n",
    "                                temperature_end=300*BSS.Units.Temperature.kelvin,\n",
    "                                restraint=\"all\"\n",
    "                                )\n",
    "equil1 = runProcess(minimised, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq1_free_lig_s\", lig_p_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NVT equilibration for {runtime_nvt} ps without restraints..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_nvt*BSS.Units.Time.picosecond, \n",
    "                                temperature=300*BSS.Units.Temperature.kelvin,\n",
    "                                )\n",
    "\n",
    "equil2 = runProcess(equil1, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq2_free_lig_s\", lig_p_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NPT equilibration for {runtime_npt} ps while restraining non-solvent heavy atoms..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_npt*BSS.Units.Time.picosecond, \n",
    "                                pressure=1*BSS.Units.Pressure.atm,\n",
    "                                temperature=300*BSS.Units.Temperature.kelvin,\n",
    "                                restraint=\"heavy\",\n",
    "                                )\n",
    "equil3 = runProcess(equil2, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq3_free_lig_s\", lig_p_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NPT equilibration for {runtime_npt} ps without restraints..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_npt*BSS.Units.Time.picosecond, \n",
    "                                pressure=1*BSS.Units.Pressure.atm,\n",
    "                                temperature=300*BSS.Units.Temperature.kelvin,\n",
    "                                )\n",
    "lig_equil_fin = runProcess(equil3, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq4_free_lig_s\", lig_p_solvated, [\"PRM7\", \"RST7\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy equilibration and minimisation is then run for the Ligand 1-PreQ1 system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############# repeat for ligand + protein system. Include extra restrain=\"backbone\" step.\n",
    "print(\"\\n#### Working on solvated ligand+protein.\")\n",
    "print(f\"Minimising in {minim_steps} steps..\")\n",
    "protocol = BSS.Protocol.Minimisation(steps=minim_steps)\n",
    "minimised = runProcess(system_solvated, protocol)\n",
    "\n",
    "\n",
    "print(f\"PMEMD NVT equilibration for {runtime_short_nvt} ps while restraining all non-solvent atoms..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_short_nvt*BSS.Units.Time.picosecond, \n",
    "                                temperature_start=0*BSS.Units.Temperature.kelvin, \n",
    "                                temperature_end=300*BSS.Units.Temperature.kelvin,\n",
    "                                restraint=\"all\"\n",
    "                                )\n",
    "equil1 = runProcess(minimised, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq1_mol_sys_s\", system_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NVT equilibration for {runtime_nvt} ps while restraining all backbone atoms..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_nvt*BSS.Units.Time.picosecond, \n",
    "                                temperature=300*BSS.Units.Temperature.kelvin, \n",
    "                                #restraint=\"backbone\"\n",
    "                                )\n",
    "equil2 = runProcess(equil1, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq2_mol_sys_s\", system_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NVT equilibration for {runtime_nvt} ps without restraints..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_nvt*BSS.Units.Time.picosecond, \n",
    "                                temperature_end=300*BSS.Units.Temperature.kelvin,\n",
    "                                )\n",
    "\n",
    "equil3 = runProcess(equil2, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq3_mol_sys_s\", system_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NPT equilibration for {runtime_npt} ps while restraining non-solvent heavy atoms..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_npt*BSS.Units.Time.picosecond, \n",
    "                                pressure=1*BSS.Units.Pressure.atm,\n",
    "                                temperature=300*BSS.Units.Temperature.kelvin,\n",
    "                                restraint=\"heavy\",\n",
    "                                )\n",
    "equil4 = runProcess(equil3, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq4_mol_sys_s\", system_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(f\"PMEMD NPT equilibration for {runtime_npt} ps without restraints..\")\n",
    "protocol = BSS.Protocol.Equilibration(\n",
    "                                runtime=runtime_npt*BSS.Units.Time.picosecond, \n",
    "                                pressure=1*BSS.Units.Pressure.atm,\n",
    "                                temperature=300*BSS.Units.Temperature.kelvin,\n",
    "                                )\n",
    "sys_equil_fin = runProcess(equil4, protocol, pmemd=True)\n",
    "\n",
    "BSS.IO.saveMolecules(f\"{tmp_dir}/eq5_mol_sys_s\", system_solvated, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "# finally, save last snapshot of both equilibrated objects.\n",
    "os.system(\"mkdir -p prep/ligands\")\n",
    "os.system(\"mkdir -p prep/protein\")\n",
    "\n",
    "print(\"Saving solvated/equilibrated systems.\")\n",
    "print(\"\\n Ligand:\")\n",
    "print(lig_equil_fin)\n",
    "BSS.IO.saveMolecules(f\"prep_bases/ligands/free_lig_equil_solv\", lig_equil_fin, [\"PRM7\", \"RST7\"])\n",
    "\n",
    "print(\"\\n Ligand + protein:\")\n",
    "print(sys_equil_fin)\n",
    "BSS.IO.saveMolecules(f\"prep_bases/protein/mol_sys_equil_solv\", sys_equil_fin, [\"PRM7\", \"RST7\"])\n",
    "print(\"First 20 molecules in ligand + protein system:\")\n",
    "for mol in sys_equil_fin.getMolecules()[:20]:\n",
    "    print(mol)\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4: Running Molecular Dynamics Simulations using AMBER, GROMACS and OpenMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Molecular dynamics simulations of the Ligand 1 in the binding site of the PreQ1 riboswitch were initiated using the BSS.Protocol and BSS.Process functions. 100 ns MD simulations of the system  were run using AMBER, GROMACS and OpenMM. The MD simulations were initiated using the following [Python Script](https://github.com/michellab/RNA_AFE.) and queued using the following [Slurm Script](https://github.com/michellab/RNA_AFE.). An example python script is shown below:\n",
    "\n",
    "```python\n",
    "import BioSimSpace as BSS\n",
    "import os\n",
    "\n",
    "print(\"Loading system...\")\n",
    "system= BSS.IO.readMolecules(\"../prep/protein/mol_sys_equil_solv.*\")\n",
    "protocol=BSS.Protocol.Production(runtime= 100 * BSS.Units.Time.nanosecond, report_interval=10_000)\n",
    "process = BSS.Process.Amber(system, protocol,exe = \"/usr/local/amber22/bin/pmemd.cuda\", work_dir=\"output_amber\")\n",
    "process.start()\n",
    "process.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5: Molecular Dynamics Simulation Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The trajectories from the MD simulations were initially visualised using PyMOL and subsequently analysed with MDAnalysis. Firstly, the root mean square deviation (RMSD) of Ligand 1 throughout the simulations was calculated relative to the ligand’s initial position using the equation below. This calculation was performed to identify the level of deviation of Ligand 1 from its starting position during the simulations.\n",
    "\\begin{equation}\n",
    "RMSD(X, X^{ref}) = \\sqrt{\\frac{1}{n} \\sum_{i=1}^{n} |X_i - X_i^{ref}|^2 }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Possiblie selections include:\n",
    "* for the PreQ1 Backbone specify `nucleic`.\n",
    "* for the Ligand specifiy `resname LIG`.\n",
    "* for both `nucleic or resname LIG`.\n",
    "* for specific residues (e.g residue 13) specify `resid 13`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmsd_plot (topology, trajectory, selection, name):\n",
    "    import MDAnalysis as mda\n",
    "    from MDAnalysis.analysis import rms, align\n",
    "    import matplotlib.pyplot as plt\n",
    "    import numpy as np\n",
    "\n",
    "    # Load the trajectory using DCDReader\n",
    "    u = mda.Universe(topology, trajectory)\n",
    "\n",
    "    # Select backbone atoms\n",
    "    backbone_1 = selection\n",
    "\n",
    "    # Calculate RMSD\n",
    "    R = rms.RMSD(u, u, select= backbone_1)\n",
    "    R.run()\n",
    "\n",
    "    # Access RMSD data\n",
    "    results = R.results.rmsd.T\n",
    "    rmsds= results[2]\n",
    "    time= results[1]\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(time, rmsds)\n",
    "    ax.set_xlabel(\"Time / ps\")\n",
    "    ax.set_ylabel(\"RMSD / $\\\\AA$\")\n",
    "    ax.grid(True)\n",
    "    plt.savefig(name)\n",
    "\n",
    "rmsd_plot ('openmm.parm7', 'thinned_openmm.nc', 'nucleic', 'rmsd_openmm_nucleic')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following this, the root mean square fluctuation (RMSF) for each atom in the PreQ1 riboswitch backbone was calculated relative to its average position during the simulation using the equation below. This was done to identify the areas of flexibility within the PreQ1 riboswitch backbone and hence highlight any potential sources of error for the RBFE calculations.\n",
    "\\begin{equation}\n",
    "\\rho_i = \\sqrt{\\langle (X_i - \\langle X_i \\rangle)^2 \\rangle}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Possiblie selections include:\n",
    "* for the PreQ1 Backbone specify `nucleic`.\n",
    "* for the Ligand specifiy `resname LIG`.\n",
    "* for both `nucleic or resname LIG`.\n",
    "* for specific residues (e.g residue 13) specify `resid 13`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmsf_plot (topology_1, trajectory_1, selection_1, name_1):\n",
    "     u_1 = mda.Universe(topology_1, trajectory_1)\n",
    "     average = align.AverageStructure(u_1,u_1, select= selection_1, ref_frame=0).run()\n",
    "     ref = average.results.universe\n",
    "\n",
    "     aligner = align.AlignTraj(u_1, ref, select= selection_1, in_memory=True).run()\n",
    "     rna = u_1.select_atoms(selection_1)\n",
    "     R = rms.RMSF(rna).run()\n",
    "     fig, ax = plt.subplots()\n",
    "     ax.plot(R.rmsf)\n",
    "     ax.set_xlabel(\"Relative Atom Index\")\n",
    "     ax.set_ylabel(\"RMSF / $\\\\AA$\")\n",
    "     tick_step = np.arange(0,1001, 100)\n",
    "     ax.set_xticks(tick_step)\n",
    "     plt.savefig(name_1)\n",
    "\n",
    "rmsf_plot ('openmm.parm7', 'thinned_openmm.nc', 'nucleic', 'rmsf_openmm_nucleic')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12.0"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}