数理生物学に入門してみる（３）振動と波

　今回は教科書の第3, 4章、授業では4~6, 10回目で扱った振動のモデル・時空間パターンについてまとめます。ホジキンハクスレー方程式やチューリングパターンについて触れます。

前回

cake-by-the-river.hatenablog.jp

振動性・興奮性
時空間パターン
- 拡散方程式と濃度勾配
- チューリングパターン
進行波のモデル
- 化学進行波と界面
- ケーブル方程式と進行波
補足（末尾）
- 拡散だけで進行波が起きないことの説明

振動性・興奮性

　細胞の中では周期的に量が変動する場合もあります。ここでは神経の活動電位を表すホジキン-ハクスレー方程式などを例に振動について考えます。

振動の条件

　化学反応で振動する例として有名なものに、Belousov-Zhabotinsky反応（通称BZ反応）があります。BZ反応はマロン酸をブロモ化する複雑な反応ですが、振動のエッセンスだけ取り出してきたモデルとして、ブラッセレータがあります。ブラッセレータは分子 X, Y が中間体として A, B (>0) という外部の成分をもとに時間変化する系です。

$A \to X$

$2X + Y \to 3X$

$B + X \to Y + C$

$X \to D$

の関係であり、時間変化は全反応速度定数を 1 とおくと

$\displaystyle \frac{d X}{dt} = A - (B+1)X + X^2Y = f(X, Y)$

$\displaystyle \frac{d Y}{dt} = BX - X^2Y = g(X, Y)$

となります。固定点は

$\displaystyle A - (B+1)X + X^2Y = 0, BX - AX^2Y = 0 \Rightarrow (X^*, Y^*) = (A, \frac{B}{A})$

と求まります。固定点でのヤコビ行列

$\displaystyle J = \begin{pmatrix} \frac{\partial \dot{X}}{\partial X} & \frac{\partial \dot{X}}{\partial Y} \\ \frac{\partial \dot{Y}}{\partial X} & \frac{\partial \dot{Y}}{\partial Y} \end{pmatrix} = \begin{pmatrix} B-1 & A^2 \\ -B & -A^2 \end{pmatrix}$

の固有値 $\lambda_1, \lambda_2$ は

$\left| \begin{array}{ccc} B-1-\lambda & A^2 \\ -B & -A^2-\lambda \end{array} \right| = \lambda^2 - (B - A^2 - 1)\lambda + A^2B = 0$

を満たすので

$\lambda_1 \lambda_2 = A^2B > 0$

$\lambda_1 + \lambda_2 = B - A^2 - 1$

という関係が求まります。すなわち、2固有値は同符号で、 $B > A^2 + 1$ のときに不安定、 $B < A^2 + 1$ のときは安定になることがわかります。

　　固定点が不安定な時、上記は振動反応になっていると考えられます。振動反応は各構成成分の量が周期的に変化する状態であり、それはリミットサイクルにあると考えられます。リミットサイクルとは、十分な時間がたった時に、（漸近的に）ある閉じた軌道に状態が落ち着く、そのときの軌道（振動解）を指します。実際にブラッセレータを相平面に表してみると、軌道に収束する様子が分かります。

f:id:cakkby6:20201128111554p:plain — ブラッセレータ：(A, B) = 左 (1, 1), 右 (1, 3)。ピンクが $f(X, Y)$ の、緑が $g(X, Y)$ のヌルクライン。右側ではリミットサイクルが生じている

ブラッセレータ：(A, B) = 左 (1, 1), 右 (1, 3)。ピンクが $f(X, Y)$ の、緑が $g(X, Y)$ のヌルクライン。右側ではリミットサイクルが生じている

　固定点は2つのヌルクラインの交点の位置に存在しており、リミットサイクルになっている場合、固定点の左側では $g(X, Y)$ のヌルクラインが、右側では $f(X, Y)$ のヌルクラインが上にくる形をとっています。ここで、各状態は２つのヌルクラインとの関係によって速度方向が切り替わることに注目します。つまり、ピンクの線と緑の線に囲まれた4つの領域で速度ベクトルの向きが切り替わり、リミットサイクルではそれが回転するように決まっているということになります。

