Python：用矩陣編寫一個小的神經網絡

Question

我想模擬python中的兩個神經元網絡。爲每個神經元編寫單獨的方程很簡單，但是由於我想將代碼更加概括一些，以便在不重寫方程的情況下很容易增加神經元的數量。兩個神經網絡方程如下：

基本上，我具有通過在電壓方程中的最後一項耦合在一起的兩個霍奇金赫胥黎神經元。所以我想要做的是以這種方式編寫代碼，以便我可以輕鬆地擴展網絡。爲此，我爲神經元電壓創建了一個向量V：[V1，V2]，並創建一個向量X，其中X對門控變量m，h和n進行建模。所以我會有X = [[m1，h1，n1]，[m2，h2，n2]]。然而，目前的代碼並沒有產生尖峯，而是看起來電壓剛剛升高到無窮大。這表明門控變量X存在問題。門控變量m，h和n應始終在0和1之間，因此看起來好像門控變量剛好達到1並停留在那裏會導致電壓下降向上。我不確定是什麼導致他們停留在1.代碼正在運行，沒有產生任何錯誤。

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

我故意我不使用odeint整合我常微分方程，因爲我想隨機性添加到更新的版本，方程

要使用上面的歐拉法。無論如何，如果任何人都可以幫助我弄清楚如何修正這段代碼，以便發生預期的尖峯行爲，那就太棒了。謝謝！

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。具體而言，您的門控變量集成不起作用。你沒有正確更新它們。檢查缺少時間步數學。