Aritalab:Lecture/NetworkBiology/BirthDeath/Continuous

出生死亡過程の一般形

確率過程における状態変化を ΔX(t) で表します。つまり

ΔX(t) = X (t + Δt) - X (t)

一般的な記法では X (t) = i のとき、出生率と死亡率はそれぞれ

λ_i (出生), μ_i (死亡)

で表します。非常に短い時間 Δt の間には、人数 i は一人生まれて i + 1 になるか、一人死んで i − 1 になるはずです。これを数式で書くと増減分を j で表して

$\begin{align} p_{i,i+j}(\Delta t) &= \mbox{Prob} \{ \Delta X(t) = j | X(t) = i \}\\ &= \begin{cases} \lambda_i \Delta t + o(\Delta t) & j = 1 \\ \mu_i \Delta t + o(\Delta t) & j = -1 \\ 1 - (\lambda_i + \mu_i) \Delta t + o(\Delta t) & j = 0 \\ o(\Delta t) & j \neq -1, 0, 1 \end{cases} \end{align}$

と書けます。この式を用いると、時刻 t から t + Δt の間に人口が i から j まで変化する確率は

$\begin{align} Eq.1: p_{i,j}(t + \Delta t) &= p_{i,j-1}(t) \lambda_{j-1} \Delta t \\ &+ p_{i,j+1}(t) \mu_{j+1} \Delta t \\ &+ p_{i,j}(t) [1 - (\lambda_j + \mu_j)\Delta t] + o(\Delta t) \end{align}$

と展開できます。

j = 0 のとき人口はマイナスにならないため j − 1 という状態は考えません。また μ₀ = 0 ですから

$\textstyle Eq.2: p_{i,0}(t + \Delta t) = p_{i,1}(t) \mu_{1} \Delta t + p_{i,0}(t) [1 - \lambda_0 \Delta t] + o(\Delta t)$

と簡略化されます。

j = N という人口の最大値が存在する場合、j + 1 という状態は考えません。また λ_N = 0 ですから

$\textstyle Eq.3: p_{i,N}(t + \Delta t) = p_{i,N-1}(t) \lambda_{N-1} \Delta t + p_{i,N}(t) [1 - \mu_N \Delta t] + o(\Delta t)$

となります。

コルモゴロフの微分方程式

Eq.1 - 3 から $\textstyle p_{i,j}(t)$ を引いて Δ t で割りましょう。

$\begin{align} \frac{d p_{i,j}(t)}{d t} &= \lambda_{j-1} p_{i,j-1}(t) + \mu_{j+1} p_{i,j+1}(t) - (\lambda_j + \mu_j) p_{i,j}(t) \\ \frac{d p_{i,0}(t)}{d t} &= - \lambda_0 p_{i,0}(t) + \mu_{1} p_{i,1}(t) \\ \frac{d p_{i,N}(t)}{d t} &= \lambda_{N-1} p_{i,N-1}(t) - \mu_{N} p_{i,N}(t) \end{align}$

この関係は $d{\mathbf p}/dt = {\mathbf p}Q$ という行列表記に書けます。p は各状態における存在確率です。

$\textstyle p_i (t) = \mbox{Prob}\{ X(t) = i \}$

$Q = \begin{bmatrix} - \lambda_0 & \lambda_0 & 0 & 0 & \cdots & 0 \\ \mu_1 & - \lambda_1 - \mu_1 & \lambda_1 & 0 & \cdots & 0 \\ 0 & \mu_2 & - \lambda_2 - \mu_2 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ldots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & \cdots & -\mu_N \end{bmatrix}$

この過程は、状態 0, 1, 2, ... N の間を動き、λ₀ > 0 、μ_N > 0 なら状態 0, N が反射壁のランダムウォークとみなせます。人口が増える確率は $\textstyle t_{i,i+1} = \lambda_i / (\lambda_i + \mu_i)$ 、減る確率は $\textstyle t_{i,i-1} = \mu_i / (\lambda_i + \mu_i)$ です。遷移行列 T は以下になります。

