読者です 読者をやめる 読者になる 読者になる

Life is Beautiful

主に進化生物学の理論のブログです。不定期更新予定。

横着者のための、モンテカルロ進化シミュレーション

少し前まで、確率シミュレーションが好きではなくて、数理モデルを前にして力学系の手法だけで料理を試みていた僕ですが、友人からの熱烈な進言をうけて、C++でシミュレーションを回す練習を始めてみました。

集団ベースモデルによるHamilton and May (1977) シミュレーション

Hamilton and May modelとは、パッチ状の構造をもったメタ個体群における移動分散の進化を論じたモデルです。

たくさんのパッチを考えます。各パッチは1個体だけ収容できるとします。各パッチ内の個体は子供をたくさん産み、死にます。 そして子どもたちは確率  d で移動して他のパッチへ移動しますが、残りの  1-d は生まれたパッチに留まるとします。 ただし移動した子どもたちは確率  1-c でしか生き残れないものとします。逆に言うと、移動した個体のうち、割合  c は死んでしまうというものです。

さて、進化的に安定な移住率  d^* はいくらでしょう?

移動分散率  d_\mathrm{m} をもつ突然変異個体が、移動分散率  d をもつ集団で稼げる適応度は \begin{align} W = \frac{1-d_\mathrm{m}}{1-d_\mathrm{m}+(1-c)d} + \frac{(1-c) d_\mathrm{m}}{1-cd} \end{align} ですので  W \leq 1 が任意の  d_\mathrm{m} で成立するための  d の条件を求めると、 \begin{align} d^* = \frac{1}{1+c} \end{align} と決まります。

個体ベースモデルはいささか時間がかかりますので、次のようなアルゴリズムにしたがって、集団が単型であるという仮定のもとで進化の方向を予測*1してみましょう。

  • メルセンヌ・ツイスタによって実一様乱数  d\in (0,1) をとってくる。この値を、集団における移動分散率の初期値とする。
  • メルセンヌ・ツイスタによって、[0,1] の実一様乱数 r_\mathrm{m}\in (0,1)をとってくる。 p=0.025 という小さい値にたいして、  r_\mathrm{m} > pなら、突然変異は起こらないが、そうでない場合、在来型が一個体だけ変異型になる *2
  • 在来型  d は突然変異の結果、  d_\mathrm{m} = d + \delta に変異する。ただし  \delta は、メルセンヌ・ツイスタにしたがって、実一様乱数として選び( \delta \in (-0.025,0.025))、もちろん  d_\mathrm{m}\in [ 0,1 ]となるように制限する。
  • 変異型の適応度  W を計算し、これが1より大きいなら( W>1)変異型がのっとる。そうでないなら在来型がとどまる。*3
  • これを1000世代くりかえし、収束先を調べる。
  • なお、パッチの数は無限大とし、移動のコスト  c =0.6 とした(したがって、ESSの予測値は 5/8=0.625になる)。

その結果がこちら。10回ほど走らせた結果をのせます。横軸は世代、縦軸は移動分散率です(つまり移動分散率の変化の時間経過をトラックしている)。初期値がバラバラなのは、ランごとに乱数を振っているためです。

f:id:lambtani:20160809213210j:plain

ほぼすべて、0.625に収束してくれていますね。めでたしめでたし。なお、10回のランにかかった時間は、数秒でした。

後記: なお、このアルゴリズムには、収束安定性が保証されないという問題がある。すなわちESSに本当に到達できるかどうかは考えずに、シミュレーションをまわしている。そもそも、収束安定性が成り立たないような場合、「集団は常に単型」であるという仮定が、(実証的にはもちろんのこと)理論的にも極めてアヤシイものになってしまう。

*1:正確には、有利な遺伝子が広まるスピードは非常に早く、有利なタイプは一世代で、在来型を駆逐してしまう、という仮定のもとでの、形質値の進化の予測

*2:実は、このシミュレーションにおいては、このような場合分けステップは不要。だが個体ベースのシミュレーションをまわす場合は、集団中のすべての個体に対して、変異するかどうかの問いかけを行なう必要がある

*3:個体ベースの場合、各タイプの各島での適応度成分を計算し、それに比例した確率でパッチを占有できる、といったような計算を行なう