Python：行列を使った小さなニューラルネットワークの作成

Question

私はPythonで2つのニューロンネットワークをシミュレートしようとしています。方程式を何度も書き直すことなくニューロンの数を増やすのが簡単であるように、コードをもう少し一般化したいので、各ニューロンについて別々の方程式を書くのは簡単です。 2つのニューラルネットワークの方程式は以下のとおりである。

基本的に、私は、電圧方程式の最後の項を介して互いに結合された2つのホジキン・ハクスリーのニューロンを持っています。だから私がしたいのは、簡単にネットワークを拡張できるようにコードを書くことです。これを行うために、ニューロン電圧[V1、V2]のベクトルVを作成し、Xがゲーティング変数m、h、nをモデル化するベクトルXを作成します。ですから、私はX = [[m1、h1、n1]、[m2、h2、n2]]を持つでしょう。しかし、現在のところ、コードはスパイキングを生成していませんが、むしろ無限に電圧が上がるように見えます。これはゲーティング変数Xに問題があることを示唆しています。ゲーティング変数m、h、nは常に0と1の間でなければなりません。ゲーティング変数がちょうど1に達していて、アップ。私は何が原因で1にとどまるのかは分かりません。コードは実行されており、エラーは発生していません。私は意図的に私は、後に、したがって、方程式にstochasticityを追加したいと思いますので、私のODEを統合するodeintを使用していない

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

NN=2 #Number of Neurons in Model 

dt=0.01 
T = sp.arange(0.0, 1000.0, dt) 
nt = len(T) # total number of time steps 

    # Constants 
gNa = 120.0 # maximum conducances, in mS/cm^2 
gK = 36.0 
gL = 0.3 
ENa = 50.0 # Nernst reversal potentials, in mV 
EK = -77 
EL = -54.387 

#Coupling Terms 
Vr = 20 
w = 1 
e11 = e22 = 0 
e12 = e21 = 0.1 
E = np.array([[e11, e12], [e21, e22]]) 

#Gating Variable Transition Rates 
def alpham(V): return (0.1*V+4.0)/(1.0 - sp.exp(-0.1*V-4.0)) 
def betam(V): return 4.0*sp.exp(-(V+65.0)/18.0) 
def alphah(V): return 0.07*sp.exp(-(V+65.0)/20.0) 
def betah(V): return 1.0/(1.0 + sp.exp(-0.1*V-3.5)) 
def alphan(V): return (0.01*V+0.55)/(1.0 - sp.exp(-0.1*V-5.5)) 
def betan(V): return 0.125*sp.exp(-(V+65.0)/80.0) 
def psp(V,s): return ((5*(1-s))/(1+sp.exp(-(V+3)/8)))-s 

#Current Terms 
def I_Na(V,x): return gNa * (x[:,0]**3) * x[:,1] * (V - ENa) #x0=m, x1=h, x2=n 
def I_K(V,x): return gK * (x[:,2]**4) * (V - EK) 
def I_L(V): return gL * (V - EL) 
def I_inj(t): return 10.0 

#Initial Conditions 
V = np.zeros((nt,NN)) #Voltage vector 
X = np.zeros((nt,NN,3)) #Gating Variables m,h,n (NN neurons x 3 gating variables) 
S = np.zeros((nt,NN)) #Coupling term 

dmdt = np.zeros((nt,NN)) 
dhdt = np.zeros((nt,NN)) 
dndt = np.zeros((nt,NN)) 

V[0,:] = -65.0 
X[0,:,0] = alpham(V[0,:])/(alpham(V[0,:])+betam(V[0,:])) #m 
X[0,:,1] = alphah(V[0,:])/(alphah(V[0,:])+betah(V[0,:])) #h 
X[0,:,2] = alphan(V[0,:])/(alphan(V[0,:])+betan(V[0,:])) #n 

alef = 5.0/(1+sp.exp(-(V[0,:]+3)/8.0)) 
S[0,:] = alef/(alef+1) 