$T = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 & 0 \\ \mu_1/(\lambda_1 + \mu_1) & 0 & \lambda_1/(\lambda_1 + \mu_1) & \cdots & 0 & 0 \\ 0 & \mu_2/(\lambda_2 + \mu_2) & 0 & \cdots & 0 & 0 \\ 0 & 0 & \mu_3/(\lambda_3 + \mu_3) & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 1 & 0 \\ \end{bmatrix}$

定常分布

生成行列がQで表される連続時間マルコフ連鎖の定常分布 π は

$\pi Q = 0,\ \sum^{\infty}_{i=0} \pi_i = 1 \ \ (\pi_i \geq 0)$

で表されます。

定常分布の一意性

Eq.1 で遷移が定義される連続時間出生死亡過程は、状態数が無限のとき、以下の条件で定常状態 π が一つに定まります。

$\mu_i > 0 \mbox{ and } \lambda_{i-1} > 0 \mbox{ and } \sum^{\infty}_{i=1} \frac{\lambda_0\lambda_1 \cdots \lambda_{i-1}}{\mu_1\mu_2\cdots \mu_i} < \infty$

このとき定常分布は

$\begin{align} \pi_i &= \frac{\lambda_0\lambda_1\ \cdots \lambda_{i-1}}{\mu_1\mu_2\cdots \mu_i} \pi_0 \ (i = 1, 2, \cdots) \\ \pi_0 &= \frac{1}{1 + \sum^{\infty}_{i=1}\frac{\lambda_0\lambda_1 \cdots \lambda_{i-1}}{\mu_1\mu_2\cdots \mu_i} } \end{align}$

証明

定常分布が満たす式 π Q = 0 を直接解けばよい。

$0 = - \lambda_0\pi_0 + \mu_1 \pi_1$ から $\pi_1 = (\lambda_0 / \mu_1) \pi_0$

が導かれる。更に、

$\begin{align} \mu_2\pi_2 &= (\lambda_1 + \mu_1) \pi_1 - \lambda_0\pi_0 \\ &= \Big[ \frac{(\lambda_1 + \mu_1)\lambda_0}{\mu_1} - \lambda_0 \Big]\pi_0 \\ \pi_2 &= \frac{\lambda_0\lambda_1}{\mu_1\mu_2} \pi_0 \end{align}$

これを繰り返せば数学的帰納法によって

$\begin{align} \pi_{i+1} &= \frac{\lambda_0\lambda_1\cdots\lambda_i}{\mu_1\mu_2\cdots\mu_{i+1}} \pi_0 \end{align}$

状態数が無限大の場合は

$\begin{align} \sum^{\infty}_{i=0} \pi_i = 1 &\Leftrightarrow \pi_0(1 + \sum^{\infty}_{i=1} \pi_i/\pi_0 ) = 1 \\ &\Leftrightarrow \pi_0 = \frac{1}{1 + \sum^{\infty}_{i=1}\frac{\lambda_0\lambda_1 \cdots \lambda_{i-1}}{\mu_1\mu_2\cdots \mu_i} } \end{align}$

ただし

$\textstyle \sum^{\infty}_{i=1}\frac{\lambda_0\lambda_1 \cdots \lambda_{i-1}}{\mu_1\mu_2\cdots \mu_i} < \infty$

の成立が必要で、これは出生率が死亡率よりも少ないことを意味します。

λ と μ が定数の場合

全ての状態において λ, μ が固定値の場合、定常分布は

$\pi_i = \Big( 1 - \frac{\lambda}{\mu} \Big)\Big( \frac{\lambda}{\mu} \Big)^i \ \ (i \geq 0)$

となります（幾何分布）。人口が発散しないためには λ < μ であることが重要です。定常分布の平均値が平均人口にあたります。

$C = \frac{\lambda/\mu}{1 - \lambda/\mu} = \frac{\lambda}{\mu - \lambda}$

人口の増加率 (λ) で平均人口を割れば平均寿命になります。

