Notebook

Izhikevich model with DifferentialEquations.jl¶

Reference¶

ゼロから作る Spiking Neural Networks (第2版)
https://booth.pm/ja/items/1585421

DifferentialEquations.jl 公式ドキュメント
https://docs.sciml.ai/stable/

In [1]:

using DifferentialEquations
using Plots

まずは微分方程式をあらわす関数を定義するところから始める。

チュートリアルの Defining Parameterized Functions と同じことをやる。 https://docs.sciml.ai/stable/tutorials/ode_example/#Defining-Parameterized-Functions-1

時刻 t での変数の値が u という配列でまとめて与えられる。ここでは、

膜電位
回復電流
直近の発火時刻

の3つの変数を u[1], u[2], u[3] に対応させる。発火時刻は計算には必要ないが、後の解析で便利なので付け加えておいた。そして、それぞれの変数の時刻 t における時間微分を du という配列の要素に書き込めばよい。

必要なパラメータは p で受け渡す。基本的に何でも入れることができて、ここでは入力電流 I(t) を関数として渡している。

In [2]:

function iz!(du, u, p, t)
    (C, a, b, k, vᵣ, vₜ), I = p
    v, i, _ = u
    
    du[1] = C \ (k * (v - vᵣ) * (v - vₜ) - i + I(t))
    du[2] = a * (b * (v - vᵣ) - i)
    du[3] = 0
end

params = (
    C = 100,
    a = 0.03,
    b = -2.0,
    k = 0.7,
    vᵣ = -60,
    vₜ = -40
)

Out[2]:

(C = 100, a = 0.03, b = -2.0, k = 0.7, vᵣ = -60, vₜ = -40)

次に、膜電位のリセットを実装する。

DiffEq にはコールバックという機能があり、条件が満たされるとイベントを発生させて、変数の値を書き換えたりすることができる。ボールが床で跳ね返るサンプルがわかりやすい。 https://docs.sciml.ai/stable/features/callback_functions/#Example-1:-Bouncing-Ball-1

コールバックには2つの関数が必要で、イベントの発生条件と発生時の動作を記述する。

condition(u, t, integrator) がゼロを返すときにイベントが発生する。Iz モデルの場合はピーク電圧への到達である。

affect!(integrator) でイベントを記述する。integrator.u や integrator.t で変数にアクセスできる。Iz モデルでやるべきことは、

電圧のリセット
発火による電流増加
発火時刻の記録（任意）

である。

In [3]:

condition(u, t, integrator) = u[1] - 35

function affect!(integrator)
    integrator.u[1] = -50
    integrator.u[2] += 100
    integrator.u[3] = integrator.t
end

cb = ContinuousCallback(condition, affect!)

Out[3]:

ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}(condition, affect!, affect!, DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0)

アルゴリズムは、とりあえずおすすめされている Tsit5() を使う。
https://docs.sciml.ai/stable/solvers/ode_solve/#Non-Stiff-Problems-1

ただし、後で示すプロットでの t=350 の微妙な膨らみは Tsit5() 特有のもので、アルゴリズムを Vern7() などに変更すると消滅する。

計算結果は sol を操作することでアクセスできる。
https://docs.sciml.ai/stable/basics/solution/

例えば、電圧 sol[1, :]、時刻 sol.t など。

In [5]:

I(t) = ifelse(50 ≤ t ≤ 350, 100, 0)
tspan = (0.0, 400.0)
u0 = [-60.0, 0.0, 0.0]
alg = Tsit5() # Vern7()

prob = ODEProblem(iz!, u0, tspan, (params, I))
@time sol = solve(prob, alg; callback=cb)
spiketimings = unique(sol[3, :])[2:end]

  3.030972 seconds (4.53 M allocations: 209.335 MiB, 2.95% gc time)

Out[5]:

4-element Array{Float64,1}:
  98.18504046749854
 171.63601907655283
 247.75605552228132
 323.78565753254725

DiffEq は Plots と連携していて、sol をそのままプロットできる。
https://docs.sciml.ai/stable/basics/plot/

固有のプロットオプションもいくつか存在する。

vars でプロットしたい変数を指定する。0が t、1はu[1] に対応する。

In [6]:

plot(sol;
    vars=(0, 1),
    plotdensity=20000,
    xlabel="Time (ms)",
    ylabel="Membrane potential (mV)",
    legend=false
)

plot!(t -> I(t) / 20 - 70, range(tspan...; length=1000))

scatter!(_ -> 40, spiketimings; shape=:dtriangle)

Out[6]:

In [8]:

function firerate(current)
    prob = ODEProblem(iz!, [-60.0, 0.0, 0.0], (0.0, 100000.0), (params, _ -> current))
    sol = solve(prob, Tsit5(); callback=cb)
    ts = unique(sol[3, :])[2:end]
    length(ts) ≥ 2 ? 1000 / (ts[end] - ts[end-1]) : 0.0
end

currents = range(50, 55; length=200)
@time frequencies = firerate.(currents)

scatter(currents, frequencies;
    xlabel="Input current (pA)",
    ylabel="Firing rate (Hz)",
    xrange=extrema(currents),
    yrange=extrema(frequencies),
    markersize=2,
    legend=false
)

  7.849069 seconds (17.99 M allocations: 1.274 GiB, 4.86% gc time)

Out[8]: