Aritalab:Lecture/Bioinformatics/Logistic

１変数基本モデルの歴史

数理生物学は、マクロからミクロまでの様々な生命現象を数学的にモデルし解釈する分野です。手法の多くは微分方程式や離散系などの力学系 (dynamical system) を用い、時間や空間に沿った状態変化を解析します。歴史的にみると、数理生物学は人口の数理モデルに端を発しています。

マルサスモデル

人口は集団サイズが大きいので連続量として表現でき、古くから解析の対象でした。資源に制約がない居住環境では、時刻 t における人口を P(t) としてその成長を以下のように書けるでしょう。

$\frac{d P(t)}{dt} = r P(t)$

ここで λ は人口の成長率で出生率から死亡率を引いたものにあたります。これを解けば

$P(t) = P(0) e^{r t}$

となり指数関数的に増加します。これをマルサス型の人口モデル (Malthusian population model) といいます。

この人口モデルを離散時間に対して作れば

$P(t+1) = r P(t)$

となり、これを解くと

$P(t) = P(0) r^t$

になります。連続モデルと離散モデルはそう違わないといえるでしょう。これから P(t) の代わりに x を用います。明示的には書きませんが、時間依存の関数であることに注意してください。

ロジスティックモデル

19世紀に入ると、フェアハルスト (Pierre F Verhulst) は上限のある人口増加を定式化しました。

$\frac{dx}{dt} = x r \Big( 1 - \frac{x}{K} \Big)$

マルサス型の増加 xr と密度効果 (1 - x/K) をかけた形になっています。この式では成長率が定数ではなく一定値 K (>0) に単調に近づきます。この K を環境容量 (carrying capacity) と呼びます。つまり資源の上限に近づくと増加率に負の効果を与えます（ロジスティック効果）。この理由で、ロジスティックモデルと呼ばれます。 これを解くと

$x= \frac{K}{1 + (\frac{K}{x_0}-1) e^{-rt}}$

ここで x₀ は初期値（定数）です。上限 K に漸近することがよくわかります。この形状はシグモイドとも呼ばれます。実際の人口増加はモデルできないことが知られていますが、微生物の細胞数などにはよく当てはまります。

ロジスティック式において x = 0 と x = K は個体数が変化しない平衡状態 (equilibrium) にあたります。このうち x = K はそこからずれても元に戻る安定点 (stable node)、x = 0 は不安定点 (unstable node) です。

局所安定と大域安定

平衡から少しずらしても元に戻る平衡点を局所安定点と呼びます。局所安定な平衡点は一つとは限らないので、収束先は初期値に依存します。もし、初期条件に関わらず収束する安定点がある場合、それを大域安定な平衡点と呼びます。ロジスティック式の場合、x = K は大域安定点になっています。

変数分離法

非線形の微分方程式は常に解けるとは限らないのですが、ロジスティック式は変数分離法を使って解けます。 まず式を変形します。

$r dt = \Big( \frac{1}{x} + \frac{1}{K - x} \Big) dx$

それから積分します。

$\begin{align} \int \Big( \frac{1}{x} + \frac{1}{K - x} \Big) dx &= \int r dt \\ \log |x| - \log |K - x| &= rt + C \\ \log \big| \frac{x}{K-x} \big| &= rt + C \\ \big| \frac{x}{K-x} \big| &= \exp(rt + C) \\ x &= \frac{K e^C}{e^C + e^{-rt}} \end{align}$

ここで初期条件 $x(t=0) = x_0$ を使うと上記の式が得られます。

リッカーモデル

このモデルの離散時間版はカオスを含むとても複雑な挙動を取ります。

$x_{t+1} = x_t r \Big( 1 - \frac{x_t}{K} \Big)$

離散時間のモデルで密度の効果が指数的に効くものをリッカー (Ricker) モデルといいます。

$x_{t+1} = x_t \exp \Big[ r \Big( 1 - \frac{x_t}{K} \Big) \Big]$

この一次近似をとると $\exp \Big[ r \Big( 1 - \frac{x_t}{K} \Big) \Big] = 1 + r \Big(1 - \frac{x_t}{K} \Big)$ になるのでロジスティック式になっています。

安定人口モデル

20世紀になると、年齢別の出生率を考慮したモデルがロトカ (Alfred Lotka) などにより研究されました。ロトカはアメリカの統計学者で20世紀初頭に活躍した人です。

$B(t) = G(t) + \int^t_0 \beta(a)\gamma(a) B(t-a) da$

上式で B(t) は単位時間あたりの出生数、β は年齢別の出生率、γ は a 歳までの生存率、G(t) は時刻 0 における初期人口による出生数になります。これを安定人口モデル (stable population model) と呼びます。

レスリー行列

安定人口モデルの離散版は行列を用いて表現できます。年齢を A グループに分け、A x A 行列で単位時間あたりの人口遷移を表現します。

$L := \begin{bmatrix} b_{1} & b_{2} & \cdots & \cdots & b_{A} \\ p_{1} & 0 & \cdots & \cdots & 0 \\ 0 & p_{2} & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 &p_{A-1}& 0 \end{bmatrix}$

この行列を考案者 Patrick Leslie の名前をとってレスリー行列と呼びます。 最初の行は各年齢グループによる単位時間あたりの出生数で、次行から下の p_i は年齢 i の個体が単位時間生存して i+1 になる確率です。時刻 t における年齢 i の人口を p_i と書いた人口ベクトルを

$\mathbf{p}(t) := (p_1(t), \cdots p_A(t))^T$

とすると単位時間後の人口は $\mathbf{p}(t+1) = L \mathbf{p}(t)$ と表されます。時間 t 後の人口は $L^{t} \mathbf{p}(0)$ となります。 このモデルでは、全個体数に対する各年齢の割合が時間が経つにつれて一定に落ち着きます。また個体数の増加率も一定値 λ に落ち着き、この値はレスリー行列の最大固有値です。

時間遅れと周期解

ロジスティック式は K が大域安定点でしたが、密度効果の部分に時間遅れが伴うと周期解を取る場合があります。これは直感的にわかります。ある時刻 t に環境容量 K に達してもその効果が出てくるのは時刻 t + τ (τは時間遅れ）となる場合、フィードバックが τ だけ遅れてかかります。そのため K を挟んで量が振動する周期になりえます。

$\frac{dx_t}{dt} = x_t r \Big( 1 - \frac{x_{t-\tau}}{K} \Big)$

こうした時間遅れを効果的に生み出す仕組みが、離散的なロジスティック式です。

$x_{t+1} = x_t r \Big( 1 - \frac{x_t}{K})$

両辺を K で割って x_t / K を新たに x_t と書きましょう。このように定性的な影響を変えずにパラメータを減らすことを無次元化 (nondimensionalization) といいます。

$x_{t+1} = r x_t (1 - x_t)$

この差分方程式の平衡点は $x^* = rx^*(1-x^*)$ を解いて $x^* = 0, \frac{(r-1)}{r}$ となります。つまり図でいうと y = rx(1-x) と y = x の交点に相当します。

安定性解析

平衡点からのずれを h とおいて解析します。まず差分方程式を f(x) と書き、一般論から始めましょう。平衡点の近くでテイラー展開して 2 次以上の項を無視すると

$x^* + h_{t+1} = f(x^* + h_t) \simeq f(x^*) + f'(x^*)h_t + \cdots$

すなわち

$h_{t+1} = f'(x^*)h_t$

$h_t \propto \big( f'(x^*) \big)^t$

という、ずれの時間発展式が得られます。この式から

• 0 < f'(x^*) < 1 のとき、ずれはそのまま減少して平衡点に収束
• -1 < f'(x^*) < 0 のとき、ずれは符号を変えながら減少して平衡点に収束
• f'(x^*) > 1 のとき、ずれは発散
• f'(x^*) < -1 のとき、ずれは符号を変えながら発散

することがわかります。

ロジスティック式の安定性

ロジスティック式に話を戻します。

$x_{t+1} = f( x_t ) = r x_t (1 - x_t)$

$f'(x_t) = r(1 - 2x_t)$

この式で x が正値であることを保つには x の初期値として 0 と 1 の間に限ります。また r も 0 と 4 の間に制限する必要があります。さらに r < 1 では必ず $x_{t+1} < x_{t}$ となるので局所安定です。そのため 4 > r > 1 の場合を考えましょう。２つの平衡点 x^* = 0, (r - 1)/r に分けて考えます。

x^* = 0 のとき

• 1 < r のときに必ず不安定

x^* = (r - 1)/r のとき

f'( (r - 1)/r ) = 2 - r になります。 上の一般論から

• 0 < 2-r < 1 のとき (つまり 1 < r < 2)、短調収束（局所安定）
• -1 < 2-r < 0 のとき (つまり 2 < r < 3)、減衰しながら振動収束（局所安定）
• 2 - r < -1 のとき (つまり 3 < r < 4)、振動して収束しない（不安定）

となります。さらに 3 と 4 の間の解析をするため周期 2 の場合を考えましょう。

$x_{t+2} = f(f(x_t)) = F(x_t)$

ロジスティック式では F が 4 次式になります。これを調べていくと

• 3 < r < 3.45 のとき周期 2 の振動解
• 3.45 < r < 3.57 のとき周期が 4, 8, 16 と周期が増えてゆく振動解
• r = 3.6786 で初めて奇数周期が登場
• r = 3.82.. 以上であらゆる周期が登場、カオス的振る舞い

となることが知られています。