Simulazione di un treno di impulsi di neurone in python

Question

Il modello su cui sto lavorando ha un neurone (modellato dalle equazioni di Hodgkin-Huxley), e il neurone riceve una serie di input sinaptici da altri neuroni perché si trova in una rete. Il modo standard per modellare gli input è con un treno spike costituito da un gruppo di impulsi di funzione delta che arrivano a una velocità specificata, come un processo di Poisson. Alcuni degli impulsi forniscono una reazione eccitatoria al neurone e alcuni forniscono un impulso inibitorio. Così la corrente sinaptica dovrebbe essere così:

Qui, Ne è il numero di neuroni eccitatori, Ni è inibitorio, le h di sono 0 o 1 (1 con probabilità p) che rappresenta o meno di un spike è stato trasmesso con successo, e $ t_k^l $ nella funzione delta è il tempo di scarica del picco del neurone kth (lo stesso per $ t_m^n $). Quindi l'idea alla base di come abbiamo provato a codificarlo era supporre che prima avessi 100 neuroni che fornivano impulsi nel mio neurone HH (80 dei quali sono eccitatori, 20 dei quali sono inibitori). Abbiamo quindi formato un array in cui una colonna enumerava i neuroni (quindi i neuroni # 0-79 erano eccitatori e # 80-99 erano inibitori). Abbiamo quindi controllato per vedere se c'è un picco in un intervallo di tempo, e se c'era, scegliere un numero casuale tra 0-1 e se è inferiore alla mia probabilità specificata p, quindi assegnarlo al numero 1, altrimenti farlo 0. Noi quindi traccia il voltaggio in funzione del tempo per guardare quando i picchi del neurone.

Penso che il codice funzioni, MA il problema è che non appena aggiungo altri neuroni nella rete (un documento afferma che hanno usato 5000 neuroni totali), ci vuole un tempo per funzionare, il che è semplicemente irrealizzabile per fare numeri simulazioni. La mia domanda è: c'è un modo migliore per simulare un treno di impulsi che pulsa in un neurone in modo che il calcolo sia sostanzialmente più veloce per un gran numero di neuroni nella rete? Ecco il codice che abbiamo provato: (è un po 'lungo perché le equazioni HH sono abbastanza dettagliato):

import scipy as sp 
import numpy as np 
import pylab as plt 

#Constants 
C_m = 1.0 #membrane capacitance, in uF/cm^2""" 
g_Na = 120.0 #Sodium (Na) maximum conductances, in mS/cm^2"" 
g_K = 36.0 #Postassium (K) maximum conductances, in mS/cm^2""" 
g_L = 0.3 #Leak maximum conductances, in mS/cm^2""" 
E_Na = 50.0 #Sodium (Na) Nernst reversal potentials, in mV""" 
E_K = -77.0 #Postassium (K) Nernst reversal potentials, in mV""" 
E_L = -54.387 #Leak Nernst reversal potentials, in mV""" 

def poisson_spikes(t, N=100, rate=1.0): 
    spks = [] 
    dt = t[1] - t[0] 
    for n in range(N): 
     spkt = t[np.random.rand(len(t)) < rate*dt/1000.] #Determine list of times of spikes 
     idx = [n]*len(spkt) #Create vector for neuron ID number the same length as time 
     spkn = np.concatenate([[idx], [spkt]], axis=0).T #Combine tw lists 
     if len(spkn)>0:   
      spks.append(spkn) 
    spks = np.concatenate(spks, axis=0) 
    return spks 

N = 100 
N_ex = 80 #(0..79) 
N_in = 20 #(80..99) 
G_ex = 1.0 
K = 4 

dt = 0.01 
t = sp.arange(0.0, 300.0, dt) #The time to integrate over """ 
ic = [-65, 0.05, 0.6, 0.32] 

spks = poisson_spikes(t, N, rate=10.) 

def alpha_m(V): 
     return 0.1*(V+40.0)/(1.0 - sp.exp(-(V+40.0)/10.0)) 

def beta_m(V): 
     return 4.0*sp.exp(-(V+65.0)/18.0) 

def alpha_h(V): 
     return 0.07*sp.exp(-(V+65.0)/20.0) 

def beta_h(V): 
     return 1.0/(1.0 + sp.exp(-(V+35.0)/10.0)) 

