Pythonでニューロンスパイク列をシミュレート

Question

私が取り組んでいるモデルにはニューロン（Hodgkin-Huxley方程式でモデル化されています）があり、ニューロン自体がネットワーク内にあるため他のニューロンからシナプス入力を受け取ります。入力をモデル化する標準的な方法は、ポアソン過程として指定された速度で到達する一連のデルタ関数パルスからなるスパイクトレインである。いくつかのパルスはニューロンへの興奮反応を提供し、いくつかは抑制パルスを提供する。そうシナプス電流は次のようになります

ここ

を、Neが興奮性ニューロンの数であり、Niが阻害され、Hさんの（確率pで1）0または1のいずれかであるか否かを示しますデルタ関数の$ t_k^l $は、k番目のニューロンの1番目のスパイクの放電時間です（$ t_m^n $と同じです）。だから私たちがこれをコーディングしようとした背景にある基本的な考え方は、最初に私のHHニューロンに80個のパルスを提供する100個のニューロンがあると仮定することでした（そのうちの80個は興奮性であり、20個は抑制性です）。次に、ある列がニューロンを列挙する配列を形成した（その結果、ニューロン＃0〜79は興奮性のものであり、＃80〜99は阻害性であった）。次に、ある時間間隔にスパイクが存在するかどうかを確認し、0-1の間の乱数を選択し、指定された確率pよりも小さい場合は番号1を割り当て、それ以外の場合は0にします。時間の関数として電圧をプロットして、ニューロンがいつスパイクするかを調べる。

私はコードがうまくいくと思いますが、問題は、ネットワークにいくつかのニューロンを追加すると（1つの論文では5000個のニューロンが使用されていると主張しています）、実行するには永遠にかかるということです。シミュレーション。私の質問は、ネットワーク内の多数のニューロンの計算が大幅に高速になるように、スパイク列をニューロンにパルス状に流すようシミュレートするより良い方法があるかどうかです。 （HH方程式は非常に詳述されているので、それは少し長いです）：ここでは試したコードがある

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

私は、具体的ニューラルネットワークのために設計された他のパッケージに慣れていないんだけど、私は自分自身を書きたかったです、私はかなりの数学的な詳細を必要とする確率論的分析を行うことを計画しているので、そのようなパッケージがそのような詳細を提供するかどうかはわかりません。

Answer 1

プロファイリングは、あなたの時間のほとんどは、この2行に費やされていることを示しています。すべてをスピードアップ

if np.any(...) 

：に

if sum(delta_k) > 0: 

と

これらのそれぞれを変更する

if sum(delta_m) > 0: 

線のプロファイリングをもっとやりたいのなら、kernprofを見てみましょう： https://github.com/rkern/line_profiler

Answer 2

ウェルチの答えを補完では、統合を加速するscipy.integrate.odeintを使用することができます。自分のコンピュータ上で10以上で

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

によって

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

計算速度置き換えます。