Hodgkin-Huxley 方程式

　神経回路における興奮は細胞膜の電位の変化によってあらわされています。神経細胞には主に3種類のイオンチャネルが存在しており、それらによってカリウムイオンやナトリウムイオンなどの濃度差を生じさせ、活動電位＝興奮を生じます。これらは、コンデンサーや抵抗などの電気回路のモデルで考えることができます。具体的には、次のような回路（等価回路）で表現することが出来ます。

f:id:cakkby6:20201217113220p:plain — 等価回路：電流は黒丸の点に入る方向を正として考える。青矢印はイオンの流れる方向で、起電力の向きが対応している。

キルヒホッフの法則から

$I = I_{memb} + I_{Na} + I_K + I_L$

コンデンサーと抵抗に対してそれぞれ

$\displaystyle I = C\frac{dV}{dt}, I = g(V-E)$

が成り立つことから、膜全体の電圧と電流を用いて

$\displaystyle C_{memb}\frac{dV}{dt} = -g_{Na}(V - E_{Na}) - g_K(V - E_K) - g_L(V-E_L) + I$

と表せます。 $g$ はコンダクタンスで、イオンチャネルの開口度合を反映しているものといえます。

　このコンダクタンスはイオンチャネルに固有であり、電圧の値に応じて変化します。HodgkinとHuxleyはこれをボルテージクランプによる実験でカーブフィッティングし確定させることで、うまく神経回路のシミュレーションができるようになったそうです。具体的には、コンダクタンスを活性化パラメータ $m(V)$ と不活性化パラメータ $h(V)$ を用いて

$\displaystyle g_i (V) = G_i m_i^{p_i}(V) h_i^{q_i} (V)$

として考え（ $G_i$ は最大伝導率に相当し、チャネルの数などに対応しています）、それぞれのパラメータをフィッティングしたことになります。例えばNa+の活性化パラメータは

$\displaystyle \frac{dm_{Na}}{dt} = \alpha_{Na}^{(m)}(1-m_{Na}) - \beta_{Na}^{(m)}m_{Na}$

$\displaystyle \alpha_{Na}^{(m)} = \frac{0.1(25-V)}{e^{(25-V)/10}-1}$

$\displaystyle \beta_{Na}^{(m)} = 4e^{(-V/18)}$

として決まるそうです。

　K+チャネルでは $m_K$ のみ存在し、 Na+チャネルでは $m_{Na}, h_{Na}$ の二つが存在しているとして、結局 $V, m_K, m_{Na}, h_{Na}$ の4変数が絡む Hodgkin-Huxley 方程式が得られます（連立方程式で）。この方程式のそれぞれの変数の動きを見てみることで
①膜電位の変化の上昇からNa+チャネルが即座に開きだす（ $m_{Na}$ が上昇、 $h_{Na}$ が減少）
②遅れてK+チャネルが開きだす（ $m_K$ が上昇）
③Na+チャネルが閉じてくる（ $h_Na$ が上昇）
という動きが活動電位を生じていることが分かります。

Fitzhugh-Nagumo 方程式

　Hodgkin-Huxley方程式は4変数もあるので一目で見ても動きの理解が難しいです。そこで、なんとか性質を残しつつ変数を削ってみようと思います。実際には、 $m_{Na} = const.$ と $h_{Na} + m_{K} = const.$ という近似を導入して

$\displaystyle C_{memb}\frac{dV}{dt} = -G_{Na}m_{Na, \infty}^3(0.8-m_K)(V - E_{Na}) - G_Km_K^4(V - E_K) - g_L(V-E_L) + I$

$\frac{dm_K}{dt} = \alpha_K^{(m)}(1-m_K) - \beta_K^{(m)}m_K$

と表します。これで $m_K, V$ の2変数のみになり、これを相平面で確認すると、リミットサイクルが生じていることが分かります。

f:id:cakkby6:20201219160549p:plain — 縮約された相平面：Pattern Formation in Spatially Extended Systems: Interplay of Variability and Noise (E. Glatt 2008) より引用

　この2変数に縮約されたモデルを、さらに分かりやすいように簡単にしたモデルとして、 Fitzhugh-Nagumo方程式というものがあります。

$\displaystyle \epsilon \frac{dv}{dt} = v(a-v)(v-1) - m + I$

$\frac{dm}{dt} = bv-cm$

これを実際に相平面にしてみると、興奮性が入力 $I$ に応じて生じている様子が分かります。

f:id:cakkby6:20201219161039p:plain — Fitzhugh-Nagumo方程式：左から $I = 1.0, 1.3, 1.8$ の場合。ピンクが前者のヌルクラインで、

Fitzhugh-Nagumo方程式：左から $I = 1.0, 1.3, 1.8$ の場合。ピンクが前者のヌルクラインで、

　リミットサイクルが生じている場合、固定点まわりでの２つのヌルクラインの傾きがそろっており、またブラッセレータとは逆となっていることが分かります。これをヤコビ行列の符号の組で考えると

$\displaystyle J = \begin{pmatrix} + & - \\ + & - \end{pmatrix} or \begin{pmatrix} - & + \\ - & + \end{pmatrix}$

$\displaystyle J = \begin{pmatrix} + & + \\ - & - \end{pmatrix} or \begin{pmatrix} - & - \\ + & + \end{pmatrix}$

となります。Fitzhugh-Nagumoモデルは上の場合で、ブラッセレータは下の場合に相当しており、両方とも振動性を示す際の相平面の特徴として考えることができます。前者は activator-inhibitor , 後者は substrate-inhibition などと呼ばれるそうです。

　また、Fitzhugh-Nagumoモデルでも、Hodgkin-Huxleyモデルにみられるような「興奮性」が現れていると言えます。「興奮性」というのは、サイクル的な動きが得られるために閾値が存在していること、そして入力に対して過渡的な変化が生じる性質のことです。神経の活動電位も「上昇⇒下降⇒不応」と変化しますが、これは興奮系にみられる特徴だと言えます。一方、入力を上げるとリミットサイクルに落ち着いている様子は「自励振動」と呼ばれるようです。

時空間パターン

拡散方程式と濃度勾配

　ここからは空間内での物質の濃度勾配などの話をするために、拡散方程式という偏微分方程式を導いておきます（お気持ち程度に）。

　簡単のために、考える物質が $x$ 方向のみに動く（濃度を $c(x, t)$ ）として、 $x$ 方向で $[x, x+\Delta x]$ の範囲内（直方体型の領域とします）での物質量の時間的変化について定式化します。領域内の物質量の時間的変化（ $\Delta t$ ）は濃度の時間的な違いと流入の違いという2つの方向から次のように表せます。

$\displaystyle (c(x, t+\Delta t) - c(x, t))\Delta x = (J_x (x) - J_x(x +\Delta x)) \Delta t$

$J_x$ は $x$ 方向への物質の流れ（ $x$ 方向を正として）になります。これをそれぞれを微小量にして微分型で表すと

$\displaystyle \frac{\partial c(x, t)}{\partial t} = - \frac{\partial J_x (x)}{\partial x}$

となります。さらにここで、「濃度が高い場所から低い場所へと物質が流れ、それは勾配 $grad$ に比例する」というフィックの第1法則を用いることで、

$\displaystyle J_x(x) = - D \frac{\partial c(x, t)}{\partial x}$

と書けます。これにより先ほどの式は、

$\displaystyle \frac{\partial c(x, t)}{\partial t} = D \frac{\partial^2 c(x,t)}{\partial x^2}$

となります。これを3次元などに拡張した場合も

$\displaystyle \frac{\partial c(x, t)}{\partial t} = D \nabla^2 c(x,t)$

と書けます。この式が拡散方程式と呼ばれる式です。拡散方程式では時間反転に対して非対称であり、それは拡散が不可逆な現象であることに対応していると言えます。

　この拡散方程式の解（1次元）としては、ガウス型の分布が当てはまります（ $c_0 = c(0, 0)$ で原点にデルタ関数の形で局在する状態を初期条件として）。実際に拡散方程式から導くにはフーリエ解析を行えばよいです。

