2017-03-16 87 views
3

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

enter image description here

基本上,我具有通過在電壓方程中的最後一項耦合在一起的兩個霍奇金赫胥黎神經元。所以我想要做的是以這種方式編寫代碼,以便我可以輕鬆地擴展網絡。爲此,我爲神經元電壓創建了一個向量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整合我常微分方程,因爲我想隨機性添加到更新的版本,方程

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

+0

如果你正在談論HH模型,這是一個生物物理學模型,我會毫不猶豫地用'neural-network'來標記這個。 –

+0

@ juanpa.arrivillaga好吧,我改變了它 – Brenton

+0

我建議讓代碼的各個部分先工作,然後逐步增加複雜性。例如,考慮只爲一個神經元獲得被動電流,然後獲得一個神經元的尖峯電流,然後添加第二個神經元的無源電流,然後添加兩個神經元的尖峯行爲。現在問題可能出現在方程式或集成代碼中,並且以步驟式方式工作應該使調試變得易於處理。 – Justas

回答

0

既然你有

\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]]) 

我相信這可以得到改善。特別是如果你想要一個更一般的耦合。

1

檢查您的輸入電流和您的電導。如果你根據現代化的形式化編寫你的代碼,它也會讓你的生活更輕鬆,那就是dm/dt =(m_inf - m)/ tau。具體而言,您的門控變量集成不起作用。你沒有正確更新它們。檢查缺少時間步數學。

+0

感謝您的幫助!事實證明,我用於dmdt,dhdt和dndt的等式沒有正確更新。它現在已經修復了! – Brenton