HH.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Hodgkin-Huxley model\n",
    "\n",
    "In this notebook we will use Python to simulate the Hodgkin-Huxley (HH) neuron model.  This model is arguably the *most* important computational model in neuroscience.  We'll focus here on simulating this model and understanding its pieces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background information about the HH model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a video that describes some of the biophysical details of the HH model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://player.vimeo.com/video/140084450\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.VimeoVideo at 0x118650210>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.lib.display import VimeoVideo\n",
    "VimeoVideo('140084450')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some additional usual videos and references:\n",
    "\n",
    "- <a href=\"http://klewel.com/conferences/epfl-neural-networks/index.php?talkID=4\" rel=\"external\">Lecture by Prof. Gerstner, *Detailed Neuron Model (a)*</a>\n",
    "\n",
    "- <a href=\"http://klewel.com/conferences/epfl-neural-networks/index.php?talkID=5\" rel=\"external\">Lecture by Prof. Gerstner, *Detailed Neuron Model (b)*</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Preliminaries \n",
    "\n",
    "Before beginning, let's load in the Python packages we'll need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: MacOSX\n"
     ]
    }
   ],
   "source": [
    "from pylab import *\n",
    "%matplotlib\n",
    "rcParams['figure.figsize']=(12,3)                # Change the default figure size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, let's import the functions we'll need to simulate the HH model, which are available on this repository:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from HH_functions import HH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Part 1:   The Hodgkin-Huxley (HH) equation code.\n",
    "\n",
    "To start, let's examine the code for the HH model. We can do so in (at least) two ways.\n",
    "\n",
    "- Go to the Case Studies repository, and examine the Python file\n",
    "`HH_functions.py`\n",
    "- Examine the code inline with `inspect`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(['def HH(I0,T0):\\n',\n",
       "  '    dt = 0.01;\\n',\n",
       "  '    T  = math.ceil(T0/dt)  # [ms]\\n',\n",
       "  '    gNa0 = 120   # [mS/cm^2]\\n',\n",
       "  '    ENa  = 115;  # [mV]\\n',\n",
       "  '    gK0  = 36;   # [mS/cm^2]\\n',\n",
       "  '    EK   = -12;  # [mV]\\n',\n",
       "  '    gL0  = 0.3;  # [mS/cm^2]\\n',\n",
       "  '    EL   = 10.6; # [mV]\\n',\n",
       "  '\\n',\n",
       "  '    t = np.arange(0,T)*dt\\n',\n",
       "  '    V = np.zeros([T,1])\\n',\n",
       "  '    m = np.zeros([T,1])\\n',\n",
       "  '    h = np.zeros([T,1])\\n',\n",
       "  '    n = np.zeros([T,1])\\n',\n",
       "  '\\n',\n",
       "  '    V[0]=-70.0\\n',\n",
       "  '    m[0]=0.05\\n',\n",
       "  '    h[0]=0.54\\n',\n",
       "  '    n[0]=0.34\\n',\n",
       "  '\\n',\n",
       "  '    for i in range(0,T-1):\\n',\n",
       "  '        V[i+1] = V[i] + dt*(gNa0*m[i]**3*h[i]*(ENa-(V[i]+65)) + gK0*n[i]**4*(EK-(V[i]+65)) + gL0*(EL-(V[i]+65)) + I0);\\n',\n",
       "  '        m[i+1] = m[i] + dt*(alphaM(V[i])*(1-m[i]) - betaM(V[i])*m[i]);\\n',\n",
       "  '        h[i+1] = h[i] + dt*(alphaH(V[i])*(1-h[i]) - betaH(V[i])*h[i]);\\n',\n",
       "  '        n[i+1] = n[i] + dt*(alphaN(V[i])*(1-n[i]) - betaN(V[i])*n[i]);\\n',\n",
       "  '    return V,m,h,n,t\\n'],\n",
       " 22)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import inspect\n",
    "inspect.getsourcelines(HH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  Examine this code.  Can you make sense of it?  Can you identify the\n",
    "gating variables?  The rate functions?  The equations that define the dynamics?\n",
    "We'll answer these questions in this in notebook, but try so on your own first.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whenever examining code, it's useful to consider the *inputs* to the code, and the *outputs* produced by the code.  There are two inputs to `HH0`:\n",
    "\n",
    "- `I0` = the current we inject to the neuron.\n",
    "- `T0` = the total time of the simulation in [ms].\n",
    "\n",
    "And there are five outputs:\n",
    "\n",
    "- `V` = the voltage of neuron.\n",
    "- `m` = activation variable for Na-current.\n",
    "- `h` = inactivation variable for Na-current.\n",
    "- `n` = activation variable for K-current.\n",
    "- `t` = the time axis of the simulation (useful for plotting)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2:  At low input current (`I0`), examine the HH dynamics.\n",
    "\n",
    "  To understand how the HH model works, we'll start by focusing on the\n",
    "  case when `I0` is small. Let's fix the input current to zero,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "I0 = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and let's simulate the model for 100 ms,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "T0 = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've now defined both inputs to the `HH` function, and can execute it, as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "[V,m,h,n,t]=HH(I0,T0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the function returns five outputs, which we assign to the variables `V`, `m`, `h`, `n`, and `t`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  What are the dynamics of the voltage (variable `V`) resulting\n",
    "from this simulation?<br>\n",
    "HINT:  Plot `V` vs `t`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**   What are the dynamics of the gating variables (`m`, `h`, `n`)\n",
    "resulting from this simulation?<br>\n",
    "HINT:  Plot them!\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  What are the final values (after the 100 ms of simulation) of\n",
    "`V`, `m`, `h`, and `n`?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observation for Part 2\n",
    "At this value of input current (`I0=0`), the model dynamics\n",
    "approach a \"fixed point\", whose location we can identify as a point in four dimensional space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3:  At high input current (`I0`), examine the HH dynamics of a spike.\n",
    "  Let's now increase the input current to the HH model and get this model\n",
    "  to generate repeated spiking activity.  To do so, let's set,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "I0 = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now simulate this model,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "[V,m,h,n,t] = HH(I0,T0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  What happens to the dynamics?<br>\n",
    "\n",
    "HINT:  Plot V vs t.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  ### Observation for Part 3\n",
    "  You should have found that, at this value of input current, the model\n",
    "  **generates repeated spikes**.\n",
    "  \n",
    "  Let's now explore how the combined gates\n",
    "  and dynamics evolve.  To do so, let's start by focusing our plot on a\n",
    "  single spike.  As a first step, we'll make a new figure with a seperate subfigure to plot\n",
    "  the voltage,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x11866dd10>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure()\n",
    "subplot(211)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `subplot` command divides the figure into two rows, and one column, and tells Python we'll start in the first row. See Python Help for more details:\n",
    "\n",
    "`subplot??`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's plot the voltage, and choose the time axis to focus on a single spike,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(t,V,'k')\n",
    "xlim([42, 56])\n",
    "ylabel('V [mV]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  Okay, we've now plotted the voltage dynamics for a single spike (and\n",
    "  colored the curve black).  Let's now plot the three gating variables.\n",
    "  To do so, we'll move to the next subplot,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot(212);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(the next row in the figure).  Within this subplot, let's start by displaying the gating variable `m` over the same x-limits,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(t,m,'r', label='m')\n",
    "xlim([42, 56]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  Notice that, in the call to `plot` we included the input `label`. This will be useful when we create a legend ... <br><br>Within this subplot, we can also simultaneously show the gating\n",
    "  variables `h` and `n`,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(t,h,'b', label='h')\n",
    "plot(t,n,'g', label='n');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Label the x-axis,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlabel('Time [ms]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's add a legend to help us keep track of the different curves,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** Using the figure you created above, describe how the gates swing open and closed during a spike.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ASIDE: \n",
    "Here's a nice plotting trick, to link the x-axes of our two subfigures.  Linking the axes is useful so that, when we zoom or move one subfigure, the other subfigure will match the x-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "ax1 = subplot(211);                 # Define axis for 1st subplot,\n",
    "ax2 = subplot(212, sharex=ax1);     # ... and link axis of 2nd subplot to the 1st.\n",
    "ax1.plot(t,V,'k')                   # Plot the voltage in the first subplot,\n",
    "xlim([42, 56]);\n",
    "ax2.plot(t,m,'r', label='m')        # ... and the gating variables in the other subplot.\n",
    "ax2.plot(t,h,'b', label='h')\n",
    "ax2.plot(t,n,'g', label='n');\n",
    "xlabel('Time [ms]');\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, in the figure, you may use the pan/zoom tool to adjust the linked subplots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4:  At high input current (`I0`), describe the dynamics of the conductances.\n",
    "  In Part 3, we explored how the three gates `m`, `h`, and `n` evolve\n",
    "  during a spike.  By combining these terms, we can visualize how the\n",
    "  *conductances* evolve during a spike.  To do so, let's stick with the\n",
    "  simulation results we generated in Part 3, and focus our plot on a\n",
    "  single spike,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "ax1=subplot(311)                # Make a subplot,\n",
    "ax1.plot(t,V,'k')               #... and plot the voltage,\n",
    "xlim([42, 56])                  #... focused on a single spike,\n",
    "ylabel('V [mV]');               #... with y-axis labeled."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to plot the conductances, let's define three new variables,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "gNa0 = 120\n",
    "gNa  = gNa0*m**3*h              # Sodium conductance\n",
    "gK0  = 36\n",
    "gK   = gK0*n**4                 # Potassium conductance\n",
    "gL0  = 0.3\n",
    "gL   = gL0*ones(shape(gK))      # Leak conductance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** Where do these terms come from?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, let's plot these conductances,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax2 = subplot(312, sharex=ax1)  #Make a second subplot,\n",
    "ax2.plot(t,gNa,'m', label='gNa')#... and plot the sodium conductance,\n",
    "ax2.plot(t,gK, 'g', label='gK') #... and plot the potassium conductance,\n",
    "ax2.plot(t,gL, 'k', label='gL') #... and plot the leak conductance.\n",
    "xlim([42, 56])                  #... focused on a single spike,\n",
    "xlabel('Time [ms]')             #... label the x-axis.\n",
    "ylabel('mS/cm^2')               #... and label the y-axis.\n",
    "legend();                       #... make a legend."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** How do the conductances evolve during a spike?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 5:  At high input current (`I0`), describe the dynamics of the *currents*.\n",
    "  In Part 4, we explored how the three conductances (`gNa`, `gK`, `gL`) evolve\n",
    "  during a spike.  Let's now visualize how the *ionic currents* evolve\n",
    "  during a spike.  To do so, let's stick with the same settings used in\n",
    "  Part 4 and examine the same simulation result.  Again, we'll focus our plot\n",
    "  on a single spike.\n",
    "  \n",
    "  \n",
    "  Now, to plot the *current*, let's define the new variables,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "gNa0 = 120\n",
    "ENa  = 115\n",
    "INa  = gNa0*m**3*h*(ENa-(V+65))    # Sodium current.\n",
    "gK0  = 36\n",
    "EK   =-12\n",
    "IK   = gK0*n**4*(EK-(V+65))        # Potassium current.\n",
    "gL0  = 0.3\n",
    "EL   = 10.6;\n",
    "IL   = gL0*(EL-(V+65))             # Leak current.\n",
    "\n",
    "ax3=subplot(313, sharex=ax1)   # Make a third subplot,\n",
    "ax3.plot(t,INa,'m', label='INa')   #... and plot the sodium current,\n",
    "ax3.plot(t,IK, 'g', label='IK')    #... and plot the potassium current,\n",
    "ax3.plot(t,IL, 'k', label='IL')    #... and plot the leak current.\n",
    "xlim([42, 56])                 #... focus on a single spike,\n",
    "xlabel('Time [ms]')            #... label the x-axis.\n",
    "ylabel('mA/cm^2')              #... and label the y-axis.\n",
    "legend();                      #... make a legend."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** How do the conductances evolve during a spike?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** You may notice a small, transient decrease in the sodium current `INa` near 47 ms. What causes this?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"donate\"></a>\n",
    "## Donate\n",
    "If you enjoy Case-Studies-Python, and would like to share your enjoyment with us, sponsor our coffee consuption <a href=\"https://www.paypal.com/donate/?hosted_button_id=DL8P5ZGS9962U\">here</a>."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md:myst"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}