HH.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "title: Hodgkin-Huxley Neuron\n",
    "project:\n",
    "  type: website\n",
    "format:\n",
    "  html:\n",
    "    code-fold: true\n",
    "    code-tools: true\n",
    "jupyter: python 3\n",
    "number-sections: true\n",
    "filters:\n",
    "    - pyodide\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": [
    "#  Preliminaries \n",
    "\n",
    "Before beginning, let's load in the Python packages we'll need:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "```"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "import requests\n",
    "url = \"https://raw.githubusercontent.com/Mark-Kramer/BU-MA665-MA666/master/HH_functions.py\"\n",
    "response = requests.get(url)\n",
    "exec(response.text)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Hodgkin-Huxley (HH) equation code.\n",
    "\n",
    "To start, let's examine the code for the HH model. To do so, visit the raw code avaiable [here](https://raw.githubusercontent.com/Mark-Kramer/BU-MA665-MA666/master/HH_functions.py).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\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",
    ":::"
   ]
  },
  {
   "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": [
    "# 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": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "I0 = 0\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and let's simulate the model for 100 ms,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "T0 = 100\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've now defined both inputs to the `HH` function, and can execute it, as follows,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "[V,m,h,n,t]=HH(I0,T0)\n",
    "```"
   ]
  },
  {
   "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": [
    "::: {.callout}\n",
    "**Q:**  What are the dynamics of the voltage (variable `V`) resulting\n",
    "from this simulation?<br>\n",
    "HINT:  Plot `V` vs `t`.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**   What are the dynamics of the gating variables (`m`, `h`, `n`)\n",
    "resulting from this simulation?<br>\n",
    "HINT:  Plot them!\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  What are the final values (after the 100 ms of simulation) of\n",
    "`V`, `m`, `h`, and `n`?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 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": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "I0 = 10\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now simulate this model,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "[V,m,h,n,t] = HH(I0,T0)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  What happens to the dynamics?<br>\n",
    "HINT:  Plot V vs t.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observations:** You should have found that, at this value of input current, the model **generates repeated spikes**.\n",
    "\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, let's plot the voltage, and choose the time axis to focus on a single spike,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(t,V,'k')\n",
    "plt.xlim([35, 55])\n",
    "plt.ylabel('V [mV]')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've now plotted the voltage dynamics for a single spike (and\n",
    "  colored the curve black).\n",
    "  \n",
    "Let's now plot the three gating variables. To do so, we'll make another plot. \n",
    "\n",
    "  Let's start by displaying the gating variable `m` over the same x-limits,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(t,m,'r', label='m')\n",
    "plt.xlim([35, 55])\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "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`, with the x-axis labeled.\n",
    "\n",
    "Let's also add a legend to help us keep track of the different curves:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(t,m,'r', label='m')\n",
    "plt.plot(t,h,'b', label='h')\n",
    "plt.plot(t,n,'g', label='n')\n",
    "plt.xlim([35, 55])\n",
    "plt.xlabel('Time [ms]');\n",
    "plt.legend();\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** Using the figure you created above, describe how the gates swing open and closed during a spike.\n",
    ":::"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "ax1 = plt.subplot(211);                 # Define axis for 1st subplot,\n",
    "ax2 = plt.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",
    "plt.xlim([35, 55])\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",
    "plt.xlabel('Time [ms]')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "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": [
    "# At high input current (`I0`), describe the dynamics of the conductances.\n",
    "  We have so far 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 the previous section, and focus our plot on a\n",
    "  single spike,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(t,V,'k')                   #Plot the voltage,\n",
    "plt.xlim([35, 55])\n",
    "plt.ylabel('V [mV]')                #... with y-axis labeled.\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to plot the conductances, let's define three new variables,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "gNa0 = 120\n",
    "# gNa  = ???                    #Sodium conductance\n",
    "gK0  = 36\n",
    "# gK   = ???                    #Potassium conductance\n",
    "gL0  = 0.3\n",
    "# gL   = ???                    #Leak conductance\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** Where do these terms come from?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, let's plot these conductances,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(t,gNa,'m', label='gNa')#... and plot the sodium conductance,\n",
    "plt.plot(t,gK, 'g', label='gK') #... and plot the potassium conductance,\n",
    "plt.plot(t,gL, 'k', label='gL') #... and plot the leak conductance.\n",
    "plt.xlim([35, 55])\n",
    "plt.xlabel('Time [ms]')             #... label the x-axis.\n",
    "plt.ylabel('mS/cm^2')               #... and label the y-axis.\n",
    "plt.legend()                        #... make a legend.\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** How do the conductances evolve during a spike?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At high input current (`I0`), describe the dynamics of the *currents*.\n",
    "  In the previous section, we explored how the three conductances (`gNa`, `gK`, `gL`) evolve\n",
    "  during a spike.\n",
    "  \n",
    "  Let's now visualize how the *ionic currents* evolve\n",
    "  during a spike.\n",
    "  \n",
    "  To do so, let's stick with the same settings used in the previous section 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": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "gNa0 = 120\n",
    "ENa  = 115\n",
    "# INa  = ??? Sodium current.\n",
    "gK0  = 36\n",
    "EK   =-12\n",
    "# IK   = ??? Potassium current.\n",
    "gL0  = 0.3\n",
    "EL   = 10.6;\n",
    "# IL   = ??? Leak current.\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(t,INa,'m', label='INa')   #... and plot the sodium current,\n",
    "plt.plot(t,IK, 'g', label='IK')    #... and plot the potassium current,\n",
    "plt.plot(t,IL, 'k', label='IL')    #... and plot the leak current.\n",
    "plt.xlim([35, 55])\n",
    "plt.xlabel('Time [ms]')            #... label the x-axis.\n",
    "plt.ylabel('mA/cm^2')              #... and label the y-axis.\n",
    "plt.legend()                      #... make a legend.\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** How do the conductances evolve during a spike?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** You may notice a small, transient decrease in the sodium current `INa` near 47 ms. What causes this?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion\n",
    "\n",
    "1. How is the HH model different from / similar to the IF and LIF models?\n",
    "2. What are the ion species simulated in the HH model?\n",
    "3. Which ion species is fastest / slowest?\n",
    "4. How many variables are in the HH model? Define each.\n",
    "5. Sketch the model equations, in schematic form, to capture the “essence” of the behavior (don’t worry about the detailed values for parameters).\n",
    "6. Which gating variables are depolarization activated / depolarization inactivated? \n",
    "7. Sketch the steady-state curves for the gating variables.\n",
    "8. Consider the current: `I_K=g_K n^4  (E_K-V)`\n",
    "    1. Define each term.\n",
    "    2. How does this current behave for different values of V? \n",
    "9. Consider the current: `I_Na=g_Na m^3 h (E_Na-V)`\n",
    "    1. Define each term.\n",
    "    2. How does this current behave for different values of V?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Challenges\n",
    "\n",
    "1. Describe the dynamics during an action potential in the HH model. How does the voltage change? How do the gates open and close? How do the ions flow?\n",
    "5. Determine how the firing rate of the HH model varies with input current I. Plot the firing rate versus I (i.e., plot the “f-I curve”).\n",
    "6. How does the firing rate of the HH model change as you increase the potassium conductance? Provide a \"simulation\" explanation and a \"physical\" explanation.\n",
    "7. How does the firing rate of the HH model change as you increase sodium conductance? Provide a \"simulation\" explanation and a \"physical\" explanation.\n",
    "8. **BONUS:** Use the forward and backward rate functions ($\\alpha$ and $\\beta$) in the HH model (see [code](https://raw.githubusercontent.com/Mark-Kramer/BU-MA665-MA666/master/HH_functions.py)) to plot the steady state function and time constant for each channel ($m$, $h$, $n$) versus voltage. In other words, create Figure 1 in the [Hodgkin-Huxley Cheat Sheet](https://mark-kramer.github.io/BU-MA665-MA666/Readings/Hodgkin-Huxley_Cheat_Sheet.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md:myst"
  },
  "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",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}