实现 Izhikevich 神经元模型

Question

我正在尝试实现 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 并且 u 重置为 u + d .

鉴于这样一个简单的等式，我预计不会有任何问题。但是，虽然图表应该看起来像 ，但我得到的只是：

我完全不知道自己做错了什么，因为几乎没有什么可做的。我一直在寻找其他实现，但我正在寻找的代码总是隐藏在某个 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。非常感谢您帮助我指出这一点。

Answer 1

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