Aritalab:Lecture/Bioinformatics/Lotka

ロトカ・ヴォルテラ方程式

ロトカ・ヴォルテラ (Lotka-Volterra) 方程式は、被食者 (prey) ・捕食者 (predator) の関係を表す代表例です。 このモデルでは、被食者 x は放っておくと増加し (ax)、捕食者によって減少します (bxy) 。捕食者 y は捕食によって増加し (dxy)、放っておくと減少します (cy)。ここで a, b, c, d は全て正のパラメータです。

$\begin{align} \frac{dx}{dt} &= ax - bxy = bx (a/b - y)\\ \frac{dy}{dt} &= dxy - cy = dy (x - c/d) \end{align}$

これ以降、 x , y を適当なスケーリングによって無次元化しておきます。問題の本質は変わりません。

$\begin{align} \frac{dx}{dt} &= x ( b - y)\\ \frac{dy}{dt} &= y (-a + x) \end{align}$

この解は

• $x(t) = y(t) = 0$
• $x(t) = 0, y(t) = y(0) e^{-at}$
• $y(t) = 0, x(t) = x(0) e^{bt}$

これらがxy平面において x > 0, y > 0 象限の境界をなしており、(x , y) の動く範囲を中に包含しています。平衡点は原点および (a , b) です。この点を境に dx/dt, dy/dt の符号が変化します。その符号から x = a, y = b を境に4つのエリアに分割され、半時計方向に回転する軌道の存在がわかります。

リアプノフ関数

回転する軌道で保存される量を考えてみます。(a , b) に依存する恒等式として

$(\frac{a}{x} - 1) \frac{dx}{dt} + (\frac{b}{y} - 1) \frac{dy}{dt} = 0$

が成立します。つまり微分してこの値になる式は x , y に依存しません。

$\begin{align}\frac{dH(x,y)}{dt} &= \frac{d}{dt}[ a\log x - x + b \log y - y] = 0 \\ H(x,y) &= a\log x - x + b\log y - y = const.\end{align}$

Hを x で微分すると (a/x - 1), 二階微分は - a / x² なので x は a で最大値をとります。同様に y は b で最大値をとります。つまり H は (a , b) を中心に - ∞ まで減少し、その等高線は閉曲線を成します。つまり軌道は周期的です。

このように、平衡点の安定性を示すのに使う関数をリアプノフ関数と呼びます。

周期の振幅と振動数は a , b に依存します。平均値を調べてみましょう。ちょっと技巧的ですが log x(t) の時間微分を積分します。

$\begin{align} \frac{1}{T}\int^T_0 \frac{d}{dt} \log x(t) dt &= \frac{1}{T}\int^T_0 \frac{dx}{dt} \frac{d}{dx} \log x(t) dt = \frac{1}{T}\int^T_0 x (b - y) \frac{1}{x} dt \\ \log x(T) - \log x(0) &= 0 = b - \frac{1}{T}\int^T_0 y(t) dt \end{align}$

つまり y の時間平均は b です。同様に x の時間平均は a です。

種内競争を持つ被・捕食者系

ロトカ・ヴォルテラ方程式は捕食者がいないと被食者が指数的に増えてしまうので、これをロジスティック式で補正します。

$\begin{align} \frac{dx}{dt} &= x (b - y - ex)\\ \frac{dy}{dt} &= y (-a + x - fy) \end{align}$

ここでは −ex , −fy という項を追加しました。 このとき $dx/dt = 0 , dy/dt = 0$ となるのはそれぞれ $y = -ex + b, fy = x - a$ という直線になります。

二本の直線はy切片がそれぞれ正と負です。第一象限で交わる場合もあるし、そうでない場合もあるので、まず交わる場合を考えます。

直線上の、交点の左右における時間微分値

直線	dx/dt	dy/dt
$y = -ex + b$	0	負 → 交点 → 正
$fy = x - a$	正 → 交点 → 負	0

正負の値から、反時計回りに回転することがわかります。ただしこの解析だけからは、交点に収束するのか交点から発散するかはわかりません。

