Python:用矩阵编写一个小型神经网络

Question

我正在尝试在 python 中模拟一个双神经元网络。为每个神经元编写单独的方程式非常简单，但由于我想对代码进行更多概括，这样就可以很容易地增加神经元的数量，而无需一遍又一遍地重写方程式。两个神经网络方程如下:

基本上，我有两个 Hodgkin-Huxley 神经元，它们通过电压方程中的最后一项耦合在一起。所以我想做的是以这种方式编写代码，以便我可以轻松扩展网络。为此，我为神经元电压创建了一个向量 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 来积分我的 ODE，因为我想稍后将随机性添加到方程中，因此想使用上面的 Euler 方法。不管怎样，如果有人能帮我弄清楚如何修复这段代码，以便出现预期的尖峰行为，那就太好了。谢谢!

Answer 1

检查您的输入电流和电导。如果您按照现代化的形式主义编写代码，它也会使您的生活更轻松，即 dm/dt = (m_inf - m)/tau。但是，具体来说，您的门控变量集成不起作用。您没有正确更新它们。检查缺少的时间步数学。