(45) Hodgkin-Huxley系の伝搬速度 - 読み書きプログラミング

ニューロンの活動電位の現象論的方程式Hodgkin-Huxley系を用いた軸索上の電位の伝搬速度について考えます。
出発点は、Hodgkin-Huxley系とケーブル理論になります。
Hodgikin-Huxley系はイオンチャンネルを加味した細胞膜の活動電位の方程式系です。

$\begin{align*} I_{ext}(t)&=c_m\frac{de}{dt}+g_Kn(t)^4(e-e_K)+g_{Na}m(t)^3h(t)(e-e_{Na})+g_l(e-e_l)\\ \frac{dm}{dt}&=a_m(e)(1-m(t))-b_m(e)m(t)\\ \frac{dh}{dt}&=a_h(e)(1-h(t))-b_h(e)h(t)\\ \frac{dn}{dt}&=a_n(e)(1-n(t))-b_n(e)n(t) \end{align*}$

$e(t)$ : 膜電位
$m$ : Naイオンチャネル活性化変数
$h$ : Naイオンチャネル不活性化変数
$n$ : カリウムイオンチャンネル開確率
$c_m$ : (単位長)静電容量
$e_k,e_{Na},e_l$ : カリウム、ナトリウム、漏洩平衡電位
$g_k,g_{Na},g_l$ : カリウム、ナトリウム、漏洩ピークコンダクタンス

Hodgkin-Huxleyモデルを線上に並べ、外部電流 $l_{ext}$ で相互作用させるモデルを考えることができます。このような拡張の一般論はケーブル理論として知られています。ケーブル理論では、表面電流はオームの法則から以下のように表されます。

$I_{ext}(x,t)=\frac{1}{r_a}\frac{\partial e}{\partial x}$

$r_a$ : 神経内の単位長抵抗

また、Kirchhoff'の法則から外部電流の差(位置微分)が内部電流になります。

$\frac{\partial I_{ext}}{\partial x}=I_{int}$

Hodgkin-Huxleyモデルとケーブル理論を組み合わせると、軸索上の膜電位の伝搬方程式を導出することができます。

$\frac{1}{r_a}\frac{\partial^2e}{\partial x^2}=c_m\frac{\partial e}{\partial t}+g_Kn(x, t)^4(e(x, t)-e_K)+g_{Na}m(x, t)^3h(x, t)(e(x,t)-e_{Na})+g_l(e(x,t)-e_l)$

4変数の二階の偏微分方程式が導かれましたが、今回は $e(x,t)=e(x-ut)$ という形の伝搬波形を仮定します。
ここで、 $u$ は伝搬速度を示す定数です。

すると、偏微分方程式は以下とイオンチャネルの式を合わせた、4変数常微分方程式になります。

$\frac{1}{r_au^2}\frac{d^2e}{dt^2}=c_m\frac{de}{dt}+g_Kn^4(e-e_K)+g_{Na}m^3h(e-e_{Na})+g_l(e-e_l)$

イオンチャネルの現象論的仮定を組み込んで、伝搬速度を求めましょう。膜電位eを発散させない $vi=r_au^2$ を二分木探索で調べます。

/* 7.4.4m */
cm : 1.0$
gk : 36$
gna : 120$
gl : 0.3$
ek : -12.0$
ena : 115.0$
el : 10.5988$
am(V) := 0.1*(25 - V)/(exp(0.1*(25 - V)) - 1)$
bm(V) := 4*exp(-V/18.0)$
ah(V) := 0.07*exp(-V/20)$
bh(V) := 1/(exp(0.1*(30 - V)) + 1)$
an(V) := 0.01*(10.0 - V)/(exp(0.1*(10 - V)) - 1)$
bn(V) := 0.125*exp(-V/80.0)$
data : makelist(0, i, 0, 100)$
dt : 0.001$
vlow : 9.169921875$
vhigh : 9.17041015625$
vipre : -1$
while true do (
  ct : 0,
  vi : (vlow + vhigh)/2,
  if vi = vipre then return() else vipre : vi,
  n : an(0)/(an(0) + bn(0)),
  m : am(0)/(am(0) + bm(0)),
  h : ah(0)/(ah(0) + bh(0)),
  e : 0.1,
  edot : 0,
  last_t : for t : 0 step dt thru 10 do (
    edotdot : vi*(cm*edot + gk*n^4*(e - ek) + gna*m^3*h*(e - ena) + gl*(e - el)),
    dm : am(e)*(1 - m) - bm(e)*m,
    dh : ah(e)*(1 - h) - bh(e)*h,
    dn : an(e)*(1 - n) - bn(e)*n,
    m : m + 3*dm*dt,
    h : h + 3*dh*dt,
    n : n + 3*dn*dt,
    edot : edot + edotdot*dt,
    e : e + edot*dt,
    ct : ct + 1,
    if e > 150 then (vhigh : vi, return(t)),
    if e < -50 then (vlow : vi, return(t)),
    if mod(ct, 100) = 0 then data[ct/100] : e
    ),
  print(last_t, " ", vhigh, " ", vlow),
  plot2d([discrete, makelist(i, i, 1, 50), makelist(data[i], i, 1, 50)])
  )$