書記が数学やるだけ#901 神経細胞の数理モデル
神経細胞の数理モデルのうち，ホジキン-ハクスリーモデルとフィッツフュー-南雲モデルについて中身を確認する。
問題
説明
ホジキン–ハクスリーモデルは、1952年にA.L.ホジキンとA.F.ハクスリーによって提案された、神経細胞膜における活動電位の発生を記述する生理学的モデルである。このモデルは、ナトリウム（Na⁺）およびカリウム（K⁺）イオンの流れと、それを制御する電位依存性ゲート変数を含む4つの微分方程式から構成される。これにより、実際のニューロンにおける活動電位の立ち上がり・再分極過程を高精度に再現できるが、数理的には非線形で複雑な構造をもつ。
フィッツヒュー–南雲方程式は、ホジキン–ハクスリーモデルを数理的に簡略化したものである。活動電位を表す速い変数 v（膜電位）と、回復過程を表す遅い変数 w（回復変数）の2変数系で構成される。
解答
ホジキン–ハクスリーモデルにおいて、膜電流はナトリウム（Na⁺）チャネル電流、カリウム（K⁺）チャネル電流、および漏れ電流から構成される。これらの各電流は、実験データに基づいてフィッティングされた経験的関数として表されており、チャネルの開閉過程を反応速度論的に記述している。具体的には、各イオンチャネルの開確率・閉確率を示すゲーティング変数が導入され、それらの時間発展を1次の常微分方程式で表すことにより、チャネルの動態と膜電位変化の関係が定式化されている。このようにしてホジキンとハクスリーは、膜電位とイオンチャネルの反応速度論を結びつける形で、神経興奮の電気的挙動を定量的に再現する数理モデルを導出した。
ホジキン–ハクスリーモデルの数値シミュレーションを以下に示す。外部電流をパルス状に与えると、膜電位が上昇し、ナトリウムチャネルが開口することでさらに膜電位が上昇し、脱分極が生じる。その後、カリウムチャネルの開口とナトリウムチャネルの閉口によって膜電位は低下し、過分極を経て静止電位へと戻る。
import numpy as np
import matplotlib.pyplot as plt
#------------------------------------------
# Hodgkin–Huxleyモデルの定義
#------------------------------------------
C_m = 1.0 # 膜容量 (µF/cm^2)
g_Na = 120.0 # Na+ 最大コンダクタンス (mS/cm^2)
g_K = 36.0 # K+ 最大コンダクタンス (mS/cm^2)
g_L = 0.3 # リークコンダクタンス (mS/cm^2)
E_Na = 50.0 # Na+ 平衡電位 (mV)
E_K = -77.0 # K+ 平衡電位 (mV)
E_L = -54.387 # リーク平衡電位 (mV)
#------------------------------------------
# 外部電流（パルス型）
#------------------------------------------
def I_ext(t):
# 10～40 msの間のみ10 µA/cm^2の外部電流を加える
return 10.0 if 10.0 <= t <= 11.0 else 0.0
#------------------------------------------
# α, β関数の定義
#------------------------------------------
def alpha_m(V): return 0.1 * (V + 40.0) / (1.0 - np.exp(-(V + 40.0) / 10.0))
def beta_m(V): return 4.0 * np.exp(-(V + 65.0) / 18.0)
def alpha_h(V): return 0.07 * np.exp(-(V + 65.0) / 20.0)
def beta_h(V): return 1.0 / (1.0 + np.exp(-(V + 35.0) / 10.0))
def alpha_n(V): return 0.01 * (V + 55.0) / (1.0 - np.exp(-(V + 55.0) / 10.0))
def beta_n(V): return 0.125 * np.exp(-(V + 65) / 80.0)
#------------------------------------------
# 時間ステップと初期条件
#------------------------------------------
T = 60.0 # シミュレーション時間 [ms]
dt = 0.01 # 時間刻み [ms]
time = np.arange(0, T + dt, dt)
V = np.zeros_like(time)
m = np.zeros_like(time)
h = np.zeros_like(time)
n = np.zeros_like(time)
V[0] = -65.0 # 静止膜電位 [mV]
m[0] = alpha_m(V[0]) / (alpha_m(V[0]) + beta_m(V[0]))
h[0] = alpha_h(V[0]) / (alpha_h(V[0]) + beta_h(V[0]))
n[0] = alpha_n(V[0]) / (alpha_n(V[0]) + beta_n(V[0]))
#------------------------------------------
# 数値積分（オイラー法）
#------------------------------------------
I_Na = np.zeros_like(time)
I_K = np.zeros_like(time)
I_L = np.zeros_like(time)
I_ext_trace = np.zeros_like(time)
for i in range(1, len(time)):
t = time[i]
I = I_ext(t)
I_ext_trace[i] = I
# イオン電流
I_Na[i] = g_Na * (m[i-1] ** 3) * h[i-1] * (V[i-1] - E_Na)
I_K[i] = g_K * (n[i-1] ** 4) * (V[i-1] - E_K)
I_L[i] = g_L * (V[i-1] - E_L)
# 膜電位の変化
dVdt = (I - I_Na[i] - I_K[i] - I_L[i]) / C_m
V[i] = V[i-1] + dt * dVdt
# ゲーティング変数の更新
m[i] = m[i-1] + dt * (alpha_m(V[i-1]) * (1 - m[i-1]) - beta_m(V[i-1]) * m[i-1])
h[i] = h[i-1] + dt * (alpha_h(V[i-1]) * (1 - h[i-1]) - beta_h(V[i-1]) * h[i-1])
n[i] = n[i-1] + dt * (alpha_n(V[i-1]) * (1 - n[i-1]) - beta_n(V[i-1]) * n[i-1])
#------------------------------------------
# プロット
#------------------------------------------
plt.figure(figsize=(10, 10))
# (1) 膜電位
plt.subplot(4, 1, 1)
plt.plot(time, V, 'k')
plt.ylabel('V (mV)')
plt.ylim(-100, 50)
# (2) 外部電流
plt.subplot(4, 1, 2)
plt.plot(time, I_ext_trace, 'r')
plt.ylabel('I_ext (µA/cm²)')
plt.ylim(-1, 11)
# (3) ゲーティング変数
plt.subplot(4, 1, 3)
plt.plot(time, m, label='m')
plt.plot(time, h, label='h')
plt.plot(time, n, label='n')
plt.ylabel('gating')
plt.legend()
plt.tight_layout()
plt.show()
継続的な外部電流を与えた場合、その強さによって膜電位の挙動が変化する。外部電流が小さい場合には脱分極が起こらないが、電流を増加させると脱分極が生じるようになり、さらに強い電流を与えると膜電位は周期的な振動パターンを示すようになる。
def I_ext(t):
return 1.0 if 10.0 <= t <= 50.0 else 0.0
def I_ext(t):
return 5.0 if 10.0 <= t <= 50.0 else 0.0
def I_ext(t):
return 10.0 if 10.0 <= t <= 50.0 else 0.0
フィッツヒュー–南雲方程式において、外部電流が0の場合、パラメータの値に応じて、系は興奮性をもつ1つの安定平衡点をとる場合と、1つの不安定平衡点と2つの安定平衡点をもつ場合のいずれかになる。このとき、v > a となるv軸上の点を初期値として与えると、相軌道は大きく迂回して原点へ戻る軌跡を描き、これは神経細胞における一過的な興奮応答に対応する。
一方、外部電流を加えた場合には、電流の大きさがある範囲にあるとき、系はリミットサイクル振動を示し、周期的な発火（スパイク）を再現することができる。
参考：
相平面をシミュレーションによって表示する。まず、外部電流が存在しない場合について、一過性の興奮応答を示す。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
a, b, eps = 0.7, 0.9, 0.08
Iext = 0
# 時間設定
dt = 0.01
T = 200
steps = int(T/dt)
# 変数の初期値
v = np.zeros(steps)
w = np.zeros(steps)
v[0] = 0
w[0] = 1.0
# 関数定義
def f(v, w):
dv = v - v**3 / 3 - w + Iext
dw = eps * (v + a - b * w)
return dv, dw
# 時間発展（オイラー法）
for i in range(steps-1):
dv, dw = f(v[i], w[i])
v[i+1] = v[i] + dv * dt
w[i+1] = w[i] + dw * dt
# ヌルクライン計算
v_null = np.linspace(-2, 2, 400)
w_null_v = v_null - v_null**3/3 + Iext
w_null_w = (v_null + a)/b
# 相平面のベクトル場
V, W = np.meshgrid(np.linspace(-2, 2, 30), np.linspace(-1, 3, 30))
DV, DW = f(V, W)
# ===== 図1: 相平面 =====
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.streamplot(V, W, DV, DW, color='gray', density=1.0)
plt.plot(v_null, w_null_v, 'r', label='v-nullcline')
plt.plot(v_null, w_null_w, 'b', label='w-nullcline')
plt.plot(v, w, 'k', linewidth=2, label='trajectory')
plt.xlabel('v')
plt.ylabel('w')
plt.legend()
plt.title('Phase Plane (FitzHugh–Nagumo)')
plt.grid(True)
# ===== 図2: 時系列 =====
time = np.arange(steps) * dt
plt.subplot(1,2,2)
plt.plot(time, v, label='v(t)', color='r')
plt.plot(time, w, label='w(t)', color='b')
plt.xlabel('Time')
plt.ylabel('Variables')
plt.title('Time evolution of v and w')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
次に、外部電流を加えた場合について、周期的な振動（持続的興奮）が生じる様子を示す。
a, b, eps = 0.7, 0.8, 0.08
Iext = 0.5
本記事のもくじはこちら：
