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).

Ce notebook a été élaboré lors d'un TP dans le cadre du Master 1 de sciences cognitives de l'Université d'Aix-Marseille.

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
fig_width = 15
phi = (np.sqrt(5)+1)/2
phi = phi**2

Mainen & Sejnowski, 1995¶

contexte¶

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 ou dans le cortex (V1) du chat.

Ici, nous allons essayer de répliquer la figure 1 de Mainen & Sejnowski (1995):

Mainen Sejnowski 1995

QUESTION: écrire un résumé rapide du papier (max 5 lignes) et pourquoi ce résultat est a priori surprenant:

REPONSE: Le but de l’étude était de déterminer directement la précision temporelle avec laquelle les neurones corticaux sont capables d'encoder un stimulus dans un train de pics. Lors d'un examen de la fiabilité de la génération de pics à l'aide d'enregistrements de neurones dans des tranches néocorticales de rats. Dans cet article les auteurs ont étudié la réponse d'un neurone de la couche 5 du cortex en fonction de deux types de stimuli. Le résultat principal semble a priori surprenant car on s'attendrait à obtenir une réponse reproductible pour un signal continu et une réponse non reproductible pour une sigal bruité, or, cet article semble démontrer l'inverse.

représentation du temps¶

In [2]:

dt = .5 # pas de discrétisation du temps
time = np.arange(0, 1000, dt)

Création d'une fonction temporelle (version séquentielle):

In [3]:

start = 150
end = 750
value = 200

def Inp(time=time, start=start, end=end, value=value):
    x=[]
    for t in range(len(time)):
        if start < time[t] < end :
            x.append(value)
        else:
            x.append(0)
    return x

I = Inp(time)

In [4]:

%%timeit
I = Inp(time)

1.63 ms ± 100 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

version vectorisée:

In [5]:

def Inp(time=time, start=start, end=end, value=value):
    I = np.zeros_like(time)
    I[time>start] = value
    I[time>end] = 0
    return I
  
I = Inp(time)

In [6]:

%%timeit
I = Inp(time)

21.3 µs ± 8.13 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)

QUESTION: essayer de décrire pourquoi le temps de calcul pour créer le vecteur est différent dans les deux versions:

REPONSE: Le temps de calcul est plus rapide quand le code est vectorisé. La vectorisation automatique du code, est une fonctionnalité du compilateur qui permet à certaines parties des programmes séquentiels d'être transformés en programme parallèles équivalents afin de produire du code qui sera optimisé pour un l'unité centrale de traitement (CPU).

In [7]:

fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi))
ax.plot(time, I)
ax.set_xlabel('Time (ms)')
ax.set_ylabel('I (pA)');

un modèle simple de neurone intègre-et-tire `leaky_IF`¶

Commençons avec cet équation du potentiel membranaire:

$$ \tau \cdot \frac{dV}{dt} = -(V - V_{rest}) + R*I(t) $$

avec émission d'un "spike" si $V > V_{threshold}$, et alors $V= V_{reset}$:

In [8]:

Vthreshold = -53
def leaky_IF(time=time, inp=I, tau=30, v0=-69, R=0.11, 
                Vthreshold=Vthreshold, Vreset=-80, Vspike=30, 
                VRest=-70):
    V = np.ones_like(time)*v0
    dt = time[1]
    for t in range(len(time)-1):
        dV = dt * (-(V[t] - VRest) + R*inp[t])/tau
        V[t+1] = V[t] + dV
        
        if V[t]>Vthreshold:
            V[t+1]= Vspike
        if V[t] == Vspike:
            V[t+1]=Vreset
         
    return V

QUESTION: régler le paramètre $R$ pour obtenir une dizaine de potentiels d'action - quelle est l'interprétation de ce paramètre et quelle est l'unité de mesure?

In [9]:

V = leaky_IF(time, I)
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, V)
ax.plot(time, np.ones_like(time)*Vthreshold, 'g--')
ax.set_xlabel('Time (ms)')
ax.set_ylabel('v (mV)');

REPONSE: Ce paramètre est la résistance de notre modèle, représentée en $\Omega$ (ohm). On observe qu'un $R$ grand nous fera atteindre le V(seuil) plus rapidement et fournira un profil de décharge haute fréquence contrairement à un $R$ plus faible.

QUESTION: quel est l'effet de $I_0$ sur la fréquence de décharge?

In [10]:

