IF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "title: \"Integrate and Fire 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 integrate and fire (I&F) neuron model.  We'll investigate, in particular, how the spiking activity varies as we adjust the input current $I$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Preliminaries\n",
    "Before beginning, let's load in the Python package 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": [
    "#  Numerical solutions - Introduction\n",
    "\n",
    " \n",
    "How do we compute a numerical solution to the integrate and fire model?\n",
    "  The basic idea is to rearrange the differential equation to get $V(t+1)$ on\n",
    "  the left hand side, and $V(t)$ on the right hand side.  Then, if we know\n",
    "  what's happening at time $t$, we can solve for what's happening at time $t+1$.\n",
    "\n",
    "For example, consider the differential equation:\n",
    "\n",
    "$$\n",
    "  \\dfrac{dV}{dt} = \\dfrac{I}{C}\n",
    "$$\n",
    "\n",
    "In words, we can think of:\n",
    "\n",
    "$dV$ as the \"change in voltage V\",\n",
    "\n",
    "$dt$ as the \"change in time t\".\n",
    "\n",
    "Let's consider the case that we record the voltage $V$ in discrete time steps. So\n",
    "we observe:\n",
    "\n",
    "$V[0], V[1], V[2], \\ldots$\n",
    "\n",
    "at times:\n",
    "\n",
    "$dt, \\, 2*dt, \\, 3*dt, \\ldots$\n",
    "\n",
    "where $dt$ is the time between our samples of $V$.\n",
    "\n",
    "We can now write the \"change in voltage V\" as:\n",
    "\n",
    "$$\n",
    "  dV = V(t+1) - V(t)\n",
    "$$\n",
    "\n",
    "Notice that the change in voltage is the difference in V between two\n",
    "sequential time samples.  Now, let's rewrite $\\dfrac{dV}{dt}$ as,\n",
    "\n",
    "$$\n",
    "  \\dfrac{dV}{dt} = \\dfrac{ V(t+1) - V(t) }{ dt } \n",
    "$$\n",
    "\n",
    "where we've replaced $dV$.  Now, let's substitute this expression into the equation at the top of this file:\n",
    "\n",
    "$$\n",
    "  \\dfrac{ V(t+1) - V(t) }{ dt } = \\dfrac{I}{C}.  \n",
    "$$\n",
    "\n",
    "Solving this equation for $V(t+1)$ you'll find  that:\n",
    "\n",
    "$$\n",
    "  V(t+1) = V(t) + dt*(I/C)\n",
    "$$\n",
    "\n",
    "Notice that, in this expression, we use our current value of the voltage V(t) and the model (I/C) to determine the next value of the voltage V(t+1).\n",
    "\n",
    "Now, let's program this equation in Python.  First, let's set the values\n",
    "for the parameters $I$ and $C$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "C=1.0\n",
    "I=1.0\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  We also need to set the value for $dt$.  This defines the time step for our\n",
    "  model.  We must choose it small enough so that we don't miss anything\n",
    "  interesting.  We'll choose:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "dt=0.01\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assume the units of time are seconds.  So, we step forward in time by $0.01$ s.\n",
    "\n",
    "The right hand side of our equation is nearly defined, but we're still missing one thing, $V(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  What value do we assign to $V(t)$?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So here's an easier question:  what *initial* value do we assign to $V(t)$?\n",
    " \n",
    " To start, we'll create an array of zeros to hold our results for $V$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "V = np.zeros([1000,1])\n",
    "np.shape(V)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This array `V` consists of 1000 rows and 1 column. We can think of each row entry as corresponding to a discrete step in time. Our goal is to fill-in the values of `V` (i.e., step forward in time), in a way consistent with our model.\n",
    "\n",
    "Let's choose an initial value for `V` of 0.2, which in our simple model we'll assume represents the rest state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "V[0]=0.2\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  Given the initial state `V[0]=0.2`, calculate `V[1]`.  Then calcualte `V[2]`.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  After the two calculations above, we've moved forward two time steps into\n",
    "  the future, from $t=0$ s to $t=0.01$ s, and then from $t=0.01$ s to $t=0.02$ s.  But what\n",
    "  if we want to know $V$ at $t=10$ s?  Then, this iteration-by-hand procedure becomes\n",
    "  much too boring and error-prone.  So, what do we do?  Let's make the\n",
    "  computer do it ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical solutions - Implementation\n",
    "\n",
    "Let's computerize this iteration-by-hand procedure to find `V[999]`.\n",
    "\n",
    "To do so, we'll use a [for-loop](https://wiki.python.org/moin/ForLoop).\n",
    "\n",
    "Here's what it looks like:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "for k in range(1,999):\n",
    "    V[k+1] = V[k] + dt*(I/C)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  Does this loop make sense?  Describe what's happening here.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  Why does the `range` command end at `999`?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute this for-loop and examine the results in vector `V`.  To do so, let's plot `V`:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(V)\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  What happens to the voltage after 1000 steps?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  This plot is informative, but not great.  Really, we'd like to plot the\n",
    "  voltage as a function of *time*, not steps or indices.  To do so, we\n",
    "  need to define a time axis:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "t = np.arange(0,len(V))*dt\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  What's happening in the command above?  Does it make sense?  (If not, trying printing or plotting `t`.)\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, with *time* defined, let's redo the plot of the voltage with the axes labeled appropriately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "plt.figure()\n",
    "plt.plot(t,V)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('V')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's put it all together . . ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I&F CODE (version 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  So far, we constructed parts of the I&F model in bits-and-pieces.\n",
    "  Let's now collect all the code, compute a numerical solution to\n",
    "  the I&F model, and plot the results (with appropriate axes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "I   = 1                           #Set the parameter I.\n",
    "C   = 1                           #Set the parameter C.\n",
    "dt  = 0.01                        #Set the timestep.\n",
    "V   = np.zeros([1000,1])             #Initialize V.\n",
    "V[0]= 0.2;                        #Set the initial value of V.\n",
    "\n",
    "for k in range(1,999):            #March forward in time,\n",
    "    V[k+1] = V[k] + dt*(I/C)      #... updating V along the way.\n",
    "\n",
    "t = np.arange(0,len(V))*dt        #Define the time axis.\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(t,V)                         #Plot the results.\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Voltage [mV]')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  Adjust the parameter `I`.  What happens to `V` if `I=0`?  Can you set `I` so that `V` > 20 within 10 s?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voltage threshold\n",
    "\n",
    "  Notice, our model is missing something important:  **the reset**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  Without the reset, how does the voltage behave as $t\\rightarrow\\infty$ (if $I>0$)?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Now, let's update our model to include the reset.  To do so, we'll add two things to our code.\n",
    "  \n",
    "  - First, we'll define the voltage threshold `Vth`, and\n",
    "  reset voltage `Vreset`.\n",
    "  - Second, we'll check to see if `V` exceeds\n",
    "  `Vth` using an [if-statement](https://docs.python.org/3/tutorial/controlflow.html);  if it does, then we'll set `V` equal to\n",
    "  `Vreset`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:**  How will you update the code below to include the reset?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "Vth = 1;        #Define the voltage threshold.\n",
    "Vreset = 0;     #Define the reset voltage.\n",
    "\n",
    "for k in range(1,999):            #March forward in time,\n",
    "    V[k+1] = V[k] + dt*(I/C)      #Update the voltage,\n",
    "    ### ADD SOMETHING HERE??? --------------------------\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I&F CODE (version 2)\n",
    "\n",
    "Now, let's put it all together to make a complete I&F model (with a thershold and reset), simulate it, and plot the result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{pyodide-python}\n",
    "### ADD YOUR CODE!\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** Adjust the parameter `I`.  What happens to `V` if `I=10`?  If `I=100`?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** Adjust the parameter `C`.  What happens to `V` if `C=0.1`?  If `C=10`?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "::: {.callout}\n",
    "**Q:** What is \"spiking\" in this I&F model?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion\n",
    "\n",
    "1. Describe the main components of the **IF** model. Draw it (in some way).\n",
    "   \n",
    "3. Describe the main components of the **LIF** model. Draw it (in some way).\n",
    "4. Describe the differences and similarities between the IF and LIF models.\n",
    "5. The IF model is meant to mimic a neurons activity. What is realistic about the IF model? What is unrealistic?\n",
    "6. Describe the roles of the IF model parameters `Vreset` and `Vthreshold`.\n",
    "7. Consider the IF model. Sketch voltage (V) versus time (t) for a small input current, for a large input current.\n",
    "    1. How does an increase in capacitance (C) impact the dynamics?\n",
    "    2. Can you interpret this physically?\n",
    "8. Plot the f-I curve for the IF model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Challenges\n",
    "\n",
    "1. Consider the **LIF** model. Plot (in Python) voltage (V) versus time (t) for a small input current, for a large input current.\n",
    "\n",
    "3. Perform numerical simulations to determine how doubling the capacitance (C) impacts the firing rate of the LIF model. Based on your numerical simulations, does the firing rate increase or decrease? Provide a physical explanation for your results (i.e., how does the charge move around the circuit to impact the firing rate?) HINT: Capacitance represents the *capacity* of a cell to hold charge.\n",
    "4. Perform numerical simulations to determine how doubling the resistance (R) impacts the firing rate of the LIF model. Based on your numerical simulations, does the firing rate increase or decrease? Provide a physical explanation for your results (i.e., how does the charge move around the circuit to impact the firing rate?) HINT: Higher resistance makes it harder for charge to flow in or out of the cell.\n",
    "8. Determine how the firing rate of the LIF model varies with input I. Plot the firing rate versus I (i.e., plot the “f-I curve”).\n",
    "9. **ADVANCED**. Simulate the [Izhikevich neuron](/Readings/Izhikevich_2003.pdf). Provide examples of two different types of spiking activity."
   ]
  },
  {
   "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
}