$W = C / \lambda = 1 / (\mu - \lambda)$

様々な出生死亡過程

死亡率 μ が 0 の場合を純出生過程と呼びます。 出生率 λ が 0 の場合を純死亡過程と呼びます。 様々な場合に対して考えてみます。

ポアソン過程 λ_i =λ, μ = 0

人数に関係なく出生数 λ が一定の場合を考えます。（それぞれの個体はだんだん子供を産まなくなる。）

$\begin{align} p_{i,i+j}(\Delta t) &= \begin{cases} \lambda \Delta t + o(\Delta t) & j = 1\\ 1 - \lambda\Delta t + o(\Delta t) & j = 0 \\ o(\Delta t) & j \geq 2\\ 0 & j < 0 \end{cases} \end{align}$

コルモゴロフの微分方程式は

$\begin{align} \frac{d p_0(t)}{dt} &= - \lambda p_0(t) \\ \frac{d p_i(t)}{dt} &= - \lambda p_i(t) + \lambda p_{i-1}(t) \ \ (i\leq 1) \end{align}$

X(0) = 0 とします。時刻 0 で状態 0 の確率が 1 なので$p_0(0) = 1$。これを解いて $\textstyle p_0(t) = e^{-\lambda t}$

ここから順次計算して

$\begin{align} p_0(t) &= e^{-\lambda t}\\ p_1(t) &= \lambda t e^{-\lambda t}\\ p_2(t) &= (\lambda t)^2 \frac{e^{-\lambda t}}{2!}\\ \vdots &= \vdots \\ p_i(t) &= (\lambda t)^i \frac{e^{-\lambda t}}{i!} \end{align}$

すなわちポアソン分布になります。ポアソン分布の平均と分散は以下になります。

$\textstyle m(t) = \sigma^2(t) = \lambda t$

Yule過程: λ_i = iλ, μ = 0

出生数が人口に比例する場合を考えます。 最初は 1 からスタートし X(0) = 1 とします。

$\begin{align} p_{i,i+j}(\Delta t) &= \begin{cases} \lambda i \Delta t + o(\Delta t) & j = 1\\ 1 - \lambda i \Delta t + o(\Delta t) & j = 0 \\ o(\Delta t) & j \geq 2\\ 0 & j < 0 \end{cases} \end{align}$

コルモゴロフの微分方程式は

$\begin{align} \frac{d p_0(t)}{dt} &= 0 \\ \frac{d p_i(t)}{dt} &= - \lambda i p_i(t) + \lambda (i-1) p_{i-1}(t) \ \ (1\leq i) \end{align}$

ここから順次計算して

$\begin{align} p_0(t) &= e^{\lambda t}/(1 - e^{-\lambda t})\\ p_1(t) &= e^{\lambda t} \\ p_2(t) &= e^{\lambda t} (1 - e^{-\lambda t}) \\ \vdots &= \vdots \\ p_i(t) &= e^{\lambda t} (1 - e^{-\lambda t})^{i-1} \end{align}$

すなわち幾何分布になります。平均と分散は

$\begin{align} m(t) &= e^{\lambda t}\\ \sigma^2(t) &= e^{2\lambda t}(1 - e^{\lambda t}) \end{align}$

平均と分散は指数的に増大します。幾何分布になるのは初期値 X(0) = 1 としているためで、これを N に一般化すると負の二項分布になります。

純死亡過程: λ = 0, μ_i = i μ

死亡数が人口に比例する場合を考えます。 最初は N からスタートし X(0) = N とします。

$\begin{align} p_{i,i+j}(\Delta t) &= \begin{cases} \mu i \Delta t + o(\Delta t) & j = -1\\ 1 - \mu i \Delta t + o(\Delta t) & j = 0 \\ o(\Delta t) & j \leq -2\\ 0 & j > 0 \end{cases} \end{align}$

コルモゴロフの微分方程式は

$\begin{align} \frac{d p_i(t)}{dt} &= \mu (i+1) p_{i+1}(t) - \mu i p_i(t)\\ \frac{d p_N(t)}{dt} &= - \mu N p_N(t)\\ \end{align}$

