はじめに
Wolfram Mathematicaは世界的に人気な数式処理ソフトの一つです。特に教育機関で広く普及しており、理系出身の方には大学から配られたライセンスでMathematicaを使ったことがあるという方も多いのではないでしょうか。かくいう私も学生時代にたびたび使用していました。MathematicaではEuler法やRunge–Kutta法などの数値解析知識がなくても直感的な記述で微分方程式を解いたり、可視化することができます。微分方程式を解くには単にDSolveと記述するだけです。
これは二階線形微分方程式$y''=y$の一般解を導くMathematicaのコード（Wolfram Language）です。ここでは条件を与えずに二階の微分方程式を解いてもらったので未知数$C_1$と$C_2$が登場していますが、初期値や境界条件を与えることもできます。つまり、Mathematicaではもっと複雑なことができます。
そこで、今回は脳細胞の発火（脳細胞が隣の脳細胞に電気信号を送る現象を発火といいます）モデルであるFitzHugh-Nagumoモデルを題材にWolfram Mathematicaによる非線形ダイナミクスの可視化に取り組みます。
Hodgkin–Huxleyモデルとは
FitzHugh-Nagumoモデルへ取り掛かる前に、その原型となったHodgkin–Huxleyモデル（J Physiol. 1952 Aug 28;117(4):500–544.）を紹介だけします。次に示すようにHodgkin–Huxleyモデルはかなり煩雑です。
\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*}
ここでは詳細に立ち入りませんが、これは実験的に検証され電気生理学の礎ともなった歴史的にも意味のある式です。記号$V$であらわされているのは電圧です。
生物の細胞表面にはチャネルたんぱく質というものがあり、それがイオンを取り込むことで細胞の内外の物質輸送を制御しています。特に神経細胞ではチャネルたんぱく質によってイオンが流入すると細胞内の電気伝導度が変化するという点が重要です。この式の一行目ではカリウムイオンの取り込みを行う$\mathrm{K}^+$チャネルと、ナトリウムイオンの取り込みを行う$\mathrm{Na}^+$チャネルの神経細胞における働きが記述されています。
このモデルの本質を失わない形で2次元に簡略化されたのが今回扱うFitzHugh-Nagumoモデルです。
FitzHugh-Nagumoモデルとは
FitzHugh-Nagumoモデルは神経細胞などの電気的興奮性細胞の活動電位を表現したモデルの一種です。定義は次の微分方程式で与えられます。
\begin{align*} \frac{dv}{dt} &= c\left(v-\frac{v^3}{3}-w+I\right)\\
\frac{dw}{dt} &= v-bw+a \end{align*}
ここで、$v$は膜電位です。膜電位を可視化するには、Mathematicaのノートブックに
sol = NDSolve[
{v'[t] == v[t] - v[t]^3/3 - w[t] + 0.5, (* dv/dt *)
w'[t] == 0.08 (v[t] + 0.7 - 0.8 w[t]), (* dw/dt *)
v[0] == -1, (* 初期条件 v(0) *)
w[0] == 1 (* 初期条件 w(0) *)},
{v, w}, {t, 0, 200}];
Plot[Evaluate[v[t] /. sol], {t, 0, 100}, AxesLabel -> {"t", "v(t)"}]
と書き、実行します。NDSolve[]は微分方程式を数値解析的に解くという意味で、その結果はここではsolに代入されています。あとはPlot[]でsolの内容を表示すればよいのです。次のような図が出力されたら成功です。
図を見ると周期的に電位が上昇（発火）していることがわかります。ただし、ここでは変数や初期条件を適当に設定しています。パラメータを変化させればこの形も変化します。
このような微分方程式はC言語でRunge–Kutta法をコーディングして微分方程式を数値解析的に解き、gnuplotで可視化するという方法もありますが複雑です。一方でMathematicaであれば非常に簡単なため、即時可視化したい、実験の当たりをあらかじめつけておきたいというときには便利です。
ヌルクライン
ここまでFitzHugh-Nagumoモデルの$v$にだけ着目してきましたが、ここからは$w$にも着目し、全体の軌道を考えていきましょう。
ダイナミクスを考えるうえで特に重要なのはヌルクラインです。ヌルクラインは微分が0になる領域です。
FitzHugh-Nagumoモデルのヌルクラインは手計算でもすぐに求めることができるので先に確かめておきましょう。モデルの二式の左辺にそれぞれ0を代入すれば、
\begin{align*}
w &= v-\frac{v^3}{3}+I\\
w &= \frac{v}{b}+\frac{a}{b}
\end{align*}
が出てきます。これがヌルクラインです。三次関数と一次関数です。
では、ここから実際に可視化の手続きをしていきましょう。先ほどのsolに格納された微分方程式の解をそのまま使うので、同じノートブックに続けて次のように書いて実行してください。
traj = ParametricPlot[
Evaluate[{v[t], w[t]} /. sol], {t, 0, 100},
PlotStyle -> Blue, PlotRange -> {{-3, 3}, {-3, 3}}, AspectRatio -> 1];
(* ヌルクライン *)
vNull = Plot[v - v^3/3 + 0.5, {v, -3, 3},
PlotStyle -> Red, PlotRange -> {{-3, 3}, {-3, 3}}, AspectRatio -> 1];
wNull = Plot[(v + 0.7)/0.8, {v, -3, 3},
PlotStyle -> Green, PlotRange -> {{-3, 3}, {-3, 3}}, AspectRatio -> 1];
Show[traj, vNull, wNull,
PlotRange -> {{-3, 3}, {-3, 3}},
AspectRatio -> 1,
AxesLabel -> {"v", "w"}]
するとヌルクラインが表示されます。
青線は実際のニューロンの発火軌道を意味します。そして、赤線$dv/dt=0$のヌルクライン、緑線は$dw/dt=0$のヌルクラインです。途切れているように見える青線の端の部分は初期状態です。軌道が次第に周期的なアトラクタに吸引されていく様子が分かります。
もっと知りたい人へ
ヌルクラインだけではなく、より詳細なダイナミクスの特徴を知りたいという場合はベクトル場を見るというもの一つの方法です。これまでの結果をそのまま使うので同じノートブックに次のように打ち込んでみてください。
(* ベクトル場の支配方程式 *)
f[v_, w_] := v - v^3/3 - w + 0.5;
g[v_, w_] := 0.08 (v + 0.7 - 0.8 w);
(* ベクトル場 *)
vec = StreamPlot[{f[v, w], g[v, w]}, {v, -3, 3}, {w, -3, 3},
StreamStyle -> Gray, AspectRatio -> 1];
Show[vec, traj, vNull, wNull,
PlotRange -> {{-3, 3}, {-3, 3}},
AspectRatio -> 1,
AxesLabel -> {"v", "w"}]
美しいベクトル場を可視化することができました。
参考文献
[1]ストロガッツ『非線形ダイナミクスとカオス』
[2]The FitzHugh-Nagumo (FHN) model
[3]川平友規『レクチャーズ オン Mathematica』