Aritalab:Lecture/NetworkBiology/Markov Chains/Birth-death Process

参考文献

• Lieberman E,Hauert C, Nowak MA (2005) Evolutionary Dynamics on Graphs. Nature 433:312-316

出生死亡過程

モラン過程

個体数を n とし、状態 i から i+1, i-1 の状態へそれぞれ確率$\alpha, \beta \ (\alpha + \beta < 1)$で遷移する場合を考える。 また状態 0 と n は吸収状態とする。（従ってユニークな定常分布を考えるわけではない。） このマルコフ連鎖を1958年にモデルを発表した遺伝学者PAP MoranにちなんでMoran過程という。

状態 i から出発して状態 n に到達する確率を $x_i$ と書く。

$\begin{cases} x_0 &= 0 \\ x_i &= \beta x_{i-1} + (1 - \alpha - \beta) x_i + \alpha x_{i+1} \quad (i = 1, \ldots N-1) \\ &= x_i + \alpha (x_{i+1} - x_i) - \beta (x_i - x_{i-1}) \\ x_n &= 1 \\ \end{cases}$

ここで記法

$\begin{align} y_i &= x_i - x_{i-1} \quad (i = 1, \ldots, N) \\ \gamma &= \beta/\alpha \end{align}$

を導入すると

$\textstyle y_{i+1} = \gamma y_i, \quad \sum^n_{i=1} y_i = 1$

を得る。Moran過程では通常、少ない個体数からn個体にまで増殖するケースを考えるため、$\gamma$は1より小さい値と考えてよい。

$y_1 = x_1, \ y_2 = \gamma x_1, \ y_3 = \gamma^2 x_1, \cdots, y_n = \gamma^{n-1} x_1$

を足し合わせると

$x_1 = \frac{1}{\sum^{n-1}_{j=0} \gamma^j} = \frac{1 - \gamma}{1 - \gamma^n}$

また

$\begin{align} x_i &= y_i + x_{i-1} = \gamma^{i-1} x_1 + y_{i-1} + x_{i-2} \\ &= \gamma^{i-1} x_1 + \gamma^{i-2} x_1 + x_{i-3}\\ &= \ldots\\ &= \textstyle x_1 \sum^{i-1}_{j=0} \gamma^j \end{align}$

とあわせて

$x_i = \frac{\sum^{i-1}_{j=0} \gamma^j}{\sum^{n-1}_{j=0} \gamma^j} = \frac{1 - \gamma^i}{1 - \gamma^n}$

式の解釈

n 個体の集団において、全ての個体が最初タイプAであるとする。 このときタイプBという突然変異が1個生じたとする。 タイプBが $\gamma < 1$を満たす、すなわち$\alpha$が$\beta$よりも大きくて個体数を増やす傾向にあると仮定する。 タイプBが何らかの優性種であると考えてもよい。 初めに1個生じたタイプBが集団 n の中に固定される（集団全体をカバーする）確率は$x_1 \xrightarrow{n \rightarrow \infty} 1 - \gamma$。 つまり、Aより10倍優位な変異体であっても ($\gamma = 1/10$)、集団全体にその変異が保持される確率は 9/10 でしかない。 通常の遺伝子変異で生存戦略に極めて優性となるケースは非常に稀だとすれば、殆どの変異は集団全体に保持されないことになる。

木村資生の中立進化説

木村資生（きむらもとお）は1970年代に活躍した国立遺伝学研究所の研究者である。 ほとんどの遺伝子の突然変異は、自然淘汰において有利でも不利でもないとする中立進化説は当時大きな反響を呼んだ。

突然変異が完全に中立な場合、$\gamma = 1$　から $x_1 = 1/n$となる。 各個体の子孫がそれぞれ集団全体をカバーする可能性があり、各個体が均等な機会を持つから、これは当然の値ともいえる。 さて、確率 m でタイプBという突然変異が生じるとき、変異体Bの生じる個体数は mn になる。 各個体がそれぞれ集団をカバーしうるから、全てタイプAの集団からはじまって全てタイプBの集団に進化する速度は、$mn/n = m$ となり、集団のサイズに依存しない。 ここから「分子時計」仮説が導かれる。

ネットワーク上の出生死亡過程

モラン過程ではA個体とB個体の個数のみを考え、個体間のネットワークという概念が無かった。 ここでは n 個体が遷移行列 ${\mathbf P}$ であらわされるネットワークを考えよう。 各頂点は状態AまたはBをとり、単位時間にある頂点 i が隣接する頂点 j に確率 $p_{ij}$ で自分の状態をコピーする。 また、BはAに比べ $1/r$ の割合で選ばれる確率が高いとする。ここで扱うのは$0 \leq r \leq 1$の場合である。

モラン過程と同等なネットワーク

初めに1個体であったBが全体に広がる確率がMoran過程と同じ $(1-r)/(1-r^n)$ であるための条件を考えよう。 現在 m 頂点が状態B、残りが状態Aと仮定し、頂点 i が状態Bのとき $v_i = 1$, 状態Aのとき $v_i = 0$と書こう。 B個体が m+1 に増える確率は、状態B (頂点 i とする)が選ばれて隣接する状態A (頂点 j とする) を上書きすればよいから