$\displaystyle c(x,t) = \frac {c_0}{\sqrt{4 \pi D t}} exp(- \frac {x^2}{4Dt})$

例えば、実際に起きている拡散現象は主にブラウン運動と呼ばれる形であることが多く、この解が意味するところは、最初局所的に集まっていた物質も、ブラウン運動に沿って全体的にバラついていく（全体の分散が大きくなっていく）ことに相当します。 $c(x, t) / c_0$ は初期の濃度で割っていることから確率密度として考えて良く、このような無次元化によりガウス分布としての分散（平均二乗変位）が

$\langle x^2 \rangle = 2Dt$

として求まります。確かに分散の値は $t$ に応じて上昇していくことが分かります。3次元の場合は $6Dt$ となり、より速くなると言えます。

　一つの例として、モルフォゲンの拡散について考えてみます。モルフォゲンとは、細胞内でのそれの濃度の違いがその後の細胞の運命（どのような組織になるかなど）を変えるというもので、細胞の集まり全体で拡散的な濃度勾配を形成しているという特徴があります（フレンチフラッグモデルなど）。

　 $x \in [0, L]$ の範囲で $c(0, t) = c_max$ として $x=0$ で常にモルフォゲンが産生されているとして、またモルフォゲンが常に分解を受ける状態だとすると、次のような定式化になります。

$\displaystyle \frac{\partial c(x, t)}{\partial t} = -k c(x,t) + D \frac{\partial^2 c(x,t)}{\partial x^2}$

これを解くと、モルフォゲンは定常状態で指数関数的減少を示します。実際にこのような仕組みが Bicoid といったモルフォゲンで成り立っているようです。

チューリングパターン

　空間的な濃度変化の他の例として、Alan Turing のチューリングパターンがあります。これは、拡散方程式の線形安定性解析という側面から考えることが出来ます。

　今考える拡散方程式の一般形を

$\displaystyle \frac{\partial c_i(x, t)}{\partial t} = f_i (\{c\}) + D \frac{\partial^2 c_i(x,t)}{\partial x^2}$

とします。 $c_i$ は物質 $i$ の濃度で、 $\{c\} = \{c_1, c_2, \dots , c_n\}$ としています。これの固定点 $c^*$ からの摂動 $\delta c_i$ に対してテイラー展開の一次までで考えてみます。

$c_i (x, t) = c_i^* + \delta c_i(x, t)$

$\displaystyle \therefore \frac{\partial \delta c_i}{\partial t} = \sum_{j=1}^N J_{ij}\delta c_j + D \frac{\partial \delta c_i}{\partial x^2}$

$\displaystyle J_{ij} = \left.\frac{\partial f_i}{\partial c_j} \right|_{c^*}$

これをフーリエ変換してみる（ $\delta c_i \to \delta \hat{c}_i$ ）と

$\displaystyle \frac{\partial \delta \hat{c}_i}{\partial t} = \sum_{j=1}^N (J_{ij} - D_i k^2 \delta_{ij}) \delta \hat{c}_j(k,t)$

となります。ここで、和の中にある行列部分は合わせて

$\displaystyle M_{ij} = \begin{pmatrix} J_{11} - D_1 k^2 & J_{12} & \cdots & J_{1N} \\ J_{21} & J_{22} - D_2 k^2 & \cdots & J_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ J_{N1} & J_{N2} & \cdots & J_{NN} - D_N k^2\end{pmatrix}$

と書けます。ここで得られた行列に対して第1回同様の議論を経ることで、安定固定点かどうかなどの判別が出来ます。この時、特にどの固有値の虚部もゼロの条件で、ある波数での固有値の実部が正の場合、その波数の変動が大きくなっていくことで空間周期的な解が得られ、これをチューリング不安定性と呼びます。

　チューリング不安定性の条件を、2種類の系で求めてみます。

$\displaystyle M = \begin{pmatrix} J_{11} - D_1 k^2 & J_{12} \\ J_{21} & J_{22} - D_2 k^2\end{pmatrix}$

