D-NEMD_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><font color=\"green\">Running dynamical-nonequilibrium simulations in <a href=\"http://www.gromacs.org/\">Gromacs</a></font></h1>     \n",
    "<br />\n",
    "The notebook forms part of the CCPBioSim training.\n",
    "\n",
    "<h2><font color=\"green\">Aim:</font></h2> \n",
    "<p>\n",
    "This workshop will introduce you how to set up and analyse dynamical-nonequilibrium simulations using <a href=\"http://www.gromacs.org/\">Gromacs</a>. It will also provide a quick overview of dynamical-nonequilibrium simulations and how such an approach can be used to study allostery in proteins. \n",
    "<hr />\n",
    "</p>\n",
    "\n",
    "<h2><font color=\"green\">Requirements:</font></h2>\n",
    "<p>\n",
    "This workshop assumes that you are familiar with the <a href=\"http://www.gromacs.org/\">Gromacs</a> software, including its inputs and outputs file formats and have a basic knowledge of the Linux command line. \n",
    "<hr />\n",
    "</p>\n",
    "\n",
    "<h2><font color=\"green\">Author (2021):</font></h2> \n",
    "<p><a href=\"https://asfosite.wordpress.com\">A. Sofia F. Oliveira</a><hr /></p>\n",
    "\n",
    "<h2><font color=\"green\">How to use this notebook</font></h2>\n",
    "    \n",
    "<p><a href=\"http://www.gromacs.org/\">Gromacs</a> is a command-line software, and as such, you will need to open the terminal from the Jupyter home page click on New and then Terminal. To work with this notebook, you will have to switch between the terminal tab and the notebook tab. It might be worthwhile separating them into two windows and arranging them across the screen in such a way that you can look at both simultaneously. Throughout this tutorial, lines beginning with $ will be commands that should be entered into the terminal..</p> \n",
    "<p>You will also need to view some structure files. This can be done using whatever visualization software works for you (e.g. <a href=\"https://pymol.org/2/\">PyMOL</a>, <a href=\"https://www.ks.uiuc.edu/Research/vmd/\">VMD</a>, <a href=\"https://www.cgl.ucsf.edu/chimera/index.html\">Chimera</a>).</p>\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    "<h2><font color=\"green\">Introduction to dynamical-nonequilibrium simulations</font></h2>\n",
    "  \n",
    "<p>\n",
    "Proteins have the ability to respond to changes in their environment. External perturbations, such as voltage, pH and temperature changes, light absorption, or binding of ligands, ions or other proteins, force the systems out of equilibrium and can shift the populations of conformational states, thus favouring specific molecular changes over others. The response to external perturbations can result in structural changes and form the basis for many complex dynamic biological processes, such as allostery.\n",
    "</p>\n",
    "\n",
    "<p>\n",
    "By combining MD simulations in equilibrium and nonequilibrium conditions,  the dynamical-nonequilibrium approach allows for the determination of the time-dependent structural response of a system (for more details, see \n",
    "<a href=\"https://www.tandfonline.com/doi/full/10.1080/08927022.2015.1121543\">Ciccotti, 2016</a> and <a href=\"https://epjb.epj.org/articles/epjb/abs/2021/07/10051_2021_Article_157/10051_2021_Article_157.html\">Oliveira, 2021</a>).\n",
    "</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "  <img src=\"pics/Fig1_300dpi.png\" width=500 height=179>\n",
    "  <figcaption><center>Fig 1 - Schematic representation of the dynamical-nonequilibrium MD approach. The equilibrium MD trajectory (blue line) provides the initial distribution for the nonequilibrium simulations at time t= 0 ns (red lines).</center></figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "The rationale for the approach is simple: if we introduce an external perturbation to a simulation sampling an equilibrium state and, by doing so, start a parallel nonequilibrium simulation, the structural response of the protein to the perturbation can be directly measured by comparing the equilibrium and nonequilibrium trajectories at equivalent points in time (as long as enough sampling is gathered). Note that the nature of the perturbation introduced in the system will depend on the question to be answered. For proteins, it can be, for example, a chemical reaction or the binding/unbinding of ligands or/and ions. \n",
    "</p>\n",
    "\n",
    "<h3>Now, let the fun begin!</h3>\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2><font color=\"green\">Visualise the protein-ligand complex</font></h2>\n",
    "\n",
    "<p>In this workshop, we will utilise as an example a nicotinic acetylcholine receptor. These receptors are cation-selective channels that belong to the superfamily of pentameric ligand-gated ion channels. These receptors are expressed throughout the nervous system, where they contribute to many functions, including cognition and nicotine addiction. The binding of a ligand, such as acetylcholine or nicotine, to the receptor triggers channel opening, allowing for a flow of positive ions to diffuse across the membrane.</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "  <img src=\"pics/Fig2_300dpi.png\" width=300 height=403>\n",
    "  <figcaption><center>Fig 2 - 3D structure of a nicotinic acetylcholine receptor.</center></figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Lets start by having a look into the pdbs directory.</p>\n",
    "\n",
    "<code> $ cd pdbs/</code>\n",
    "\n",
    "<code> $ ls </code>\n",
    "\n",
    "<p>Inside this directory, you will find a file called nicotinic_receptor.pdb containing a 3D structure of the complex formed by the receptor and nicotine inserted into a lipid membrane. In the cell below you can use the NGLview widget to look at the complex (hint: nct is the residue name of nicotine).</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d67858b3ff649b0936eaab877cb9a6c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>NGLWidget</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "NGLWidget()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import nglview\n",
    "\n",
    "view = nglview.show_file(\"pdbs/nicotinic_receptor.pdb\")\n",
    "view\n",
    "view.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2><font color=\"green\">Setup of the dynamical-nonequilibrium simulations</font></h2>    \n",
    "<p>\n",
    "The procedure to set up the dynamical-nonequilibrium simulations is easy and consists of, first, running long equilibrium MD simulations for the equilibrium reference state (unperturbed simulations). The equilibrium simulation(s) will generate a distribution of configurations that will be used as the starting points for an ensemble of short nonequilibrium trajectories under the effect of the external perturbation (perturbed simulations).</p>\n",
    "<p>Note that the perturbation to be introduced in the system depends on whatever question the user wants to address, and it can be either instantaneous or gradual <a href=\"https://www.tandfonline.com/doi/full/10.1080/08927022.2015.1121543\">Ciccotti, 2016</a> and <a href=\"https://epjb.epj.org/articles/epjb/abs/2021/07/10051_2021_Article_157/10051_2021_Article_157.html\">Oliveira, 2021</a>). Note also that in order to obtain statistically significant responses, tens/hundreds of short nonequilibrium trajectories may be required.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>In this workshop, you will use the removal of nicotine from the receptor's orthosteric binding site as the external perturbation (similarly to <a href=\"https://www.cell.com/structure/fulltext/S0969-2126(19)30130-3?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0969212619301303%3Fshowall%3Dtrue\n",
    "    \">Oliveira, 2019a</a> and <a href=\"https://pubs.acs.org/doi/abs/10.1021/jacs.9b09055\">Oliveira, 2019b</a>), and you will learn all the steps needed to set up and analyse a single dynamical-nonequilibrium simulation. However, remember that you will need to prepare many more nonequilibrium simulations to obtain statistically significant results (for some examples, see <a href=\"https://epjb.epj.org/articles/epjb/abs/2021/07/10051_2021_Article_157/10051_2021_Article_157.html\">Oliveira, 2021</a>).</p>\n",
    "<p>\n",
    "</p>    \n",
    "<figure>\n",
    "  <img src=\"pics/Fig3_300dpi.png\" width=500 height=403>\n",
    "  <figcaption><center>Fig 3 - Scheme for the dynamical-nonequilibrium simulations.</center></figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Let us first create a new folder to place all the necessary files to run the dynamical-nonequilibrium simulation:</p>\n",
    "\n",
    "<code> $ cd ../</code>\n",
    "\n",
    "<code> $ mkdir d-nemd </code>\n",
    "\n",
    "<code> $ cd d-nemd/</code>\n",
    "\n",
    "<p>Next, you can use the <a href=\"https://manual.gromacs.org/2021.3/onlinehelp/gmx-trjconv.html\">trjconv</a> tool available in gromacs with the <code>-dump</code> and <code>-tu</code> options to extract from the equilibrium simulation the initial configurations for the short nonequilibrium trajectories. These configuration files need to include all the atoms in the system except the ones belonging to nicotine. Note that in the following command, you will extract the system's configuration at t= 50 ns.</p>  \n",
    "\n",
    "<code> $ gmx trjconv  -f ../data/example1/traj.xtc -o frame.pdb -s ../data/example1/md.tpr -n ../data/example1/index.ndx -dump 50 -tu ns </code>\n",
    "\n",
    "<p>and select the <code>System_without_nicotine</code> group for output. If you saved the atom velocities from the equilibrium simulations, please use .gro or .g96 file format for the output instead of the .pdb format.</p>   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>The next step is to remove nicotine from the topology file. This can  be easilydone by commenting (adding `;` in the beginning of the line) or removing the following lines from the topol.top file:</p>\n",
    "\n",
    "<code> ; #include \"nct.itp\" </code>\n",
    "\n",
    "<code> ; nct                 2 </code>\n",
    "\n",
    "<p>To make the necessary changes to the topol.top file you can either use whatever text editor works for you or the <code>sed</code> command:</p>\n",
    "\n",
    "<code> $ sed 's/#include \"nct.itp\"/;#include \"nct.itp\"/;s/nct                 2/;nct                 2/' ../data/example1/topol.top > topol.top </code>\n",
    "\n",
    "<p>Don't forget that you also need to copy all the itp files describing your system and, if needed, the force-field to the d-nemd folder</p>\n",
    "\n",
    "<code>$ cp ../data/example1/*itp . </code>\n",
    "\n",
    "<code>$ cp -r ../data/example1/AMBER99SB-ILDN_SLIPIDS.ff . </code>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Next, you will need to remove any reference to nicotine from the mdp parameters files (including from the temperature coupling and energy groups) while keeping everything else as it was for the equilibrium simulations. If needed, you can also change the number of step (<code>nsteps</code> option) in the mdp file. To do all the necessary changes to the mdp file, you can either use a text editor or the <code>sed</code> command:</p>\n",
    "\n",
    "<code> $ sed 's/Protein_nct/Protein/;s/ nct / /' ../data/example1/md.mdp > md.mdp </code>\n",
    "\n",
    "<p><font color=\"red\"><strong>Important note:</strong></font> if you saved the atom velocities from the equilibrium simulations, please use them and do not generate new velocities in the beginning of the nonequilibrium simulations (<code>gen_vel = no </code>)</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Finally, you will need to create a new index file using the <a href=\"https://manual.gromacs.org/2021.3/onlinehelp/gmx-make_ndx.html?highlight=make_ndx\">make_ndx</a> tool available in gromacs:</p>\n",
    "\n",
    "<code> $ gmx make_ndx  -f frame.pdb -o index.ndx </code>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><h3>Congratulations! You are now ready to run your dynamical-nonequilibrium simulation.</h3></p> \n",
    "\n",
    "<p>For that, you just need to use <a href=\"https://manual.gromacs.org/current/onlinehelp/gmx-grompp.html\">grompp</a> and <a href=\"https://manual.gromacs.org/current/onlinehelp/gmx-mdrun.htmlmdrun\">mdrun</a> as in any equilibrium MD simulation. Note that due to the system size, the dynamical-nonequilibrium simulation that you just prepared will take a while to run. As such you will use for the next section of the tutorial, previously calculated trajectory files.</p> \n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2><font color=\"green\">Analysis of the dynamical-nonequilibrium simulations</font></h2>    \n",
    "\n",
    "<p>The response of the system to the perturbation introduced can be extracted using the Kubo-Onsager relation through the direct comparison of the observable of interest between each pair of unperturbed equilibrium, and \"branching\" perturbed nonequilibrium trajectories at equivalent time t, followed by the averaging of the difference over all pairs of simulations.</p> \n",
    "<p>For a description of the theoretical framework underlying the Kubo-Onsager relation, see <a href=\"https://www.tandfonline.com/doi/full/10.1080/08927022.2015.1121543\">Ciccotti, 2016</a> and <a href=\"https://epjb.epj.org/articles/epjb/abs/2021/07/10051_2021_Article_157/10051_2021_Article_157.html\">Oliveira, 2021</a>.</p>\n",
    "\n",
    "<figure>\n",
    "  <img src=\"pics/Fig4_300dpi.png\" width=500 height=179>\n",
    "  <figcaption><center>Fig 4 - Scheme for the analysis of the dynamical-nonequilibrium simulations.</center></figcaption>\n",
    "</figure>\n",
    " \n",
    "<p>Note how straightforward it is with this approach to measure the system's response as it simply requires the subtraction of the observable between two trajectories. In particular, for the structural response of proteins, the difference in the position of each individual Cα atom between pairs of equilibrium and nonequilibrium trajectories at a given time point is determined and afterwards averaged over all pairs of simulations. The time evolution of the average difference in position will highlight the cascade of events associated with the protein’s structural response. <strong>Now, let's start!</strong></p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>In the following example, you will analyse the structural response of the receptor at 0, 100 and 5000 ps after nicotine removal. All necessary files to extract the receptor's response are located inside the data/example2 directory. Let us have a look into the data/example2 directory.</p>\n",
    "\n",
    "<code>$ cd ../data/example2 </code>\n",
    "\n",
    "<code>$ ls </code>\n",
    "\n",
    "<p>There you will find two directories called equilibrium_simulations and nonequilibrium_simulations. Each folder contains an index file, the trajectory file and its corresponding tpr file.</p> \n",
    "\n",
    "<code>$ ls equilibrium_simulations/ </code> \n",
    "\n",
    "<code>$ ls nonequilibrium_simulations/ </code>    \n",
    "    \n",
    "<p>The nonequilibrium trajectory inside the nonequilibrium_simulations directory started from the equilibrium configuration at t= 50 ns (the configuration you have extracted in the <font color=\"green\">“Setup of the dynamical-nonequilibrium simulations”</font> section above).</p>\n",
    "\n",
    "<figure>\n",
    "  <img src=\"pics/Fig5_300dpi.png\" width=400 height=179>\n",
    "  <figcaption><center>Fig 5 - Scheme for the analysis of the dynamical-nonequilibrium simulations that have started from the equilibrium configuration at t= 50 ns.</center></figcaption>\n",
    "</figure>\n",
    "\n",
    "<p>Next, let us create a new folder for the analysis of the data:</p>\n",
    "\n",
    "<code> $ cd ../../ </code>\n",
    "\n",
    "<code> $ mkdir -p d-nemd_analysis/equilibrium d-nemd_analysis/nonequilibrium d-nemd_analysis/deviations</code>\n",
    "\n",
    "<code> $ cd d-nemd_analysis/ </code>\n",
    "\n",
    "<p>Now extract the position of the individual Cα atom from the dynamical-nonequilibrium simulation at 0, 100 and 5000 ps. Please don't forget to correct the periodic boundary conditions before extracting the Cα atom xyz positions using the -fit option.</p>\n",
    "\n",
    "<code> $ cd nonequilibrium </code>\n",
    "\n",
    "<code> $ gmx trjconv -f ../../data/example2/nonequilibrium_simulations/traj.xtc -o conf_0ps.pdb -s ../../data/example2/nonequilibrium_simulations/md.tpr -n ../../data/example2/nonequilibrium_simulations/index.ndx -fit rot+trans -dump 0 -tu ps </code>\n",
    "\n",
    "<p>Please choose the C-alpha group for both least squares fit and output.</p>\n",
    "\n",
    "<code> $ gmx trjconv -f ../../data/example2/nonequilibrium_simulations/traj.xtc -o conf_100ps.pdb -s ../../data/example2/nonequilibrium_simulations/md.tpr -n ../../data/example2/nonequilibrium_simulations/index.ndx -fit rot+trans -dump 100 -tu ps </code>\n",
    "\n",
    "<code> $ gmx trjconv -f ../../data/example2/nonequilibrium_simulations/traj.xtc -o conf_5000ps.pdb -s ../../data/example2/nonequilibrium_simulations/md.tpr -n ../../data/example2/nonequilibrium_simulations/index.ndx -fit rot+trans -dump 5000 -tu ps </code>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Now let us extract the position of the Cα atoms from the equilibrium simulation at equivalent points in time. Again do not forget to correct the periodic boundary conditions before extracting the atomic positions.</p>\n",
    "\n",
    "<code> $ cd ../equilibrium </code>\n",
    "\n",
    "<code> $ gmx trjconv -f ../../data/example2/equilibrium_simulations/traj.xtc -o conf_50000ps.pdb -s ../../data/example2/equilibrium_simulations/md.tpr -n ../../data/example2/equilibrium_simulations/index.ndx -fit rot+trans -dump 50000 -tu ps </code>\n",
    "\n",
    "<p>As before, please choose the C-alpha group for both least squares fit and output.</p>\n",
    "\n",
    "<code> $ gmx trjconv -f ../../data/example2/equilibrium_simulations/traj.xtc -o conf_50100ps.pdb -s ../../data/example2/equilibrium_simulations/md.tpr -n ../../data/example2/equilibrium_simulations/index.ndx -fit rot+trans -dump 50100 -tu ps </code>\n",
    "\n",
    "<code> $ gmx trjconv -f ../../data/example2/equilibrium_simulations/traj.xtc -o conf_55000ps.pdb -s ../../data/example2/equilibrium_simulations/md.tpr -n ../../data/example2/equilibrium_simulations/index.ndx -fit rot+trans -dump 55000 -tu ps </code>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Next, you will need to calculate the deviation of each Cα atom between each pair of equilibrium and nonequilibrium simulations.</p> \n",
    " \n",
    "<code> $ cd ../deviations </code>\n",
    "    \n",
    "<p>This calculation can be done by directly comparing the positions of individual Cα atoms at equivalent points in time (namely after 0, 100 and 5000 ps of simulation) for every pair of unperturbed equilibrium and perturbed nonequilibrium simulations.</p> \n",
    "<p>Note that as we are only calculating the difference in position between 2 conformations, no special code is needed and, any tool able to calculate the distance between two atoms can be used for this analysis. <a href=\"http://www.gromacs.org/\">Gromacs</a>, in particular, has several tools that can be used for this calculation. In the following example, you will use the <a href=\"https://manual.gromacs.org/current/onlinehelp/gmx-rmsf.html\">rmsf</a> tool, but both the <a href=\"https://manual.gromacs.org/current/onlinehelp/gmx-rms.html\">rms</a> or <a href=\"https://manual.gromacs.org/current/onlinehelp/gmx-distance.html\">distance</a> tools can also be used for this calculation. To get help about the options available in each tool, please either check the <a href=\"https://manual.gromacs.org/documentation/current/reference-manual/index.html\">Gromacs reference manual</a> or use the <code>-h</code> option.</p>\n",
    "\n",
    "<code>$ gmx rmsf -f ../nonequilibrium/conf_0ps.pdb -s ../equilibrium/conf_0ps.pdb -o rms_0ps.xvg -od rmsdev_0ps.xvg -fit </code>\n",
    "\n",
    "Please choose the C-alpha group for the root mean square calculation.\n",
    "\n",
    "<code> $ gmx rmsf -f ../nonequilibrium/conf_100ps.pdb -s ../equilibrium/conf_100ps.pdb -o rms_100ps.xvg -od rmsdev_100ps.xvg -fit </code>\n",
    "\n",
    "<code> $ gmx rmsf -f ../nonequilibrium/conf_5000ps.pdb -s ../equilibrium/conf_5000ps.pdb -o rms_5000ps.xvg -od rmsdev_5000ps.xvg -fit </code> \n",
    "\n",
    "<code> rm rms_*.xvg </code>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>You can then use your favourite plotting software (for example, <a href=\"http://www.gnuplot.info/\">Gnuplot</a>) to plot the response over time. <font color=\"red\">Question: why are the deviations zero at t= 0 ns?</font><p/> \n",
    "\n",
    "<figure>\n",
    "  <img src=\"pics/Fig6_300dpi.png\" width=600 height=179>\n",
    "  <figcaption><center>Fig 6 - Time evolution of the response of the Cα atoms of Chain A and B.</center></figcaption>\n",
    "</figure>\n",
    "\n",
    "<p>It is very important to note that the previous plots represent the protein's response from a single dynamical-nonequilibrium replicate, and as such, these results are not statistically significant. Many more (usually tens/hundreds) dynamical-nonequilibrium simulations are needed to achieve statistical significance (<a href=\"https://www.tandfonline.com/doi/full/10.1080/08927022.2015.1121543\">Ciccotti, 2016</a> and <a href=\"https://epjb.epj.org/articles/epjb/abs/2021/07/10051_2021_Article_157/10051_2021_Article_157.html\">Oliveira, 2021</a>).<p/> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Once you have extracted the response from the ensemble (tens/hundreds) of dynamical-nonequilibrium simulations performed, you can average the results and compute its associated statistical errors to assess the significance of the results (Fig 7). If needed, the statistical error associated with the response can be made as small as necessary just by increasing the number of nonequilibrium trajectories. Note in Fig 7, how the error associated with the response decreases with an increasing number of dynamical-nonequilibrium trajectories.<p/> \n",
    "\n",
    "<figure>\n",
    "  <img src=\"pics/Fig7_300dpi.png\" width=600 height=179>\n",
    "  <figcaption><center>Fig 7 - Average Cα-positional deviations (black line) and corresponding standard errors (red vertical lines) of chain A in the 5 ns after nicotine removal. The average deviations were calculated using the Kubo-Onsager relation for the pairwise comparison between the nonequilibrium and equilibrium simulations.</center></figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><font color=\"red\"><strong>It is also very important to highlight that it is not possible to <em>a priori</em> know the amount of sampling needed to obtain statistically significant structural responses, as this will depend not only on the system but also on the perturbation introduced.</strong></font> Determining the statistical errors associated with the structural response is important to test if the sampling gathered is sufficient.</p>\n",
    "\n",
    "Additionaly, you can also map the average Cα positional deviations (at different point in time) into a 3D structure in order to identify the spatial location of the residues that respond to the perturbation introduced. This can be done, for example, by using the last column (B-factor colum) in a pdb file (Fig 8). \n",
    "\n",
    "<figure>\n",
    "  <img src=\"pics/Fig8_300dpi.png\" width=700 height=179>\n",
    "  <figcaption><center>Fig 8 - Signal propagation pathway from the agonist binding site to the transmembrane domain. Average Cα-positional deviation at times 0, 0.05, 0.5, 1, 3 and 5 ns following nicotine removal from the binding pockets.</center></figcaption>\n",
    "</figure>\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2><font color=\"green\">Summary</font></h2> \n",
    "\n",
    "<p>In this tutorial, you have learned the basic concepts and how to set up and analyse dynamical-nonequilibrium simulations. It is important to bear in mind that in this type of simulations, the perturbation(s) introduced in the system depend on the question and user's choice, and that different perturbations may give different results and might identify different communication networks. As such, careful thought needs to be given to the choice of the perturbation in order to figure out which is the most (biologically) relevant trigger(s) behind the functionally relevant conformational changes you are interested in studying.</p>\n",
    "<p>Before running any dynamical-nonequilibrium simulations, I suggest you to review the literature (e.g. <a href=\"https://link.springer.com/article/10.1007/BF01011477\">Ciccotti, 1979</a>, <a href=\"https://www.mdpi.com/1099-4300/16/1/233\">Ciccotti, 2014</a>, <a href=\"https://www.tandfonline.com/doi/full/10.1080/08927022.2015.1121543\">Ciccotti, 2016</a> and <a href=\"https://epjb.epj.org/articles/epjb/abs/2021/07/10051_2021_Article_157/10051_2021_Article_157.html\">Oliveira, 2021</a>) and check the <a href=\"https://manual.gromacs.org/documentation/current/reference-manual/index.html\">Gromacs reference manual</a> for any changes to the tools needed to set up and analyse these simulations. This tutorial should not be viewed as comprehensive.</p>\n",
    "\n",
    "<p>If you have any suggestions on how to improve this tutorial or if you notice a mistake, feel free to <a href=\"sofia.oliveira@bristol.ac.uk\">email me</a>.</p>\n",
    "<hr />"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}