第一象限で交わらない場合は、第一象限に交点があるときの挙動を切り取った形になっています。x 軸に相当する部分 (y = 0) を考えると $y = -ex + b$ の x 切片 b/e に収束します。 y 軸に相当する部分 (x = 0) を考えると、常に減少です。すなわち、二本の直線が y 座標が負の部分で交わるとき、 (0, b/e) が安定不動点になります。

リアプノフ関数

交点の周りでの挙動を調べるため、ロトカ・ヴォルテラ方程式と同様のリアプノフ関数を考えてみましょう。

$H(x,y) = u \log x - x + v \log y - y$

この関数は二本の直線の交点（x座標は正）

$(u , v) = \Big(\frac{a+bf}{ef+1}, \frac{b-a}{ef+1}\Big)$

で最大値をとり、そこから離れるほど−∞に近づきます。ここまではロトカ・ヴォルテラ方程式とそう変わりません。

交点では$v = -eu + b, fv = u - a$が成立するので

$\begin{align} \frac{H(x,y)}{dt} &= \frac{\partial H}{\partial x} \frac{dx}{dt} + \frac{\partial H}{\partial y} \frac{dy}{dt} \\ &= (\frac{u}{x} - 1) \frac{dx}{dt} + (\frac{v}{y} - 1) \frac{dy}{dt} \\ &= (\frac{u}{x} - 1) x (b - y - ex) + (\frac{v}{y} - 1)y (-a + x - fy)\\ &= (u-x) [(v + eu) - y - ex] + (v - y) [(-u + fv) + x - fy]\\ &= e(u-x)^2 + f(v-y)^2 \geq 0\\ \end{align}$

つまり $d H(x,y)/dt$ は (u , v) においてのみ 0 となり（安定）、それ以外の場合は安定点に向かって上り勾配です。H(x , y) の値と考え合わせると (u , v) は大域安定でもあります。

定性的解釈としては、如何に種内競争の項 e や f が小さくても、存在しさえすれば、それらが大域的安定点を作り出すことを示しています。つまり、ロトカ・ヴォルテラ方程式の周期性は構造的に安定ではありません。

安定性解析

ここで話を一般化するため、一般的な線形微分方程式系を考えましょう。

$\begin{align} \frac{dx}{dt} &= f(x, y)\\ \frac{dy}{dt} &= g(x, y) \end{align}$

定常状態を (x^*, y^*) とし、そこから少しずれた位置の挙動を調べます。それには (x^* + u, y^* + v) とした差分 u , v をテイラー展開します。ずれが微小であることから2次以上の部分は無視します。

$\begin{align} \frac{du}{dt} &= f(x^* + u, y^* + v) - f(x^*, y^*) = f_x(x^*, y^*)u + f_y(x^*, y^*)v + \cdots \\ \frac{dv}{dt} &= g(x^* + u, y^* + v) - g(x^*, y^*) = g_x(x^*, y^*)u + g_y(x^*, y^*)v + \cdots \end{align}$

これを行列の形に書くと

$\frac{d}{dt}\begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} f_x & f_y \\ g_x & g_y \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix} = J \begin{bmatrix} u \\ v \end{bmatrix}$

行列 J をヤコビ行列 (Jacobian) といいます。線形微分方程式であれば J の固有値 λ₁, λ₂ を用いた解があります。

λ₁ ≠ λ₂ のとき

$\begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} c_{11} \\ c_{21} \end{bmatrix} \exp(\lambda_1 t) + \begin{bmatrix} c_{12} \\ c_{22} \end{bmatrix} \exp(\lambda_2 t)$

λ₁ ＝ λ₂ ＝ λ のとき

$\begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} c_{11} \\ c_{21} \end{bmatrix} \exp(\lambda t) + \begin{bmatrix} c_{12} \\ c_{22} \end{bmatrix} t \exp(\lambda t)$

