神経細胞モデルを実装してみよう #Python3

ここでいう神経細胞モデルとは、神経細胞の発火の時間変化を表現可能なモデリングのことを指します。このページは
[1] E.M.Izhikevich (2004) "Which Model to Use For COrtical Spiking Neurons?", IEEE transaction on neural networks,Vol.15, No.5,pp.1063-1070
[2] https://neuronaldynamics.epfl.ch/index.html
に基づきます．

以下のように文字を定めます。
$u,v,w$：微分方程式で表される変数,膜電位(参考にした文献によって微妙に文字が違います)
$I$：入力電流
$R$：抵抗値
$\tau$ ：時定数

Hodgkin-Huxley model

神経細胞のモデルとして、最も精度が高くて複雑なものが、このHHモデルです。神経科学の世界では最も重要な内容の１つです。膜電位・やイオンチャンネルについての４つの微分方程式で表現されます。
式の詳細はこちらを参照。実装は長いので略。
https://bsd.neuroinf.jp/wiki/Hodgkin-Huxley%E6%96%B9%E7%A8%8B%E5%BC%8F

Fitzhugh-Nagumo Model

HHモデルは4つの微分方程式で構成されており複雑です。そこで，実際の細胞のダイナミクスに基づいて，無理のない範囲で近似することで，2変数に減らしました。その1例がFNモデルです。2変数になったことで，2次元平面での力学系として解析することが可能になりました。

\tau\frac{dv}{dt} = v-\frac{1}{3}v^3 - w + RI(t) \\
\tau_w\frac{dw}{dt} = a + bv - w

import numpy as np
import matplotlib.pyplot as plt

def fitzhugh_nagumo(tau,tau_w,a,b,v_init,u_init,I):
    dt = 0.01
    v = np.zeros(np.size(I))
    u = np.zeros(np.size(I))
    R=1
    v[0] = v_init
    u[0] = u_init
    for i in range(np.size(I)-1):
        v[i+1] = v[i] + (v[i] - v[i]*v[i]*v[i]/3 -u[i]+R*I[i])*dt/tau
        u[i+1] = u[i] + (v[i] - b*u[i] +a)*dt/tau_w
    return v,u

v1,u1 = fitzhugh_nagumo(0.1,0.1,0.7,0.8,-5,-1,np.concatenate([np.ones(10)*2,np.zeros(490)])) #large pulse input
v2,u2 = fitzhugh_nagumo(0.1,0.2,0.7,0.8,-5,-1,np.ones(500)) #constant input

dt = 0.01
x = np.linspace(0,500,500)*dt
plt.subplot(2,1,1)
plt.plot(x,v1,label="large pulse input")
plt.plot(x,v2,label="constant input")
plt.legend()

plt.subplot(2,1,2)
plt.plot(v1,u1)
plt.plot(v2,u2)
plt.show()

上記コードは，定常入力時とパルス入力時の応答を示します。
u-v平面の図も表示でき，定常入力時はリミットサイクルが発生することがわかります。

Integrate-and-Fire model(積分発火モデル)

HNモデルよりも更に変数が減って１つしかないため，計算コストが非常に小さいというメリットがありますが，あまりちゃんとした発火パターンをモデリングすることができません。

\begin{align}
\tau \frac{dv}{dt} &= -(v-v_{rest}) + RI(t)\\
&if\quad v \geq v_{thres} \quad v=c
\end{align}

import numpy as np
import matplotlib.pyplot as plt

def IandF(tau,v_init,v_thresh,v_reset,I):
    dt = 0.001
    v = np.zeros(np.size(I))
    v[0] = v_init
    R=5
    for i in range(0,np.size(v)-1):
        v[i+1] = v[i] + ((-(v[i] -v_init) + R*I[i]) * dt)/tau
        if v[i+1]>= v_thresh:
            v[i+1] = v_reset
    return v

dt = 0.01
x = np.linspace(0,1000,1000)*dt
v = IandF(1,-70,-50,-65,np.ones(1000)*10);

plt.plot(x,v)
plt.legend()
plt.show()

Izhikevich model

比較的少ない計算量で，かつ高精度な結果を出すことができるのがこのモデルとされています。u,vに関わる２つの微分方程式で表現されます。

\begin{align}
\dot{v} &= 0.04v^2 + 5v + 140 -u + I \\
\dot{u} &= a(bv-u) \\
\\
if &\quad v \geq 30(mV) \\
 v &= c \\
 u &= u+ d \\
\end{align}

import numpy as np
import matplotlib.pyplot as plt

def izhikevich(u_init,v_init,I):
    dt = 0.01
    u = np.zeros(np.size(I))
    v = np.zeros(np.size(I))
    u[0]= u_init;
    v[0]= v_init;

    a = 0.02*10
    b = 0.2
    c = -65
    d = 2

    for i in range(0,np.size(I)-1):
        v[i+1] = v[i] + (0.04*v[i]*v[i] + 5*v[i] + 140 -u[i] +I[i])*dt*10
        u[i+1] = u[i] + a*(b*v[i]-u[i])*dt

        if v[i+1]>=30:
            v[i+1] = c
            u[i+1] = u[i+1]+d
    return u,v

u,v = izhikevich(30,-80,np.ones(1000)*10)

dt = 0.01
x = np.linspace(0,1000,1000)*dt
plt.plot(x,v)
plt.show()