（著）山たー
FitzHugh-Nagumo モデル
FitzHugh-Nagumo(フィッツフュー・南雲)モデルは神経細胞の活動電位のダイナミクスをモデル化したものです。活動電位のモデル化としてはHodgkin–Huxleyモデルが最初ですが、Hodgkin–Huxleyモデルは4変数から成る非線形微分方程式であり、少し複雑です（次式）。
\begin{align*} I&=C_{m}{\frac {{\mathrm {d} }V_{m}}{{\mathrm {d} }t}}+{\bar {g}}_{\text{K}}n^{4}(V_{m}-V_{K})+{\bar {g}}_{\text{Na}}m^{3}h(V_{m}-V_{Na})+{\bar {g}}_{l}(V_{m}-V_{l})\\
\frac{dn}{dt}&=\alpha _{n}(V_{m})(1-n)-\beta _{n}(V_{m})n\\ \frac{dm}{dt}&=\alpha _{m}(V_{m})(1-m)-\beta _{m}(V_{m})m\\ \frac{dh}{dt}&=\alpha _{h}(V_{m})(1-h)-\beta _{h}(V_{m})h
\end{align*}
FitzHugh-NagumoモデルはHodgkin–Huxleyモデルを簡易にしたものとも言え、2変数の非線形微分方程式で記述されています。
\begin{align*} \frac{du}{dt} &= c\left(u-\frac{u^3}{3}-v+I\right)\\ \frac{dv}{dt} &= u-bv+a \end{align*} ここで$a,b,c$は定数であり、$a=0.7, b=0.8, c=10$がよく使われます。$u$は膜電位(membrane
potential)で、$v$はrecovery variableという隠れ変数です。 $I$は外部刺激電流に対応します。
FitzHugh-Nagumoモデルの詳細はScholarpediaなどに譲るとして、今回はFitzHugh-Nagumoモデルをmatplotlibを用いてアニメーションとして可視化することを目標とします。
膜電位変化を見る
コード
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate
I = 0.34 #external stimulus
a = 0.7
b = 0.8
c = 10
def FHN(state, t):
"""
FitzHugh-Nagumo Equations
u : the membrane potential
v : a recovery variable
"""
u, v = state
dot_u = c * (-v + u - pow(u,3)/3 + I)
dot_v = u - b * v + a
return dot_u, dot_v
#initial state
u0 = 2.0
v0 = 1.0
fig = plt.figure()
t = np.arange(0.0, 10, 0.01)
len_t = len(t)
dt = 5 #time steps
#animationの1step
def update(i):
global y, y0
#y0の初期値の設定
if i ==0:
y0 = [u0, v0]
plt.cla() #現在描写されているグラフを消去
#微分方程式を解く
y = integrate.odeint(FHN, y0, t)
#1Step(=dt)後のyの値を次のステップでのy0の値に更新する
y0 = (y[dt,0], y[dt,1])
#配列としてu,vを取得
u = y[:,0]
v = y[:,1]
#描画
plt.plot(t, u, label="u : membrane potential", color="#ff7f0e")
plt.plot(t, v, label="v : recovery variable", color="#1f77b4")
plt.plot(t[len_t-1], u[len_t-1],'o--', color="#ff7f0e") #uのmarker
plt.plot(t[len_t-1], v[len_t-1],'o--', color="#1f77b4") #vのmarker
plt.title("Membrane Potential / Volt")
plt.grid()
plt.legend(bbox_to_anchor=(0, 1),
loc='upper left',
borderaxespad=0)
ani = animation.FuncAnimation(fig, update, interval=100,
frames=300)
plt.show() #表示
ani.save("FitzHugh-Nagumo.mp4") #保存
結果
出力は次のようになりました。
出力はmp4なのですが、それをgifに変換しています。
相空間での軌道を見る
コード
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate
I = 0.34 #external stimulus
a = 0.7
b = 0.8
c = 10
def FHN(state, t):
"""
FitzHugh-Nagumo Equations
u : the membrane potential
v : a recovery variable
"""
u, v = state
dot_u = c * (-v + u - pow(u,3)/3 + I)
dot_v = u - b * v + a
return dot_u, dot_v
#initial state
u0 = 2.0
v0 = 1.0
fig = plt.figure()
t = np.arange(0.0, 1, 0.01)
t0 = np.arange(0.0, 20, 0.01)
#全ての軌道を表示するため
y_all = integrate.odeint(FHN, [u0, v0], t0)
u_all = y_all[:,0]
v_all = y_all[:,1]
len_t = len(t)
dt = 5 #time steps
#animationの1step
def update(i):
global y, y0
#y0の初期値の設定
if i ==0:
y0 = [u0, v0]
plt.cla() #現在描写されているグラフを消去
#微分方程式を解く
y = integrate.odeint(FHN, y0, t)
#1Step(=dt)後のyの値を次のステップでのy0の値に更新する
y0 = (y[dt,0], y[dt,1])
#配列としてu,vを取得
u = y[:,0]
v = y[:,1]
#描画
plt.plot(u_all, v_all, color="k", dashes=[1, 5])
plt.plot(u, v,color="r")
plt.plot(u[len_t-1],v[len_t-1],'o--', color="r") #uのmarker
plt.xlabel("u : membrane potential / Volt")
plt.ylabel("v : recovery variable")
plt.xlim([-2.5,2.5])
plt.ylim([-0.5,1.5])
plt.title("FitzHugh-Nagumo Phase Space")
plt.grid()
ani = animation.FuncAnimation(fig, update, interval=100,
frames=300)
plt.show() #表示
ani.save("FitzHugh-Nagumo_phase_space.mp4") #保存
y_allなど"_all"と付けた変数は下の結果における黒い点線を表示するためのものです。また、値域を固定するため、xlim,ylimを記述しています。
結果
これはリミットサイクルです。安定した軌道になるかは外部刺激電流Iの値によって変化します。閾値はI=0.34です。
膜電位変化と相空間を並べて見る
最後に2つを並べたアニメーションを作ってみます。
コード
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate
I = 0.34 #external stimulus
a = 0.7
b = 0.8
c = 10
def FHN(state, t):
"""
FitzHugh-Nagumo Equations
u : the membrane potential
v : a recovery variable
"""
u, v = state
dot_u = c * (-v + u - pow(u,3)/3 + I)
dot_v = u - b * v + a
return dot_u, dot_v
#initial state
u0 = 2.0
v0 = 1.0
t = np.arange(0.0, 5, 0.01)
#全軌道の表示用
t0 = np.arange(0.0, 20, 0.01)
y_all = integrate.odeint(FHN, [u0, v0], t0)
u_all = y_all[:,0]
v_all = y_all[:,1]
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10,4))
fig.suptitle("FitzHugh-Nagumo Model")
len_t = len(t)
dt = 5 #time steps
#animationの1step
def update(i):
global y, y0
#y0の初期値の設定
if i ==0:
y0 = [u0, v0]
#現在描写されているグラフを消去
ax1.cla()
ax2.cla()
#微分方程式を解く
y = integrate.odeint(FHN, y0, t)
#1Step(=dt)後のyの値を次のステップでのy0の値に更新する
y0 = (y[dt,0], y[dt,1])
#配列としてu,vを取得
u = y[:,0]
v = y[:,1]
#Phase Space
ax1.plot(u_all, v_all, color="k", dashes=[1, 6])
ax1.plot(u[len_t-20:len_t-1], v[len_t-20:len_t-1],color="r")
ax1.plot(u[len_t-1],v[len_t-1],'o--', color="r") #uのmarker
ax1.set_xlabel("u : membrane potential / Volt")
ax1.set_ylabel("v : recovery variable")
ax1.set_xlim([-2.2,2.2])
ax1.set_ylim([-0.5,1.5])
ax1.set_title("Phase Space")
ax1.grid()
#Membrane Potential
ax2.plot(t, u, label="u : membrane potential", color="#ff7f0e")
ax2.plot(t, v, label="v : recovery variable", color="#1f77b4")
ax2.plot(t[len_t-1], u[len_t-1],'o--', color="#ff7f0e") #uのmarker
ax2.plot(t[len_t-1], v[len_t-1],'o--', color="#1f77b4") #vのmarker
ax2.set_title("Membrane Potential / Volt")
ax2.set_ylim([-2.2,2.0])
ax2.grid()
ax2.legend(bbox_to_anchor=(0, 1),
loc='upper left',
borderaxespad=0)
ani = animation.FuncAnimation(fig, update, interval=100,
frames=300)
#plt.show() #表示
ani.save("FitzHugh-Nagumo_all.mp4") #保存
結果
・I=0.34の場合
・I=0.33の場合
刺激電流Iが閾値0.34を下回る時はリミットサイクルとならないことが分かります。今回は刺激電流Iを定数にしましたが、ステップ関数など変数にしても面白そうです。
参考にしたサイト
・Hodgkin–Huxley model - Wikipedia
・pythonを使ってFitzHugh南雲方程式のnullclineを描く - 技術メモ
・animation — Matplotlib 2.2.2 documentation
非線形物理は神経科学において重要なのは分かっているのですが、今一つ自分の研究と結びつかない感があります。
コメントをお書きください