見やすくするため

$J_{11}=a, J_{12} = b, J_{21} = c, J_{22} = d, J_{11} - D_1 k^2 = \tilde{a}, J_{22} - D_2 k^2 = \tilde{d}$

とすると、固有方程式は

$\lambda ^ 2 - (\tilde{a} + \tilde{d}) \lambda + \tilde{a}\tilde{d} - bc = 0$

となります。まずここで、 $k = 0$ , つまり空間全体で一様な成分が安定である（不安定な場合、解が無限に発散することになります）という条件を考えると、

$\tilde{a} + \tilde{d} = a + d < 0$

$\tilde{a}\tilde{d} - bc = ad - bc > 0$

となり、これを整理すると、ヤコビ行列 $J_{ij}$ の符号が

$\displaystyle J = \begin{pmatrix} + & - \\ + & - \end{pmatrix} or \begin{pmatrix} - & + \\ - & + \end{pmatrix}$

$\displaystyle J = \begin{pmatrix} + & + \\ - & - \end{pmatrix} or \begin{pmatrix} - & - \\ + & + \end{pmatrix}$

に分類されることが分かります。これは先ほどの activator-inhibitor , substrate-inhibition の例と同じになります。

　一方で、ある $k \neq 0$ において不安定な条件になると仮定すると、

$\tilde{a}\tilde{d} - bc <= 0 \Rightarrow -bc <= (a - D_1 k^2)×-(d - D_2 k^2)$

となります（第1式は $k=0$ が安定という条件下で必ず満たされるので）。この式を満たしやすい $D_1, D_2$ の条件として、

$D_2 \gg D_1$

が考えられます。これは実際の系について考えてみれば、inhibitor の拡散が activator の拡散に比べて速いということになります。このような場合、activator の先が inhibitor によって抑制されるものの、さらにその先では activator の減少による inhibitor の減少がみられることで再び activator が増えて…といった周期的な構造が見られ、これがチューリングパターンになります。

　なお、この時に増幅される波長 $k_c$ は、

$bc = (a - D_1 k_c^2)(d - D_2 k_c^2)$

を満たすもので、解くことにより

$\displaystyle k_c = \left(\frac{ad-bc}{D_1D_2}\right)^\frac{1}{4}$

と求まります。

進行波のモデル

化学進行波と界面

　BZ反応と呼ばれる有名な化学反応があります（ggってみてください）。特徴的なこととして、スパイラル（らせん波）という形状で繰り返し波が生成されるという現象があります。このようなスパイラルは、細胞分裂面の決定や、心臓における心室細動などで見られるそうで、普遍的な現象だと言えます。ここからは、このようなスパイラルについて進行波の数学モデルを用いて考えていきます。

　一次元の進行波から考えてみると、次の形のものが進行波となります。

$u(t, x) = u(x - ct)$

このとき、波動方程式を考えると $c$ が進行波の速度となります。波は $+c$ と $-c$ の二方向の波の重ね合わせとなっており、エネルギーの保存が成り立っています。一方、拡散方程式のみからなる系では進行波は現れません（補足で説明）。進行波のためには反応項を必要とすると言えます。

　反応項の一例として、Allen-Cahn方程式を取り上げます。

$\displaystyle \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + ku(1-u)(u-a)$

反応項である最後の項について、

$\displaystyle ku(1-u)(u-a) = - \frac{d}{du}[\frac{ka}{2}u^2 - \frac{k(a+1)}{3}u^3 + \frac{k}{4}u^4 ]$

が成り立ちます。このとき微分の中身の項はポテンシャルと呼ばれており、4次関数で2つの安定固定点で1つの不安定固定点を挟んだホップ分岐後のような形をしています。統計力学におけるIsingモデルも同様のポテンシャルを持ち、Isingモデルの連続版が Allen-Cahn方程式に相当するようです。以下では、4次関数における2つの"谷"のうち、 $u = 1$ の方が $u = 0$ よりも深い場合について考えます。

