{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/laurentperrinet/2021_UE-neurosciences-computationnelles/blob/master/C_MainenSejnowski1995.ipynb)\n",
    "\n",
    "Le but ici de cette première tache est de créer un \"raster plot\" qui montre la reproducibilité d'un train de spike avec des répétitions du même stimulus. En particulier, nous allons essayer de répliquer la figure 1 de [Mainen & Sejnowski (1995)](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.299.8560&rep=rep1&type=pdf). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "fig_width = 15\n",
    "phi = (np.sqrt(5)+1)/2\n",
    "phi = phi**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mainen & Sejnowski, 1995\n",
    "\n",
    "## contexte\n",
    "Le but de cette première tache est de créer un \"raster plot\" qui montre la reproducibilité d'un train de spike avec des répétitions du même stimulus, comme dans ce travail dans la [rétine de rongeurs](https://laurentperrinet.github.io/2019-04-03_a_course_on_vision_and_modelization/#/1/3) ou dans le [cortex (V1) du chat](https://laurentperrinet.github.io/2019-04-03_a_course_on_vision_and_modelization/#/1/6).\n",
    "\n",
    "Ici, nous allons essayer de répliquer la figure 1 de [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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: écrire un résumé rapide du papier (max 5 lignes) et pourquoi ce résultat est *a priori* surprenant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# représentation du temps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = .5 # pas de discrétisation du temps\n",
    "time = np.arange(0, 1000, dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Création d'une fonction temporelle (version séquentielle):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "start = 150\n",
    "end = 750\n",
    "value = 350\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)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "I = Inp(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "version vectorisée:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "I = Inp(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: essayer de décrire pourquoi le temps de calcul pour créer le vecteur est différent dans les deux versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, I)\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I (pA)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: relancer ce calcul en ajustant les paramètres pour correspondre à la figure 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "## un modèle simple de neurone intègre-et-tire `leaky_IF`\n",
    "\n",
    "Commençons avec cet équation du potentiel membranaire:\n",
    "\n",
    "$$\n",
    "\\tau \\cdot \\frac{dV}{dt} = -(V - V_{rest}) + R*I(t)\n",
    "$$\n",
    "\n",
    "avec émission d'un \"spike\" si $V > V_{rest}$, et alors $V= V_{rest}$ pour $3 ms$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "Vthreshold = -53\n",
    "def leaky_IF(time=time, inp=I, tau=30, v0=-69, R=0.12, \n",
    "                Vthreshold=Vthreshold, Vreset=-80, Vspike=30, \n",
    "                VRest=-70):\n",
    "    V = np.ones_like(time)*v0\n",
    "    dt = time[1]\n",
    "    for t in range(len(time)-1):\n",
    "        dV = dt * (-(V[t] - VRest) + R*inp[t])/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": {},
   "source": [
    "QUESTION: régler le paramètre $R$ pour obtenir une dizaine de potentiels d'action - quel est l'interprétation de ce paramètre et quelles est l'unité de mesure?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.plot(time, np.ones_like(time)*Vthreshold, 'g--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: quel est l'effet de $I_0$ sur la fréquence de décharge?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.5, 1.2, 5):\n",
    "    I_0_ = rho*250\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.set_ylim(-83, 40)\n",
    "    ax.set_ylabel(\"potentiel de membrane (mV)\")\n",
    "    ax.set_xlabel('Time (ms)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plusieurs essais montrent que c'est parfaitement reproductible, contrairement à la figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "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.plot(time, np.ones_like(time)*Vthreshold, 'g--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: ce modèle semble ne pas reproduire les résultats, une explication?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Création d'un input bruité\n",
    "\n",
    "Un modèle linéaire de diffusion permet de créer simplement un bruit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def Bruit(time=time, tau_n=30, I_n=900, I_0=20, start=start, end=end):\n",
    "    dt = time[1]\n",
    "    x=np.ones_like(time)\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, Bruit())\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.plot(time, Bruit())\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I_b (pA)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: ce modèle représente-t-il bien celui dans le papier? régler $I_n$ et $I_0$ pour obtenir quelque chose qui corresponde mieux."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neurone LIF avec Input bruité\n",
    "\n",
    "Observons maintenant la réponse de notre neurone LIF à cette entrée:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "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, Bruit())\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.plot(time, V1.T)\n",
    "ax.plot(time, np.ones_like(time)*Vthreshold, 'g--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.eventplot([dt*np.where(V1.T[:, i] == 30)[0] for i in range(0, n_trials)], \n",
    "              colors='black', lineoffsets=1, linelengths=0.9);\n",
    "ax.set_ylabel('numéro essai')\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": {},
   "source": [
    "QUESTION: régler $I_n$ et $I_0$ pour obtenir quelque chose qui corresponde mieux à la sortie observée:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.5, 1.2, 5):\n",
    "    I_0_ = rho*250\n",
    "    print('I_0=', I_0_)\n",
    "    V= leaky_IF(time, Bruit(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(\"potentiel de membrane (mV)\")\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.6, 1.4, 5):\n",
    "    I_n_ = rho*500\n",
    "    print('I_n=', I_n_)\n",
    "    V= leaky_IF(time, Bruit(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(\"potentiel de membrane (mV)\")\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: obtient-on bien quelque chose de reproductible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "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, Bruit(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.plot(time, np.ones_like(time)*Vthreshold, 'g--')\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bruit \"gelé\" ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: quel est la nature du bruit utilisé dans l'article? pourquoi peut-on le décrire comme un [bruit gelé](https://www.oxfordreference.com/view/10.1093/oi/authority.20110803095836900) ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: comment implanter un tel bruit? que savez-vous des générateurs de bruit utilisés dans un ordinateur?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(np.random.seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def Bruit(time=time, tau_n=20, I_n=400, I_0=200, seed=42, start=start, end=end):\n",
    "    np.random.seed(seed)\n",
    "    dt = time[1]\n",
    "    x=np.ones_like(time)\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, Bruit())\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.plot(time, Bruit())\n",
    "ax.set_xlabel('Time (ms)')\n",
    "ax.set_ylabel('I_b (pA)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plusieurs trials\n",
    "Ici on montre le maintien du temps des spikes en utilisant un input bruité (bruit gelé)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: régler le paramètre $I_0$ et $I_n$ pour obtenir une dizaine de potentiels d'action:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "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, Bruit(I_n=400, I_0=300))\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('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "On reproduit le paneau B: avec un bruit gelé, les traces des neurones sont reproductibles.\n",
    "\n",
    "Celà prouve aussi que l'on a \"oublié\" d'inclure un bruit intrinsèque à la dynamique du neurone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def leaky_IF(time=time, inp=I, tau=30, v0=-65, R=0.12, \n",
    "                Vthreshold =-53, Vreset = -80, Vspike=30, \n",
    "                VRest=-70, b=100, seed=None):\n",
    "    np.random.seed(seed)\n",
    "    V = np.ones_like(time)*v0\n",
    "    dt = time[1]\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": {},
   "source": [
    "Plusieurs essais montrent qu'avec un créneau les temps des spikes perdent leur reproducibilité, comme sur la figure:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: régler $I_0$ et $I_n$ pour obtenir un nombre qualitativement similaire de spikes en sortie du neurone. Pour celà, essayer de controler le nombre de spikes :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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, Bruit(I_n=I_n_, I_0=175))\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('v (mV)');\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: voir l'influence de $I_0$ sur le comportement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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, Bruit(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('v (mV)');\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: voir l'influence de $I_0$ sur le comportement, *quand le bruit est nul* :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for rho in np.linspace(0.9, 1.1, 5):\n",
    "    I_0_ = rho*42\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, Bruit(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('v (mV)');\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QUESTION: reproduire le paneau A: quand le bruit est nul, les traces des neurones ne sont pas reproductibles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "seed = 2020\n",
    "VA = np.zeros((n_trials,len(time)))\n",
    "b_A = Bruit(I_n=0, I_0=42, 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('v (mV)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))\n",
    "ax.eventplot([dt*np.where(VA.T[:, i] == 30)[0] for i in range(0, n_trials)], \n",
    "              colors='black', lineoffsets=1, linelengths=0.9);\n",
    "ax.set_ylabel('numéro essai')\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": {},
   "source": [
    "QUESTION: reproduire le paneau B: avec un bruit gelé, les traces des neurones sont reproductibles, même quand le neurone possède un bruit intrinsèque:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "inputHidden": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "VB = np.zeros((n_trials, len(time)))\n",
    "b_B = Bruit(I_n=42, I_0=42, 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('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour résumer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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)\n",
    "axs[2][1].pcolor(time, range(n_trials), VB, vmax=Vthreshold)\n",
    "for ax in axs.ravel(): \n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    ax.set_ylabel('v (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": {},
   "source": [
    "QUESTION: conclure rapidement: à quel point a-t-on expliqué le phénomène? quelle est la conclusion sur la préférence des neurones à des signaux dynamiques?"
   ]
  }
 ],
 "metadata": {
  "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.9.2"
  },
  "nteract": {
   "version": "0.14.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}