实施 Izhikevich 神经元模型

Implementing the Izhikevich neuron model

我正在尝试实现 Izhikevich 模型的脉冲神经元。这种神经元的公式非常简单:

v[n+1] = 0.04*v[n]^2 + 5*v[n] + 140 - u[n] + I

u[n+1] = a*(b*v[n] - u[n])

其中 v 是膜电位,u 是恢复变量。

如果 v 超过 30,则重置为 c 并重置 uu + d.

鉴于如此简单的方程式,我预计不会有任何问题。但是虽然图表应该看起来像 , all I'm getting is this:

我完全不知道自己做错了什么,因为几乎没有什么可做错的。我一直在寻找其他实现,但我正在寻找的代码总是隐藏在某处的 dll 中。但是我很确定我正在做的正是作者 (2) 的 Matlab 代码所做的。这是我的完整 R 代码:

v = -70
u = 0
a = 0.02
b = 0.2
c = -65
d = 6
history <- c()

for (i in 1:100) {
  if (v >= 30) {
    v = c
    u = u + d
  }
  v = 0.04*v^2 + 5*v + 140 - u + 0
  u=a*(b*v-u); 
  history <- c(history, v)
}

plot(history, type = "l")

对于任何曾经实施过 Izhikevich 模型的人,我错过了什么?

有用的链接: (1) http://www.opensourcebrain.org/projects/izhikevichmodel/wiki (2) http://www.izhikevich.org/publications/spikes.pdf

回答

原来是我看错公式了。显然 v' 表示新的 v = v + 0.04*v^2 + 5*v + 140 - u + I。我的老师会写这是 v' = 0.04*v^2 + 6*v + 140 - u + I。我非常感谢你的帮助向我指出这一点。

看看下面在 R 中实现 Izhikevich 模型的代码。它导致以下 R 图:

常规尖峰细胞:

喋喋不休的细胞:

R 代码:

# Simulation parameters
dt = 0.01 # ms
simtime = 500 # ms
t = 0

# Injection current
I = 15
delay = 100 # ms

# Model parameters (RS)
a = 0.02
b = 0.2
c = -65
d = 8

# Params for chattering cell (CH)
# c = -50
# d = 2

# Initial conditions
v = -80 # mv
u = 0

# Input current equation
current = function()
{
  if(t >= delay)
  {
    return(I)
  }

  return (0)
}

# Model state equations
deltaV = function()
{
  return (0.04*v*v+5*v+140-u+current())
}

deltaU = function()
{
  return (a*(b*v-u))
}

updateState = function()
{
  v <<- v + deltaV()*dt
  u <<- u + deltaU()*dt

  if(v >= 30) 
  {
    v <<- c
    u <<- u + d
  }
}

# Simulation code
runsim = function()
{
  steps = simtime / dt
  resultT = rep(NA, steps)
  resultV = rep(NA, steps)

  for (i in seq(steps))
  {
    updateState()
    t <<- dt*(i-1)

    resultT[i] = t
    resultV[i] = v
  }

  plot(resultT, resultV, 
       type="l", xlab = "Time (ms)", ylab = "Membrane Potential (mV)")
}

runsim()

一些注意事项:

  • 我从 Izhikevich's site 中为 "Regular Spiking (RS)" 单元格选择了参数。您可以从该页面上的两个 upper-right 图中选择其他参数。取消注释 CH 参数以获得 "Chattering" 类型单元格的绘图。
  • 正如评论者所建议的,问题中的前两个方程是错误实现的微分方程。实现第一个的正确方法类似于:"v[n+1] = v[n] + (0.04*v[n]^2 + 5*v[n] + 140 - u[n] + I) * dt"。例如,请参见上面的代码。 dt 是指用户指定的时间步积分变量,通常 dt << 1 ms.
  • 题中的for循环,应该先更新状态变量u和v,然后再检查条件。
  • 正如其他人所指出的,这两种细胞类型都需要电流源。我使用了 this page on the author's site 中的 15 个(我相信这些是微微放大器)(屏幕截图中 I 的底部值)。我还实现了当前开始的延迟(100 毫秒参数)。
  • 模拟代码应该实现某种时间跟踪,以便更容易知道结果图中出现尖峰的时间。上面的代码实现了这一点,并运行了 500 毫秒的模拟。