　さて、この方程式の解として次のようなものを考えます。 $x \to -\infty$ で $u = 0$ , $x \to \infty$ で $u = 1$ であり、ある $x$ でシグモイド関数様に $0$ から $1$ へと変位が移る解を考えます。このとき、 $u=1$ の解の方がポテンシャルが深くより安定であるため、物理的直観によると波は $u=1$ の領域が増加する左側へと移動すると考えられます。

　そこで、進行波解 $u(x-ct) = u(z)$ を仮定してこの方程式を解いてみます。まずは（次元成分で両辺を割って）無次元化し、さらにもとの方程式を $z$ の2階常微分方程式に帰着することで、

$\displaystyle \frac{d^2 u}{d z^2} + c \frac{d u}{d z} + u(1-u)(u-a) = 0$

この式の厳密解を求めるために、テクニカルですが $\displaystyle \frac{d u}{d z}$ を両辺にかけて $z$ で積分します。

$\displaystyle \int \frac{d^2 u}{d z^2}\frac{d u}{d z} dz = \frac{1}{2}\left (\frac{d u}{d z}\right )^2 + C$

$\displaystyle \int u(1-u)(u-a) \frac{d u}{d z} dz = \int u(1-u)(u-a) du = \frac{z^2}{12} \left (2a(2z-3) - z(3z-4) \right ) + C'$

ここで、2式を実際に $z = - \infty$ から $z = \infty$ まで積分すれば、それぞれ

$\displaystyle 0, \frac{1}{6}(a - \frac{1}{2})$

となります。つまり、 $a = \frac{1}{2}$ のときに

$\displaystyle \frac{1}{2}\left (\frac{d u}{d z}\right )^2 + c \int \left (\frac{d u}{d z}\right )^2 dz - \frac{z^2}{4} \left (z+1 \right )^2 = 0$

が成り立つことになります。もし $c = 0$ ならば上の方程式は単なる常微分方程式の問題に帰着され、その解は

$\displaystyle u = \frac{1}{1 + exp(\frac{-z}{\sqrt{2}})}$

となります。 $c = 0$ の時の解は $t = 0$ の時の解とも考えられ、それは $z = x$ の時に相当します。いま進行波を考えていたために、 $c \neq 0$ の場合であっても上の式の形が解になるということが考えられます。一方、 $a \neq \frac{1}{2}$ の場合は

$\displaystyle \frac{d^2 u}{d z^2} + c \frac{d u}{d z} + u(1-u)(u-\frac{1}{2}) + u(1-u)(\frac{1}{2} - a) = 0$

と分解してあげれば前半部分が今までの議論で消えることから

$c \frac{d u}{d z} = c \frac{u(1-u)}{\sqrt{2}} = u(1-u)(\frac{1}{2} - a)$

$\therefore c =\sqrt{2} (a - \frac{1}{2})$

と求まりました。重要なのは、最後の波の速度が一定であり、無次元化を元に戻すことで

$c = \sqrt{2} (a - \frac{1}{2}) \sqrt{kD}$

と $\sqrt{kD}$ に比例するようになるということです。今求まった解はいわゆるフロント（界面）をもつ波の形をしており、この方程式での界面の移動速度（伝播速度）が $\sqrt{kD}$ に比例することが分かりました。

ケーブル方程式と進行波

　神経の膜電位が伝わっていく過程は、ケーブル理論という形で定式化できます。このケーブル理論は最終的に上記の拡散の理論と似た式の形へと落とし込むことが出来ます。

　神経を次のように考えます。細胞膜において、細胞側（内側）と外側とで電流が流れ（電線が存在）、その間を膜やチャネルがコンデンサーや導電体（抵抗）として繋いでいるというモデルです。下にモデルの詳細な図を引用してきました（脳科学辞典様より）。

　さて、細胞外側を流れる電流・電圧を $I_o, V_o$ 、内側は $I_i, V_i$ 、繋いでいる部分（接点とします）の電流は $I_t$ とします。ある接点 $x$ とその隣の節点 $x+dx$ の間 $[x, x+dx]$ の電荷の変化（考える上では $x$ での内側と外側の電線の交点での電流のやり取り）は

