View source for Aritalab:Lecture/NetworkBiology/Markov Chains/Birth-death Process

{{Lecture/Header}}

==参考文献==
* Lieberman E, Hauert C, Nowak MA (2005) Evolutionary Dynamics on Graphs. [http://www.ped.fas.harvard.edu/people/faculty/publications_nowak/nature05.pdf Nature 433:312-316]
* Nowak MA (2006) [http://www.amazon.com/Evolutionary-Dynamics-Exploring-Equations-Life/dp/0674023382 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 [http://rspa.royalsocietypublishing.org/content/464/2098/2609.abstract?ijkey=685441775235bf5fdd5a71bb66ef1743749f7f92&keytype2=tf_ipsecsha ''Proc R Soc A'' 464:2609-2627]

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

状態 ''i'' から出発して状態 ''N'' に到達する確率を <math>x_i</math> と書く。これを固定確率と呼ぶ。

<math>
\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}
</math>

ここで記法

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

を導入すると

<math>\textstyle y_{i+1} = \gamma y_i, \quad \sum^N_{i=1} y_i = 1</math>

を得る。Moran過程では通常、少ない個体数から ''N'' 個体にまで増殖するケースを考えるため、<math>\gamma</math>は 1 より小さい値 ( つまり ''&beta;'' > ''&alpha;'' ) と考えてよい。

<math> y_1 = x_1, \ y_2 = \gamma x_1, \ y_3 = \gamma^2 x_1, \cdots, y_n = \gamma^{N-1} x_1 </math>

を足し合わせると

<math>
x_1 = \frac{1}{\sum^{N-1}_{j=0} \gamma^j} = \frac{1 - \gamma}{1 - \gamma^N}
</math>

また

<math>
\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}
</math>

とあわせて

<math>
x_i = \frac{\sum^{i-1}_{j=0} \gamma^j}{\sum^{N-1}_{j=0} \gamma^j}
=  \frac{1 - \gamma^i}{1 - \gamma^N} 
</math>

となる。

===式の解釈===
''N'' 個体の集団において、全ての個体が最初タイプ　A　であるとする。
このときタイプ　B　という突然変異が　1　個生じたとする。
タイプ　B　が <math>\gamma < 1</math>　を満たす、すなわち　<math>\alpha</math>　が　<math>\beta</math>　よりも大きくて個体数を増やす傾向にあると仮定する。
タイプ　B　が何らかの優性種であると考えてもよい。

初めに　1　個生じたタイプBが集団 ''N'' の中に固定される（集団全体をカバーする）確率はモラン過程より　<math>x_1 \xrightarrow{N \rightarrow \infty} 1 - \gamma</math> になる。
つまり、A より10倍優位な変異体であっても (<math>\gamma = 1/10</math>)、集団全体にその変異が保持される確率は 9/10 でしかない。
通常の遺伝子変異で生存戦略に極めて優性となるケースは非常に稀だと考えられる。そうすると、優性な変異が集団全体に保持される可能性はそう大きくない。

===木村資生の中立進化説===
木村資生（きむらもとお）は1970年代に活躍した国立遺伝学研究所の研究者である。
ほとんどの遺伝子の突然変異は、自然淘汰において有利でも不利でもないとする[http://www.nig.ac.jp/museum_080501/evolution/C/bunsi-02.html 中立進化説]は当時大きな反響を呼んだ。

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

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

===モラン過程と同等なネットワーク===
初めに1個体であった B が全体に広がる確率がMoran過程と同じ <math>\textstyle \frac{1-r}{1-r^N}</math> であるための条件を考えよう。
現在 ''m'' 頂点が状態 B、残りが状態 A と仮定する。
B のほうが A より <math>1/r</math> 倍選択されやすい場合、''N'' 個の頂点のうち、状態Aのものが選択される確率は <math>m/r</math> だけ個体数が増えていることに等しいから <math>1 /(N-m + m/r)</math> となる。
状態 B の頂点が選択される確率は<math>1/r(N-m+ m/r)</math>である。

頂点 ''i'' が状態Bのとき <math>v_i = 1</math>, 状態Aのとき <math>v_i = 0</math>と書こう。

B 個体が ''m'' + 1 に増える確率は、状態 B の頂点 ''i'' (<math> v_i = 1</math>) が選ばれて、隣接する状態 A の頂点 ''j'' (<math> v_j = 0</math>) を上書きすればよい。

<math>\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} </math>

同様に、B 個体が''m'' &minus; 1 に減少する確率は

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

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

<math>\textstyle \frac{B_{m\rightarrow (m-1)}}{B_{m\rightarrow (m+1)}} = r</math>

が成立せねばならず、すべての ''i'', ''j'' について

<math>\textstyle \sum_i\sum_j p_{ij} v_i (1-v_j) = \sum_i\sum_j p_{ij} (1-v_i) v_j</math>

となる。

個体 ''k'' ただ一つが B である特殊例 <math>v_k = 1,\ v_i = 0\, (\forall i \neq k)</math> を考えると <math>\textstyle \sum_j p_{jk} = \sum_j p_{kj} = 1</math> が必要条件となる。
ただし右側の等号は、考えているネットワークがマルコフ連鎖であるため、頂点から出る辺確率の和が1であることから導かれる。
これがすべての頂点について成立するため、定常分布での存在確率が全頂点において等しいときには（別の言葉でいうと、全頂点が対等な関係にある場合は）、Moran過程と同じプロセスになる。

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

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

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

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

ここでは星型 (star) のネットワークはMoran過程より固定確率が高くなることを説明する。
結果から先に言うと、固定確率は <math>\textstyle x_i = \frac{1-r^{2i}}{1-r^{2n}}</math> であらわされる。

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

つまり

<math>\textstyle \frac{B_{m\rightarrow (m-1)}}{B_{m\rightarrow (m+1)}}
= r^2 </math>

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

==ネットワーク上のゲーム==

各頂点は 協力 (C: Cooperate) か裏切り (D: Defect) のいずれかの戦略を取るとする。
協力者は隣接頂点それぞれにコスト (cost) ''c'' を支払い、隣接する頂点は利得 (benefit) ''b'' を受け取る。
裏切る場合はコストを払わず、隣からの利得のみを受け取る。
頂点の適応度 (fitness) を支払いと利得の合計で判断し、戦略を更新する方法として以下の2通りを考えよう。
# 出生死亡過程: 単位時間毎に適応度に応じて頂点を選択し、隣接頂点をランダムに選んで同じ戦略へとコピーする。
# 死亡出生過程: 単位時間毎に頂点がランダムに選ばれて死亡する。隣接する頂点がそれぞれの適応度に応じて自分の戦略をその頂点にコピーする。

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

;出生死亡過程の場合
:C戦略とD戦略をとる頂点の境界部分にある2頂点を考慮すればよい。
C戦略の端にいる頂点は、隣接する2頂点分のコストを払って利得を片側からしか得ないので合計は <math>b-2c</math> になるが、D戦略の端にいる頂点は利得 ''b'' を得られるので、いかなる ''b, c'' であってもD戦略が優位になる。
;死亡出生過程の場合
:C戦略とD戦略をとる頂点の境界部分から2つ離れたところまでの合計4頂点を考慮する。
D戦略を取る側は、C戦略と隣接する頂点が ''b''、そのまた隣は ''0'' になる。
C戦略を取る側は、D戦略と隣接する頂点が <math>b-2c</math>、そのまた隣は <math>2(b-c)</math>になる。
従って <math>b/c > 2</math> の場合はC戦略が優位。
一般化すると、次数 ''k'' の格子であれば <math> b/c > k</math>の場合にC戦略が優位になる。