<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Guia1</span> por <a xmlns:cc="http://creativecommons.org/ns#" href="http://nbviewer.ipython.org/urls/raw.githubusercontent.com/elsuizo/Redes_neuronales_Fuzzy/master/guia1.ipynb?create=1" property="cc:attributionName" rel="cc:attributionURL">Martin Noblía Marcos Panizzo Damián Presti</a> se distribuye bajo una Licencia Creative Commons Atribución-CompartirIgual 4.0 Internacional.
En este ejercicio trabajaremos con una neurona I&F con los siguientes parámetros:
tau = 20. # time constant(ms)
R = 10. # input Resistance(MOhms)
C = tau / R # capacidad
v_thres = 40 # thresold voltage(mV)
a) Programar la Neurona cuando se aplica un pulso de corrient constante ( $I_{max}$ )
b) Usar el método de Euler forward y backward. Mostrar que el método forward se vuelve inestable si $dt > 2 \tau$
c) Calcular la curva teórica que predice la frecuencia en función de la corriente aplicada.
d) Verificar la predicción por medio de simulaciones. O sea, simular distintos valores de corrientes y medir la frecuencia de emisión de spikes. Solapar sus mediciones computacionales con la curva teórica en c).
En este modelo suponemos que debajo de un potencial umbral $v_{thres}$ de disparo la neurona se comporta como un circuito RC y cuando sobrepasa $v_{thres}$ dispara un spike.
Como sabemos la ecuación diferencial que modela un circuito RC es:
$$C\frac{dv}{dt}=\frac{-v}{R}+I$$o equivalentemente:
$$\frac{dv}{dt}=\frac{-v}{RC}+\frac{I}{C}$$Para resolver numéricamente esta ecuación diferencial lineal se nos proponen los métodos de Euler(Forward y Backward) que consisten en aproximar la derivada como $\frac{V^{j}-V^{j-1}}{dt}$, donde $dt$ es el paso de discretización y en el caso forward la iteración de referencia es la $(j-1)$ por ello la ecuación diferencial discretizada queda:
$$\frac{V^{j}-V^{j-1}}{dt} + \frac{V^{j-1}}{\tau}=f^{j-1}$$o equivalentemente:
$$V^{j}=(1-\frac{dt}{\tau}) V^{j-1}+dt f^{j-1}$$Para el caso del método Backward la iteración de referencia es la $(j)$ por ello la ecuación diferencial discretizada queda:
$$\frac{V^{j}-V^{j-1}}{dt} + \frac{V^{j}}{\tau}=f^{j}$$o equivalentemente:
$$V^{j}=\frac{V^{j-1}+dt \, f^{j}}{1+\frac{dt}{\tau}}$$Modelamos a la neurona de acuerdo a la siguiente clase que tiene un metodo de simulación pedido
In [1]:
#comando para empotrar los gráficos dentro de la página
%matplotlib inline
In [2]:
import numpy as np
class NeuronaIF:
"""
Clase NeuronaIF que modela una neurona Integrate and fire
"""
#Constructor de la clase
def __init__(self, tau, R, v_thres):
self.tau = tau
self.R = R
self.v_thres = v_thres
self.C = tau / R
# Metodo de simulacion
def simulate(self, dt=1, method='forward', t_max=1000, I_max=.2, n_sims=1):
"""Metodo para simular la neurona"""
self.n_sims = n_sims
self.I_max = I_max
self.dt = dt
self.t_max = t_max
self.method = method
t_v = np.arange(0,self.t_max - self.dt + 1, dt,dtype=float)
nt_v = len(t_v)
dt_tau = self.dt / self.tau
one_dt_tau = 1 - dt_tau # forward Euler
m1_dt_tau = 1 / (1 + dt_tau) # backward Euler
Nt = len(t_v)
i_v_i = np.zeros_like(t_v) # nA
#i_v_i[(Nt*3/4.):(Nt*3/4.))
i_v_i[0:(Nt*3/4)] = I_max # nA
for j in xrange(0,n_sims):
i_v_tot = i_v_i
v_v = np.zeros_like(t_v)
s_v = np.zeros_like(t_v)
for i in xrange(1,nt_v):
# forward Euler
if (self.method == 'forward'):
v_v[i] =((1 - self.dt / self.tau) * v_v[i-1] )+ (self.dt * i_v_tot[i-1] / self.C)
# backward Euler
elif (self.method == 'backward'):
v_v[i] = (m1_dt_tau * v_v[i-1] )+ (self.dt * (i_v_tot[i] / self.C))
if (v_v[i] >= self.v_thres):
v_v[i] = 0.
s_v[i] = 1.
vteo = self.I_max * self.R * (1 - np.exp(-t_v / self.tau))
return t_v, v_v, i_v_i, vteo
def freq(self, t_ref):
self.t_ref = t_ref
# Corriente minima
I_min = self.v_thres / self.R
I = np.linspace(I_min,7)
f = 1 / (self.t_ref - self.tau * np.log(1 - self.v_thres/ (I * self.R)))
return I, f
# Metodo para imprimir los objetos
def __str__(self):
return '[NeuronaIF: tau:%s R:%s v_thres:%s]' % (self.tau, self.R, self.v_thres)
In [3]:
#Objeto de la clase NeuronaIF
neurona_1 = NeuronaIF(20., 10., 40.)
In [4]:
#Simulamos la neurona con los parametros default
t_v, v_v, i_v_i, vteo = neurona_1.simulate()
In [5]:
#Configuración de gráficos
import matplotlib.pyplot as plt #módulo para gráficos
plt.rcParams['figure.figsize'] = 8,6 #parámetros de tamaño
In [6]:
#Gráfico de neurona
fig, axes = plt.subplots()
axes.plot(t_v, v_v,'b')
axes.plot(t_v, i_v_i, 'r')
axes.plot(t_v, vteo, 'k')
axes.set_xlabel('$t$', fontsize=20)
axes.set_title('Simulacion de una Neurona IF', fontsize=20)
#fig.savefig('guia1_2_files/fig1.eps',format='eps',dpi=1000)
Out[6]:
De acuerdo a la ecuación de discretización del método forward si llamamos $a = 1-\frac{dt}{\tau}$ y evaluamos vemos que se puede generalizar la formula a la siguiente expresión:
$$V^{j}=a^{j-1} V^{1}+dt\sum_{i=1}^{j-1}a^{j-i-1}f^{i}$$Al ser una serie geométrica $|a|< 1$ y además para asegurar convergencia $dt < 2\tau$
Vamos a simular nuestra neurona con un $dt = 41$ que no cumpliría la condición de convergencia
In [7]:
# Simulacion de una neurona con un dt=41
t_v, v_v, i_v_i, vteo = neurona_1.simulate(41)
In [8]:
#Gráfico de neurona
fig, axes = plt.subplots()
axes.plot(t_v, v_v,'b')
axes.plot(t_v, i_v_i, 'r')
axes.plot(t_v, vteo, 'k')
axes.set_xlabel('$t$', fontsize=20)
axes.set_title('Simulacion de una Neurona IF', fontsize=20)
#fig.savefig('guia1_2_files/fig2.eps',format='eps',dpi=1000)
Out[8]:
Como sabemos la solución a a la ecuación diferencial es:
$$v(t)=IR(1-exp(\frac{-t}{\tau}))$$cuando $v(t)=v_{thres}$ tenemos:
$$v_{thres}=IR(1-exp(\frac{-t_{thres}}{\tau}))$$Por ello si despejamos $t_{thres}$ nos queda:
$$t_{thres}=-\tau \, log(\frac{1-v_{thres}}{IR})$$Luego como la frecuencia es la inversa del tiempo y además teniendo en cuenta que existe un tiempo muerto($t_{ref}$) donde la neurona no puede ser exitada
$$f=\frac{1}{t_{ref}+t_{thres}}=\frac{1}{t_{ref}-\tau \, log(1-\frac{v_{thres}}{IR})}$$
In [9]:
I,f = neurona_1.freq(.01) # el parametro que utilizamos es t_ref = 0.1 ms
In [10]:
# Plots de las frecuencias de spiking
plt.plot(I,f,'o')
plt.xlabel('I[nA]', fontsize=20)
plt.ylabel('frequency', fontsize=20)
#plt.savefig('guia1_2_files/fig3.eps',format='eps',dpi=1000)
Out[10]:
In [11]:
# Cálculo de la frecuencia real
frecuencia = np.zeros(len(I)) #vector de frecuencias reales
for i in xrange(1,len(I)):
aux = neurona_1.simulate(1, 'forward', 5000, I[i], 1) #simulación de la neurona para cada valor de corriente I[i]
tiempo = np.where(aux[1]==0.)
frecuencia[i] = 1./(tiempo[0][1]) #cálculo de frecuencia con el delta de tiempo
In [12]:
# Plots de las frecuencias de spiking
plt.plot(I,f,'o')
plt.plot(I,frecuencia,'r*')
plt.xlabel('I[nA]', fontsize=20)
plt.ylabel('frequency', fontsize=20)
plt.legend(['Frecuencia teorica', 'Frecuencia real'])
#plt.savefig('guia1_2_files/fig4.eps',format='eps',dpi=1000)
Out[12]:
In [12]: