Life is Beautiful

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

真剣に取り組みたい人のための、モンテカルロ進化シミュレーション

C++勉強が進みます。今日はタカハトゲームの進化をシミュレートをしてみましょう。

タカハトゲームとは、強気に振る舞うか、弱気に振る舞うかというゲームだと考えて構いません。

  • タカ同士が出会うと、死闘を繰り広げ、資源もコストも分け合います。利得は  \frac{V-C}{2}
  • タカとハトと出会うと、ハトは問答無用で資源を得 (利得は  V)、ハトは資源をえることができず、利得は0です。
  • ハト同士が出会うと、仲良く資源をわけあい、利得は \frac{V}{2}です。

利得行列は次のようにかけます:

\begin{align} \begin{pmatrix} \frac{V-C}{2} & V\\ 0 & \frac{V}{2} \end{pmatrix} \end{align}

個人の戦略(形質)  p は、タカとして振る舞う確率だとします(よって  1-p はハトとして振る舞う確率)。このような  p は混合戦略と呼ばれます。 いつもどおり、集団が戦略  p で占められているとします。このとき、変異体  p_\mathrm{m} が侵入できるかを計算してみましょう。

これに基いて利得を計算しますが、ここでカンタンのために次のように仮定します:

  • 産子数は非常に大きい;
  • ゲームをプレイするのは  N 人の大人たちである;
  • 利得は大人の産子数に寄与する; *1

これにより、突然変異個体のこども一個体あたりの利得  \pi_\mathrm{m} は、

\begin{align} \pi_\mathrm{m} = p_\mathrm{m} \left(\frac{V-C}{2} p_{-\mathrm{m}}+ V (1-p_{-\mathrm{m}}) \right) +(1-p_\mathrm{m}) \left( 0+ \frac{V}{2}(1-p_{-\mathrm{m}}) \right) \end{align}

となります。ここで、 \begin{align} p_{-\mathrm{m}} = \frac{\overline{p} - \frac{1}{N} p_\mathrm{m} }{1-\frac{1}{N}} \end{align} は突然変異個体を除いた、集団における平均値です。

利得が適応度に及ぼす効果をスケールするため、淘汰の強さを  s として、変異体の適応度  W を、 W = 1+s\pi_\mathrm{m} とします。 s が大きくなればなるほど、利得が適応度に及ぼす効果は強くなります(強淘汰)。  s=0 の場合は、利得が一切進化に影響を及ぼさないという、中立的な状況に相当します*2

さらに、有限集団を仮定しているために、在来型pの適応度が突然変異型の戦略に依存してしまいます。在来型の利得は

\begin{align} \pi_\mathrm{w} = p \left(\frac{V-C}{2} q+ V (1-q) \right) +(1-p) \left( 0+ \frac{V}{2}(1-q) \right) \end{align}

です。ここで、 \begin{align} q = \frac{1}{N-1}p_\mathrm{m} +\frac{N-2}{N-1}p \end{align} は在来型が対戦する相手の平均戦略です。よって突然変異型の相対適応度は、

\begin{align} w = \frac{N(1+s\pi_{\mathrm{m}})}{1+s\pi_{\mathrm{m}}+(N-1)(1+s\pi_{\mathrm{w}})} \end{align} となります。これにもとづきESSを計算すると、*3

\begin{align} p^* = \min \left( \frac{VN}{C(N-2)}, 1 \right) \end{align} が得られます。つまり V>C は(つまりタカ同士の喧嘩でコストを上回る資源が得られることは)、常にタカ派として振る舞うのがESSであるための十分条件ということになります。

これを個体ベースモデルで確認してみましょう。アルゴリズムは以下。まずは、集団サイズを決めて、初期の形質値をランダムに振ってやることから始まります。今回は、個体  i の形質 p[i] を配列ではなくベクトルで準備してやりました。

  •  N 個体で集団がきっちり埋まっているとする。
  •  i に対してメルセンヌ・ツイスタによって一様乱数  r \in(0,1) をひっぱってきて、これが  r \leq R を満たすなら、 p_i の値は  p_i + \delta に更新される(ただし今回は  R=0.25 と大きめに設定した)。ここで  \delta は、 \delta\in (-0.07,0.07) を満たす、メルセンヌ・ツイスタによる一様乱数。変異後の値はもちろん0以上1以下になるように制限をする。
  •  p の集団平均  \overline{p} を計算する。
  •  \overline{p} から、 i 番目の個体の利得  \pi_i を計算し、これから適応度  W_i = 1+s\pi_iを計算する。
  •  W_i から、平均適応度  \overline{W} を計算してやる。
  • 2個体がランダムに選ばれ、適応度の大きい個体が、小さい個体を排除する。
  • 元のサイズ  N に戻る。これを1024世代繰り返す。
  • パラメータは  V=0.1, C=0.2, N=50, s=0.01 とした。よって予測される ESS は、 0.52 ということになる。

後記:ミスを発見。。修正しました。

結果が下。8回まわした結果をのせます。点線はESS予測値。一つ一つの色は、ラン。戦略平均の時間変化をトラックしています。

f:id:lambtani:20160818160222j:plain

やはりシミュレーションの回数が少ないと、ばらつきが大きいですねー。そしてMoranプロセスを採用しているので、完全に世代交代するまでにかかる時間の期待値は、Wright-Fisher の  N 倍。100個体も準備しているので、2000世代まわしても、実質的にWright-Fisherで20世代まわした程度の結果ということになります。 しかしコードはMoranプロセスのほうがだいぶんラクですね。どこでラクするか。

ちゃんと収束してくれて何より。なお、適応度が最大の個体が、適応度最小の個体を排除する、というアルゴリズムを組んでも、同じ結果がえられます。*4

*1:子どもたちがゲームをプレイし、利得は子どもの生存率に寄与する 、と仮定しても、このモデルにおいては同一の結果が得られる。 場合には、無限集団極限  N\to\infty に相当する結果が得られる(後述)。ただし一般には、survival effect と fecundity effect とは、厳密に区別をせねばならない

*2:ここで、 W = 1-t+t\pi_\mathrm{m} とする流儀もある。この場合、  s t のオッズであり、 t \in [0,1) である。Sigmund (2010)など。

*3:本当はこれは狭義ESSではなく、Nash均衡である。なぜなら、ESSはあくまでも十分小さい規模  \epsilon での侵入を想定しているのだが、ここでは侵入規模は \epsilon=1/N となっているからである。もちろん無限集団を仮定すれば、すみやかに狭義条件が得られ、その場合のESSは V/C になることが確かめられる。これは、産仔数が非常に大きいという極限においては、子供たちがゲームをプレイするという仮定からも得られる結果である。産子数が有限の場合はさらなる計算が要求され、最終的によりタカ戦略にバイアスしたESSが得られる

*4:戦略値の分散も十分小さく保たれていること、すなわち進化的分岐は起こっていないことも、予備シミュレーションで示してあります。

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

少し前まで、確率シミュレーションが好きではなくて、数理モデルを前にして力学系の手法だけで料理を試みていた僕ですが、友人からの熱烈な進言をうけて、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:個体ベースの場合、各タイプの各島での適応度成分を計算し、それに比例した確率でパッチを占有できる、といったような計算を行なう

定常状態という仮定

定常状態。平衡状態とも言う。理論・数理生物学の研究者はこれが本当に好きだ。批判ではない。ただの事実である。

その事実には原因がある。定常状態の性質を調べないと、パタンが見いだせないからである。平衡状態にないような、何かの時間変化のことを、transient dynamicsと呼ぶ。たとえば漸近的な平衡状態に到達していない系はそうである。周期的な挙動を示すシステムについては、ここでは深く考えない。

Transient dynamicsでのある時刻 tにおける振る舞いは、連続的に動かせるパラメータに関して一般的には非可算無限個、存在し得る。非可算無限とは、1つ、2つ、…10億4391つ、という整数での数え上げが不可能なものであり、無限の中でも濃度の強い無限である。このような場合、人の頭の中では振る舞いをパターン化できない。ところが平衡状態に徹すると、パターン化は明解であることが多い。たとえば集団遺伝学的な意味で、定常状態の予測として「(1)突然変異が侵入できずに絶滅する;(2)野生型が駆逐され、突然変異が固定する;(3)突然変異と野生型とがある割合で混在・共存する;(4)突然変異と野生型の頻度が振動する」といった有限個のパターンにわけられ、かつその条件を網羅的シミュレーション・網羅的な数値計算・解析的な計算、のどれかで解明できれば、系を包括的に理解できたような気がする。これが、「力学系」の基礎となる考え方である。

transient dynamicsはずばり「時(たとえば時刻t)と場合(パラメータ)による」という結論にしかなっていない。それは、「動態」に関するあらゆる議論の前提である。

ではそのように「時と場合による」のだとすると、あるtransient dynamicsにおいて状態Aが達成される確率を導入したくなる。あるいは、そのような状態Aが、「あらゆるとりうる状態」の集合{ \mathcal{S}}= \{ A,B,C,D,...\}において、どの程度の割合を占めているのかを考えると良さそうである。すると、{ \mathcal{S}}に含まれる要素(つまり状態)が無限個である場合には、状態 Aの占める割合は0である。そのような性質は、generic*1ではないという。ちょっと状況が変われば変化するような性質は非常にデリケートなバランスの上に成立しているはずであり、実現される望みは薄い、という考え方である。*2

しかしいっぽう、(僕が興味のあるような)野外でのダイナミクスは、ほぼすべてtransientなはずだ。つまり平衡状態にはない。ということは、やはりtransientなダイナミクスを考えねばなるまいのではないか!!?という気持ちが湧く。

それでも定常状態を仮定・考察する理由を、ポジティブにいくつか検討してみよう。

  • まずそもそも僕らは、分類したい。脳みそのキャパシティーは有限である。有限の特徴で、何かを理解したい。これにより、分類されて共通のカテゴリに入れられたものを(完全に同一視、まではしなくとも)アナロジーで以って理解して、脳みそに収納しやすくなる。たとえば12×12=144を覚えたら、102×102=10404であることはすぐに覚えられる。古文単語が覚えにくいなら、ゴロあわせで覚えてしまう。これは、ゴロ合わせで覚えることで、イメージされるものに現実(単語)を逆に結びつけるという、高等テクニックである。褒め過ぎか。

  • また、定常状態は再現性が得やすい。たとえば野外で得られたパターンを、微生物などを用いて再現する場合に、そのミニ・コスモでは素早く定常状態に収束するかもしれない。操作実験の結果を野外でのパタンと比較するとなると、やはり定常状態を仮定することで比較が容易になる。そうでないと、比較することができない。

  • 更に、定常状態にもクラスがある。その中に、準定常状態というものがある。他のダイナミクスに比べて速やかに進むような性質については、はじめから「定常状態」という静的なものを考えてしまうという近似である。これは、力学系の変数を減らすことができる。おお、便利。だがこれらは便宜的な理由であり、弱いかもしれない。

  • そもそも、野外でのダイナミクスにも、似たような傾向が見られている。つまり種をまたいで*3、あるいは地域をまたいで、似たような傾向が見られていることがあるのである。これは「パターン」と言ってよさそうである。平衡状態を仮定するとパターンを認識しやすくなるのだから、ただの横着ではなく、「分類する」という目的にかなった仮定である!

  • 最後に、野外でのダイナミクスはそもそも定常状態にあるかもしれない。たとえば自然選択は、長い時間をかけて、生物の性質を形作る。その自然選択が何世代にも何万世代にもわたって作用してきたとなると、その結果は、なんらかの定常状態に収束している可能性がある。ただ、たとえば遺伝率の高い形質や、可塑性の顕著な形質には、注意を払わねばなるまい。

このような理由で、数理生物学では定常状態を仮定することがそうそう馬鹿げているようには思えない。もちろんtransient dynamicsは、知りたい究極的な性質になりうる。そうであるなら、その理由を明確にしたうえで、適切に取り組むべきである。そうでないなら、分類思考に基いて、パターンをはっきりさせたほうが、明確なメッセージを提示できるのだと思う。

*1:通有的、という言葉が充てられている。国府ら、2000:

https://www.amazon.co.jp/dp/4254126727www.amazon.co.jp

*2:僕らが知っているような「通有的な平衡点」、すなわちパラメタの摂動に対してロバストに局所安定な平衡点とか、逆に不安定な平衡点とかは、双曲型と言われる。厳密な定義は、固有値によって、すなわち局所的な性質によって与えられる。

*3:これはこれで、生物種の性質の間に無視できぬ系統的な似通りが存在することがあって、単純に「種をまたいだ共通のパターン」を定性的に受け入れるのがまずいことがある。Felsenstein 1985: http://www.journals.uchicago.edu/doi/10.1086/284325

Dropboxの有料プランに加入した

ついにDropboxの有料プランに加入してしまいました。一括で払うと一年で12000円。いままではあらゆる手を尽くして容量拡大の努力を行なっていたのですが、

  • 引用論文の管理ーーPDFリンクのはられたBibdeskによる管理
  • 発表資料の管理ーー年月と学会名でディレクトリわけ
  • 投稿論文の管理ーー執筆順(これはいつ10に到達するだろう・・・)

Dropboxで一括で行なうことにしたため、限界に達しました。

これまでは容量の90%をすでに使用していて、アラートで「限界が近いよ!」と注意されていたのですが、今後のことを考えると、まあ悪く無いかなと。バックアップのサービスも手厚く受けられるようだし。

問題は、Evernoteとのすみ分け。いまやチャットルームのためだけに契約している感じすらある。。便利なのだけど、Evernoteはそろそろ見切りをつけたほうがいいのかも知れない。

引用文献リストでDOIを非表示にする

ISSNはDOIは、文献のグローバルな個体番号のようなもので、例えばDOIを

dx.doi.org/

に続けて打ち込むだけで論文URLに飛べる、便利なものです。

しかしDOIやISSNを論文の引用文献から消せないかと思っていましたら、方法が見つかりました:

tex.stackexchange.com

たとえばURLを隠す場合。bstファイルを新しく作ります: plainnat.bstファイルを作業ディレクトリにコピーし、それをテキストエディタで開きます。 そして

FUNCTION {format.url}
{ url empty$
    { "" }
    { URLのフォーマットを宣言する部分 }
  if$
}

となっている部分を

FUNCTION {format.url}
{ url empty$
    { "" }
    { "" }
  if$
}

と書き換えればおしまい。DOIについても、FUNCTIONの部分からDOIのフォーマット定義を消してしまえばよいだけです。

もしものときのため、この作業は、あくまでbstファイルをコピーしたうえで行なうようにしましょう。

水素水で営利を追求する企業への不買

非常に良いエントリーを見つけました。

金儲けのために水素水販売に手を出した8つの大手企業 | netgeek

私は、伊藤園PanasonicSharpなど、消費者を騙して利潤を得る企業の商品を購入しません。また、他者から何か物品の購入アドバイスを受けたら、これらの企業の商品を決して勸めないことにします。

複数の数式を引用したい

複数の数式をLaTeX内で引用・表示したい場合はどうすればよいか。

いま、4つの式A,B,C,Dを別個の

{equation}

環境で表示していて、A,B,C,Dをいっきに引用しようと思ったら、詰みました。

で、調べてみると素晴らしいパッケージが。

その名もCleveref CTAN: Package cleveref

\usepackage{cleveref}

でインクルードして

\Crefrange{A}{D}

とすると、(たとえば)こんな感じ(数式Aが、数式7に相当します): f:id:lambtani:20160426194605p:plain

かしこい。マニュアルを読んでみると、もう少し(いやもっと)いろいろできそう。