Symulacja pociągu spike'a neuronu w pytonie

Question

Model, nad którym pracuję, ma neuron (wzorowany na równaniach Hodgkina-Huxleya), a sam neuron otrzymuje wiele wejść synaptycznych z innych neuronów, ponieważ znajduje się w sieci. Standardowym sposobem modelowania danych wejściowych jest kolec złożony z puli impulsów funkcji delta, które osiągają określoną prędkość, jako proces Poissona. Niektóre impulsy wywołują pobudzającą reakcję na neuron, a niektóre zapewniają impuls hamujący. Tak więc obecny synaptyczne powinna wyglądać następująco:

Tu Ne jest liczba neuronów pobudzających, Ni jest hamujący, że H są 0 lub 1 (1 z prawdopodobieństwem p) reprezentowanie czy nie skok został pomyślnie przesłany, a funkcja $ t_k^l $ w funkcji delta jest czasem rozładowania 1-tego przyrostu k-tego neuronu (to samo dla $ t_m^n $). Tak więc podstawową ideą tego, jak próbowaliśmy kodować to było założenie, że najpierw miałem 100 neuronów dostarczających impulsy do mojego neuronu HH (80 z nich jest pobudzających, z których 20 jest hamujących). Następnie utworzyliśmy tablicę, w której jedna kolumna wymieniła neurony (tak, że neurony # 0-79 były pobudzające, a # 80-99 były hamujące). Następnie sprawdziliśmy, czy w pewnym przedziale czasowym występuje skok, a jeśli tak, wybierz liczbę losową z przedziału od 0 do 1, a jeśli jest ona poniżej mojego określonego prawdopodobieństwa p, to przypisz mu numer 1, w przeciwnym razie wybierz 0. następnie wykreśl napięcie w funkcji czasu, aby zobaczyć, kiedy neuron skacze.

Myślę, że kod działa, ALE problem polega na tym, że jak tylko dodaję więcej neuronów do sieci (jeden papier twierdzi, że użył 5000 neuronów), trwa to na zawsze, co jest niemożliwe do wykonania symulacje. Moje pytanie brzmi: czy istnieje lepszy sposób symulacji pociągu spajkowego pulsującego w neuronie, tak że obliczenia są znacznie szybsze dla dużej liczby neuronów w sieci? Oto kod staraliśmy: (nie jest to trochę długo, ponieważ równania HH są dość szczegółowe):

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

nie jestem zaznajomiony z innymi pakietami zaprojektowanych specjalnie dla sieci neuronowych, ale chciałem napisać własną rękę, głównie dlatego, że planuję przeprowadzić analizę stochastyczną, która wymaga sporo szczegółów matematycznych i nie wiem, czy te pakiety dostarczają takich szczegółów.

Answer 1

Profilowanie pokazuje, że większość czasu są wydawane w tych dwóch wierszach:

if sum(delta_k) > 0: 

i

if sum(delta_m) > 0: 

Zmiana każdego z nich do:

if np.any(...) 

przyspiesza wszystkiego przez współczynnik 10. Rzuć okiem na kernprof, jeśli chcesz zrobić więcej profilowania linii po linii: https://github.com/rkern/line_profiler

Answer 2

W uzupełnieniu odpowiedzi Welch, można użyć scipy.integrate.odeint aby przyspieszyć proces integracji: zastąpienie

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

przez

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

przyspiesza obliczanie przez więcej niż 10 na moim komputerze.