05_analysis/02_physicochemical_properties.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physicochemical properties\n",
    "\n",
    "This notebook demonstrates some of the higher-level analysis tools from MDTraj that focus on the physicochemical properties of a macromolecule. These can be useful to understand (de)stabilizing factors. A complete overview of analysis functions in MDTraj can be found here: http://mdtraj.org/latest/analysis.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "before-import"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import mdtraj\n",
    "import nglview\n",
    "import numpy as np\n",
    "import pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = mdtraj.load(\"traj.dcd\", top=\"init.pdb\")\n",
    "traj.restrict_atoms(traj.topology.select(\"protein\"))\n",
    "df = pandas.read_csv(\"scalars.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Hydrogen bonding\n",
    "\n",
    "(Internal) hydrogen bonds are one of the important factors that stabilize protein folds and also in other macromolecular systems they often play an important role. We will use the [Wernet-Nillson criterion](http://mdtraj.org/latest/api/generated/mdtraj.wernet_nilsson.html) implemented in MDTraj to count the number of hydrogen bonds in each frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hbonds = mdtraj.wernet_nilsson(traj)\n",
    "numhbonds = [len(pairs) for pairs in hbonds]\n",
    "plt.close(1)  # This is needed to rerun the code cell correctly\n",
    "fig, ax1 = plt.subplots(num=1)\n",
    "ax1.plot(df[\"Time (ps)\"], numhbonds, label=\"Num. HBonds\")\n",
    "ax1.set_xlabel(\"Time [ps]\")\n",
    "ax1.set_ylabel(\"Number of hydrogen bonds\", color=\"C0\")\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(df[\"Time (ps)\"], df[\"Potential Energy (kJ/mole)\"], color=\"C1\", label=\"Pot. Energy\")\n",
    "ax2.set_ylabel(\"Potential energy [kJ/mol]\", color=\"C1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The potential energy finds its equilibrium value twice as fast compared to the number of hydrogen bonds. This shows (again) that one should not just trust a single property to determine the equilibration phase. Most of the increase of the potential energy is simply the thermal activation of the solvent, which is not necessarily representative for the solute."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Secondary structure assignment\n",
    "\n",
    "The secondary structure assignment of every residue (at each time step) can be obtained with the [compute_dssp](http://mdtraj.org/latest/api/generated/mdtraj.compute_dssp.html) function of MDTraj. For small proteins like the Villin headpiece, the results can be plotted conveniently with Matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dssp = mdtraj.compute_dssp(traj)\n",
    "# Convert string characters to numbers\n",
    "dssp_num = np.vectorize(lambda char: \"HEC\".index(char))(dssp)\n",
    "plt.close(2)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=2)\n",
    "im = ax.pcolormesh(df[\"Time (ps)\"], np.arange(dssp.shape[1]), dssp_num.T, shading=\"nearest\")\n",
    "cbar = fig.colorbar(im, values=[0, 1, 2], ticks=[0, 1, 2])\n",
    "cbar.ax.set_yticklabels([\"H\", \"E\", \"C\"])\n",
    "ax.set_xlabel(\"Time [ps]\")\n",
    "ax.set_ylabel(\"Residue index\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example already shows nicely the presence of the three alpha helices and how they change somewhat over time.\n",
    "\n",
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> Modify the above example to plot the 8-level DSSP assignment, instead of the simplified 3-level assignment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Radial distribution function\n",
    "\n",
    "The coordination of a functional group by another chemical moiety can be characterized with a [radial distribution function](https://en.wikipedia.org/wiki/Radial_distribution_function) (RDF). In general, a radial distribution function is a ratio of two probability densities as function of an inter-atomic distance:\n",
    "\n",
    "- In the numerator, one puts the empirically observed probability density for finding a given pair of atoms.\n",
    "- The denominator contains the probability in case the two atoms would be uniformly distributed over space.\n",
    "\n",
    "In case of dense liquids, the RDF converges to one for large distances because there is no long-range structure in a normal liquid. In case of a single isolated system (like a protein), the RDF converges to zero because all empirical inter-atomic distances are bound by the size of the molecule.\n",
    "\n",
    "We can use [compute_rdf](http://mdtraj.org/latest/api/generated/mdtraj.compute_rdf.html) in combination with [select_pairs](http://mdtraj.org/latest/api/generated/mdtraj.Topology.html#mdtraj.Topology.select_pairs) to easily generate RDFs.\n",
    "\n",
    "In the example below, we compute the RDF for all amine nitrogens with amide oxygens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r, gr = mdtraj.compute_rdf(\n",
    "    traj,\n",
    "    traj.topology.select_pairs(\"protein and name N\", \"protein and name O\"),\n",
    "    r_range=(0.1, 0.5),\n",
    "    bin_width=0.0005,\n",
    ")\n",
    "plt.close(3)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=3)\n",
    "ax.plot(r, gr)\n",
    "ax.set_xlabel(\"Radius [nm]\")\n",
    "ax.set_ylabel(\"Radial distribution function [gr]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The peak at approx $0.225 \\textrm{nm}$ is due to internal hydrogen bonds, mostly in alpha helices. Subsequent peaks are likely caused by amide groups in the alpha helices forming cooperative hydrogen bonds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Solvent accessible surface area\n",
    "\n",
    "The Shrake-Rupley method can be used to measure the solvent-accessible surface area of each residue, showing which are exposed to the solvent water molecules. For this, we use the [shrake_rupley](http://mdtraj.org/latest/api/generated/mdtraj.shrake_rupley.html) function from MDTraj."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sasa = mdtraj.shrake_rupley(traj[::10], mode=\"residue\")\n",
    "plt.close(4)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=4)\n",
    "im = ax.pcolormesh(df[\"Time (ps)\"][::10], np.arange(dssp.shape[1]), sasa.T, shading=\"nearest\")\n",
    "cbar = fig.colorbar(im)\n",
    "cbar.set_label(\"Solvent accessible surface areas [nm^2]\")\n",
    "ax.set_xlabel(\"Time [ps]\")\n",
    "ax.set_ylabel(\"Residue index\")"
   ]
  }
 ],
 "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
}