04_protein/01_villin_headpiece.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular dynamics of an explicitly solvated small protein\n",
    "\n",
    "For the simulations in this notebook, we will make use of the crystal structure of the Villin Headpiece subdomain [1YRF](https://www.rcsb.org/structure/1YRF). This is relatively small protein (35 residues) consisting of three alpha helices. This is a prototypical fast folding protein, see [10.1016/j.jmb.2006.03.034](https://doi.org/10.1016/j.jmb.2006.03.034) and [10.1073/pnas.0502495102](https://doi.org/10.1073/pnas.0502495102), and it is therefore a popular benchmark for protein folding molecular dynamics simulations, e.g. [10.1073/pnas.1800690115](https://doi.org/10.1073/pnas.1800690115) (force field accuracy) and [10.1109/SC.2014.9](https://doi.org/10.1109/SC.2014.9) (computational performance).\n",
    "\n",
    "If you want to run a longer simulation in this notebook, after increasing the number of time steps below, it may be useful to submit this notebook as a non-interactive job.\n",
    "This is explained in section `03_noninteractive` of the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shutil\n",
    "from sys import stdout\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import mdtraj\n",
    "import nglview\n",
    "import pandas\n",
    "from openmm import *\n",
    "from openmm.app import *\n",
    "from openmm.unit import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PDB file for 1YRF contains more information than we need. It holds several possible states of some residues mixed into one file and OpenMM cannot handle this. Most atomic positions are the same for all these states, but a few have an extra `A`, `B`, `C`, ... just after the atom name. The following function splits such a PDB file into multiples ones. Each of these files is suitable for starting an OpenMM simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_pdb(fn_pdb):\n",
    "    pos = 16\n",
    "    groups = {}\n",
    "    counter = 0\n",
    "    with open(fn_pdb) as f:\n",
    "        for line in f:\n",
    "            if line.startswith(\"ATOM\"):\n",
    "                state = line[pos]\n",
    "                line = line[:pos] + \" \" + line[pos + 1 :]\n",
    "                groups.setdefault(state, []).append((counter, line))\n",
    "                counter += 1\n",
    "    print(\"Group keys:\", groups.keys())\n",
    "    for key, lines_group in groups.items():\n",
    "        if key == \" \":\n",
    "            continue\n",
    "        lines_both = lines_group + groups[\" \"]\n",
    "        lines_both.sort()\n",
    "        new_fn_pdb = f\"{fn_pdb[:-4]}_{key.lower()}.pdb\"\n",
    "        with open(new_fn_pdb, \"w\") as f:\n",
    "            for _counter, line in lines_both:\n",
    "                f.write(line)\n",
    "        print(\"Written\", new_fn_pdb)\n",
    "\n",
    "\n",
    "split_pdb(\"1yrf.pdb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "> Inspect the contents of the files `1yrf_a.pdb` and `1yrf_b.pdb`. Verify for at least one atom that these files contain the correct information taken from `1yrf.pdb`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdb = PDBFile(\"1yrf_a.pdb\")\n",
    "modeller = Modeller(pdb.topology, pdb.positions)\n",
    "forcefield = ForceField(\"amber14-all.xml\", \"amber14/tip3pfb.xml\")\n",
    "modeller.addHydrogens(forcefield)\n",
    "modeller.addSolvent(forcefield, model=\"tip3p\", padding=1 * nanometer)\n",
    "system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, constraints=HBonds)\n",
    "temperature = 300 * kelvin\n",
    "pressure = 1 * bar\n",
    "integrator = LangevinIntegrator(temperature, 1 / picosecond, 2 * femtoseconds)\n",
    "system.addForce(MonteCarloBarostat(pressure, temperature))\n",
    "simulation = Simulation(modeller.topology, system, integrator)\n",
    "simulation.context.setPositions(modeller.positions)\n",
    "simulation.minimizeEnergy(maxIterations=100)\n",
    "positions = simulation.context.getState(getPositions=True).getPositions()\n",
    "with open(\"init.pdb\", \"w\") as f:\n",
    "    PDBFile.writeFile(simulation.topology, positions, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> Which changes to the input structure can be made, for which the force field can still be applied? Try the following:\n",
    ">\n",
    "> - Remove an atom from the PDB file.\n",
    "> - Remove an entire residue from the PDB file.\n",
    "> - Use an engineered form of the Villin Headpiece, e.g. the Lys24Nle mutant (PDB 1WY3) or the Lys24Nle/Lys29Nle double mutant (PDB 2F4K).\n",
    ">\n",
    "> Finally, restore the input to the original one and rerun the above code cell, to have a good starting point for the next cell.\n",
    ">\n",
    "> With the methodology shown here, it is not possible to define custom mutations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.reporters = []\n",
    "simulation.reporters.append(DCDReporter(\"traj.dcd\", 10))\n",
    "simulation.reporters.append(\n",
    "    StateDataReporter(stdout, 100, step=True, temperature=True, elapsedTime=True)\n",
    ")\n",
    "simulation.reporters.append(\n",
    "    StateDataReporter(\n",
    "        \"scalars.csv\",\n",
    "        10,\n",
    "        step=True,\n",
    "        time=True,\n",
    "        potentialEnergy=True,\n",
    "        totalEnergy=True,\n",
    "        temperature=True,\n",
    "    )\n",
    ")\n",
    "simulation.step(30000)\n",
    "\n",
    "# The last line is only needed for Windows users,\n",
    "# to close the DCD file before it can be opened by nglview.\n",
    "del simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will need the output files of this simulation for the following notebook. Copy over the files `init.pdb`, `scalars.csv` and `traj.dcd` to the directory `../05_analysis`. This can be done with the following code, or with any file manager that comes with the operating system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shutil.copy(\"init.pdb\", \"../05_analysis\")\n",
    "shutil.copy(\"scalars.csv\", \"../05_analysis\")\n",
    "shutil.copy(\"traj.dcd\", \"../05_analysis\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> Make an estimate of the computational cost (wall time) on your current hardware to run a sufficiently long simulation to observe spontaneously a protein folding event for this fast folder. The required simulation time is approximately $4 \\mathrm{\\mu s}$. Would such a calculations be feasible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pandas.read_csv(\"scalars.csv\")\n",
    "plt.close(\"ener\")\n",
    "fig, ax = plt.subplots(num=\"ener\")\n",
    "df.plot(kind=\"line\", x=\"Time (ps)\", y=\"Potential Energy (kJ/mole)\", ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = mdtraj.load(\"traj.dcd\", top=\"init.pdb\")\n",
    "view = nglview.show_mdtraj(traj)\n",
    "view"
   ]
  }
 ],
 "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}