4

我正在尝试实现 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 . 我非常感谢您帮助我指出这一点。

4

1 回答 1

2

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

常规加标单元

Izhikevich 正则尖峰 RS 细胞图

颤抖的细胞Izhikevich Chattering CH 细胞图

和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 的网站上选择了“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 << 1 ms。
  • 在问题的 for 循环中,应先更新状态变量 u 和 v,然后检查条件。
  • 正如其他人所指出的,这两种细胞类型都需要电流源。我在作者网站的这个页面上使用了 15 个(我相信这些是 pico amps) (截图中 I 的底部值)。我还实现了当前开始的延迟(100 ms 参数)。
  • 模拟代码应该实现某种时间跟踪,以便更容易知道结果图中何时出现尖峰。上面的代码实现了这一点,并运行了 500 毫秒的模拟。
于 2017-01-20T18:33:07.427 回答