$\textstyle B_{m\rightarrow (m+1)} = \frac{1}{r} \frac{\sum_i\sum_j p_{ij} v_i (1-v_j)}{m/r + N - m}$

同様に、m-1 に減少する確率は

$\textstyle B_{m\rightarrow (m-1)} = \frac{\sum_i\sum_j p_{ij} (1-v_i) v_j}{m/r + N - m}$

Moran過程と同等になるには各頂点の状態に非依存で

$\textstyle \frac{B_{m\rightarrow (m-1)}}{B_{m\rightarrow (m+1)}} = r$

が成立せねばならず

$\textstyle \sum_i\sum_j p_{ij} v_i (1-v_j) = \sum_i\sum_j p_{ij} (1-v_i) v_j$

個体 k ただ一つがBである特殊例 $v_k = 1,\ v_i = 0 (\forall i \neq k)$ を考えると $\textstyle \sum_j p_{jk} = \sum_j p_{kj} = 1$ となる。 ただし右側の等号は、考えているネットワークがマルコフ連鎖であるため、頂点から出る辺確率の和が1であることに対応する。 つまり、定常分布での存在確率が全頂点において等しい場合は、Moran過程と全く同じプロセスになる。

定常分布での存在確率が全頂点で等しくなるのは接続関係が対称なネットワークだけに限らず、非対称な場合も存在する。 例えば一方向にのみ遷移する環状ネットワークは方向性を持つので非対称だがMoran過程と同等になる。

固定確率を低めるネットワーク

いかに優性な個体Bが生じても、固定確率が中立進化の場合に同じ $1/n$ となるネットワークは簡単に構成できる。 例えば既約なネットワークに根 (root) となる頂点を一つ付け足した場合を考えれば、根にあたる頂点に優性な変異Bが生じないと固定されないので固定確率が $1/n$ となる。同じ理由で、根にあたる頂点が複数あるネットワークは決して固定されることがない。

固定確率を高めるネットワーク

星型 (star) のネットワークはMoran過程より固定確率が高くなる。 現在 m 頂点が状態B、残りが状態Aと仮定する。 星型の枝が非常に多いとき、最初のB個体は非常に高い確率で周囲に位置する。 ここではB個体が増えていくプロセスを考えてみよう。

Bが増える場合 (最低3ステップ必要)

中央がBでないときに周囲のBを選んで中央をBにする (m+1)

中央のBを選んで、周囲のBをひとつ増やす (m+2)

周囲のAを選んで中央をAにする (m+1)

というステップの繰り返しになるので $\textstyle B_{m\rightarrow (m+1)} = \frac{m}{r(N-m+m/r)} \frac{1}{r(N-m+m/r)} \frac{N-m}{N-m+m/r}$

Bが減る場合

中央のAを選ぶ (m)

周囲のBをひとつ減らす (m-1)

周囲のAを選んで中央をAにする (m-1)

なので $\textstyle B_{m\rightarrow (m-1)} = \frac{1}{N-m+m/r} \frac{m}{r(N-m+m/r)} \frac{N-m}{N-m+m/r}$

したがって

$\textstyle \frac{B_{m\rightarrow (m-1)}}{B_{m\rightarrow (m+1)}} = \frac{m(N-m)}{m(N-m)/r^2} = r^2$

優性率 r がMoran過程に比べて $r^2$ と増幅されていることがわかる。

ネットワーク上のゲーム

各頂点は 協力 (C: Cooperate) か裏切り (D: Defect) のいずれかの戦略を取るとする。 協力者は隣接頂点それぞれにコスト (cost) c を支払い、隣接する頂点は利得 (benefit) b を受け取る。 裏切る場合はコストを払わず、隣からの利得のみを受け取る。 頂点の適応度 (fitness) を支払いと利得の合計で判断し、戦略を更新する方法として以下の2通りを考えよう。

1. 出生死亡過程: 単位時間毎に適応度に応じて頂点を選択し、隣接頂点をランダムに選んで同じ戦略へとコピーする。
2. 死亡出生過程: 単位時間毎に頂点がランダムに選ばれて死亡する。隣接する頂点がそれぞれの適応度に応じて自分の戦略をその頂点にコピーする。

もっとも単純なケースとして単純サイクルを考える。

出生死亡過程の場合

C戦略とD戦略をとる頂点の境界部分にある2頂点を考慮すればよい。

C戦略の端にいる頂点は、隣接する2頂点分のコストを払って利得を片側からしか得ないので合計は $b-2c$ になるが、D戦略の端にいる頂点は利得 b を得られるので、いかなる b, c であってもD戦略が優位になる。

死亡出生過程の場合

C戦略とD戦略をとる頂点の境界部分から2つ離れたところまでの合計4頂点を考慮する。

D戦略を取る側は、C戦略と隣接する頂点が b、そのまた隣は 0 になる。 C戦略を取る側は、D戦略と隣接する頂点が $b-2c$、そのまた隣は $2(b-c)$になる。 従って $b/c > 2$ の場合はC戦略が優位。 一般化すると、次数 k の格子であれば $b/c > k$の場合にC戦略が優位になる。