dmdt[0,:] = alpham(V[0,:])*(1-X[0,:,0])-betam(V[0,:])*X[0,:,0] 
dhdt[0,:] = alphah(V[0,:])*(1-X[0,:,1])-betah(V[0,:])*X[0,:,1] 
dndt[0,:] = alphan(V[0,:])*(1-X[0,:,2])-betan(V[0,:])*X[0,:,2] 

#Euler-Maruyama Integration 
for i in xrange(1,nt): 
    V[i,:]= V[i-1,:]+dt*(I_inj(i-1)-I_Na(V[i-1,:],X[i-1,:])-I_K(V[i-1,:],X[i-1,:])-I_L(V[i-1,:]))+dt*((Vr-V[i-1,:])/w * np.dot(E,S[i-1,:]))  

    #Gating Variable 
    dmdt[i,:] = dmdt[i-1,:] + alpham(V[i-1,:])*(1-X[i-1,:,0])-betam(V[i-1,:])*X[i-1,:,0] 
    dhdt[i,:] = dhdt[i-1,:] + alphah(V[i-1,:])*(1-X[i-1,:,1])-betah(V[i-1,:])*X[i-1,:,1] 
    dndt[i,:] = dndt[i-1,:] + alphan(V[i-1,:])*(1-X[i-1,:,2])-betan(V[i-1,:])*X[i-1,:,2] 
    z = np.array([dmdt[i-1,:],dhdt[i-1,:],dndt[i-1,:]]).T 

    #Gating Variable Constraints (0<m,h,n<1) 
    X[i,:,0] = max(0,min(X[i,:,0].all(),1)) 
    X[i,:,1] = max(0,min(X[i,:,1].all(),1)) 
    X[i,:,2] = max(0,min(X[i,:,2].all(),1)) 

    #Update Gating Variables 
    X[i,:,:]= X[i-1,:,:]+dt*(z) 

    #Coupling Term 
    S[i,:] = S[i-1,:]+dt*psp(V[i-i,:],S[i-1,:]) 

V1 = V[:,0] 
V2 = V[:,1] 

plt.plot(T,V1, 'red') 
plt.plot(T,V2, 'blue') 

plt.show() 

は、上記オイラー法を使用します。とにかく、誰かがこのコードを修正して、予想されるスパイキング動作が起こるようにする手助けをすることができれば、それは素晴らしいことでしょう。ありがとうございました！

Answer 1

あなたは

\dot(x)_j = f(x_j) + \sum_j C_ij g(x_j, x_i)

を持っているので、あなたもそのようなあなたの乗車右側を書くことができます。 x_jはベクトルである。 （v、m、h、n s）である。

この発振器の例では、出発点のようになります。

x.shapeはn_onsite_variables、n_oscillators

def f(x): 
    return np.array([x[1, :], -x[0, :]]) 


def g(xi, xj): 
    return xi[0] - xj[0] 

def rhs(x, t, c): 
    coupling = np.array([sum(cij*g(xi,xj) for cij, xj in zip(ci, x.T)) 
              for ci, xi in zip(c, x.T)]) 
    coupling = np.outer(np.arange(2), coupling) # coupling in x'' 

    return f(x) + coupling 

x0 = np.random.random(size=(2, 3)) 
>>> array([[ 0.74386362, 0.85799588, 0.70501992], 
    [ 0.65903405, 0.41575505, 0.93166973]]) 

def ring(n): 
    c = np.eye(n, k=1) + np.eye(n, k=-1) 
    c[0, -1] = 1 
    c[-1, 0] = 1 
    return c 

c = ring(3) 
x1 = rhs(x0, 0, c) 
>>> array([[ 0.65903405, 0.41575505, 0.93166973], 
    [-0.81915217, -0.59088767, -0.89683958]]) 

私はこれを向上させることができると確信しています。特に、より一般的なカップリングが必要な場合。

Answer 2

入力電流とコンダクタンスを確認してください。近代化された形式に従ってコードを書くと、dm/dt =（m_inf - m）/ tauというあなたの人生も楽になります。具体的には、ゲーティング変数の統合は機能しません。あなたはそれらを正しく更新していません。タイムステップの計算がないか確認してください。