在 python 中模拟神经元脉冲序列

Posted 2023-03-12

【中文标题】在 python 中模拟神经元脉冲序列

【英文标题】：Simulating a neuron spike train in python

【发布时间】：2016-07-03 04:26:55

【问题描述】：

我正在研究的模型有一个神经元（由 Hodgkin-Huxley 方程建模），并且神经元本身接收来自其他神经元的一堆突触输入，因为它位于网络中。对输入进行建模的标准方法是使用由一组以指定速率到达的 delta 函数脉冲组成的尖峰序列，作为泊松过程。一些脉冲提供对神经元的兴奋反应，而一些提供抑制脉冲。所以突触电流应该是这样的：

这里，Ne 是兴奋性神经元的数量，Ni 是抑制性神经元，h 是 0 或 1（概率 p 为 1）表示是否成功传输了一个尖峰，而 delta 中的 $t_k^l$函数是第 k 个神经元的第 l^th 个尖峰的放电时间（$t_m^n$ 相同）。所以我们尝试编码的基本想法是假设首先我有 100 个神经元向我的 HH 神经元提供脉冲（其中 80 个是兴奋性的，其中 20 个是抑制性的）。然后我们形成一个阵列，其中一列列举了神经元（因此神经元#0-79 是兴奋性的，#80-99 是抑制性的）。然后我们检查在某个时间间隔内是否有尖峰，如果有，则在 0-1 之间选择一个随机数，如果它低于我指定的概率 p，则将其分配为数字 1，否则将其设为 0。我们然后将电压绘制为时间的函数，以查看神经元何时出现峰值。

我认为代码有效，但问题是，一旦我在网络中添加更多神经元（一篇论文声称他们总共使用了 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()

我不熟悉其他专门为神经网络设计的包，但我想自己写，主要是因为我打算做随机分析，这需要相当多的数学细节，我不知道那些包提供了这样的细节。

【问题讨论】：

你考虑过使用带有python接口的NEURON模拟器吗？它处理数值积分性能问题，并让您专注于方程式。可以绘制任何变量，并且它带有内置的 HH 通道。当然，您可以使用自己的 DE 编写自己的通道机制。 neuron.yale.edu/neuron

您可能还对NEST Simulator 感兴趣，它被明确构建为在点神经元动力学方面速度很快，并且可以很好地扩展到最大的网络。输入是神经元和突触的方程，连接性，输出是脉冲序列。 Python 界面，HH 已经可用，可定制，...对于分析，您可能对 Elephant 感兴趣，其中包含许多用于统计分析的构建块。

【参考方案1】：

分析表明您的大部分时间都花在了这两行：

    if sum(delta_k) > 0:

和

    if sum(delta_m) > 0:

将这些更改为：

    if np.any(...)

将一切速度提高 10 倍。如果您想逐行进行更多分析，请查看 kernprof： https://github.com/rkern/line_profiler

【讨论】：

谢谢你。这无疑有助于加快计划！

【参考方案2】：

作为对 welch 回答的补充，您可以使用 scipy.integrate.odeint 来加速集成：替换

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

通过

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

在我的计算机上将计算速度提高了 10 倍以上。

【讨论】：

所以通常我使用 odeint，但是当我有一堆 delta 函数输入时，我对使用它有些担心（我是数学家，所以我们从不使用 delta 函数，但物理学家会这样做我被困在使用它们）。当我通过打印时间来查看时间积分时，它没有使用 dt = 0.01 我不确定这是否会导致问题（实际上，我什至不确定“尖峰”是否出现了...似乎通过将成功概率设置为 1 表示兴奋和 0 表示抑制，没有尖峰......这很奇怪）

是的，我担心可能存在集成问题。 odeint 中有一个 hmaxparameter 可能可以使用，但也可能导致一些奇怪的行为。

【参考方案3】：

如果您有 NVidia grpahics 板，您可以使用 numba/numbapro 来加速您的 Python 代码并达到实时 4K 神经元，每个神经元有 128 个突触前神经元。

【讨论】：

这将有助于扩展这个答案，也许有一个融合 numba 库的小例子。