固有値を実数とすると、λ₁ < 0, λ₂ < 0 であれば、平衡点からのズレ u , v が時間とともに e^{λ t} → 0 に収束します。どちらかでも正であれば、ズレは発散します。 固有値が複素数の時、２つの固有値は複素共役になり λ₁, λ₂ = a ± i b です。虚数の exp は絶対値が 1 になることを利用して

$\lim_{t \rightarrow \infty}| \exp(a \pm ib)t | = \lim_{t \rightarrow \infty} \exp(at) \times | \exp(\pm ibt) | = \lim_{t \rightarrow \infty} \exp(at)$

です。つまり虚数の部分は影響を与えません。固有値の実部が共に負であれば、平衡点からのズレはやはり 0 に収束することがわかります。 さらにこの議論は多元連立の場合もそのまま成立します。

まとめると、

• 固有値の実部が全て負であれば、平衡点は局所安定
• 固有値の実部が一つでも正であれば、平衡点は不安定

となります。固有値が複素数の時は近傍で振動がおこります。

2 x 2 行列 $A = \begin{bmatrix} a & b\\ c & d\end{bmatrix}$ の場合は固有方程式 x² - (a + d) x + (ad - bc) = 0 の解が固有値です。 固有値実部が負であるための必要十分条件は tr(A) = a + d < 0 かつ det(A) = ad - bc > 0 です。

競争方程式

安定性解析をベースに、被・捕食者系を捉え直してみます。 ここで係数の符号が異なる競争系を考えましょう。

$\begin{align} \frac{dx}{dt} &= x (1 - ax - by)\\ \frac{dy}{dt} &= y (1 - cx - dy) \end{align}$

種 x と y の相互作用の関係は対称形で、被・捕食者系とは y の単独の増殖率および c の符号が違っています。安定点は

$(u , v) = \Big(\frac{d-b}{ad-bc}, \frac{a-c}{ad-bc}\Big)$

ヤコビ行列は

$J = \begin{bmatrix} -a & -b \\ -c & -d \end{bmatrix}$

です。被・捕食者系と異なり、二本の直線はいずれも傾きが負です。直線どうしは第一象限で交わる場合もあるし、交わらない場合もあります。

交わる場合から考えます。ad − bc > 0 の場合、交点 (u , v) の分母は正になり、交点が第一象限にあることから d > b > 0, a > c > 0 です。さらに ad - bc > 0, (-a + -d) < 0 を満たすので、ヤコビ行列の固有値は負になります。つまり交点は大域安定な不動点です。 同じことは dx/dt, dy/dt の符号を考えても (u , v) に近づくのでわかります。 ここまでは被・捕食者系と全く同じです。

次に ad − bc < 0 の場合を考えます（被・捕食者系は c の符号が負のため、この場合が生じませんでした）。交点の分母が負になるので b > d > 0, c > a > 0 です。このとき、x 軸上 (y = 0) の切片 1/a および y 軸上 (x = 0) の切片 1/d に収束します。初期条件に依存してどちらかの種が絶滅することがわかります。交点は不安定不動点です。

リアプノフ関数

競争方程式のリアプノフ関数を考えましょう。

$H(x,y) = u \log x - x + v \log y - y$

二本の直線の交点は

$(u , v) = \Big(\frac{d-b}{ad-bc}, \frac{a-c}{ad-bc}\Big)$

です。 交点では$1 = au + bv, 1 = cu + dv$が成立するので

$\begin{align} \frac{H(x,y)}{dt} &= \frac{\partial H}{\partial x} \frac{dx}{dt} + \frac{\partial H}{\partial y} \frac{dy}{dt} \\ &= (\frac{u}{x} - 1) \frac{dx}{dt} + (\frac{v}{y} - 1) \frac{dy}{dt} \\ &= (\frac{u}{x} - 1) x (1-ax-by) + (\frac{v}{y} - 1)y (1-cx-dy)\\ &= (u-x) [a(u-x) +b(v-y)] + (v - y) [c(u-x)+d(v-y)]\\ &= a(u-x)^2 + (b+c)(u-x)(v-y) + d(v-y)^2\\ \end{align}$

二次形式の形をとり (u , v) では 0 になっています。