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

Wiki Top

Up one level

レポートの書き方

Arita Laboratory

参考文献

• Lieberman E, Hauert C, Nowak MA (2005) Evolutionary Dynamics on Graphs. Nature 433:312-316
• Nowak MA (2006) Evolutionary Dynamics: Exploring the Equations of Life, Belknap Press of Harvard University Press
• Broom M, Rychtar J (2008) An analysis of the fixation probability of a mutant on special classes of non-directed graphs Proc R Soc A 464:2609-2627

出生死亡過程

モラン過程

個体数を 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と仮定する。 BのほうがAより$1/r$倍選択されやすい場合、N 個の頂点のうち、状態Aのものが選択される確率は$1 /(N-m + m/r)$となる。 状態Bの頂点が選択される確率は$1/r(N-m+ m/r)$である。

頂点 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過程より固定確率が高くなることを説明する。 結果から先に言うと、固定確率は $\textstyle x_i = \frac{1-r^{2i}}{1-r^{2n}}$ であらわされる。

星型の枝が非常に多い場合を考え、周囲の頂点には状態Bが m 個、状態Aが N-m 個あると仮定しよう。 周囲の頂点で状態Bを増やすには、中央の頂点が状態Bの際に選択されねばならず、 周囲の頂点で状態Aを増やすには、中央の頂点が状態Aの際に選択されねばならない。 中央の頂点は状態BとAの間を行き来するが、状態Bのほうが $1/r$ 倍選択されやすい。 中央の頂点が状態Bにいるほうが $1/r$ 倍選択されやすく、更に中央の頂点の状態にかかわらず周囲の頂点のうちBが選択されるほうが $1/r$ 倍多いため、結果としてBが $1/r^2$ 倍選択されることになる。

つまり

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

適応度がMoran過程に比べて $r^2$ に増幅されていることがわかる。この議論の詳細な解析は、冒頭の参考文献 (Broom & Rychtar)を参照してもらいたい。

ネットワーク上のゲーム

各頂点は 協力 (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戦略が優位になる。