def alpha_n(V): 
     return 0.01*(V+55.0)/(1.0 - sp.exp(-(V+55.0)/10.0)) 

def beta_n(V): 
     return 0.125*sp.exp(-(V+65)/80.0) 

def I_Na(V, m, h): 
     return g_Na * m**3 * h * (V - E_Na) 

def I_K(V, n): 
     return g_K * n**4 * (V - E_K) 

def I_L(V): 
     return g_L * (V - E_L) 

def I_app(t): 
     return 3 

def I_syn(spks, t): 
    """ 
    Synaptic current 
    spks = [[synid, t],] 
    """ 
    exspk = spks[spks[:,0]<N_ex] # Check for all excitatory spikes 
    delta_k = exspk[:,1] == t # Delta function 
    if sum(delta_k) > 0: 
     h_k = np.random.rand(len(delta_k)) < 0.5 # p = 0.5 
    else: 
     h_k = 0 

    inspk = spks[spks[:,0] >= N_ex] #Check remaining neurons for inhibitory spikes 
    delta_m = inspk[:,1] == t #Delta function for inhibitory neurons 
    if sum(delta_m) > 0: 
     h_m = np.random.rand(len(delta_m)) < 0.5 #p =0.5 
    else: 
     h_m = 0 

    isyn = C_m*G_ex*(np.sum(h_k*delta_k) - K*np.sum(h_m*delta_m)) 

    return isyn 


def dALLdt(X, t): 
     V, m, h, n = X 
     dVdt = (I_app(t)+I_syn(spks,t)-I_Na(V, m, h) - I_K(V, n) - I_L(V))/C_m 
     dmdt = alpha_m(V)*(1.0-m) - beta_m(V)*m 
     dhdt = alpha_h(V)*(1.0-h) - beta_h(V)*h 
     dndt = alpha_n(V)*(1.0-n) - beta_n(V)*n 
     return np.array([dVdt, dmdt, dhdt, dndt]) 

X = [ic] 
for i in t[1:]: 
    dx = (dALLdt(X[-1],i)) 
    x = X[-1]+dt*dx 
    X.append(x) 

X = np.array(X)  
V = X[:,0]   
m = X[:,1] 
h = X[:,2] 
n = X[:,3] 
ina = I_Na(V, m, h) 
ik = I_K(V, n) 
il = I_L(V) 

plt.figure() 
plt.subplot(3,1,1) 
plt.title('Hodgkin-Huxley Neuron') 
plt.plot(t, V, 'k') 
plt.ylabel('V (mV)') 

plt.subplot(3,1,2) 
plt.plot(t, ina, 'c', label='$I_{Na}$') 
plt.plot(t, ik, 'y', label='$I_{K}$') 
plt.plot(t, il, 'm', label='$I_{L}$') 
plt.ylabel('Current') 
plt.legend() 

plt.subplot(3,1,3) 
plt.plot(t, m, 'r', label='m') 
plt.plot(t, h, 'g', label='h') 
plt.plot(t, n, 'b', label='n') 
plt.ylabel('Gating Value') 
plt.legend() 

plt.show() 

io non sono a conoscenza di altri pacchetti progettati specificamente per le reti neurali, ma ho voluto scrivere il mio, principalmente perché ho intenzione di fare un'analisi stocastica che richiede un bel po 'di dettagli matematici, e non so se quei pacchetti forniscano tali dettagli.

Answer 1

Profiling mostra che la maggior parte del vostro tempo viene speso in queste due righe:

if sum(delta_k) > 0: 

e

if sum(delta_m) > 0: 

Modifica ciascuno di questi per:

if np.any(...) 

accelera tutto in su di un fattore di 10. Dai un'occhiata a kernprof se vuoi fare più profiling line by line: https://github.com/rkern/line_profiler

Answer 2

In complemento alla risposta di Welch, è possibile utilizzare scipy.integrate.odeint per accelerare l'integrazione: la sostituzione

X = [ic] 
for i in t[1:]: 
    dx = (dALLdt(X[-1],i)) 
    x = X[-1]+dt*dx 
    X.append(x) 

da

from scipy.integrate import odeint 
X=odeint(dALLdt,ic,t) 

costi il ​​calcolo di oltre il 10 sul mio computer.