cs291i-izzy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0. Intro to Jupyter Notebook\n",
    "\n",
    "Jupyter Notebook is an open-source web application that you can use to create and share documents that contain live code, equations, visualizations, and text.\n",
    "To get started, go to `Help > User Interface Tour`. Also check out `Help > Keyboard Shortcuts`.\n",
    "\n",
    "Jupyter Notebook has two different keyboard input modes.\n",
    "- `Edit mode` allows you to type code or text into a cell and is indicated by a green cell border. Click into a cell to reach Edit mode.\n",
    "- `Command mode` binds the keyboard to notebook level commands and is indicated by a grey cell border with a blue left margin. Hit Escape to reach Command mode.\n",
    "\n",
    "Some useful keyboard shortcuts in command mode:\n",
    "\n",
    "- `A`: insert cell above\n",
    "- `B`: insert cell below\n",
    "\n",
    "Some useful keyboard shortcuts in edit mode:\n",
    "\n",
    "- `Tab`: Code completion or indent\n",
    "- `Shift+Tab`: tooltip\n",
    "- `Shift+Tab+Tab`: expand tooltip\n",
    "- `Ctrl+Enter`: Execute a cell\n",
    "- `Shift+Enter`: Execute a cell and select the cell below it\n",
    "\n",
    "Also: `Kernel > Interrupt`, `Kernel > Restart & Clear Output`, `Kernel > Restart & Run All`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import NumPy for math:\n",
    "import numpy as np\n",
    "\n",
    "# Import Matplotlib for plotting:\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Izhikevich spiking neuron\n",
    "\n",
    "Below you will implement an Izhikevich spiking neuron and see how it reacts to injected current.\n",
    "\n",
    "**TODO**: Throughout the notebook, fill in the empty lines marked with **...**. Then execute the cell to see the resulting output.\n",
    "\n",
    "\n",
    "## 1.1 Generating stimuli\n",
    "\n",
    "Thus, the first step is to generate a current stimulus.\n",
    "We want a function that returns a current pulse just like in the Izhikevich paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_stim(Imax, tonset, tmax, dt=1):\n",
    "    \"\"\"A current pulse starting at ``tonset`` with amplitude ``Imax``\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    Imax : float\n",
    "        the amplitude of the current pulse (micro amps)\n",
    "    tonset : float\n",
    "        the time in milliseconds at which the pulse starts\n",
    "        (amplitude is 0 up to that point)\n",
    "    tmax : float\n",
    "        the end point of the current pulse\n",
    "    dt : float, optional, default: 1\n",
    "        the time step to use (in milliseconds)\n",
    "            \n",
    "    Returns\n",
    "    -------\n",
    "    A dictionary with field 't' that contains the time steps and\n",
    "    'I' that contains the current amplitude at each time step.\n",
    "    \n",
    "    Examples\n",
    "    --------\n",
    "    # A pulse in the range t=[0, 100) with amp=10 starting at 15ms\n",
    "    >>> stim = get_stim(10, 15, 100, dt=1)\n",
    "\n",
    "    \"\"\"\n",
    "    # Generate the time array in [0, tmax] with time step dt:\n",
    "    tstim = np.arange(0, tmax, dt)\n",
    "    # Allocate an array for the current values:\n",
    "    Istim = np.zeros_like(tstim)\n",
    "    # TODO Set current to Imax after tonset:\n",
    "    Istim[tstim >= tonset] = ...\n",
    "    # TODO Set up a dictionary according to the documentation above:\n",
    "    stim = {...}\n",
    "    return stim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO Generate a pulse in the range t=[0, 1000) with amp=14 starting\n",
    "# at 100ms, with a 0.5 ms time step:\n",
    "stim = ...\n",
    "\n",
    "plt.plot(stim['t'], stim['I'], linewidth=3)\n",
    "plt.xlabel('time (ms)')\n",
    "plt.ylabel('current (uA)')\n",
    "plt.title('Current injection');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 The Izhikevich model\n",
    "\n",
    "The Izhikevich neuron is a dynamical systems model that can be described by a two-dimensional system of ordinary differential equations:\n",
    "\n",
    "$$\\begin{align*} \\frac{dv}{dt} = &  ~ 0.04v^2 + 5v + 140 - u + I & \\text{(1)} \\\\ \\frac{du}{dt} = & ~ a (bv - u) & \\text{(2)} \\end{align*}$$\n",
    "\n",
    "Here, (1) describes the membrane potential $v$ for a given current $I$, where $I$ is the sum of all synaptic and external currents; that is, $I = i_{syn} + I_{ext}$ (see 3.2 Synapses). (2) describes a recovery variable $u$; the parameter $a$ is the rate constant of the recovery variable, and the parameter $b$ describes the sensitivity of the recovery variable to the subthreshold fluctuations of the membrane potential. All parameters in (1) and (2) are dimensionless; however, the right-hand side of (1) is in a form such that the membrane potential $v$ has mV scale and the time $t$ has ms scale. $a, b, c, d$ are open parameters that have different values for different neuron types.\n",
    "\n",
    "In contrast to other simple models such as the leaky integrate-and-fire (LIF) neuron, the Izhikevich neuron is able to generate the upstroke of the spike itself. Thus the voltage reset occurs not at the threshold, but at the peak ( $v_{cutoff}=+30$) of the spike. The action potential downstroke is modeled using an instantaneous reset of the membrane potential whenever $v$ reaches the spike cutoff, plus a stepping of the recovery variable:\n",
    "\n",
    "$$\\begin{align*} v(v>30) = & ~ c & \\text{(3)}\\\\ u(v>30) = & ~ u + d & \\text{(4)}\\end{align*}$$\n",
    "\n",
    "The inclusion of $u$ in the model allows for the simulation of typical spike patterns observed in biological neurons. The four parameters $a, b, c, d$ can be set to simulate different types of neurons. For example, regular spiking (RS) neurons (class 1 excitable) have $a=0.02, b=0.2, c=-65, d=8$. Fast spiking (FS) neurons (class 2 excitable) have $a=0.1, b=0.2, c=-65, d=2)$. For more information on different neuron types, see Izhikevich (2003, 2004).\n",
    "\n",
    "In the below code, we want to encapsulate the Izhikevich model in a function that takes as input the stimulus dictionary created above (`stim`), as well as the neuron parameters `a`, `b`, `c`, and `d`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_izhikevich(stim, a=0.02, b=0.2, c=-65, d=6, V0=-70):\n",
    "    \"\"\"Izhikevich neuron model\n",
    "        \n",
    "    Parameters\n",
    "    ----------\n",
    "    stim : dict\n",
    "        A stimulus dictionary with fields 't' and 'I'\n",
    "    a : float\n",
    "        rate constant of the recovery variable\n",
    "    b : float\n",
    "        sensitivity of the recovery variable to the subthreshold\n",
    "        fluctuations of the membrane potential\n",
    "    c : float\n",
    "        reset potential after spike\n",
    "    d : float\n",
    "        step of the recovery variable after spike\n",
    "    V0 : float\n",
    "        the initial value (at t=0) for the membrane potential\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    VV, uu : np.array, np.array\n",
    "        An array containing the membrane potential at each time step (``VV``)\n",
    "        as well as an array containing the value of the recovery variable at\n",
    "        each time step (``u``).\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    # A regular-spiking neuron:\n",
    "    >>> V, u = run_izhikevich(stim, a=0.02, b=0.2, c=-65, d=6)\n",
    "    \n",
    "    # A phasic bursting neuron:\n",
    "    >>> V, u = run_izhikevich(stim, a=0.02, b=0.25,  c=-55,  d=0.05, V0=-64)\n",
    "    \n",
    "    \"\"\"\n",
    "    # TODO Initialize model parameters:\n",
    "    V = ...\n",
    "    u = b * V\n",
    "    \n",
    "    # Initialize empty arrays for the output variables:\n",
    "    VV = np.zeros_like(stim['t'])\n",
    "    uu = np.zeros_like(stim['t'])\n",
    "    \n",
    "    # Infer the time step from the 't' array:\n",
    "    tau = np.diff(stim['t'])[0]\n",
    "    \n",
    "    # The main loop - this has to be executed at every time step.\n",
    "    # The enumerate statement returns two things: the index and value of\n",
    "    # each element in stim['I']:\n",
    "    for t, I in enumerate(stim['I']):\n",
    "        # TODO Update the two state variables:\n",
    "        V += tau * ...\n",
    "        u += tau * ...\n",
    "        \n",
    "        # Event handling:\n",
    "        if V > 30:\n",
    "            # TODO Spike!\n",
    "            VV[t] = 30\n",
    "            V = ...\n",
    "            u += ...\n",
    "        else:\n",
    "            # No spike...\n",
    "            VV[t] = V\n",
    "        uu[t] = u\n",
    "    return VV, uu\n",
    "\n",
    "def plot_izhikevich(stim, V):\n",
    "    plt.plot(stim['t'], V, label='membrane potential (mV)', linewidth=3)\n",
    "    plt.plot(stim['t'], stim['I'] * 10 / np.max(stim['I']) - 120, label='current (uA)')\n",
    "    plt.ylim(-150, 50)\n",
    "    plt.yticks(np.linspace(-90, 30, 5))\n",
    "    plt.xlabel('time (ms)')\n",
    "    plt.legend(loc='lower right');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Tonic Spiking\n",
    "\n",
    "If you did everything right above, you should see a tonic-spiking neuron fire ~9 spikes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iext = get_stim(14, 10, 200, dt=0.1)\n",
    "V, u = run_izhikevich(Iext, a=0.02, b=0.2, c=-65, d=6)\n",
    "plot_izhikevich(Iext, V)\n",
    "plt.title('tonic spiking');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 Phasic bursting\n",
    "\n",
    "A phasic-bursting neuron spikes totally differently.\n",
    "You can simulate it with `a=0.02`, `b=0.25`, `c=-55`, `d=0.05`, and `V0=-64`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iext = get_stim(0.6, 10, 200, dt=0.1)\n",
    "# TODO\n",
    "V, u = ...\n",
    "plot_izhikevich(Iext, V)\n",
    "plt.title('phasic bursting');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.5 Calculating Firing Rates\n",
    "\n",
    "Real neurons spike, but artificial neurons often do not.\n",
    "To establish a correspondence between spikes and firing rate, we need to count spikes in a long-enough window:\n",
    "\n",
    "$$\\textrm{rate} = \\frac{\\textrm{number of spikes}}{\\textrm{time window}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_num_spikes(V):\n",
    "    \"\"\"Counts spikes\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    V : np.array\n",
    "        An array containing the voltage trace of an Izhikevich neuron\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    spikes : int\n",
    "        Number of spikes\n",
    "    \"\"\"\n",
    "    # TODO\n",
    "    n_spikes = ...\n",
    "    return n_spikes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's put our function to the test, and create a plot showing the number of spikes as a function of the injected current:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a range of injected currents:\n",
    "in_amp = np.arange(-10, 50, 5)\n",
    "# Allocate an empty list for the number of spikes:\n",
    "out_spikes = []\n",
    "\n",
    "# For every element in the injected current array:\n",
    "for I in in_amp:\n",
    "    # Generate the stimulus:\n",
    "    Iext = get_stim(I, 10, 1000, dt=0.1)\n",
    "    # Run the Izhikevich model:\n",
    "    V, u = run_izhikevich(Iext, a=0.02, b=0.2, c=-65, d=6)\n",
    "    # TODO Calculate the number of spikes\n",
    "    n_spikes = ...\n",
    "    # Append the number of spikes to a list:\n",
    "    out_spikes.append(n_spikes)\n",
    "    \n",
    "# Plot the result:\n",
    "plt.plot(in_amp, out_spikes, linewidth=3)\n",
    "plt.xlabel('input current (uA)')\n",
    "plt.ylabel('number of spikes')\n",
    "plt.title('tonic spiking')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is your ReLU!\n",
    "\n",
    "What about a tonic-bursting neuron?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO repeat for tonic bursting\n",
    "out_spikes = ...\n",
    "...\n",
    "...\n",
    "\n",
    "# Plot the result:\n",
    "plt.plot(in_amp, out_spikes, linewidth=3)\n",
    "plt.xlabel('input current (uA)')\n",
    "plt.ylabel('number of spikes')\n",
    "plt.title('tonic bursting');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the two compare?\n",
    "\n",
    "One more time: What about a phasic-bursting neuron? `a=0.02, b=0.25,  c=-55,  d=0.05, V0=-64`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO repeat for phasic bursting\n",
    "in_amp = np.arange(-10, 50, 5)\n",
    "out_spikes = ...\n",
    "\n",
    "...\n",
    "...\n",
    "\n",
    "# Plot the result:\n",
    "plt.plot(in_amp, out_spikes, linewidth=3)\n",
    "plt.xlabel('input current (uA)')\n",
    "plt.ylabel('number of spikes')\n",
    "plt.title('phasic bursting');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does this one compare?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Electrical stimulation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Idealized point source\n",
    "\n",
    "The electric field $V(r)$ of a point source in an infinite homogeneous medium is given by:\n",
    "\n",
    "$$V(r) = \\frac{\\sigma I}{4 \\pi r},$$\n",
    "\n",
    "where $\\sigma$ is the resistivity of the extracellular solution (typically Ames medium), $I$ is the amplitude of the constant current pulse, and $r$ is the distance from the stimulating electrode to the point at which the voltage is being computed.\n",
    "\n",
    "We can write a function for that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def electric_potential(I, r, sigma=110):\n",
    "    \"\"\"Electric potential of a point source\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    I : float\n",
    "        amplitude of the constant current pulse\n",
    "    r : float\n",
    "        distance from the stimulation electrode to thte point at which\n",
    "        the voltage is being computed\n",
    "    sigma : float\n",
    "        resistivity of the extracellular solution (Ohm cm)\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    V : float\n",
    "        the electric potential as a function of r\n",
    "    \"\"\"\n",
    "    if r == 0:\n",
    "        # Special case: avoid division by zero\n",
    "        V = sigma * I\n",
    "    else:\n",
    "        # TODO implement the equation above\n",
    "        ...\n",
    "    return V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A list of radii to use as distance:\n",
    "radii = np.arange(0, 100, 5)\n",
    "# Allocate an empty list for the number of spikes:\n",
    "n_spikes = []\n",
    "# Assume the current pulse amplitude is 100uA (at the electrode):\n",
    "I = 100\n",
    "# Assume a constant resistance so we can convert Vel to Iext:\n",
    "R = 1\n",
    "for r in radii:\n",
    "    # TODO Calculate the electric potential at radius r using your function:\n",
    "    Vel = ...\n",
    "    # Generate the stimulus, convert Vel to a current\n",
    "    Iext = get_stim(Vel / R, 10, 1000, dt=0.1)\n",
    "    # TODO Run the Izhikevich model:\n",
    "    V, u = ...\n",
    "    # TODO Calculate the number of spikes and append to list:\n",
    "    n_spikes.append(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the result:\n",
    "plt.plot(radii, n_spikes, linewidth=3)\n",
    "plt.xlabel('distance from the electrode (um)')\n",
    "plt.ylabel('number of spikes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 The activation function\n",
    "\n",
    "Rattay describes the regions where electrical stimulation leads to spiking with an \"activation function\".\n",
    "He found that the activation function is proportional to the second spatial derivative of the electric potential.\n",
    "\n",
    "We can approximate the second derivative with finite differences, using NumPy's `gradient` function. Consider this Gaussian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 150\n",
    "r = np.arange(-1000, 1000)\n",
    "g = np.exp(-r ** 2 / (2 * sigma ** 2))\n",
    "plt.plot(r, g, linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the first derivative is given as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gprime = np.gradient(g)\n",
    "plt.plot(r, gprime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second derivative looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gprime2 = np.gradient(np.gradient(g))\n",
    "plt.plot(r, gprime2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sign of the current amplitude is usually subject to convention.\n",
    "For example, \"excitatory\" currents usually have a negative sign.\n",
    "So the region that leads to activation of neurons is in the center, surrounded by two flankers causing inhibition.\n",
    "\n",
    "We can put all this in a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def activation_function(Vel):\n",
    "    \"\"\"Rattay's activation function\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    Vel : np.array\n",
    "        Array of electric potentials\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    act : np.array\n",
    "        Activation function\n",
    "    \"\"\"\n",
    "    # TODO implement the 2nd derivative\n",
    "    act = ...\n",
    "    return act"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assume the current pulse amplitude is 100uA (at the electrode):\n",
    "I = 100\n",
    "# Assume a constant resistance so we can convert Vel to Iext:\n",
    "R = 1\n",
    "\n",
    "# All electric potentials:\n",
    "Vel = np.array([electric_potential(I, r) for r in radii])\n",
    "# Rattay's activation function:\n",
    "act = activation_function(Vel)\n",
    "\n",
    "# Allocate an empty list for the number of spikes:\n",
    "n_spikes_rattay = []\n",
    "for a in act:\n",
    "    # TODO Generate the stimulus using Rattay's activation function as\n",
    "    # injected current instead of Vel/R:\n",
    "    Iext = get_stim(..., 10, 1000, dt=0.1)\n",
    "    # Run the Izhikevich model:\n",
    "    V, u = ...\n",
    "    # TODO Calculate the number of spikes and append to list:\n",
    "    n_spikes_rattay.append(...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the result to the earlier one:\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(15, 5), sharey=True)\n",
    "\n",
    "ax[0].plot(radii, n_spikes, linewidth=3)\n",
    "ax[0].set_xlabel('distance from the electrode (um)')\n",
    "ax[0].set_ylabel('number of spikes')\n",
    "ax[0].set_title('Spikes based on Vel')\n",
    "\n",
    "ax[1].plot(radii, n_spikes_rattay, linewidth=3)\n",
    "ax[1].set_xlabel('distance from the electrode (um)')\n",
    "ax[1].set_ylabel('number of spikes')\n",
    "ax[1].set_title('Spikes based on activation function')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}