これを解くと二項分布になります。

$\begin{align} p_i(t) = \binom{N}{i}e^{-\mu t i}(1 - e^{- \mu t})^{N - i} \end{align}$

平均と分散は以下になります。

$\begin{align} m(t) &= N e^{-\mu t} \\ \sigma^2(t) &= N e^{-\mu t}(1 - e^{- \mu t}) \end{align}$

平均と分散はともに指数的に減少します。

出生死亡過程: λ_i = i λ, μ_i = i μ

$\begin{align} p_{i,i+j}(\Delta t) &= \begin{cases} \mu i \Delta t + o(\Delta t) & j = -1\\ \lambda i \Delta t + o(\Delta t) & j = 1\\ 1 - (\lambda + \mu) i \Delta t + o(\Delta t) & j = 0 \\ o(\Delta t) & j \neq -1, 0, 1\\ \end{cases} \end{align}$

コルモゴロフの微分方程式は

$\begin{align} \frac{d p_i(t)}{dt} &= \lambda(i-1) p_{i-1}(t) + \mu (i+1) p_{i+1}(t) - (\mu + \lambda) i p_i(t)\\ \frac{d p_0(t)}{dt} &= \mu p_1(t) \end{align}$

計算は大変そうですがこれを解くと

$p_0(t) = \begin{cases} \Big(\frac{\mu - \mu e^{(\mu - \lambda)t} }{\lambda - \mu e^{(\mu - \lambda)t}}\Big)^{N} & \mbox{ if } \lambda \neq \mu \\ \Big( \frac{\lambda t}{1 + \lambda t} \Big)^N & \mbox{ if } \lambda = \mu \\ \end{cases}$

無限時間過ぎた後を考えると

$p_0(\infty) = \begin{cases} 1 & \mbox{ if } \lambda \leq \mu \\ \Big( \frac{\mu}{\lambda} \Big)^N & \mbox{ if } \lambda > \mu \\ \end{cases}$

平均と分散は以下になります。

$\lambda \neq \mu$ のとき

$\begin{align} m(t) &= N e^{(\lambda - \mu) t} \\ \sigma^2(t) &= N \frac{\lambda + \mu}{\lambda - \mu} e^{(\lambda -\mu) t}(e^{(\lambda - \mu) t} - 1) \end{align}$

$\lambda = \mu$ のとき

$\begin{align} m(t) &= N \\ \sigma^2(t) &= 2 N \lambda t \end{align}$

移入のある過程: λ_i = i λ + ν, μ_i = i μ

$\begin{align} p_{i,i+j}(\Delta t) &= \begin{cases} \mu i \Delta t + o(\Delta t) & j = -1\\ (\nu + \lambda i) \Delta t + o(\Delta t) & j = 1\\ 1 - [\nu + (\lambda + \mu) i ]\Delta t + o(\Delta t) & j = 0 \\ o(\Delta t) & j \neq -1, 0, 1\\ \end{cases} \end{align}$

コルモゴロフの微分方程式は

$\begin{align} \frac{d p_i(t)}{dt} &= [\lambda(i-1) + \nu] p_{i-1}(t) + \mu (i+1) p_{i+1}(t) - (\mu i + \lambda i + \nu) p_i(t)\\ \frac{d p_0(t)}{dt} &= - \nu p_0(t) + \mu p_1(t) \end{align}$

これも頑張ると解けて平均が以下になります。

$m(t) = \begin{cases} \frac{e^{(\lambda - \mu)t} (N \mu - N \lambda - \nu) + \nu}{\mu - \lambda} & \mbox{ if } \lambda \neq \mu \\ \nu t + N & \mbox{ if } \lambda = \mu \\ \end{cases}$

つまり $\lambda > \mu$ であれば指数的に増加し、等しい時は線形に増加します。しかし $\lambda < \mu$ の場合は最終的に $\textstyle \frac{\nu}{\mu - \lambda}$ に落ち込みます。