$- I_o(x+dx) + I_o(x) - I_t(x)dx = 0$

$- I_i(x+dx) + I_i(x) + I_t(x)dx = 0$

と書けます。よって $dx \to 0$ の極限を取ると

$\displaystyle \frac{\partial I_o}{\partial x} = -I_t$

$\displaystyle \frac{\partial I_i}{\partial x} = I_t$

という式が得られます。一方で、電圧はオームの法則に従って下がることから

$\displaystyle \frac{\partial V_o}{\partial x} = -I_o r_o$

$\displaystyle \frac{\partial V_i}{\partial x} = -I_i r_i$

と表されます。いつも考える電位差 $V = V_o - V_i$ と、電流の総和 [I_{tot} = I_o + I_i] の関係は、以上の関係式から

$\displaystyle \begin{eqnarray} I_{tot} &=& -\frac{1}{r_o}\frac{\partial V_o}{\partial x} -\frac{1}{r_i}\frac{\partial V_i}{\partial x} &=& -\frac{1}{r_o}\frac{\partial V}{\partial x} + \frac{r_o + r_i}{r_o}I_i \end{eqnarray}$

と求まります。一方で、 $I_{tot}$ はどの場所でも保存されるため $\displaystyle \frac{\partial I_{tot}}{\partial x} = 0$ となり、 $\displaystyle I_t(x) = \frac{\partial I_i}{\partial x}$ は

$\displaystyle I_t(x) = \frac{\partial}{\partial x}\left ( \frac{1}{r_i + r_o} \frac{\partial V}{\partial x} \right )$

と書けます。最後に、 $\displaystyle I_t = I_{ion} + C_m \frac{\partial V}{\partial t}$ と書けること（ $I_{ion}$ はイオンチャネルを通る電流）から、次に示すケーブル方程式が導かれます。

$\displaystyle C_m \frac{\partial V}{\partial t} = - I_{ion} + \frac{\partial}{\partial x}\left ( \frac{1}{r_i + r_o} \frac{\partial V}{\partial x} \right )$

これは、 $\displaystyle D = \frac{1}{C_m(r_i+r_o)}$ である拡散方程式の形の式となります。つまり、記事の最初の方で扱った神経系の方程式（ $I_{ion}$ 部分）に拡散方程式の項を追加することで考えることが出来ると言えます。このことと直前で扱った進行波的発想とを組み合わせてみましょう。

　Fitzhugh-Nagumo方程式に今のケーブル方程式から導かれる拡散方程式的項を追加すると、

$\displaystyle \epsilon \frac{du}{dt} = u(a-u)(u-1) - v + I + D \frac{\partial^2 u}{\partial x^2}$

$\frac{dv}{dt} = bu-\gamma v$

となります。よく見ると、最初の式はAllen-Cahn方程式とほぼ同じになっています。つまり、最初の項の解析にAllen-Cahn方程式の時の結果を考慮できると言えます。例えば、波の速度に関して $\sqrt(D)$ に比例することなどが考えられます。

　最後に、スパイラルについて考えてみます。Fitzhugh-Nagumo方程式は興奮系の方程式ですが、スパイラルを起こすことがあります。このようなスパイラルは所謂セルオートマトンというモデルで簡単に考えることが出来ます。セルオートマトンは興奮系の「安定状態」「興奮状態」「不応期」の３状態をとる格子状のモデルであり、周りの格子の状態に応じて次の自分の状態が変化します。そのようなモデルでは、不応期を越すような大きな周期での界面の相互作用（回り込む動き）があればらせんになることが分かります。

　なお、スパイラルの回転を生み出す点では空間的な位相の周期積分が $2\pi$ になり、特異点とも呼ばれています。興奮系ではこのようなスパイラルを生み出す性質があり、また興奮系でなく振動系であってもスパイラルを起こすこともあるようです。

　今回で力学系からパターン形成にかかわる内容までを扱いました。次回以降は（先生が変わったこともあり）内容が生物に近い進化シミュレーションを扱います。

補足（末尾）

拡散だけで進行波が起きないことの説明

　そのうち書きます。