for rho in np.linspace(0.5, 2.0, 5):
    I_0_ = rho*250
    print('I_0=', I_0_)
    V= leaky_IF(time, Inp(value=I_0_))

    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
    ax.plot(time, V)
    ax.set_ylim(-83, 40)
    ax.set_ylabel("potentiel de membrane (mV)")
    ax.set_xlabel('Time (ms)')

    plt.show()

I_0= 125.0

I_0= 218.75

I_0= 312.5

I_0= 406.25

I_0= 500.0

REPONSE: La valeur I0 mesure l'intensité du courant entrant injecté, correspond à l'injection d'un courant constant. Plus la valeur est haute, plus la fréquence de décharge augmente.

Plusieurs essais montrent que c'est parfaitement reproductible, contrairement à la figure:

In [11]:

n_trials = 15
V1 = np.zeros((n_trials,len(time)))
for i in range(n_trials):
    V1[i, :] = leaky_IF()

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, V1.T)
ax.plot(time, np.ones_like(time)*Vthreshold, 'g--')
ax.set_xlabel('Time (ms)')
ax.set_ylabel('v (mV)');

QUESTION: ce modèle semble ne pas reproduire les résultats, une explication?

In [12]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.eventplot([dt*np.where(V1.T[:, i] == 30)[0] for i in range(0, n_trials)], 
              colors='black', lineoffsets=1, linelengths=0.9)
ax.set_ylabel('numéro essai')
ax.set_xlabel('Time (ms)')
ax.set_xlim(time.min(), time.max())
ax.set_ylim(-.5, n_trials-.5);

REPONSE: Ce modèle ne semble pas reproduire les résultats, car il n'y a pas de variabilité dans le signal d'entrée. Il faut ajouter du bruit dans ce signal pour imiter l'entrée de nature variable qu'aurait un neurone in-vivo.

Création d'un input bruité¶

Un modèle linéaire de diffusion permet de créer simplement un bruit:

In [13]:

def Bruit(time=time, tau_n=15, I_n=1000, I_0=200, start=start, end=end):
    dt = time[1]
    x=np.ones_like(time)
    for t in range(len(x)-1):
        n = np.random.randn()*I_n
        x[t+1]=(1-dt/tau_n)*x[t]+ (dt*n/tau_n)
    
    x += I_0
    x[time<start], x[time>end] = 0, 0
    
    return x

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, Bruit())
ax.set_xlabel('Time (ms)')
ax.set_ylabel('I_b (pA)');

In [14]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, Bruit())
ax.set_xlabel('Time (ms)')
ax.set_ylabel('I_b (pA)');

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.

REPONSE: Les paramètres In = 1000pA et I0 = 200pA donnent un ordre de grandeur de ce qui est présenté dans le papier.

Neurone LIF avec Input bruité¶

Observons maintenant la réponse de notre neurone LIF à cette entrée:

In [15]:

n_trials = 5
V1 = np.zeros((n_trials,len(time)))

for i in range(n_trials):
    V1[i, :] = leaky_IF(time, Bruit())

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, V1.T)
ax.plot(time, np.ones_like(time)*Vthreshold, 'g--')
ax.set_xlabel('Time (ms)')
ax.set_ylabel('v (mV)');

In [16]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.eventplot([dt*np.where(V1.T[:, i] == 30)[0] for i in range(0, n_trials)], 
              colors='black', lineoffsets=1, linelengths=0.9);
ax.set_ylabel('numéro essai')
ax.set_xlabel('Time (ms)')
ax.set_xlim(time.min(), time.max())
ax.set_ylim(-.5, n_trials-.5);

QUESTION: régler $I_n$ et $I_0$ pour obtenir quelque chose qui corresponde mieux à la sortie observée:

In [17]:

for rho in np.linspace(0.5, 1.2, 5):
    I_0_ = rho*200
    print('I_0=', I_0_)
    V= leaky_IF(time, Bruit(time, I_n=1000, I_0=I_0_))

    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
    ax.plot(time, V)
    ax.set_ylim(-83, 40)
    ax.set_ylabel("potentiel de membrane (mV)")
    ax.set_xlabel('Time (ms)')

    plt.show()

I_0= 100.0

I_0= 135.0

I_0= 170.0

I_0= 204.99999999999997

I_0= 240.0

REPONSE: I0 = 200 semble qualitativement correct.

In [18]:

for rho in np.linspace(0.6, 1.4, 5):
    I_n_ = rho*1000
    print('I_n=', I_n_)
    V= leaky_IF(time, Bruit(time, I_n=I_n_, I_0=200))

    fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
    ax.plot(time, V)
    ax.set_ylim(-83, 40)
    ax.set_ylabel("potentiel de membrane (mV)")
    ax.set_xlabel('Time (ms)')

    plt.show()

I_n= 600.0

I_n= 799.9999999999999

I_n= 1000.0

I_n= 1200.0

I_n= 1400.0

QUESTION: obtient-on bien quelque chose de reproductible?

REPONSE: Non, à chaque nouvelle exécution, les graphes seront différents, et ce quelques soient les valeurs de I_0, et dans la mesure ou la valeur de I_n n'approche pas de 0 (dans ce cas il n'y aurait plus de variabilité et on aurait plus de bruit). Cela a pour origine le bruit utilisé dans le modèle. Celui ci change à chaque nouvelle simulation du modèle, simulant un stimulus différent, la réponse du neurone n'est donc jamais la même:

In [19]:

n_trials = 5
V1 = np.zeros((n_trials,len(time)))

for i in range(n_trials):
    V1[i, :] = leaky_IF(time, Bruit(time, I_n=1000, I_0=200))

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, V1.T)
ax.plot(time, np.ones_like(time)*Vthreshold, 'g--')
ax.set_xlabel('Time (ms)')
ax.set_ylabel('v (mV)');

In [20]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.eventplot([dt*np.where(V1.T[:, i] == 30)[0] for i in range(0, n_trials)], 
              colors='black', lineoffsets=1, linelengths=0.9);
ax.set_ylabel('numéro essai');
ax.set_xlabel('Time (ms)')
ax.set_xlim(time.min(), time.max())
ax.set_ylim(-.5, n_trials-.5);

Bruit "gelé" ?¶

QUESTION : quel est la nature du bruit utilisé dans l'article? pourquoi peut-on le décrire comme un bruit gelé ?

REPONSE: La nature du bruit utilisé dans l'article est un bruit blanc Gaussien filtré généré par l'ordinateur. Le bruit blanc gaussien est un bruit blanc qui suit une loi normale de moyenne et variance données. On peut le décrire comme un "bruit gelé", c'est un bruit qui se répète à l'identique au cours des essais.

QUESTION : comment implanter un tel bruit? que savez-vous des générateurs de bruit utilisés dans un ordinateur?

REPONSE: On peut implanter ce bruit par informatique en "gelant" les inputs. Les générateurs de bruit utilisés dans un ordinateur, servent à générer un signal de bruit précis tout en maintenant une haute précision du ratio entre le bruit fourni par l'utilisateur et le bruit interne, au dessus de la gamme étendue du niveau de puissance et de fréquence du signal. Pour générer un bruit gelé, il suffit de fixer l'état de notre machine dans son processus aléatoire, cela peut se faire via la fonction random.seed de la bibliothèque numpy par exemple:

In [21]:

help(np.random.seed)

Help on built-in function seed:

seed(...) method of numpy.random.mtrand.RandomState instance
    seed(self, seed=None)
    
    Reseed a legacy MT19937 BitGenerator
    
    Notes
    -----
    This is a convenience, legacy function.
    
    The best practice is to **not** reseed a BitGenerator, rather to
    recreate a new one. This method is here for legacy reasons.
    This example demonstrates best practice.
    
    >>> from numpy.random import MT19937
    >>> from numpy.random import RandomState, SeedSequence
    >>> rs = RandomState(MT19937(SeedSequence(123456789)))
    # Later, you want to restart the stream
    >>> rs = RandomState(MT19937(SeedSequence(987654321)))

In [22]:

def Bruit(time=time, tau_n=15, I_n=1000, I_0=200, seed=42, start=start, end=end):
    np.random.seed(seed)
    dt = time[1]
    x=np.ones_like(time)
    for t in range(len(x)-1):
        n = np.random.randn()*I_n
        x[t+1] = (1-dt/tau_n)*x[t]+ (dt*n/tau_n)
    
    x += I_0
    x[time<start], x[time>end] = 0,0
    
    return x

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, Bruit())
ax.set_xlabel('Time (ms)')
ax.set_ylabel('I_b (pA)');

On vérifie que si on relance la fonction générant le bruit aléatoire, on obtient le même signal:

In [23]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/phi))
ax.plot(time, Bruit())
ax.set_xlabel('Time (ms)')
ax.set_ylabel('I_b (pA)');