MainenSejnowski1995.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TcHOdlymE03a"
   },
   "source": [
    "The goal here is to create a \"raster plot\" that shows the reproducibility of a spike train to different repetitions of different stimulus: a DC current versus noise. \n",
    "\n",
    "In particular, we will attempt to replicate Figure 1 of [Mainen & Sejnowski (1995)](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.299.8560&rep=rep1&type=pdf) and use computational Neuroscience as a tool to better understanfd neurophysiological observations. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7TpiNqRaE03g"
   },
   "source": [
    "# Mainen & Sejnowski, 1995"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7TpiNqRaE03g"
   },
   "source": [
    "## context"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7TpiNqRaE03g"
   },
   "source": [
    "The goal of this first task is to create a \"raster plot\" that shows the reproducibility of a spike train with repetitions of the same stimulus, as in this work in the [rodent retina](https://laurentperrinet.github.io/2019-04-03_a_course_on_vision_and_modelization/#/1/3) or in the [cat cortex (V1)](https://laurentperrinet.github.io/2019-04-03_a_course_on_vision_and_modelization/#/1/6)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7TpiNqRaE03g"
   },
   "source": [
    "## figure 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7TpiNqRaE03g"
   },
   "source": [
    "Here, we will attempt to replicate Figure 1 of [Mainen & Sejnowski (1995)](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.299.8560&rep=rep1&type=pdf):\n",
    "\n",
    "![Mainen Sejnowski 1995](http://i.stack.imgur.com/ixnrz.png \"figure 1\")\n",
    "\n",
    "\n",
    "QUESTION: write a quick summary of the paper (max 5 lines) and why this result is *a priori* surprising."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "py_kpyj7E03h"
   },
   "source": [
    "# representation of time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "py_kpyj7E03h"
   },
   "source": [
    "## a gentle introduction to python and numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rfl9SgYYE03h"
   },
   "outputs": [],
   "source": [
    "dt = .5 # size of the time discretization step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# help(np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.lookfor('binary representation') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.lookfor('evenly spaced values')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rfl9SgYYE03h"
   },
   "outputs": [],
   "source": [
    "time = np.linspace(-100, 1000, int(1100/dt))\n",
    "# time = np.arange(-100, 1000, dt) # alternative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "xNFNW8J4aUfn",
    "outputId": "0c27720e-96be-46c0-88e9-23efde0c36bf"
   },
   "outputs": [],
   "source": [
    "time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9ZFnkuXkfE2P",
    "outputId": "1c2ca12a-b5a6-4129-f5b0-64a653409db5"
   },
   "outputs": [],
   "source": [
    "time.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "_laLJ91RZp7d",
    "outputId": "85372dff-2f76-445d-ed24-6fa133c79a1a"
   },
   "outputs": [],
   "source": [
    "time[0], time[1], time[-2], time[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "-tJrnSs2fKrF",
    "outputId": "6782d757-a28e-47cd-ba30-fdce177676c4"
   },
   "outputs": [],
   "source": [
    "time[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lBIsJF7VE03i"
   },
   "source": [
    "Creation of a DC current (version one):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "id": "R0SuQceOE03i",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "start = 0\n",
    "end = 900\n",
    "value = 150\n",
    "\n",
    "def Inp(time=time, start=start, end=end, value=value):\n",
    "    x=[]\n",
    "    for t in range(len(time)):\n",
    "        if start < time[t] < end :\n",
    "            x.append(value)\n",
    "        else:\n",
    "            x.append(0)\n",
    "    return x\n",
    "\n",
    "I = Inp(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "LSIctw7qE03j",
    "outputId": "02101de2-19fc-4c79-985e-89b10e9cfab2"
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "I = Inp(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "201SptVbE03k"
   },
   "source": [
    "Vectorizing the code (version two):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0QEvPCBzE03l"
   },
   "outputs": [],
   "source": [
    "def Inp(time=time, start=start, end=end, value=value):\n",
    "    I = np.zeros_like(time)\n",
    "    I[time>start] = value\n",
    "    I[time>end] = 0\n",
    "    return I\n",
    "  \n",
    "I = Inp(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "BNLAKOlGE03l",
    "outputId": "8b88f18b-3249-4bfb-cfb1-447faca935d4"
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "I = Inp(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lzl-C57IE03l"
   },
   "source": [
    "QUESTION: try to describe why the computation time to create the vector is different in the two versions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "py_kpyj7E03h"
   },
   "source": [
    "## a gentle introduction to matplotlib / pylab / pyplot\n",
    "\n",
    "Matplotlib is the major library for plotting in the python ecosystem. Import it like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WFnYFn4JE03f"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig_width = 15 # defining the size of plots\n",
    "phi = (np.sqrt(5)+1)/2 # setting up the golden ratio\n",
    "phi = phi**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.lookfor('bar plot', module='matplotlib.pyplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "x__aDnneE03l",
    "outputId": "c83a14ca-3a68-4f52-f485-0b9ba19f6f85"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, Inp(time))\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I (pA)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time.cumsum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6uLA7gLtE03m"
   },
   "source": [
    "QUESTION: re-run this calculation by adjusting the parameters to match figure 1\n",
    "\n",
    "If the code above looks scary, do not worry - Pythons are loving animals:\n",
    "\n",
    "* a nice [one-day course in python](https://github.com/NeuromatchAcademy/precourse)\n",
    "* resources abound  and [an army of geeks are ready to answer to any of your questions](https://stackoverflow.com/questions/tagged/python)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6uLA7gLtE03m"
   },
   "source": [
    "## a simple model of an integrate-and-fire neuron `leaky_IF`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6uLA7gLtE03m"
   },
   "source": [
    "Let's start with this membrane potential equation:\n",
    "\n",
    "$$\n",
    "\\tau \\cdot \\frac{dV}{dt} = -(V - V_{rest}) + R*I(t)\n",
    "$$\n",
    "\n",
    "$$\n",
    " {dV} = ( -(V - V_{rest}) + R*I(t) ) \\cdot \\frac{dt}{\\tau}\n",
    "$$\n",
    "\n",
    "with emission of a spike if $V > V_{rest}$, and then $V= V_{rest}$ for $3 ms$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "id": "xLYKnYodE03n",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "Vthreshold = -53\n",
    "Vreset = -80\n",
    "VRest = -70\n",
    "Vspike = 30\n",
    "R = 0.13\n",
    "tau = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "id": "xLYKnYodE03n",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def leaky_IF(time=time, inp=I, tau=tau, v0=VRest, R=R, \n",
    "             Vthreshold=Vthreshold, Vreset=Vreset, Vspike=Vspike, \n",
    "             VRest=VRest):\n",
    "    \n",
    "    V = np.ones_like(time)*v0\n",
    "    dt = time[1] - time[0]\n",
    "    for t in range(len(time)-1):\n",
    "        \n",
    "        dV = dt/tau * (-(V[t] - VRest) + R*inp[t])\n",
    "        V[t+1] = V[t] + dV\n",
    "        \n",
    "        if V[t]>Vthreshold:\n",
    "            V[t+1] = Vspike\n",
    "        if V[t] == Vspike:\n",
    "            V[t+1] = Vreset\n",
    "         \n",
    "    return V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SQGbZ9CPE03n"
   },
   "source": [
    "QUESTION: Set the parameter $R$ to obtain about 10 spikes - what is the interpretation of this parameter and what is the unit of measurement?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "RaOvJOSJE03n",
    "outputId": "37cc4a3e-5fc7-4af9-87b1-ac297e894be3"
   },
   "outputs": [],
   "source": [
    "V = leaky_IF(time, I)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V)\n",
    "ax.axhline(Vthreshold, c='g', ls='--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Focus: The linear part of the response"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "RaOvJOSJE03n",
    "outputId": "37cc4a3e-5fc7-4af9-87b1-ac297e894be3"
   },
   "outputs": [],
   "source": [
    "value = 150\n",
    "V = leaky_IF(time, Inp(time=time, start=100, end=500, value=value))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V)\n",
    "ax.axhline(Vthreshold, c='g', ls='--')\n",
    "ax.axvline(100, c='k', ls='--')\n",
    "ax.axvline(500, c='k', ls='--')\n",
    "ax.axhline(VRest, c='r', ls='--')\n",
    "#ax.axline((100, VRest), slope=R*value/tau,  c='r', ls='--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_xlim(80, 200)\n",
    "ax.set_ylim(-70, -50)\n",
    "\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6eXiMgxFE03n"
   },
   "source": [
    "QUESTION: What is the effect of $I_0$ on the discharge frequency?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "zsolXKMtE03n",
    "outputId": "1ea53219-f72d-44b4-d31e-20a94fa6f875"
   },
   "outputs": [],
   "source": [
    "for rho in np.geomspace(0.5, 2., 5):\n",
    "    I_0_ = rho*150\n",
    "    print('I_0 =', I_0_)\n",
    "    V = leaky_IF(time, Inp(value=I_0_))\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.plot(time, V)\n",
    "    ax.axhline(Vthreshold, c='g', ls='--')\n",
    "    ax.set_ylim(-83, 40)\n",
    "    ax.set_ylabel(\"Membrane potential (mV)\")\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "t7tc4XplE03n"
   },
   "source": [
    "Several tests show that this is perfectly reproducible, unlike figure 1A:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "collapsed": false,
    "id": "MO_9z9CiE03o",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "7771c542-3e84-4375-bae1-cca60215d427"
   },
   "outputs": [],
   "source": [
    "n_trials = 5\n",
    "V1 = np.zeros((n_trials, len(time)))\n",
    "\n",
    "for i in range(n_trials):\n",
    "    V1[i, :] = leaky_IF()\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V1.T)\n",
    "ax.axhline(Vthreshold, c='g', ls='--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LwhxhzlhE03o"
   },
   "source": [
    "QUESTION: this model seems not to reproduce the results, any explanation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LwhxhzlhE03o"
   },
   "source": [
    "# Creating a noisy input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LwhxhzlhE03o"
   },
   "source": [
    "A linear scattering model allows you to simply create a noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "collapsed": false,
    "id": "HOf58ZJtE03p",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "3003418c-ae53-46a6-f4bd-58df058d28f2"
   },
   "outputs": [],
   "source": [
    "def genNoise(time=time, tau_n=6, I_n=200, I_0=150, start=start, end=end):\n",
    "    x = np.ones_like(time)\n",
    "    dt = time[1] - time[0]\n",
    "    \n",
    "    for t in range(len(x)-1):\n",
    "        n = np.random.randn()*I_n\n",
    "        x[t+1] = (1-dt/tau_n)*x[t] + (dt*n/tau_n)\n",
    "    \n",
    "    x += I_0\n",
    "    x[time<start], x[time>end] = 0, 0\n",
    "    \n",
    "    return x\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, genNoise())\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I_b (pA)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.hist(genNoise(I_n=15, I_0=-0, start=-100, end=1000), bins=np.linspace(-100, 100, 100), density=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "g7jqFVbuE03p",
    "outputId": "a2c12aee-ac1b-4ba4-9646-23e597a04a3d"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, genNoise())\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I_b (pA)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6gUQtDZoE03p"
   },
   "source": [
    "QUESTION: does this model represent well the one in the paper? adjust $I_n$ and $I_0$ to get something that fits better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6gUQtDZoE03p"
   },
   "source": [
    "## LIF neuron with noisy input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6gUQtDZoE03p"
   },
   "source": [
    "Let's now observe the response of our LIF neuron to this input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "collapsed": false,
    "id": "w1RkAIIKE03q",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "6b044ecf-d0a1-448d-93cb-3ca3da8999af"
   },
   "outputs": [],
   "source": [
    "n_trials = 5\n",
    "V1 = np.zeros((n_trials,len(time)))\n",
    "\n",
    "for i in range(n_trials):\n",
    "    V1[i, :] = leaky_IF(time, genNoise())\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V1.T)\n",
    "ax.axhline(Vthreshold, c='g', ls='--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "ErJV7AllE03q",
    "outputId": "41d516a2-53c4-4abf-8e4c-6d4d60f6a3c2"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.eventplot([dt*np.where(V1.T[:, i] == Vspike)[0] for i in range(0, n_trials)], \n",
    "              colors='black', lineoffsets=1, linelengths=0.9);\n",
    "ax.set_ylabel('Trial number')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_xlim(time.min(), time.max())\n",
    "ax.set_ylim(-.5, n_trials-.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6EOFKlevE03q"
   },
   "source": [
    "QUESTION: adjust $I_n$ and $I_0$ to get something that better matches the observed output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "Ykp9623aE03q",
    "outputId": "695e5e90-da08-4564-acdb-8775e0e4501e"
   },
   "outputs": [],
   "source": [
    "for rho in np.geomspace(0.5, 2., 5):\n",
    "    I_0_ = rho*250\n",
    "    print('I_0 =', I_0_)\n",
    "    V= leaky_IF(time, genNoise(time, I_n=1000, I_0=I_0_))\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.plot(time, V)\n",
    "    ax.set_ylim(-83, 40)\n",
    "    ax.set_ylabel(\"Membrane potential (mV)\")\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "Vn2tCDWzE03q",
    "outputId": "d32dd877-3639-48a4-8a2e-bf99bec9220e"
   },
   "outputs": [],
   "source": [
    "for rho in np.geomspace(0.5, 2., 5):\n",
    "    I_n_ = rho*250\n",
    "    print('I_n=', I_n_)\n",
    "    V= leaky_IF(time, genNoise(time, I_n=I_n_, I_0=150))\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.plot(time, V)\n",
    "    ax.set_ylim(-83, 40)\n",
    "    ax.set_ylabel(\"Membrane potential (mV)\")\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "piLohkOPE03q"
   },
   "source": [
    "QUESTION: Do we obtain something reproducible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "collapsed": false,
    "id": "Oso0_M-tE03q",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "90d0378a-dc50-4e72-b4c4-45ea82357b48"
   },
   "outputs": [],
   "source": [
    "n_trials = 5\n",
    "V1 = np.zeros((n_trials,len(time)))\n",
    "\n",
    "for i in range(n_trials):\n",
    "    V1[i, :] = leaky_IF(time, genNoise(time, I_n=500, I_0=150))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V1.T)\n",
    "ax.axhline(Vthreshold, c='g', ls='--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tFzW7iAEE03r"
   },
   "source": [
    "## \"Frozen\" noise?\n",
    "\n",
    "QUESTION: what is the nature of the noise used in the article? why can it be described as [frozen noise](https://www.oxfordreference.com/view/10.1093/oi/authority.20110803095836900)?\n",
    "\n",
    "QUESTION: how to implement such a noise? what do you know about noise generators used in a computer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "riFZy5P5E03s",
    "outputId": "dee032e1-60d1-4c01-f73f-bc09d38ec72f"
   },
   "outputs": [],
   "source": [
    "help(np.random.seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "aV7qF_7qwEP5",
    "outputId": "02208c00-71c3-4d6d-ab7c-7c21ce759fbc"
   },
   "outputs": [],
   "source": [
    "np.random.rand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "aV7qF_7qwEP5",
    "outputId": "02208c00-71c3-4d6d-ab7c-7c21ce759fbc"
   },
   "outputs": [],
   "source": [
    "print('Tail') if np.random.rand() >.5 else print('Head')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "aV7qF_7qwEP5",
    "outputId": "02208c00-71c3-4d6d-ab7c-7c21ce759fbc"
   },
   "outputs": [],
   "source": [
    "np.random.randn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "a5TyWSutwElr",
    "outputId": "008ad58e-7696-440f-baa9-18497fabc3f4"
   },
   "outputs": [],
   "source": [
    "np.random.seed(1973) # douglas adams\n",
    "np.random.randn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "JG9RSUxywE-a",
    "outputId": "56a1cb94-4071-4206-e970-deea76c24322"
   },
   "outputs": [],
   "source": [
    "np.random.randn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "collapsed": false,
    "id": "fkySuKi_E03s",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "dbfe25e8-e4cd-435e-e992-cbe7175252f1"
   },
   "outputs": [],
   "source": [
    "def genNoise(time=time, tau_n=6, I_n=200, I_0=150, seed=42, start=start, end=end):\n",
    "    x = np.ones_like(time)\n",
    "    np.random.seed(seed)\n",
    "    dt = time[1] - time[0]\n",
    "    for t in range(len(x)-1):\n",
    "        n = np.random.randn()*I_n\n",
    "        x[t+1] = (1-dt/tau_n)*x[t]+ (dt*n/tau_n)\n",
    "    \n",
    "    x += I_0\n",
    "    x[time<start], x[time>end] = 0,0\n",
    "    \n",
    "    return x\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, genNoise())\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I_b (pA)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "MD-GZho2E03s",
    "outputId": "294437bb-db1f-4f12-ce3a-f0fc61b894e8"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, genNoise())\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I_b (pA)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wnmp4sfdE03s"
   },
   "source": [
    "## Multiple trials\n",
    "Here we show the conservation of the time of the spikes using a noisy input (frozen noise)\n",
    "\n",
    "QUESTION: adjust the parameter $I_0$ and $I_n$ to obtain about ten action potentials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 408
    },
    "collapsed": false,
    "id": "QSVes4kpE03t",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "ac60dd7b-96fd-40b2-e715-62ce88a92cf9"
   },
   "outputs": [],
   "source": [
    "n_trials = 25\n",
    "V1 = np.zeros((n_trials,len(time)))\n",
    "\n",
    "for i in range(n_trials):\n",
    "    V1[i, :] = leaky_IF(time, genNoise(I_n=600))\n",
    "\n",
    "print('number of spikes per trial :', (V1>0).sum(axis=1))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V1.T)\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "CQMizgIgx_zE",
    "outputId": "2a6c5e9e-6e7f-4d44-a90e-42272c864ec2"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.eventplot([dt*np.where(V1.T[:, i] == Vspike)[0] for i in range(0, n_trials)], \n",
    "              colors='black', lineoffsets=1, linelengths=0.9);\n",
    "ax.set_ylabel('Trial number')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_xlim(time.min(), time.max())\n",
    "ax.set_ylim(-.5, n_trials-.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AGiOTVTiE03t"
   },
   "source": [
    "We reproduce panel B: with a frozen noise, the neuron traces are reproducible.\n",
    "\n",
    "This also proves that we \"forgot\" to include a noise intrinsic to the dynamics of the neuron:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "id": "7LNrOgmME03t",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def leaky_IF(time=time, inp=I, tau=30, v0=-65, R=R, \n",
    "                Vthreshold=Vthreshold, Vreset=Vreset, Vspike=Vspike, \n",
    "                VRest=VRest, b=15, seed=None):\n",
    "    np.random.seed(seed)\n",
    "    V = np.ones_like(time)*v0\n",
    "    dt = time[1] - time[0]\n",
    "    for t in range(len(time)-1):\n",
    "        n=np.random.randn()\n",
    "        dV = dt * (-(V[t] - VRest) + R*(inp[t]+b*n))/tau\n",
    "        V[t+1] = V[t] + dV\n",
    "        \n",
    "        if V[t]>Vthreshold:\n",
    "            V[t+1]= Vspike\n",
    "        if V[t] == Vspike:\n",
    "            V[t+1]=Vreset\n",
    "         \n",
    "    return V\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZOD_kKvlE03t"
   },
   "source": [
    "Several tests show that with a square wave the spike times lose their reproducibility, as shown in the figure:\n",
    "\n",
    "QUESTION: adjust $I_0$ and $I_n$ to obtain a qualitatively similar number of spikes at the neuron output. To do this, try to control the number of spikes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "5fP5D8OzE03t",
    "outputId": "50d9ee9c-e881-4671-88b0-227e4d0e48f7"
   },
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.1, 2., 5):\n",
    "    b_ = rho*15\n",
    "    print('b =', b_)\n",
    "    VA = np.zeros((n_trials,len(time)))\n",
    "\n",
    "    for i in range(n_trials):\n",
    "        VA[i, :] = leaky_IF(time, genNoise(I_n=100, I_0=150), b=b_)\n",
    "\n",
    "    print('number of spikes per trial :', (VA>0).sum(axis=1))\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.eventplot([dt*np.where(VA.T[:, i] == Vspike)[0] for i in range(0, n_trials)], \n",
    "                colors='black', lineoffsets=1, linelengths=0.9);\n",
    "    ax.set_ylabel('Trial number')\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    ax.set_xlim(time.min(), time.max())\n",
    "    ax.set_ylim(-.5, n_trials-.5);\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "NsrSOc9Y1yPu",
    "outputId": "9d38cd03-1d4b-4f53-b2b1-032ff2914e5d"
   },
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.7, 1.3, 5):\n",
    "    I_n_ = rho*200\n",
    "    print('I_n =', I_n_)\n",
    "    VA = np.zeros((n_trials,len(time)))\n",
    "\n",
    "    for i in range(n_trials):\n",
    "        VA[i, :] = leaky_IF(time, genNoise(I_n=I_n_, I_0=150))\n",
    "\n",
    "    print('number of spikes per trial :', (VA>0).sum(axis=1))\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.eventplot([dt*np.where(VA.T[:, i] == Vspike)[0] for i in range(0, n_trials)], \n",
    "                colors='black', lineoffsets=1, linelengths=0.9);\n",
    "    ax.set_ylabel('Trial number')\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    ax.set_xlim(time.min(), time.max())\n",
    "    ax.set_ylim(-.5, n_trials-.5);\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GAeQiCY1E03u"
   },
   "source": [
    "QUESTION: see the influence of $I_0$ on the behavior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "ab_PHg4ME03u",
    "outputId": "e3777687-942f-4f60-bce1-31c3decc5d7a"
   },
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.7, 1.3, 5):\n",
    "    print('rho=', rho)\n",
    "    VA = np.zeros((n_trials,len(time)))\n",
    "\n",
    "    for i in range(n_trials):\n",
    "        VA[i, :] = leaky_IF(time, genNoise(I_n=500, I_0=rho*150))\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.plot(time, VA.T)\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    ax.set_ylabel('Membrane potential (mV)');\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KfsZTyL6E03v"
   },
   "source": [
    "QUESTION: see the influence of $I_0$ on the behavior, *when the noise amplitude is zero* :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "fM3Cm8f3E03v",
    "outputId": "0ac786ab-43ac-4b83-dfd7-ea22c6a76752"
   },
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.9, 1.1, 5):\n",
    "    I_0_ = rho*150\n",
    "    print('I_0_=', I_0_)\n",
    "    VA = np.zeros((n_trials,len(time)))\n",
    "\n",
    "    for i in range(n_trials):\n",
    "        VA[i, :] = leaky_IF(time, genNoise(I_n=0, I_0=I_0_))\n",
    "\n",
    "    print('number of spikes per trial :', (VA>0).sum(axis=1))\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "    ax.plot(time, VA.T)\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    ax.set_ylabel('Membrane potential (mV)');\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wrapping up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "K08jKjkSE03v"
   },
   "source": [
    "QUESTION: reproduce panel A: when the noise is zero, the traces of the neurons are not reproducible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 408
    },
    "collapsed": false,
    "id": "rtrEqNeuE03v",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "a394205c-54a0-4c5d-b570-6e928dd392a9"
   },
   "outputs": [],
   "source": [
    "seed = 2021\n",
    "VA = np.zeros((n_trials,len(time)))\n",
    "b_A = genNoise(I_n=0, I_0=150, seed=seed)\n",
    "\n",
    "for i in range(n_trials):\n",
    "    VA[i, :] = leaky_IF(time, b_A)\n",
    "\n",
    "\n",
    "print('number of spikes per trial :', (VA>0).sum(axis=1))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, VA.T)\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 373
    },
    "id": "AFmFMCJiE03v",
    "outputId": "8d5c48a2-3b9f-40c4-9238-c277c21ee3b2"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.eventplot([dt*np.where(VA.T[:, i] == Vspike)[0] for i in range(0, n_trials)], \n",
    "              colors='black', lineoffsets=1, linelengths=0.9);\n",
    "ax.set_ylabel('Trial number')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_xlim(time.min(), time.max())\n",
    "ax.set_ylim(-.5, n_trials-.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "emVfUHg6E03v"
   },
   "source": [
    "QUESTION: reproduce panel B: with frozen noise, the neuron traces are reproducible, even when the neuron has intrinsic noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 408
    },
    "collapsed": false,
    "id": "mO_j7LAyE03v",
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false,
    "outputId": "a56f8223-eaf9-41e6-ffad-8db545d3af4b"
   },
   "outputs": [],
   "source": [
    "VB = np.zeros((n_trials, len(time)))\n",
    "b_B = genNoise(I_n=500, I_0=150, tau_n=8, seed=seed)\n",
    "for i in range(n_trials):\n",
    "    VB[i, :] = leaky_IF(time, b_B)\n",
    "\n",
    "print('number of spikes per trial :', (VB>0).sum(axis=1))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, VB.T)\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('Membrane potential (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lYD_panUE03v"
   },
   "source": [
    "To sum up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 881
    },
    "id": "nfHo5Hr9E03w",
    "outputId": "195aa64b-1841-48ef-f38d-66bfbafde378"
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(3, 2, figsize=(fig_width, fig_width))\n",
    "\n",
    "axs[0][0].plot(time, b_A)\n",
    "axs[0][1].plot(time, b_B)\n",
    "axs[1][0].plot(time, VA.T)\n",
    "axs[1][1].plot(time, VB.T)\n",
    "axs[2][0].pcolor(time, range(n_trials), VA, vmax=Vthreshold)#, shading='nearest')\n",
    "axs[2][1].pcolor(time, range(n_trials), VB, vmax=Vthreshold)#, shading='nearest')\n",
    "for ax in axs.ravel(): \n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    ax.set_ylabel('Membrane potential (mV)');\n",
    "axs[2][0].set_ylabel('trial #');\n",
    "axs[2][1].set_ylabel('trial #');\n",
    "for i in range(2):\n",
    "    axs[0][i].set_ylabel('I_n (pA)')\n",
    "    axs[0][i].set_ylim(0, 400);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nUgOm4fCE03x"
   },
   "source": [
    "QUESTION: conclude quickly: to what extent has the phenomenon been explained? What is the conclusion about the response of neurons to different dynamical signals?"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "C_MainenSejnowski1995.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernel_info": {
   "name": "python3"
  },
  "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.6.9"
  },
  "nteract": {
   "version": "0.14.3"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}