双安定進化ゲームの確率的ダイナミクスに対する空間

《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
東京大学大学院 理学系研究科 物理学専攻 修士論文
双安定進化ゲームの確率的ダイナミクスに対する
空間自由度の影響
Spatial effect on stochastic dynamics of
bistable evolutionary games
曽弘博
Kohaku So (Hongbo Zeng)
東京大学大学院 理学系研究科 物理学専攻
2014/2/3
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
概要
進化ゲーム理論は、生物や文化などの「自己複製し、進化するもの」を記述する理論である。元来提唱された
分野である生物科学はもちろんのこと、経済学をはじめとした社会科学の研究にも広く用いられている。
進化ゲーム理論は当初、無限に大きく、かつ well-mixed な集団での決定論的過程として進化を記述していたが、
この仮定は現実の集団での進化の過程を記述するには不適切な場合もあった。そこで、有限サイズでかつ構造を
持った集団での確率過程として進化を記述する研究がなされるようになり、多くの研究から、集団の構造が進化
ゲームに大きな影響を及ぼすことが示された。しかしながら、これらの研究の多くは系の定常状態などの静的な
性質に着目したものであり、動的な性質に集団の構造がどのような影響を及ぼすかについては、まだ多くは知ら
れていない。
本研究ではこのような背景を踏まえ、進化ゲームの確率的な dynamics に空間自由度が及ぼす影響を探った。具
体的には、2 つの準安定状態を持つ進化ゲームを考え、これらの準安定状態間の遷移という動的な現象に空間自由
度の有無がどのような影響を与えるかを調べた。
解析にあたっては進化の過程を確率過程として定式化し、準安定状態間の遷移確率を経路積分によって表した。
そこに半古典近似と呼ばれる近似手法を適用することにより、各準安定状態の寿命を空間自由度がない場合とあ
る場合とでそれぞれ評価した。特に、寿命が集団サイズにどのように依存性するかを詳しく評価した。空間自由
度がない場合については、他の種々のパラメーターに対する依存性の評価や、master 方程式の数値計算結果との
比較検討も行った。
結果として、空間自由度の有無が、準安定状態の寿命の集団サイズへの依存性を質的に変えることが分かった:
空間自由度がない場合は各準安定状態の寿命は集団サイズの増加に伴って指数的に増える。空間自由度がある場
合では、片方の準安定状態の寿命は集団サイズの増加に伴って指数的に増えるのに対し、もう片方の準安定状態
の寿命は (集団サイズがある一定値より大きくなれば) 集団サイズに依らない一定値をとるようになる。
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
謝辞
本研究を遂行するにあたり、丁寧に御指導をしてくださった加藤岳生准教授に深謝いたします。研究室のテーマ
とは一見関係のないテーマを選んだにもかかわらず熱心な御指導を賜ったことは、感謝の念に堪えません。
総合研究大学院大学の大槻久助教には、研究テーマの選定や結果の解釈にあたって大変示唆に富んだ有意義な議
論をしていただきました。このような議論がなければ、本論文は完成しなかったでしょう。深く感謝いたします。
また、加藤研究室のメンバーをはじめとする物性研究所・物理学専攻の皆様には、日々の議論やゼミなどの学
問的な面をはじめ、様々な面で大変お世話になりました。ここに謝意を表します。
最後になりましたが、私が落ち着いて研究に専念することができたのは、心の支えとなってくださった多くの
方々のおかげです。一人一人の名を挙げることは紙幅の都合によりできませんが、ここにお礼の言葉を述べたいと
思います。本当にありがとうございます。
1
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
目次
第 1 章 序論
4
第 2 章 進化ゲーム理論
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
6
7
replicator 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2×2 進化ゲームの分類 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
7
2.1
2.2
進化ゲーム理論とは
決定論的進化ゲーム
2.2.1
2.2.2
2.3
2.2.3 進化ゲームの例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
確率的進化ゲーム（well-mixed） . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 遷移ルール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2
2.3.3
2.4
2.5
確率的記述と決定論的記述の関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
突然変異と定常分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
確率的進化ゲーム（structured population） . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 遷移ルール . . . . . . . . . . . . . . . .
structured population における進化ゲーム . . .
2.5.1 正方格子上での囚人の dilemma . . . . .
2.5.2 格子上での snowdrift ゲーム . . . . . . .
2.5.3 Rock-Paper-Scissors ゲーム . . . . . . .
2.5.4
2.5.5
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
一般の structured population における進化ゲーム . . . . . . . . . . . . . . . . . . . . . . . 19
まとめと課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
第 3 章 確率過程と経路積分
3.1
3.2
3.3
21
確率過程の経路積分による記述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.1 簡単な例:崩壊過程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.2
3.1.3
一般の場合（1 変数） . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.1
3.2.2
半古典近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
多変数の場合（反応拡散系）
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1.4 結果のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
半古典近似とその応用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
崩壊過程の半古典近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
stochastic population の寿命：bounce 解の手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.3.1 決定論的解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3.2 確率過程の解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
第 4 章 モデルと解析
4.1
4.1.1
4.1.2
4.2
37
モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
well-mixed モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
spatial モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
well-mixed モデルの決定論的解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
4.3
4.4
well-mixed モデルにおける準安定状態間の遷移 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.3.1 経路積分表示 . . . . .
4.3.2 半古典近似 . . . . . .
spatial モデルの決定論的解析
4.4.1 rate 方程式 . . . . . .
4.4.2
4.4.3
4.5
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
連続化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
定常解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
spatial モデルにおける準安定状態間の遷移 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5.1 経路積分表示 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5.2
4.5.3
半古典近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
bounce 解の数値計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
第 5 章 結果と考察
5.1
55
well-mixed モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.1.1 集団サイズ依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.1.2
5.1.3
w 依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
利得行列依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.2
spatial モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.1 作用の集団サイズ依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.3
考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
第 6 章 まとめと展望
61
6.1
本研究のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.2
展望 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.2.1
他の系への適用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.2.2
6.2.3
より一般の構造上での進化ゲーム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
より精緻な近似手法による解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
付 録 A 確率過程の近似手法
63
A.1 master 方程式の van Kampen 展開 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
A.1.1 崩壊過程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
A.1.2 one-step process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
A.2 半古典近似の詳細 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
付 録 B 空間モデルの Hamiltonian と半古典近似
B.1
B.2
B.3
B.4
70
Hamiltonian の連続化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
運動方程式の導出 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Hamiltonian の保存 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
緩和法による数値解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
第 1 章 序論
動物、植物、細胞、ウイルス、言語、文化—我々の周りには生物・非生物を問わず、自らを「複製」し、集団中
に広まっていくような要素が数多く存在している。集団が複数種の要素からなり、異なる種の要素が異なる速さで
複製を行うとき、そこに「進化」の過程が生じる。生物は自らの子を産み、より子を多く産むものが集団中に広
がる。また、文化は人間の間で模倣を通じて増殖し、より模倣されやすいものほど集団中に広まる。
進化ゲーム理論は、このような進化の過程を数理的に説明する理論である。進化ゲーム理論においては、頻度
依存淘汰、つまり要素の複製の速さが、周りに存在する他の要素の頻度に依存するような状況を取り扱う。進化
ゲーム理論は、広範な進化プロセスを記述する一般的な理論として、これまで多くの生物学分野の研究で用いら
れてきた。また、進化ゲーム理論は経済学・政治学の分野で発達してきたゲーム理論の枠組みを応用したもので
あるが、逆に近年これらの社会科学の分野にも大きな影響を与えている。
進化ゲーム理論は 1970 年代から活発に研究が行われているが、当初の研究では、十分に要素数の大きい集団が
考えられ、かつ空間的な構造は考慮されなかった (空間的な構造を持たない集団を well-mixed な集団と呼ぶ。)。
しかし、現実の集団を考えるときには、これらの仮定は不適切な場合がある。
まず、現実の集団では空間構造が無視できない。例えば、移動能力が限られた生物は、遠くの個体と比べて近く
の個体と相互作用する可能性が高い。また、人間は赤の他人よりも知人・友人と相互作用する可能性が高いため、
人間関係のネットワーク構造が重要になる。このような空間構造・ネットワーク構造を導入した研究は 1990 年代
から盛んに行われ、これらの構造が進化のプロセスを質的に変えることが明らかになった。
また、現実の集団の要素数は有限である。要素数が十分に大きい場合には、要素の自己複製過程は決定論的な
方程式 (replicator 方程式 (2.2.1 項)) で捉えられるが、要素数が少なくなるにつれて確率性が無視できなくなる。
確率的効果は、進化の過程に大きな影響を及ぼす: 例えば、集団の要素数が十分に大きい時には安定して存在で
きた種が、要素数が少なくなると絶滅してしまうことがある。このような確率性の影響が注目を集め始めたのは、
2000 年代に入ってからのことである。
以上のように、空間構造と確率性は、進化ゲームの静的・動的性質に大きな影響を与えるため、この両者を考
察することは重要である。特に系の種の分布が大きく変化するような現象において、この 2 つの影響は顕著に現
れると期待される。そのような現象の最も典型的な例として「準安定状態間の遷移現象」が挙げられる。これは、
2 つ以上の安定状態が存在するような系において、系の状態が確率的にこれらの安定状態の間を遷移する現象で
ある。この遷移現象はごく稀にしか生じないが、系全体の性質が大きく変化するため重要な現象である。しかし、
準安定状態間の遷移に関する理論研究は、まだ多くはなされていない。
本研究では以上の背景を踏まえ、進化ゲーム理論での準安定状態間の遷移に対する空間構造の効果について、理
論的研究を行う。解析にあたっては、進化の過程を確率過程として定式化したのち、遷移確率を経路積分によって
記述し、「半古典近似」と呼ばれる近似手法を適用する。この手法によって、準安定状態間の遷移確率が、空間構
造の有無によってどのように変化するかを探る。
以下に、本論文の構成を示す:
• 2 章では進化ゲーム理論の review を行う。進化ゲーム理論の定式化、特に決定論的定式化と確率的定式化の
双方を概観したのち、空間自由度が進化ゲームにどのような影響を及ぼすかについて、先行研究をまとめる。
• 3 章では、解析に用いる手法について述べる。具体的には、経路積分による確率過程の記述の方法と、それ
に対する半古典近似について述べる。
• 4 章では、本研究で用いる進化ゲームモデルの設定を述べ、第 3 章の手法に基づき遷移確率を評価する。
4
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
• 5 章ではこれを承けて、遷移確率の種々のパラメーターに対する依存性と空間構造の影響についての解析結
果をまとめ、それらが示唆することについて論じる。
• 6 章には、本研究のまとめと今後の展望を述べる。
• 付録 A では 3 章で用いた確率過程の近似手法の詳細を説明する。
• 付録 B では 4 章の解析の詳細について説明する。
5
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
第 2 章 進化ゲーム理論
この章では、「進化ゲーム理論」とは何かを概観した後、進化ゲームの 2 つの定式化（決定論的記述、確率的記
述）について述べる。それらを踏まえて、近年の研究、特に、空間自由度が戦略の広まり方に与える影響につい
て得られている知見を紹介する。
2.1
進化ゲーム理論とは
人間社会や生態系においては、構成要素 (個人や生物個体) が他の要素との相互作用を通じて、意思決定や行動
様式の選択を行うと考えられる状況が頻繁に現れる。ゲーム理論はこのような「意思」を持った要素の相互作用
を分析する学問分野である。ゲーム理論で考察する問題は、一般に次のような形をとる。
N 人の人間（player と呼ぶ）がいて、それぞれが s 個の選択肢 (戦略と呼ぶ) を持っているとする。
個々の player がどの戦略を取るかによって、個々の player の利得 (payoff) が決まるとする。このと
き、個々の player はどの選択肢を選ぶことになるか？
ゲーム理論は、Neumann と Morgernstern らによる研究 [1] に端を発し、経済学・政治学の分野で古くから盛んに
研究が行われてきた。
古典的なゲーム理論は、その基本的な仮定として「合理的な player」の存在を基礎としている。
「合理的な player」
とは、あらゆる可能な選択肢すべての利得を知っており、かつ他のすべての player の行動を推定できる能力を持
つ player のことである。ゲーム理論は人間の意思決定に関して、ナッシュ均衡などの多くの重要な概念を生み出
した。しかし、「合理的 player」の仮定は、現実の人間の行動を考察する上で、必ずしも妥当なものではないとの
批判がある [2]。
他方、ゲーム理論は社会科学の枠を超えて、自然科学にも大きな影響を与えた。生物学の分野においては、理論
生物学者の John Maynard Smith らが、ゲーム理論の枠組みを生物進化の解析に応用した「進化ゲーム理論」と
いう理論を提案した [3]。進化ゲーム理論は player の合理性を仮定せず、「より高い利得をもたらす戦略ほど集団
中に広まる」という過程を記述する理論である。より具体的に述べると、
N 人の player がいるとし、それぞれの player は s 個の戦略のうちのいずれかをとっているとする。
個々の player が取る戦略から、個々の player の利得が決まる。player は、周りの player の戦略に応じ
て、特定のルールに従って自らの戦略を変える: 具体的には、より利得の高い他の player の戦略を模
倣する1 。
これは「戦略」を「生物の種類」に置き換えることで、次のようにも言い換えられる2 :
s 種類の生物が合計 N 個体いるとする。集団にどの種類の生物がどれだけいるかによって、それ
ぞれの個体の利得が決まる。個体は特定のルールに従って繁殖により子供を産む: 具体的には、より
利得の高い個体ほど繁殖をしやすい、とする3 。個体が 1 つ新たに生まれるときには別の個体が 1 つ死
に、集団の個体数は一定に保たれるとする。
1 戦略の更新の仕方は他にも様々なものがありえる。詳細は
[2] などを参照
2 この「種類」は生物学的な「種」でなくても良い。具体的には種内での異なるタイプなど。
3 ただしこれは、
「最も『優れた』生物（戦略）が集団中に広まる」ことや、「集団の平均的な利得が絶えず増加する」ことを意味しない。
どの種類がより「有利」かは、他の個体がどの種類であるかに依存する。2.2.2 項の “bistability” や “coexistence” の状況がこれにあたる。
また、囚人の dilemma(2.2.3 項) では、集団の平均的な利得は常に減少することが示せる。
6
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
前者での player(人間など) は、後者での生物が住める「スペース」に対応する。進化ゲーム理論は元々、生物の進
化の過程を説明する理論として提案されたが、同時に、合理性を仮定しないゲーム理論の 1 つの形態にもなって
いる4 。現在、進化ゲーム理論は、生物学や、経済学などの社会科学の分野で広く応用されている ([4] の引用文献
などを参照のこと)。
以下では、進化ゲーム理論の定式化とその応用例を、本論文に必要な範囲で概観する。より詳細な review とし
ては、生物学に近い視点から書かれたものでは [4]、物理学者向けのものでは [2, 5] などを薦める。また、教科書
としては [6, 7] などがある。
2.2
決定論的進化ゲーム
まず、進化ゲームの決定論的定式化を考える。この定式化では、考える集団の個体数が十分大きく、かつ、どの
個体も集団中の他のどの個体とも等しくゲームを行うこと (このような集団を well-mixed な集団と呼ぶ) を仮定す
る。以上の仮定の下、replicator 方程式と呼ばれる微分方程式が、集団中で各戦略をとる個体の数の時間変化を記
述する。この節ではまず、replicator 方程式の定義と、簡単な場合 (戦略が 2 種類のとき) での進化ゲームの性質に
ついて述べる。最後に、具体的な進化ゲームの例をいくつか見る。
2.2.1
replicator 方程式
s 種類の戦略 i ∈ {1, 2, · · · , s} があるとし、利得行列（payoff matrix）を
⎛
a11
⎜
⎜ a21
A=⎜
⎜ ..
⎝ .
as1
a12
...
a1s
a22
..
.
as2
...
..
.
a2s
..
.
...
ass
⎞
⎟
⎟
⎟
⎟
⎠
(2.1)
とする。ここで ai,j は、戦略 i を取る個体が戦略 j を取る個体とゲームを行ったときに得る利得である。また、集
団のサイズを一定とし、集団中で戦略 i を取る個体の割合を xi とする。
上述したように、進化ゲーム理論は「より多くの利得を得た戦略がより多くの子孫を残し、集団中に広まる様
子」を記述するものである。この様子を記述する dynamics には様々なものがあり得るが、replicator dynamics で
は、xi が次の方程式に従うことを要請する [8, 9]:
⎛
⎞
s
s
d
ai,j xj (t) −
ak,l xk (t)xl (t)⎠ i ∈ {1, 2, · · · , s}
xi (t) = xi (t) ⎝
dt
j=1
(2.2)
k,l=1
この方程式を replicator 方程式と呼ぶ。この式の意味を考えよう。右辺の (· · · ) の中は、戦略 i の増殖率を表して
おり、第 1 項は戦略 i の個体が得る平均的な利得、第 2 項は全員の平均利得をそれぞれ表す。つまり、replicator
方程式は、
「各戦略の増殖率がその戦略の利得と集団の平均利得の差に比例すること」を表している5 。なお、ここ
では replicator 方程式は進化ゲーム理論のある種の現象論として導入したが、2.3.2 項で、集団サイズが有限である
ことを陽に取り込んだ確率的定式化から (集団サイズが大きい極限で)replicator 方程式が再現されることを見る。
2.2.2
2×2 進化ゲームの分類
replicator 方程式の具体的な振る舞いを、戦略が 2 種類しかない場合 (s = 2) の場合で見ておく。この場合、系
の dynamics は実質的に 1 変数で書けるが、利得行列の値によって様々な振る舞いが見られる。
4 後者の文脈では、
「文化進化」と解釈されることも多い。
5 比例定数は時間スケールの変換で消せる。
7
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
利得行列を
A=
a
b
c
d
(2.3)
と書き、戦略 i = 1 の個体の割合を x と置くことにする（従って戦略 i = 2 の個体の割合は 1 − x）。このとき、
replicator 方程式は
d
x(t) = x(1 − x) [(a − b − c + d)x + b − d]
dt
に帰着する。α := a − c, β := d − b とおくと、
d
x(t) = x(1 − x) [(α + β)x − β]
dt
(2.4)
(2.5)
となる。x = 0, x = 1 は必ず固定点になることが分かり、また、
x∗ =
β
α+β
(2.6)
も x∗ ∈ [0, 1] なら固定点になることが分かる。
α と β の正負によって、定性的には以下の 5 通りの振る舞いが考えられる（図 2.1）[5] :
dominance 1
α > 0 かつ β < 0 のとき、x = 1 が安定固定点、x = 0 が不安定固定点となる6 。直感的な説明: α, β に対す
るこの条件は、相手が誰であっても戦略 1 を取る方が戦略 2 を取るよりも有利であることを表す。従って、
x = 1(集団中の全個体が戦略 1 をとる状態) が安定固定点になる。
dominance 2
α < 0 かつ β > 0 のとき、x = 0 が安定固定点、x = 1 が不安定固定点となる7 。直感的な説明: α, β に対す
るこの条件は、相手が誰であっても戦略 2 を取る方が戦略 1 を取るよりも有利であることを表す。従って、
x = 0(集団中の全個体が戦略 2 をとる状態) が安定固定点になる。
bistability
α > 0 かつ β > 0 のとき、x = 0, x = 1 は安定固定点、x∗ は不安定固定点となる。このようなゲームを
“coordination game ” と呼ぶ。直感的な説明: α, β に対するこの条件は、ゲームの相手と同じ戦略を取る方
が、異なる戦略を取るよりも有利であることを意味する。つまり、集団内の “多数派” に従うことが有利に
なる状況であり、双安定性が現れる。
coexistence
α < 0 かつ β < 0 のとき、x = 0, x = 1 は不安定固定点、x∗ は安定固定点となる。このようなゲームを
“anti-coordination game” と呼ぶ。直感的な説明: α, β に対する条件は、ゲームの相手と異なる戦略を取る
方が、同じ戦略を取るよりも有利であることを意味する。従って、集団内に複数の戦略が共存する状況が安
定になる。
neutrality
α = 0 かつ β = 0 のとき、∀x ∈ [0, 1] は固定点。
2.2.3
進化ゲームの例
決定論的進化ゲームについての一般論について述べ終えたので、ここでいくつか例を見ておこう。
6 より正確には、
「α
7 より正確には、
「α
> 0 かつ β ≤ 0」または「α ≥ 0 かつ β < 0」のとき
< 0 かつ β ≥ 0」または「α ≤ 0 かつ β > 0」のとき
8
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
dominance 1
dominance 2
0.0
0.2
0.4
bistability
0.6
0.8
1.0
0.0
0.2
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
x
coexistence
0.4
0.4
x
0.6
0.8
1.0
0.6
0.8
1.0
図 2.1: 2 × 2 ゲームの分類。縦軸は replicator 方程式の右辺の値を表す。黒丸は安定固定点、白丸は不安定固定点
を表す。
囚人の dilemma
囚人の dilemma は、協力行動の進化を分析するために広く用いられているゲームである。ゲームを次のように
定義する: 戦略は C(cooperator) か D(defector) のいずれかとする。C は「協力行動を行う」という戦略、D は
「協力行動を行わない」という戦略である。自らが協力行動を行うことによる費用 (cost) を −c、相手に協力行動
を行ってもらうことで得る利益 (benefit) を b とする。b > c > 0 を仮定する。
• C が D とゲームを行ったとき、C は利得 −c を得、D は利得 b を得る（ただ乗りを表す。）
• C どうしがゲームを行ったとき、お互いが協力行動を行うので、どちらも利得 b − c を得る。
• D どうしがゲームを行ったとき、どちらも協力行動を行わないので、どちらも利得 0 を得る。
これらをまとめて利得行列で表現すると（C を i = 1、D を i = 2 として）
b − c −c
b
0
(2.7)
と書ける8 。
明らかに、全員が戦略 C を取る方が全員が戦略 D を取るより利得が大きい（∵ b − c > 0）が、実際には後者が
実現してしまう9 : b > c > 0 を仮定しているので、上記の 2 × 2 ゲームの分類で考えると、α < 0, β > 0 の場合に
相当する。従って、初期条件が「全員が C」でない限り、この系は C の割合が 0 の状態に収束する。つまり、こ
のルールの下では協力行動が自発的に進化する可能性はないことが分かる。
8 囚人の
dilemma の利得行列は、より一般には
R
T
S
P
, T >R>P >S
と定義される。しかし、ここではより単純でよく用いられる (2.7) 式の形を仮定した。
9 これが、このゲームが dilemma と呼ばれる所以である。
9
(2.8)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
Hawk-Dove ゲーム
Hawk-Dove ゲームは、生物の種内における儀式的な闘争10 の進化的起源を説明するために考えられたゲームで
ある [3] 11 。
ある 1 つの生物種を考え、その中にタカ派・ハト派の 2 種類の戦略があるとする12 。この生物の個体どうしが、
資源（例えば食物）を巡って争う場合を考える。資源による利益を V 、争って負けた場合の（例えば負傷による）
利得を −C とおく。C > V > 0 を仮定する。
ゲームを次のように定義する: タカ派は相手が誰であれ全力で争い、
• 相手がハト派なら確実に勝って利得 V を得る。
• 相手がタカ派なら確率
1
2
で勝って利得 V を得て、確率
1
2
で負けて利得 −C を得る。
1
2
で負けて利得 0 を得る。
これに対してハト派は、
• 相手がタカ派なら戦わずして逃げ、利得は 0。
• 相手がハト派なら確率
1
2
で勝って利得 V を得て、確率
これらをまとめて利得行列で表現すると（タカ派を i = 1、ハト派を i = 2 として）
V
C
V
2 − 2
(2.9)
V
2
0
と書ける。
今、C > V > 0 を仮定しているので、上記の 2 × 2 ゲームの分類で考えると、α < 0, β < 0 の場合に相当する。
従って、初期条件が「全てタカ派」
「全てハト派」でない限り、この系はタカ派の割合が x∗ :=
V
C
∈ (0, 1) の状態に
収束する。つまり、集団中の個体が全てタカ派になることはなく、一定の割合のハト派と共存することが分かる。
Rock-Paper-Scissors ゲーム
今までは 2 戦略のゲームを見てきたが、ここで戦略の数が 3 つのゲームの例を見よう。次のような利得行列で
表される 3 戦略のゲームを考える（s > 0 とする）:
⎛
0
⎜
⎝1
−s
−s
0
1
⎞
1
⎟
−s⎠
0
(2.10)
これは、戦略 2 が戦略 1 に勝ち、戦略 3 が戦略 2 に勝ち、戦略 1 が戦略 3 に勝つ、という状況を表している。こ
の構造がまるで “じゃんけん” のようであることから、この種のゲームは Rock-Paper-Scissors ゲームと呼ばれて
いる [5]。
replicator 方程式は、
ẋ1
=
x1 (−sx2 + x3 − π(x))
(2.11)
ẋ2
=
x2 (x1 − sx3 − π(x))
(2.12)
ẋ3
=
x3 (−sx1 + x2 − π(x))
(2.13)
π(x)
:=
(1 − s)(x1 x2 + x2 x3 + x3 x1 )
(2.14)
10 儀式的な闘争:
生物は種間だけでなく種内でも様々な資源（食物、縄張り、繁殖相手など）を巡って争うが、どちらかが死ぬまで争うの
ではなく、比較的「穏和に」争う例が多く知られている [3]。
11 一見すると、相手が死ぬかひどく傷つくまで争いを止めない「過激」な個体が資源を独占し、有利になるとも考えられる。従って、この
ような「過激」な個体が「穏和」な個体を駆逐してしまうのではないかとも考えられる。しかし、自然界には儀式的な闘争が数多くみられる。
この起源を進化論的に説明するために考案されたのが、この Hawk-Dove ゲームである。
12 これをタカとハトの 2 種の生物と混同しないように注意。ここではあくまで 1 種の生物を考えている。
10
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
となる。自明な固定点 (1, 0, 0), (0, 1, 0), (0, 0, 1) に加え、x∗ := ( 13 , 13 , 13 ) もこの系の固定点である。
この系の dynamics を Lyapunov 函数を用いて考える [5]:
H(x) = −x1 x2 x3
(2.15)
1
とおく。H は x1 = 0 or x2 = 0 or x3 = 0 で最大値 0 をとり、x∗ で最小値 − 27
を取る13 。また、replicator 方程
式に従う軌道 x(·) に対して、
d
H(x(t))
dt
=
−x˙1 x2 x3 − x1 x˙2 x3 − x1 x2 x˙3
(2.16)
=
−(1 − s)x1 x2 x3 (1 − 3x1 x2 − 3x2 x3 − 3x3 x1 )
(2.17)
となる。x1 x2 x3 (1 − 3x1 x2 − 3x2 x3 − 3x3 x1 ) は x = x∗ または x1 x2 x3 = 0 では 0、それ以外では常に正である。
従って（図 2.2 を参照）、
• 0 < s < 1 では H は（初期状態で x1 x2 x3 = 0）なら常に減少し、x は x∗ に漸近する。
• s > 1 では H は（初期状態で x = x∗ ）なら常に増加し、x1 = 0 or x2 = 0 or x3 = 0 の軌道に漸近する。ま
た、x∗ は不安定固定点となる。
• s = 1 では H は保存量となり、系は閉軌道を持つ。
図 2.2: Rock-Paper-Scissor ゲームの相図。左: 0 < s < 1 の場合。右: s > 1 の場合。
特に s = 1 の場合には、これまで見てきた 2 × 2 ゲームの例とは違って、系が定常的に振動を続けることに注意。
2.3
確率的進化ゲーム（well-mixed）
前節では、進化ゲームの決定論的定式化に基づく記述を概観してきた。しかし、現実の生物や人間の集団は有
限サイズであり、確率的なゆらぎが存在する。また、replicator 方程式は現象論的に導入されたものであり、進化
ゲームのミクロな dynamics がどのように反映されているのかが定かでない。
この節では、有限集団における進化ゲームを、確率過程を用いて定式化する。簡単のため、まずは well-mixed
な集団（集団のどの個体も他の全ての個体と等しくゲームを行う）を考える14 。まずいくつかの遷移ルールを紹介
し、続いて確率的進化ゲームと決定論的記述の関係を見る。最後に、確率的進化ゲームにおいて戦略の広まり方
を特徴づける方法について述べる。
13 x
1 + x2 + x3 = 1 という制限に注意。
14 well-mixed でない集団については次節で述べる。
11
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
2.3.1
遷移ルール
上述したように、進化ゲームの確率的定式化に際しては、ゲームと繁殖（文化進化の文脈で言えば戦略の更新）
のミクロなルールを指定する必要がある。これらのルールには当然多種多様なものが有り得るが、ここでは代表
的なものを 2 つ紹介する15 。2 つのルールは見かけは異なるが、どちらも「高い利得を得た個体ほど、高い確率で
子孫を残す」という進化ゲームの基本的な原理を反映したものになっていることに注意してほしい。
以下、どちらのルールを考えるときも、戦略は {1, 2, · · · , s} のいずれかであり、ゲームの利得行列は
⎛
a11
⎜
⎜ a21
A=⎜
⎜ ..
⎝ .
as1
a12
...
a1s
a22
..
.
as2
...
..
.
a2s
..
.
...
ass
⎞
⎟
⎟
⎟
⎟
⎠
(2.18)
と表されるとする。ここで ai,j は、戦略 i を取る個体が戦略 j を取る個体とゲームを行ったときに得る利得であ
る。個体の総数（集団サイズとも呼ぶ）は N とし、戦略 i の個体の数を ni と書く。戦略 i の個体の得る利得 π i は
s
1 π (n) =
aij nj
N j=1
i
(n := (n1 , n2 , · · · , ns ))
(2.19)
と書けるとする16 。
Moran process
Moran process
17
では、次のような過程を考える（ここでは生物学的な解釈を用いて説明する）。
1. 戦略 i の個体に適応度（fitness）f i (n) := f (π i (n)) を割り当てる。f は正の値を取る単調増加函数であれば
よいが、具体的な形としては
• f (π) = 1 − w + wπ （ w ∈ [0, 1]、利得行列は π が非負であるようなものに限る（具体的には、a, b, c, d ≥
0）。）
• f (π) = ewπ (w ∈ R)
などが多く用いられる18 。ここで、w は利得がどれだけ繁殖の成功率（ここでは適応度）に影響するかを表
すパラメーターで、「自然選択の強さ（intensity of selection）」と呼ばれる。w の大きさは、自然選択が確
率的浮動に比べてどれほど強いかを意味する。
2. rate λ で、1 個体をランダムに選ぶ。
3. 次に、1 個体を適応度に比例した確率で選ぶ（2. で選ばれた個体が再び選ばれても良い）
：戦略 i を持つ個体
ni f i (n)
が選ばれる確率は s nj f j (n) となる。
j=1
4. 3. で選ばれた個体は自分と同じ戦略を持つ子孫を産む。2. で選ばれた個体は死に、新しく生まれた個体がそ
の空きを埋める。
15 より網羅的な定義が知りたい読者は
[2] の 4 章を参照のこと。
や考える状況にもよるが、右辺の 1/N はなくてもよい。また、ここの convention では、個体が「自分とゲームを行うこと」
は除外していない。
17 Moran process は元来は集団遺伝学の文脈で用いられてきたモデル [10] で、比較的最近になって進化ゲームに導入された [11, 12]。
18 w の範囲と利得行列につく制限は、f が非負であることを保証するためのものである。
16 convention
12
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
特に、簡単な場合として s = 2 の場合の master 方程式を示しておく: 時刻 t で戦略 1 の個体の数が n である確
率を P (n, t) とおくと、master 方程式は
∂t P (n, t)
=
W + (n)
:=
W − (n)
:=
W + (n − 1)P (n − 1, t) + W − (n + 1)P (n + 1, t) − [W + (n) + W − (n)]P (n, t)
N −n
nf 1 (n)
λ 1
nf (n) + (N − n)f 2 (n) N
n
(N − n)f 2 (n)
λ 1
2
nf (n) + (N − n)f (n) N
(2.20)
(2.21)
(2.22)
と書かれる。一般の場合も同様である。
Pairwise comparison processes
Moran process では、集団中の全ての個体それぞれが持つ利得についての情報が必要であったが、これは現実の
系には適さない場合もある。それに対して、Pairwise comparison process では、次のような過程を考える [13]：
1. 集団から、rate λ で focal 個体と role 個体をランダムに選ぶ。
2. focal 個体の利得を π f 、role 個体の利得を π r とし、p : R → [0, 1] とする。このとき、確率 p（π r − π f ）で
focal 個体は死に、そこに role 個体が自らと同じ戦略を持つ子を産む。p の函数形としては、
• p(Δπ) = 1+wΔπ
(w ∈ [0, 1]、利得行列は Δπ の絶対値が 1 以下であるようなものに限る (具体的には、
2
a − c, b − d ∈ [−1, 1]))
−1
(w ∈ R)
• p(Δπ) = 1 + e−wΔπ
などが用いられる。w は Moran process のときと同様に自然選択の強さと呼ばれ、利得が繁殖の成功率にど
の程度影響を及ぼすかを表す。
特に、簡単な場合として s = 2 の場合の master 方程式を示しておく：時刻 t で戦略 1 の個体の数が n である確
率を P (n, t) とおくと、master 方程式は
∂t P (n, t)
=
W + (n)
:=
W − (n)
:=
W + (n − 1)P (n − 1, t) + W − (n + 1)P (n + 1, t) − [W + (n) + W − (n)]P (n, t)
n N −n
λ
p(π 1 (n) − π 2 (n))
N N
n N −n
λ
p(π 2 (n) − π 1 (n))
N N
(2.23)
(2.24)
(2.25)
と書かれる。一般の場合も同様である。
2.3.2
確率的記述と決定論的記述の関係
ここまで進化ゲームの決定論的記述と確率的記述について述べてきたが、ここで両者の関係について述べる。進
化ゲームの確率的記述において、集団サイズ N が大きい極限を考えると確率的ゆらぎを無視することができ、期
待値の dynamics が replicator 方程式（またはそれに類するもの）に従うことを示す。
簡単のため戦略が 2 つの場合 (s = 2) を考え、戦略 i = 1 の個体の数を n とおく。利得行列は
a b
c
d
(2.26)
と書けるとする。このとき、master 方程式
∂t P (n, t) = W + (n − 1)P (n − 1, t) + W − (n + 1)P (n + 1, t) − [W + (n) + W − (n)]P (n, t)
13
(2.27)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
を集団サイズ N 1 で考える。この系に対して、van Kampen 展開による linear noise 近似（[31] の Chap. X）
と呼ばれる近似手法を適用し、n の期待値とその周りの揺らぎを評価する。
まず、
n(t) = N φ(t) +
√
Nξ
(2.28)
と書けると仮定する。φ(t) は時刻 t での n の期待値を表し、ξ はその周りのゆらぎを表す。また、
W± (x) :=
1 ±
W (N x) (x ∈ [0, 1])
N
(2.29)
と定義する。右辺は N に依らないと仮定する (上述の例では、λ = N とおけば成り立つ。)。更に、ξ の確率密度
関数を Π とおく：
√
1
Π(ξ, t) := √ P (N φ(t) + N ξ, t)
N
(2.30)
これらを元の master 方程式に代入して N について展開すると (計算の詳細は付録 A.1.2 を参照のこと)、期待
値 φ(·) の従う方程式として、
φ̇(t) = W+ (φ(t)) − W− (φ(t))
(2.31)
が得られる19 。
W の具体的な形を代入して、(2.31) 式がどのような形を取るかを見てみよう [13, 14]20 。
Moran process
Moran process における遷移率の定義 (2.21),(2.22) より、λ = N とすると、
W+ (φ)
=
W− (φ)
=
fi (φ)
:=
φf1 (φ)
(1 − φ)
φf1 (φ) + (1 − φ)f2 (φ)
(1 − φ)f2 (φ)
φ
φf1 (φ) + (1 − φ)f2 (φ)
f (πi (φ))
π1 (φ)
:=
aφ + b(1 − φ)
(2.35)
π2 (φ)
:=
cφ + d(1 − φ)
(2.36)
(2.32)
(2.33)
(2.34)
となる。ただし、f (π) = 1 − w + wπ または f (π) = ewπ である。
これを (2.31) 式に代入すると、
φ̇ =
φ(1 − φ) [f1 (φ) − f2 (φ)]
φf1 (φ) + (1 − φ)f2 (φ)
(2.37)
が得られる。Moran process についての仮定より、分母は非負である。このことと、π1 (φ) − π2 (φ) = (a − b − c +
d)φ + b − d であることに注意すると、(2.37) 式は定性的には replicator 方程式と同じ振る舞いをすることが分かる。
更に、w が小さい極限を考えると、
φ̇ = wφ(1 − φ) [(a − b − c + d)φ + b − d] + O(w2 )
(2.38)
となり、replicator 方程式で近似できる。
19 Π についての偏微分方程式も得られ (付録 A.1.2)、こちらは揺らぎ ξ の分布が Gaussian で与えられることを意味する（[31] Chap. X）
。
期待値の dynamics が (2.31) 式で書けるためには、本来なら Π について揺らぎ ξ が大きくならないことを示す必要がある。しかしここでは
簡単のため、そのような詳細には立ち入らない。
20 導出法はここで用いたものとは異なるが、得られる結果は同じである。
14
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
Pairwise comparison process
Pairwise comparison process における遷移率の定義 (2.24),(2.25) より、λ = N とすると、
となる。ただし、 p(Δπ) =
1+wΔπ
2
W+ (φ)
=
φ(1 − φ)p(π1 (φ) − π2 (φ))
(2.39)
W− (φ)
=
φ(1 − φ)p(π2 (φ) − π1 (φ))
(2.40)
π1 (φ) :=
aφ + b(1 − φ)
(2.41)
π2 (φ) :=
cφ + d(1 − φ)
(2.42)
−1
または p(Δπ) = 1 + e−wΔπ
である。
これを (2.31) 式に代入すると、
φ̇ = φ(1 − φ) [p(π1 (φ) − π2 (φ)) − p(π2 (φ) − π1 (φ))]
(2.43)
が得られる。より具体的に書くと、
1 + wΔπ
p(Δπ) =
2
−wΔπ −1
p(Δπ) = 1 + e
のとき： のとき： φ̇ = wφ(1 − φ) [π1 (φ) − π2 (φ)]
w
[π1 (φ) − π2 (φ)]
φ̇ = φ(1 − φ) tanh
2
(2.44)
(2.45)
となる。π1 (φ) − π2 (φ) = (a − b − c + d)φ + b − d であることに注意すると、(2.44) は（時間スケールを除いて）
replicator 方程式そのものであり、(2.45) 式も定性的には replicator 方程式と同じ振る舞いをすることが分かる。
更に、w が小さい極限を考えると (2.45) 式は
φ̇ =
w
φ(1 − φ) [(a − b − c + d)φ + b − d] + O(w3 )
2
(2.46)
となり、replicator 方程式で近似できる。
2.3.3
突然変異と定常分布
ここまで考えてきた進化ゲームの dynamics では、一度集団から消えた戦略は二度と現れず、また、元から集団
中になかった戦略が新しく生まれることもない21 。これは、どの個体の子供も、その親と同じ戦略を持つことから
来るものであった。
しかし、現実の系では、子供の戦略と親の戦略が（きわめて低い確率ではあるが）異なる場合がありうる。こ
れは生物学的には突然変異、文化進化の文脈では学習の失敗・または気まぐれと捉えられ、自然なものである。
突然変異のモデル化としては、次のようなものが主に用いられる [15, 16]:
戦略 i の個体の子は、確率 μij で戦略 j を取る。ただし、μij は μij > 0 かつ
s
j=1
μij = 1 を満
たす。
このとき、対応する master 方程式には定常解が一意に存在し、長時間極限を考えると任意の初期条件に対して
この定常解に収束することが言える ([31]V.3)。この定常解（定常確率分布）を用いて、「どの戦略がより有利か
（＝より自然選択で選ばれやすいか）」を特徴づけることができる: Pst (·) を定常分布とするとき、
ni :=
ni Pst (n)
n
とおく。ni が大きいほど、その戦略 i が自然選択で選ばれやすいと考えることができる [15, 16]。
21 このような
dynamics は “non-innovative” と呼ばれる [2]。
15
(2.47)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
2 × 2 進化ゲームでの定常分布での戦略分布
定常分布を用いた解析の例として、Antal らによる研究 [15] の結果を簡潔に紹介する:
利得行列は
a b
c d
(2.48)
と書けるとする。集団の全個体数を N とし、戦略 i = 1 の個体が n いるときの各個体の利得を
N −n
n−1
a+
b
N −1
N −1
N −n−1
n
c+
d
N −1
N −1
π 1 (n) =
π 2 (n) =
(2.49)
(2.50)
とする22 。遷移ルールは Moran process と Pairwise comparison process を含む広いクラスのものを考える。
この系に対して自然選択の強さ w が小さい極限（弱選択極限）を考えると、
戦略 i = 1 の個体が戦略 i = 2 の個体より定常分布で考えて多くなる
2
2
+b>c+d 1−
a 1−
N
N
⇐⇒
(2.51)
と表される23 。
これを囚人の dilemma
b−c
−c
b
0
(2.55)
に適用すると、
⇐⇒
戦略 C が戦略 D より（定常分布について）多くなる
2
2
c
− b> 2−
N
N
(2.56)
となる。しかし、b, c > 0 より、これは成り立たない24 。従って、囚人の dilemma(2.55) では、
（定常分布で考えれ
ば）常に D の数が C の数を上回ることが分かる。
2.4
確率的進化ゲーム（structured population）
先ほどは well-mixed な集団での進化ゲームの確率的定式化について述べた。しかし、現実の生物や社会を考え
ると、
「集団の各個体が他の全個体と等しい確率で相互作用する」という仮定は、場合によっては非現実的である。
そこで、この節では well-mixed でない、何らかの構造を持った集団（これを structured population と呼ぶ）での
22 個体が自分とゲームを行うことを除外していることに相当。
23 個体が自分とゲームを行うことを除外しなければ、利得は
π 1 (n)
=
π 2 (n)
=
n
N −n
a+
b
N
N
n
N −n
c+
d
N
N
(2.52)
(2.53)
で、結果は
戦略 i = 1 の個体が戦略 i = 2 の個体より（定常分布で考えて）多くなる
⇐⇒
a+b>c+d
(2.54)
となる。
24 個体が自分とゲームを行うことを除外しなければ、
戦略 C が戦略 D より（定常分布について）多くなる
⇐⇒
0 > 2c
(2.57)
となり、これは c > 0 である限り成り立たない。
16
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
進化ゲームの定式化について論じる。structured population の例としては、個体が空間上の各点やグラフ (ネット
ワーク) の各頂点に分布しているような状況や、個体が複数のグループのうちいずれかに属しているような状況を
考えることができる。
2.4.1
遷移ルール
まず、structured population での確率的進化ゲームの遷移ルールの定式化を概観する。
戦略は i ∈ {1, 2, · · · , s} のいずれかであり、ゲームの利得行列は
⎛
a11
⎜
⎜ a21
A=⎜
⎜ ..
⎝ .
as1
a12
...
a1s
a22
..
.
as2
...
..
.
a2s
..
.
...
ass
⎞
⎟
⎟
⎟
⎟
⎠
(2.58)
と書けるとする。ここで ai,j は、戦略 i を取る個体が戦略 j を取る個体とゲームを行ったときに得る利得である。
個体は空間（または何らかの構造を表す集合）Λ の各点に分布しており、各点 x ∈ Λ にいる個体数 Nx ∈ N と、
x の近傍 Nx ⊂ Λ が定まっているとする25 。近傍は、位置 x にいる個体が（ゲームや繁殖の際に）相互作用できる
範囲を表す26 。
位置 x にいる戦略 i の個体の数を nx,i とする。位置 x にいる戦略 i の個体の得る利得は
π x,i (n) = C
s
aij ny,j
(2.59)
y∈Nx j=1
とする27 。ただし、C は適当な定数で、通常は 1 にとることが多い。また、Nx が全ての x ∈ Λ について等しけれ
ば、C = 1/Nx と取ることもある。
遷移のルールとしては様々なものが有り得るが、主に次の 2 つを考える28 :
Moran process
1. 各個体に対して、適応度 f x,i (n) = f (π x,i (n)) を割り当てる（f の関数形は f (π) = 1 − w + wπ や f (π) = ewπ
など。）
2. rate λ で集団全体から 1 個体が random に選ばれる。この個体のいる位置を x とする。
3. x の近傍 Nx から 1 個体を適応度に比例した確率で選ぶ（このとき、前のステップで選ばれた個体が再び選
ばれても良い）。
4. 2. で選ばれた個体が死に、3. で選ばれた個体の子がその空きを埋める。
Pairwise comparison process
1. rate λ で集団全体から 1 個体がランダムに選ばれる。この個体を focal 個体と呼ぶ。focal 個体のいる位置を
x とする。
25 「近傍」N は、空間上に分布した集団なら通常の意味での近くの点たち、グラフ上の集団ならエッジ (辺) で結ばれた点たち、グループ
x
に分かれた集団なら同一のグループ、などを考えることが多い。
26 一般には、ゲームの相互作用を表す近傍と、繁殖の際の相互作用を表す近傍が同一である必要はない。しかし、ここでは簡単のため両者
が一致する場合のみを考える。一致しない場合については
[17] を参照。
27 右辺の
y∈Nx が、x にいる個体がその周り Nx にいる個体としかゲームをしないことを表す。
28 well-mixed な場合は次のような特殊なケースとしてこのルールに含まれる: Λ を 1 点から成る集合 {x} にとり、N = Λ とすればよい。
x
または、Λ を任意として、Nx = Λ(∀x ∈ Λ) としてもよい。
17
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
2. x の近傍 Nx から 1 個体をランダムに選ぶ。この個体を role 個体と呼ぶ。
3. focal 個体の利得を π f 、role 個体の利得を π r とし、p : R → [0, 1] とする。focal 個体は、確率 p(π r − π f ) で
死に、そこに role 個体は自らの子を産む。
上記の説明では突然変異は考慮していないが、これらの過程に加えて突然変異を入れることができる。このと
き、well-mixed な場合と同様に定常分布を考えることができる。
2.5
structured population における進化ゲーム
最後に、これまで述べてきた進化ゲーム理論の枠組みを踏まえ、集団の構造が進化ゲームの dynamics にどのよ
うな影響を与えるかについて、いくつかの既知の結果を紹介する。
2.5.1
正方格子上での囚人の dilemma
Szabo, Toke[18] は 2 次元正方格子上での囚人の dilemma を考えた。利得行列は
1 0
b 0
(2.60)
2
で与えられ（i = 1 が C、i = 2 が D、b > 1 とする）、Λ = {1, 2, · · · , L} , Nx = {x の第一近接点 } で、突然変
異はないとする。
Szabo らは、Monte Carlo simulation を用いて L がきわめて大きいときの定常状態29 での戦略 C の密度を求め
た。結果、b > 1 であっても C が絶滅せず D と共存することを示した ([18] の図 1-図 3 を参照)。well-mixed モデ
ルでは N → ∞ で C が絶滅してしまうのとは対照的である。このことは直感的には次のように説明できる: C が
C ばかりが集まった colony を形成する。colony の中では C ばかりなので、そこではどの C の個体も高い利得を
得ることができる。一方、D と相互作用する個体は少ないので、全体としては D による「搾取」を逃れることが
できる。
2.5.2
格子上での snowdrift ゲーム
上述の囚人の dilemma では、空間構造は協力行動の進化が促されていた。しかし、このことは必ずしも一般の
ゲームについて言えるわけではない。そのような例をここで紹介する。
Hauert,Doebeli は、様々な 2 次元格子上で、snowdrift ゲーム30 を考えた [19]。利得行列は
b − 2c b − c
b
0
(2.61)
と表されるとする (b > c > 0)。戦略 1 を C(cooperator)、戦略 2 を D(defector) と呼ぶ。突然変異はないとする。
このゲームを well-mixed な集団で考えると、2.2.2 項の分類より 2 つの戦略は共存し、戦略 1(C) の割合は
c
1 − r (r := 2b−c
) となることが分かる。
一方、Hauert, Doebeli らは格子上での snowdrift ゲームを考え、Monte Carlo simulation によって定常状態で
の C の割合を求めた。結果、パラメーター r の範囲によっては、C の割合が well-mixed な集団での割合 1 − r よ
りも小さくなることが示された ([19] の図 1, 図 2 参照)。
29 集団サイズが有限なら定常分布は「全て
C」の状態か「全て D」の状態にしか値を取らない。そのため、無限系を考える必要がある。
ゲームという名前は、次のような設定に由来する: 2 人のドライバーが吹雪の中で吹き溜まり (snowdrift) に捕まった。ドライ
バーたちには、
「C. 車を降りて除雪する。 D. 除雪せずに車の中に留まる。」という 2 つの選択肢がある。1 人で吹き溜まりを取り除くのにか
かるコスト (cost) を −c(1 人だけでも除雪はできるとする)、目的にたどり着ける利益 (benefit) を b をとする。1 人だけが除雪すれば、除雪
した人は利得 b − c を、除雪しなかった人は利得 b を得る。2 人が除雪すれば、2 人とも利得 b − 2c を得る。2 人とも除雪しなければ、2 人と
も利得 0 を得る。
30 snowdrift
18
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
このことは、直感的には次のように説明できる: 囚人の dilemma では C が colony を形成することよって C が
D による「搾取」を免れていた。しかし、snowdrift ゲームでは「相手と異なる戦略を取ると得をする」という状
況なので、colony を形成しにくい。
以上のように、snowdrift ゲームは空間構造が協力行動の進化を必ずしも促進するとは限らないことを示す興味
深い例である。
2.5.3
Rock-Paper-Scissors ゲーム
Kerr ら [20] は、3 種の生物からなる Rock-Paper-Scissors 型のゲームを well-mixed な集団と 2 次元正方格子上
で考え31 、Monte Carlo simulation と、大腸菌を用いた実験を行った。
実験にあたっては、3 種の遺伝子組換え大腸菌を用いた。これらをそれぞれ C,S,R と呼ぶ。
• C (colicinogenic): colicin という毒素を生産する。自身は colicin に対する耐性を持つ。
• S (sensitive): colicin に対する耐性を持たない。
• R (resistant): colicin に対する耐性を持つ。
C は colicin を生産して S を殺す。S は R より成長速度が速い（∵colicin 耐性は成長速度を遅くする）。R は C よ
り成長速度が速い（∵colicin を生産するコストがかからない）。この 3 つの関係から、C, S, R は「C が S に勝ち、
S が R に勝ち、R が C に勝つ」という、Rock-Paper-Scissors 型の相互作用をすることが分かる。well-mixed な状
況を撹拌されたフラスコで、空間構造がある状況をシャーレでそれぞれ再現した。
simulation にあたっては次のルールを採用した (以下は 2 次元格子上での場合。well-mixed の場合は各点の近傍
を格子全体に取ればよい。)：
• 2 次元正方格子を考え、その各点には 3 種の個体のうちいずれか 1 種の個体がいる、または空地であると
する。
• 一定の rate でランダムに 1 つの格子点を選ぶ。
• 選ばれた格子点が空地であるなら、確率 fi (i ∈ {C, S, R}) で種 i がこの点に生まれる。ただし、fi は選ば
れた格子点の近傍（8 個の格子点）にいる個体に占める種 i の個体の割合とする。
• 選ばれた格子点に種 i の個体がいれば、確率 Δi (i ∈ {C, S, R}) でこの個体は死んで、選ばれた格子点は空
地になる。ただし ΔC , ΔR は定数であるのに対して、ΔS は fC の函数とする: ΔS (fC ) = ΔS,0 + τ fC とす
る (τ は適当な定数)。
simulation、大腸菌を用いた実験いずれにおいても、well-mixed な場合ではすぐに C,S,R のうち 2 種が絶滅し
いずれか 1 種のみが生存することになったが、空間構造があると C, S, R の 3 種が共存するようになった ([20] の
図 1, 図 2 参照)。これらの結果から、空間構造は種の共存を促進することが示唆される。
2.5.4
一般の structured population における進化ゲーム
Tarnita ら [21] は 2.3.3 節の定式化を用いて、集団の構造が定常状態での各戦略の個体数の期待値に対して及ぼ
す影響を検討した:
2 × 2 ゲームを考え、利得行列は
a
b
c
d
31 ゲームのルールは本論文の形とは大きく異なっているが、定性的には似たような状況を考えている
19
(2.62)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と表されるとする。繁殖には突然変異が入るとする。考える集団は一般の structured population とする。遷移ルー
ルは Moran や Pairwise comparison を含む広いクラスのものを考える。
このとき、定常分布についての期待値で考えて戦略 1 が戦略 2 より多いことの必要十分条件は、自然選択が弱
い極限で
σa + b > c + σd
(2.63)
と表される。ただし、σ は遷移ルールの詳細と考える集団の構造に依存したパラメーターで、利得行列の要素には
依らない。利得行列の対角要素がパラメーター σ の分だけ重みづけされることに注意してほしい。σ の値はいく
つかの例について具体的に求めることができる。例えば、ある種の Moran process を k 次の正則グラフ32 上の各
頂点に 1 個体がいるというモデルで考えると、
σ=
(k + 1)N − 4k
(k − 1)N
(2.64)
となる [21]。この σ は N > 2k で σ > 1 となることに注意。
この結果を囚人の dilemma(2.7) 式に適用すると、C が D より（定常分布についての期待値で考えて）多くなる
必要十分条件は、自然選択が弱い極限で
(σ − 1)b > (σ + 1)c
(2.65)
となる。これは σ > 1 なら、b が c より十分大きければ満たされる。つまり、σ > 1 なら、well-mixed な場合に比
べて協力行動が進化しやすくなる (i.e. C の個体数が多くなる) ことが分かる（c.f. 2.3.3 項）。
2.5.5
まとめと課題
以上のように、集団の構造が進化ゲームに与える影響については様々な研究がなされている [22, 23, 24, 25]。そ
の多くでは協力行動の進化が集団の構造によって促進・阻害されたり、種の共存が促されたりと、集団の構造が戦
略の広まり方に大きな影響を与えることが示されている。
しかし、既存の研究の多くは「複数の戦略が共存しやすくなるか/しにくくなるか」や「定常分布での期待値
が集団の構造の有無によって増減するか」など、進化ゲームの静的な側面に着目している。これに対して、進化
ゲームの動的な側面に関する研究はあまり進んでいない。well-mixed な場合にはいくつかの研究が知られている
[26, 27, 28, 29, 30] が、集団の構造が dynamics にどのような影響を与えるかに関しては、まだ理論的研究は始まっ
たばかりである。
32 全ての頂点が
k 本のエッジ（線分）を持つグラフのこと。
20
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
第 3 章 確率過程と経路積分
この章では確率過程の解析に用いる手法、特に経路積分を用いた手法について論じる。
多くの確率過程はその master 方程式を厳密に解くことは不可能であり、従って具体的な問題を解くにあたって
は様々な近似が必要となる。広く利用される近似手法として van Kampen の Ω 展開による linear noise 近似 ([31]
の chap. X) が知られている。しかし、この手法は期待値から大きく外れた振る舞い (いわゆる rare event) の確率
の評価には適さない (そのような例を 3.2.2 項で見る。)。対して、rare event の確率の評価に適した手法として、
経路積分の手法に基づく半古典近似と呼ばれる手法を用いることができる ([32] の chap. 10)。
この章では、確率過程を経路積分で記述する手法を述べた後、半古典近似の手法について概観する。
3.1
確率過程の経路積分による記述
連続時間・離散状態の Markov 過程を経路積分によって記述する手法を解説する1 。一般の整数値 n を取る Markov
過程を考え、その 1 つの sample path を n(·) と書く。n(·) から決まる物理量を O とする。本節では、次のことを
示したい: O の期待値 O は経路積分を用いて
O
=
S[n(·), p(·)]
:=
DnDpO[n(·)]e−S[n(·),p(·)]
(3.1)
dt[p(t)∂t n(t) − H(n(t), p(t))]
(3.2)
のように書ける。ここで H は確率過程ごとに決まる Hamiltonian であり、具体的な決め方は後で述べる。本節で
は、与えられた確率過程のルールからどのように Hamiltonian が定まるかを、具体的な例を通して見ていく2 。
3.1.1
簡単な例:崩壊過程
この項では、崩壊過程を例にとる。粒子 A が一定の rate λ で崩壊するという過程を、時刻 t ∈ [0, T ] で考える3 :
A → ∅ (with rate λ per particle)
(3.3)
時刻 t での A 粒子の数を n(t) とし、n(0) は given とする。
step 1
まず、n(·) に依存する物理量 O[·] の期待値を経路についての和の形で書き表す。確率過程の定義より、Δt が十
分小さければ
1 本節の内容は
⎧
⎨−1 ( with probability λn(t)Δt)
n(t + Δt) − n(t) =
⎩0
( with probability 1 − λn(t)Δt)
[33] に基づく。
2 導出は長くなるので、最終的な結果だけを知りたい読者は、3.1.4
3 この確率過程の
項のまとめを認めて先に進んでも構わない。
master 方程式は、時刻 t で粒子数が n である確率を P (n, t) とすると、
∂t P (n, t) = λ(n + 1)P (n + 1, t) − λnP (n, t)
となる。
21
(3.4)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
が成り立つ。
時間 [0, T ] を M ( 1) 個に分割し、時間を離散化して考える:
Δt := T /M, tj := jΔt (j ∈ {0, 1, · · · , M }), nj := n(tj )
(3.5)
すると、n(·) に依存した物理量 O の期待値 O は
∞
O =
n1 ,n2 ,··· ,nM =0
∞
O[n0 , n1 , · · · , nM ] Prob[n0 → n1 → n2 → · · · nM ]
(3.6)
O[n0 , n1 , · · · , nM ] · Prob[n0 → n1 ] · Prob[n1 → n2 ] · · · Prob[nM −1 → nM ]
(3.7)
n1 ,n2 ,··· ,nM =0
と書ける。ただし、Prob[n0 → n1 → n2 → · · · nM ] は n0 → n1 → n2 → · · · nM という経路が実現する確率を表
し、Prob[nj → nj+1 ] は現在の状態が nj としたときに次の時間に nj+1 に移る条件付き確率を表す。2 行目に移る
ときには、確率過程の Markov 性を用いた。
step 2
ここで、新しい変数を
Jj := nj+1 − nj
と定義する。崩壊過程の例では、nj が与えられたもとで、Jj は
⎧
⎨−1 ( with probability λn Δt)
j
Jj =
⎩0
( with probability 1 − λnj Δt)
(3.8)
(3.9)
と決まる。そこで、これ以降 n(·): given に対して、J を (3.9) 式で n(·) から定まる新しい確率過程と見なすこ
とにする4 。
Pj (Δn|nj ) = (時刻 j で n = nj であったときに Jj = Δn である確率)
とおくと、この例では
⎧
⎪
⎪
λn Δt
⎪
⎨ j
Pj (Δn|nj ) =
(3.10)
(Δn = −1 の場合)
1 − λnj Δt (Δn = 0 の場合)
⎪
⎪
⎪
⎩0
(それ以外の場合)
(3.11)
となる。Pj (Δn|nj ) の定義より
Prob[nj → nj+1 ] = Pj (nj+1 − nj |nj )
(3.12)
が成り立つ。更に、Jj についての期待値を
∞
f (Jj )|nj =
f (Δnj )Pj (Δnj |nj )
(3.13)
Δnj =−∞
とおくと、
Pj (nj+1 − nj |nj ) = δnj+1 −nj ,Jj |nj
(3.14)
となる。従って、
O ∞
O[n0 , n1 , · · · , nM ]δn1 −n0 ,J0 |n0 δn2 −n1 ,J1 |n1 · · · δnM −nM −1 ,JM −1 |nM −1
(3.15)
n1 ,n2 ,··· ,nM =0
=
M
i=1
dni O[n0 , n1 , · · · , nM ]
M
−1
δ(nj+1 − nj − Jj )|nj
(3.16)
j=0
4 これは、(3.7) 式において「確率性」を J に押し込めることに対応している。(3.7) 式の和の各項においては、n(·) は確率過程ではなく、
1 つに定まっていることに注意。
22
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と表すことができる。
⎛
(∵)
⎜
∞
⎜
⎜
O[n0 , n1 , · · · , nM ]δn1 −n0 ,J0 |n0 δn2 −n1 ,J1 |n1 · · · δnM −nM −1 ,JM −1 |nM −1
⎜ O ⎜
n1 ,n2 ,··· ,nM =0
⎜
⎜
⎜
· の定義より、
⎜
⎜
∞
⎜
⎜
=
O[n0 , n1 , · · · , nM ]
⎜
⎜
n1 ,n2 ,··· ,nM =0
⎜
⎜
⎜
δn1 −n0 ,Δn0 P0 (Δn0 |n0 )δn2 −n1 ,Δn1 P1 (Δn1 |n1 )
×
⎜
⎜
Δn0 ,Δn1 ,··· ,ΔnM −1
⎜
⎜
⎜
· · · δnM −nM −1 ,ΔnM −1 PM −1 (ΔnM −1 |nM −1 )
⎜
⎜
∞
⎜
⎜
=
O[n0 , n1 , · · · , nM ]
⎜
⎜
n
,···
,n
=0
Δn
,···
,Δn
M
0
M −1 1
⎜
⎜
⎜
× δn1 −n0 ,Δn0 P0 (Δn0 |n0 )δn2 −n1 ,Δn1 P1 (Δn1 |n1 ) · · · δnM −nM −1 ,ΔnM −1 PM −1 (ΔnM −1 |nM −1 )
⎜
⎜
⎜
ここから nj についての和を nj についての積分に変換していく:
⎜
⎜
∞
∞
⎜
⎜
=
dnM O[n0 , n1 , · · · , nM ]
⎜
⎜
−∞
n
,···
,n
=0
Δn
,···
,Δn
M −1
0
M −1 1
⎜
⎜
⎜
× δn1 −n0 ,Δn0 P0 (Δn0 |n0 )δn2 −n1 ,Δn1 P1 (Δn1 |n1 ) · · · δnM −1 −nM −2 ,ΔnM −2 PM −2 (ΔnM −2 |nM −2 )
⎜
⎜
⎜
× δ(nM − nM −1 − ΔnM −1 )PM −1 (ΔnM −1 |nM −1 )
⎜
⎜
∞
∞
∞
⎜
⎜
=
dn
dnM O[n0 , n1 , · · · , nM ]
⎜
M −1
⎜
−∞
−∞
n
,···
,n
=0
Δn
,···
,Δn
M −2
0
M −1 1
⎜
⎜
⎜
× δn1 −n0 ,Δn0 P0 (Δn0 |n0 )δn2 −n1 ,Δn1 P1 (Δn1 |n1 ) · · · δ(nM −1 − nM −2 − ΔnM −2 )PM −2 (ΔnM −2 |nM −2 )
⎜
⎜
⎜
× δ(nM − nM −1 − ΔnM −1 )PM −1 (ΔnM −1 |nM −1 )
⎜
⎜
⎜
= ···
⎜
⎜
∞
∞
⎜
⎜
=
dn
·
·
·
dnM O[n0 , n1 , · · · , nM ]
1
⎜
−∞
−∞
⎜
Δn
,···
,Δn
0
M −1
⎜
⎜
⎜
× δ(n1 − n0 − Δn0 )P0 (Δn0 |n0 ) · · · δ(nM − nM −1 − ΔnM −1 )PM −1 (ΔnM −1 |nM −1 )
⎜
⎜
⎜
再び · の定義より、
⎜
∞
∞
⎜
⎜
⎜
=
dn1 · · ·
dnM O[n0 , n1 , · · · , nM ]
⎜
−∞
−∞
⎜
⎜
⎜
× δ(n1 − n0 − J0 )|n0 · · · δ(nM − nM −1 − JM −1 )|nM −1
⎜
⎜
M
M
−1
⎜
⎝
=
dni O[n0 , n1 , · · · , nM ]
δ(nj+1 − nj − Jj )|nj
i=1
j=0
step 3
δ 函数の Fourier 積分表示
δ(nj+1 − nj − Jj ) =
i∞
−i∞
dpj+1 pj+1 (Jj −nj+1 +nj )
(j ∈ {0, 1, · · · , M − 1})
e
2πi
23
(3.17)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
を用いると、
O
=
M
M
−1
dnk dpk
O[n0 , n1 , · · · , nM ]
exp[pj+1 (Jj − nj+1 + nj )] |nj
2πi
j=0
k=1
M −1
M −1
M
dnk dpk
pi+1 (ni+1 − ni )
exp(pj+1 Jj ) |nj
O[n0 , n1 , · · · , nM ] exp −
2πi
i=0
j=0
(3.18)
(3.19)
k=1
ここで、
= (e−pj+1 λnj Δt) + (1 − λnj Δt)
exp(pj+1 Jj )|nj
(3.20)
−pj+1
− 1)λnj Δt
= 1 + (e
−pj+1
− 1)nj
exp λΔt(e
と書けるので、
O (3.21)
(3.22)
⎧
⎫
M
−1
⎨M
⎬
dnk dpk
O[n0 , n1 , · · · , nM ] exp
−pj+1 (nj+1 − nj ) + λΔt(e−pj+1 − 1)nj
(3.23)
⎩
⎭
2πi
j=0
k=1
M → ∞ とすれば、（3.23）式の は = と見なしてよいと期待される:
O
=
M
dnk dpk
O[n0 , n1 , · · · , nM ] exp(−S[n, p])
M →∞
2πi
lim
(3.24)
k=1
S[n, p]
=
M
−1
[pj+1 (nj+1 − nj ) − Δt · H(nj , pj+1 )]
(3.25)
j=0
H(n, p)
=
λ(e−p − 1)n
これを
(3.26)
O
=
S[n(·), p(·)]
=
DnDp O[n(·)] exp(−S[n(·), p(·)])
T
(3.27)
dt [p(t)∂t n(t) − H(n(t), p(t))]
(3.28)
0
H(n, p)
=
λ(e−p − 1)n
(3.29)
と書く。
特に、O[·] として、ある path ngiven (·) に n(·) が等しければ 1 を、等しくなければ 0 を返す汎関数をとると、O
は path ngiven (·) が実現する確率密度 P[ngiven (·)] を与える:
P[n(·)] =
3.1.2
Dp exp(−S[n(·), p(·)])
(3.30)
一般の場合（1 変数）
一般の離散状態・連続時間 Markov 確率過程についても、同様の定式化ができる。ここでは、1 変数に限ること
にし、
n → n + r (with rate Wr (n) ) (r ∈ R ⊂ Z − {0})
なる確率過程を時間 t ∈ [0, T ] で考える5 。
5
この反応は、master 方程式で表せば (時刻 t に n である確率を P (n, t) とおいて) ∂t P (n, t) =
r∈R Wr (n)P (n, t) となる。
24
(3.31)
r∈R
Wr (n − r)P (n − r, t) −
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
時間幅 Δt が十分小さければ
⎧
⎨r (r ∈ R, with probability W (n)Δt)
r
n(t + Δt) − n(t) =
⎩0 (with probability 1 − r∈R Wr (n)Δt)
(3.32)
が成り立つ。
時間 [0, T ] を M ( 1) 個に分割し、時間を離散化して考える：
Δt := T /M, tj := jΔt (j ∈ {0, 1, · · · , M }), nj := n(tj )
(3.33)
前節と同様にして、n(·) に依存する物理量 O の期待値は
∞
O O[n0 , n1 , · · · , nM ] · Prob[n0 → n1 ] · Prob[n1 → n2 ] · · · Prob[nM −1 → nM ] (3.34)
n1 ,n2 ,··· ,nM =0
と書ける。
ここで
Jj := nj+1 − nj
と定義すると、
Jj =
(3.35)
⎧
⎨r
( with probability Wr (nj )Δt)
(r ∈ R)
⎩0 ( with probability 1 − r∈R Wr (nj )Δt)
(3.36)
となる。これ以降、n(·):given に対して、J を (3.36) 式で n(·) から定まる新しい確率過程と見なすことにする。
Pj (Δn|nj ) = (時刻 j で n = nj であったときに Jj = Δn である確率)
(3.37)
とおくと、ここでは
⎧
⎪
⎪
(r ∈ R の場合)
⎪Wr (nj )Δt
⎨
Pj (r|nj ) = 1 − r∈R Wr (nj )Δt (r = 0 の場合)
⎪
⎪
⎪
⎩0
(r ∈
/ R かつ r = 0 の場合)
(3.38)
となる。
先ほどと同様の計算により、
O
M
dni O[n0 , n1 , · · · , nM ]
i=1
=
M
−1
δ(nj+1 − nj − Jj )|nj
(3.39)
j=0
M −1
M −1
M
dnk dpk
pi+1 (ni+1 − ni )
exp(pj+1 Jj ) |nj
O[n0 , n1 , · · · , nM ] exp −
2πi
i=0
j=0
(3.40)
k=1
と表せる。
exp(pj+1 Jj )|nj
=
e
rpj+1
Wr (nj )Δt
+
1−
r∈R
=
1+
Wr (nj )Δt
(3.41)
r∈R
(erpj+1 − 1) Wr (nj )Δt
r∈R
exp Δt
(3.42)
(erpj+1 − 1) Wr (nj )
r∈R
25
(3.43)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と変形できるので、
O
M
dnk dpk
O[n0 , n1 , · · · , nM ]
2πi
k=1
⎧
⎫
−1
⎨ M
⎬
pj+1 (nj+1 − nj ) − Δt
(erpj+1 − 1)Wr (nj )
× exp −
⎩
⎭
j=0
(3.44)
r∈R
M → ∞ とすれば、3.44 式の は = と見なしてよいと期待される:
O
=
M
dnk dpk
O[n0 , n1 , · · · , nM ] exp(−S[n, p])
lim
M →∞
2πi
(3.45)
k=1
S[n, p]
=
M
−1
[pj+1 (nj+1 − nj ) − Δt · H(nj , pj+1 )]
(3.46)
j=0
H(n, p)
=
(erp − 1)Wr (n)
(3.47)
r∈R
これを
O
DnDp O[n(·)] exp(−S[n(·), p(·)])
=
S[n(·), p(·)]
H(n, p)
T
=
=
0
(3.48)
dt [p(t)∂t n(t) − H(n(t), p(t))]
(3.49)
(erp − 1)Wr (n)
(3.50)
r∈R
と書く。
特に、O[·] として、ある path ngiven (·) に n(·) が等しければ 1 を、等しくなければ 0 を返す汎関数をとると、O
は path ngiven (·) が実現する確率密度 P[ngiven (·)] を与える:
P[n(·)] =
3.1.3
Dp exp(−S[n(·), p(·)])
(3.51)
多変数の場合（反応拡散系）
ここまでは、1 変数の確率過程を考えてきた。この項では、多変数の確率過程、特に、格子上の反応拡散系を考
える: 具体的には、格子の各点に粒子が分布し、粒子が各点での反応と拡散を行うような系を考える。
簡単のため、粒子の種類は 1 種類とする。まず、反応はないとして、拡散のみを取り扱うことにする。格子点を
x ∈ Λ（Λ は有限集合と仮定する）、時刻 t での各点 x での粒子の数を nx (t) とする。時間 t から t + Δt（Δt は十
分小さいとする）の間に粒子が 1 つの点 x から点 y(y = x) に拡散で移動する確率は、位置の関数 Wx,y を用いて
nx (t)Wx,y Δt
(3.52)
と表されると仮定する。この確率過程を時間 t ∈ [0, T ] で考える。
nx (t) をまとめて n(t) と表すと、ここでの拡散過程の定義より
⎧
y
x
⎨(0, · · · , 0, −1,
ˇ 0, · · · , 0, 1̌, 0, · · · , 0) (with probability nx (t)Wx,y Δt (x = y))
n(t + Δt) − n(t) =
(3.53)
⎩(0, · · · , 0)
(with probability 1 − x,y∈Λ nx (t)Wx,y Δt)
と書ける。
26
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
時間 [0, T ] を M ( 1) 個に分割し、時間を離散化して考える：
Δt := T /M, tj := jΔt (j ∈ {0, 1, · · · , M }), nj := n(tj ), nx,j := nx (tj )
(3.54)
前項と同様にして、n(·) による物理量 O の期待値は
O
O[n0 , n1 , · · · , nM ] · Prob[n0 → n1 ] · Prob[n1 → n2 ] · · · Prob[nM −1 → nM ] (3.55)
n1 ,n2 ,··· ,nM
と書ける。
ここで、
Jj
:=
nj+1 − nj
(3.56)
Jx,j
:=
nx,j+1 − nx,j
(3.57)
と定義する。ただし、Jx,j は J j の x 成分、すなわち位置 x ∈ Λ での粒子数変化を表す。すると、元の確率過程の
定義より
⎧
y
x
⎨(0, · · · , 0, −1,
ˇ 0, · · · , 0, 1̌, 0, · · · , 0)
Jj =
⎩(0, · · · , 0)
(with probability nx,j Wx,y Δt (x = y))
(with probability 1 − x,y∈Λ nx,j Wx,y Δt)
(3.58)
となる。そこで、これ以降、n(·):given に対して、J を (3.58) 式で n(·) から定まる新しい確率過程と見なすこ
とにする。
3.1.1 項と同様の計算により
O
M
M
−1
dnz,k dpz,k
O[n0 , n1 , · · · nM ]
δ(nx,i+1 − nx,i − Jx,i ) 2πi
i=0
z∈Λ k=1
=
x∈Λ
M
dnz,k dpz,k
O[n0 , n1 , · · · nM ]
2πi
z∈Λ k=1
⎡
⎤
−1
M
−1
M
p
J
× exp ⎣−
py,j+1 (ny,j+1 − ny,j )⎦
e x,i+1 x,i y∈Λ j=0
i=0
x∈Λ
(3.59)
ni
(3.60)
ni
ここで、
x∈Λ
px,i+1 Jx,i e
⎛
=
ni
⎝1 −
⎡
exp ⎣
⎞
nx,i Wx,y Δt⎠ + ⎝
x,y∈Λ, x=y
$
⎛
⎞
epy,i+1 −px,i+1 nx,i Wx,y Δt⎠ (3.61)
x,y∈Λ, x=y
epy,i+1 −px,i+1
⎤
%
− 1 nx,i Wx,y Δt⎦
(3.62)
x,y∈Λ
なので、
O
M
dnz,k dpz,k
O[n0 , n1 , · · · nM ]
2πi
z∈Λ k=1
⎧
⎤⎫
⎡
−1 ⎨ M
⎬
$
%
⎣
× exp −
epy,j+1 −px,j+1 − 1 nx,j Wx,y Δt⎦
px,j+1 (nx,j+1 − nx,j ) −
⎩
⎭
j=0
x∈Λ
x,y∈Λ
M → ∞ とすれば (3.63) 式の は = と見なして差し支えないと期待される:
27
(3.63)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
O
=
S[n, p] =
M dnz,k dpz,k
lim
O[n0 , n1 , · · · nM ] exp(−S[n, p])
M →∞
2πi
k=1 z∈Λ
M
−1 px,j+1 (nx,j+1 − nx,j ) − H(nj , pj+1 )Δt
j=0
H(n, p) =
(3.64)
(3.65)
x∈Λ
$
%
epy −px − 1 nx Wx,y
(3.66)
x,y∈Λ
これを
O
=
S[n(·), p(·)] =
DnDp O[n(·)] exp(−S[n(·), p(·)])
T dt
px (t)∂t nx (t) − H(n(t), p(t))
0
H(n, p) =
$
(3.67)
(3.68)
x∈Λ
%
epy −px − 1 nx Wx,y
(3.69)
x,y∈Λ
と書く。
特に、O[·] として、ある path ngiven (·) に n(·) が等しければ 1 を、等しくなければ 0 を返す汎関数をとると、
O は path ngiven (·) が実現する確率密度 P[ngiven (·)] を与える:
P[n(·)] =
Dp exp(−S[n(·), p(·)])
(3.70)
次に、拡散に加えて各点での反応がある場合を考える。各点 x での反応が
nx → nx + r (with rate Wr (nx ) ) (r ∈ R ⊂ Z − {0})
(3.71)
と書けるとすると、上の計算の J を次にように改めればよい:
⎧
y
x
⎪
ˇ 0, · · · , 0, 1̌, 0, · · · , 0) (with probability nx,j Wx,y Δt (x = y))
⎪
(0,
·
·
·
,
0,
−1,
⎪
⎨
Jj =
x
(0, · · · , 0, ř, · · · , 0, · · · , 0)
⎪
⎪
⎪
⎩(0, · · · , 0)
(with probability Wr (nx,j )Δt)
(with probability 1 − x,y∈Λ, x=y nx,j Wx,y Δt − x∈Λ Wr (nx,j )Δt)
(3.72)
この J の表式を用いて同様な計算を行えば、
O =
DnDp O[n(·)] exp(−S[n(·), p(·)])
(3.73)
T
S[n(·), p(·)] =
dt
px (t)∂t nx (t) − H(n(t), p(t))
(3.74)
0
H(n, p) =
x∈Λ
(e
rpx
− 1)Wr (nx ) +
x∈Λ
$
%
epy −px − 1 nx Wx,y
(3.75)
x,y∈Λ
が得られる。
3.1.4
結果のまとめ
この節で得られた結果を以下にまとめておく。まず、1 種の粒子 A の粒子数変化に関する一般の確率過程
n → n + r (with rate Wr (n) ) (r ∈ R ⊂ Z − {0})
28
(3.76)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
に対して、時刻 0 から時刻 T までの A の粒子数 n(·) に依存した物理量 O[·] の期待値は、経路積分を用いて、
O =
DnDp O[n(·)] exp(−S[n(·), p(·)])
(3.77)
S[n(·), p(·)]
T
=
dt [p(t)∂t n(t) − H(n(t), p(t))]
(3.78)
(erp − 1)Wr (n)
(3.79)
0
H(n, p)
=
r∈R
と書ける。
さらに、粒子 A が上記の反応に加えて、格子点 Λ 上で拡散する系を考える。微小時刻間隔 [t, t + Δt] において、
格子点 x から格子点 y に 1 つの粒子 A が拡散で移動する確率を Wx,y Δt とする。時刻 t に格子点 x にある粒子 A
の粒子数を nx (t) とし、さらに粒子数の配置をまとめて n(·) とかいたとき、配置 n(·) に依存する物理量 O[·] の期
待値は
O
=
S[n(·), p(·)] =
DnDp O[n(·)] exp(−S[n(·), p(·)])
T dt
px (t)∂t nx (t) − H(n(t), p(t))
0
H(n, p) =
x∈Λ
$
%
epy −px − 1 nx Wx,y
(erpx − 1)Wr (nx ) +
x∈Λ
(3.80)
(3.81)
(3.82)
x,y∈Λ
と書ける。
本論文では、以上の経路積分表示を用いて解析を行う。
3.2
半古典近似とその応用
ここまで、種々の確率過程を経路積分によって記述する方法を見てきた。しかしながら、それらの経路積分はほ
とんどの場合、そのまま厳密に実行することはできない。従って、具体的な問題を解くにあたっては適切な近似を
行うことが必要不可欠となる。ここでは、後の解析に用いる「半古典近似」の基礎とその応用、特に、rare event
の評価精度について論じる。
3.2.1
半古典近似
時刻 0 で n = n0 であったとき、時刻 T で n = nT である確率を p(nT , T |n0 , 0) と書き、これを評価することを
考える。p(nT , T |n0 , 0) は物理量 O として「時刻 T で n = nT なら 1, そうでなければ 0 を返す」ものを取ること
により、経路積分を用いて
p(nT , T |n0 , 0)
S̃[n(·), p(·)]
DnDp e−S̃[n(·),p(·)]
=
=
(3.83)
n(0)=n0 ,n(T )=nT
dt[p(t)∂t n(t) − H̃(n(t), p(t))]
(3.84)
&
と表すことができる。ここで、 n(0)=n0 ,n(T )=nT DnDp は、n についての経路積分を n(0) = n0 , n(T ) = nT を満た
す n(·) に制限することを意味する。この経路積分で表された量を評価する。
ここで、次の仮定を置く: 系に特徴的な数 N ( 1) が存在して、
q :=
29
n
N
(3.85)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と変換することにより、
p(nT , T |n0 , 0)
S[q(·), p(·)]
DqDp e−N S[q(·),p(·)]
=
=
(3.86)
q(0)=n0 /N,q(T )=nT /N
dt[p(t)∂t q(t) − H(q(t), p(t))]
(3.87)
と書ける6 。ただし、N は積分の measure Dq に含まれているとする。このとき、p(nT , T |n0 , 0) を N 1 におい
て評価することを考える。
今、N が極めて大きいことを仮定しているので、S の実部が大きいところでは被積分汎関数は非常に小さくなっ
て経路積分に寄与せず、S の実部が小さいところからの寄与が非常に大きな割合を占めることが期待される7 。そ
こで、p(nT , T |n0 , 0) を、S が停留値を取る “点” での作用の値と、その周りでの低次のゆらぎを用いて評価する。
この近似を、半古典近似と呼ぶ8 。
S が停留値をとるような “点”9 を (qcl (·), pcl (·))10 とおき、古典解と呼ぶことにする。古典解は、「運動方程式」
∂t q(t) = ∂p H(q, p), ∂t p(t) = −∂q H(q, p)
(3.88)
q(0) = n0 /N, q(T ) = nT /N
(3.89)
と、境界条件
を満たす。
S を (qcl (·), pcl (·)) の周りで 2 次まで展開すると11 :
S[qcl (·) + δq(·), pcl (·) + δp(·)]
ただし、
A(t, t ) :=
1
dt dt δq(t)A(t, t )δq(t )
S[qcl (·), pcl (·)] +
2
1
dt dt δp(t)C(t, t )δp(t ) (3.90)
+ dt dt δq(t)B(t, t )δp(t ) +
2
δ 2 S[q, p]
δ 2 S[q, p]
δ 2 S[q, p]
)
:=
)
:=
,
B(t,
t
,
C(t,
t
δq(t)δq(t )
δq(t)δp(t )
δp(t)δp(t )
とした。これをもとの式に代入して δq(·), δp(·) についての経路積分を行うことにより、
'
N
−N S[qcl (·),pcl (·)]
dt dt δq(t)A(t, t )δq(t )
p(nT , T |n0 , 0) e
DδqDδp exp −
2
(
N
dt dt δp(t)C(t, t )δp(t )
−N dt dt δq(t)B(t, t )δp(t ) −
2
=
K e−N S[qcl (·),pcl (·)]
(3.91)
(3.92)
(3.93)
を得る。ただし、(3.92) 式の Gauss 積分の値を K とした。K を prefactor と呼ぶ。prefactor K の N 依存性は多く
の場合あまり大きくないので、その N 依存性を無視しても本論文の議論には大きな影響を与えないと期待される。
3.2.2
崩壊過程の半古典近似
簡単な例として、崩壊過程に半古典近似を適用し、厳密解に近い結果が得られることを見る。具体的には、時
刻 t = 0 での粒子数が N ( 1) としたときの、時刻 T での粒子数 nT の確率分布 p(nT , T ) を求める。この N が半
古典近似での大きいパラメーターとしての役割を果たす。0 ≤ nT ≤ N とする。
6 この後扱うような問題は全てこの仮定を満たす。
7 以下の論法は複素積分に対する鞍点法と全く同様である。
8 「半古典」という名前は、この手法が量子力学の経路積分に対する近似手法に由来することから用いられる。
9 函数空間内の
1点
は classical の cl
11 1 次の項は S が (q (·), p (·)) で停留値をとることから、0 である。
cl
cl
10 “cl”
30
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
3.1.1 節の結果 (3.29) より、q := n/N とおくと、
p(nT , T ) =
DqDp exp(−N S[q(·), p(·)])
(3.94)
q(0)=1,q(T )=nT /N
T
S[q(·), p(·)] :=
dt [p(t)∂t q(t) − H(q(t), p(t))]
(3.95)
0
λ(e−p − 1)q
H(q, p) :=
(3.96)
と書ける。
古典近似を適用して得られる運動方程式と境界条件はそれぞれ、
∂t q(t)
∂t p(t)
=
=
−λqe−p
−p
−λ(e
(3.97)
− 1)
q(0) = 1, q(T ) = nT /N =: qT
(3.98)
(3.99)
で与えられる。この方程式の解は、
q(t)
p(t)
C1
C2
C1 e−λt + C2
C2 λt
e
= log 1 +
C1
1 − qT
:=
1 − e−λT
qT − e−λT
:=
1 − e−λT
=
(3.100)
(3.101)
(3.102)
(3.103)
となる。これらの式を作用 S[·, ·] の表式 (3.95) に代入して計算すると、
−S[q, p] = −λqT T − (1 − qT ) log(1 − qT ) − qT log qT + (1 − qT ) log(1 − e−λT )
(3.104)
が得られる。従って、
p(nT , T ) e−N S[q,p] =
nnTT (N
NN
e−λT nT (1 − e−λT )N −nT
− nT )N −nT
(3.105)
となる。見易いように T を t に、nT を n に書き換えると
p(n, t) =
nn (N
NN
e−λtn (1 − e−λt )N −n
− n)N −n
(3.106)
を得る。
近似の精度：rare event の確率の評価
ここまでで半古典近似を用いて崩壊過程の確率分布を導出したが、この近似はどれくらい精度の良い近似であ
ろうか？そこでここでは、崩壊過程に対する半古典近似の結果を、厳密解、および他の近似手法で得られる結果
と比較する。
崩壊過程の master 方程式は厳密に解くことができて、
⎧
⎨ C e−λtn (1 − e−λt )N −n
N n
p(n, t) =
⎩0
31
(n ∈ {0, 1, · · · , N })
(otherwise)
(3.107)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
0
20
log p
40
60
exact
semi-classical
Gaussian
expectation value
80
1000
2
4
6
8
n
10
12
14
16
図 3.1: N = 15, λ = 1, t = 3 で 3 つの結果を比較
となる12 。半古典近似による結果 (3.106) を厳密解 (3.107) と比べると、前者は後者の N Cn を Stirling の公式で近
似したものと等しいことが分かる。
他の近似手法として、ここでは master 方程式の近似解法としてよく用いられる van Kampen 展開による linear
noise 近似を取り上げる。この近似は、求める確率分布を期待値と、その周りの Gauss 揺らぎに分解して評価する
ものである。計算の過程は付録 A.1.1 に回すが、結果として
(
'
(n − N e−λt )2
exp −
2σ 2 (t)
2πσ 2 (t)
1
p(n, t)
=
)
σ 2 (t)
:=
N (e−λt − e−2λt )
(3.108)
(3.109)
が得られる。
厳密解、半古典近似による結果、van Kampen 展開による結果を比べたのが図 3.1 である。図から、van Kampen
展開は期待値の周りでは良い近似を与えるが、期待値から外れたところでは正しい値から大きく外れることが分
かる。対して、半古典近似は期待値から離れた箇所で非常に良い近似を与えることが分かる。
このことから、半古典近似は期待値から大きく外れた領域で威力を発揮することが期待される。実際、一般に
半古典近似は linear noise 近似に比べて分布の裾、すなわち、rare な事象の確率の評価に効力を発揮することが知
られている ([32]subsection 10.4.4)。
3.3
stochastic population の寿命：bounce 解の手法
次の例として、確率的な生態系での集団の「寿命」を求める問題を考える。この問題に半古典近似を適用する
際には、「bounce 解」と呼ばれる一群の特徴的な解を考える13 。この手法は第 4 章の解析でも多用する。
n
n
:=
F についての偏微分方程式
n∈Z x P (n, t) とおき、master 方程式の両辺に x をかけて n についての和を取ると、
∂t F (x, t) = λ(1 − x)∂x F (x, t) が得られる。この方程式の解は、滑らかな函数 f を用いて F (x, t) = f (x − 1)e−λt と書ける。これと初期
N
が得られる。この F を x について展開すれば、P (n, t)
条件 F (x, 0) = xN より、f (y) = (y + 1)N なので、F (x, t) = 1 + (x − 1)e−λt
が求まる。
13 以下の議論は量子力学のトンネル問題で、metastable state の寿命を経路積分を用いて求めることと同様である。
12 F (x, t)
32
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
具体的には、次のような確率過程を考える: ある種の生物 A と「空地」∅ の合計の数は N で一定とし、A の数
は以下の反応で増減するとする:
b
遷移率は
d
A + ∅ → A + A, A → ∅ (b > d)
(3.110)
*
n+
, W − (n) = dn
W + (n) = bn 1 −
N
(3.111)
で表される14 。
この過程の定義から明らかなように、いったん系が n = 0 の状態になってしまうと、その後 n は永久に 0 に留
まる。系が最初 n = n0 であったとき、n = 0 になるまでの平均時間 τ はどれほどであろうか？換言すると、この
集団の “寿命”τ はどれほどであろうか？これが考えたい問題である。
3.3.1
決定論的解析
まず系の大まかな振る舞いを探るために、n の期待値の従う方程式を考える。
n/N の期待値を q とすると、master 方程式の van Kampen 展開（または Kramers-Moyal 展開を経由した計算）
から
b
∂t q(t) = bq(1 − q) − dq = (b − d)q 1 −
q
(3.113)
b−d
となる (付録 A.1.2)。今 b > d を仮定しているので、これは数理生物学で良く知られた logistic 方程式となる
([34]section 1.1)。qs := 1 − db はこの常微分方程式の安定固定点であり、初期条件が q = 0 であれば、t → ∞ で
q → qs となることが分かる。これより、元の確率過程に対しても、n0 = 0 であれば、n は N qs の近くに長時間留
まることが期待される。そこで、N qs =: ns とおく。
3.3.2
確率過程の解析
さて、再び確率過程に戻ろう。系の寿命を考えるために、「時刻 −T /2 に n = n0 であったときに時刻 T /2 に
n = nT である確率」p(n, T /2|n0 , −T /2) の、T が大きい（具体的には T 1/b）ところでの振る舞いを考える15 。
先ほども見たように系は (n = 0 にならなければ)ns の近くに留まることが期待されるので、p(ns , T /2|ns , −T /2)
を考えてしまっても差し支えない。この確率は T の増加につれて減衰し、T が大きいときの減衰の様子から、系
の典型的な寿命が分かると期待される。
p(ns , T /2|ns , −T /2) を経路積分で表すと、3.1 節の一般論より、
T
T
=
p ns , ns , −
DqDp e−N S[q(·),p(·)]
2
2
q(−T /2)=qs ,q(T /2)=qs
T /2
dt[p∂t q(t) − H(q(t), p(t))]
S[q(·), p(·)] :=
H(q, p)
と書ける（ただし、W± (q) :=
14 時刻
(3.114)
(3.115)
−T /2
p
:=
(e − 1)W+ (q) + (e−p − 1)W− (q)
=
b(ep − 1)q(1 − q) + d(e−p − 1)q
(3.116)
1
±
。
N W (N q)）
t で A の数が n である確率を P (n, t) とすると、master 方程式は
∂t P (n, t) = W + (n − 1)P (n − 1, t) + W − (n + 1)P (n + 1, t) − W + (n) + W − (n) P (n, t)
と書ける。
15 時刻を [0, T ] ではなく [−T /2, T /2] で考えているが、これはのちの便宜のためである。
33
(3.112)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
b= 3.0,d= 1.0
2.0
6.00
p
1.5
5.25
1.0
4.50
0.5
3.75
qs
0.0
3.00
2.25
0.5
1.50
1.0
ps
0.75
1.5
2.0
0.0
0.00
0.2
0.4
q
0.6
0.8
1.0
0.75
図 3.2: b = 3.0, d = 1.0 での相図。緑の丸は固定点、右側のカラーバーは Hamiltonian の値を表す。
半古典近似を適用すると、古典解の満たすべき運動方程式と境界条件は
∂t q(t) =
bep q(1 − q) − de−p q
∂t p(t) =
−b(e − 1)(1 − 2q) − d(e
p
(3.117)
−p
− 1)
q(−T /2) = qs , q(T /2) = qs
(3.118)
(3.119)
となる。これを厳密に解くことはできないが、相図をプロットすると解の振る舞いが分かる16 （図 3.2）。
各軌道上で H の値は一定であることに注意（energy 保存）。特に、p = 0 の軸上では常に H = 0 が成り立つ。ま
た、p = 0 の軸と交差する軌道（(qs , 0) と (0, ps ) を結ぶ軌道）は activation trajectory と呼ばれる [35]。activation
trajectory の表式を pa (q) とすると、activation trajectory 上で H = 0 が成り立つことと p = 0 より、
b
W+ (q)
= − log
(1 − q)
pa (q) = − log
W− (q)
d
(3.120)
と表されることが分かる。
さて、(3.119) 式の境界条件を満たす軌道を考えよう。厳密に T を有限に保って考えると、この境界条件を満た
すのは (qs , 0) に留まり続ける軌道しかありえない。しかし今は bT 1 なる状況を考えているので、T → ∞ な
る極限を近似的に考えることができる。この極限においては、activation trajectory も境界条件 (3.119) を満たす。
また、(qs , 0) → (0, ps ) → (0, 0) → (qs , 0) なる閉軌道も（T → ∞ で近似的に）この境界条件を満たす軌道である。
更に、この閉軌道を複数回まわるような閉軌道も条件 (3.119) を満たす。これらの軌道を “bounce 解” と呼ぶ。
16 この相図は軌道を計算して求めたものではなく、Hamiltonian
の値が等しい曲線（等 energy 曲線）をプロットしたものである。
34
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
1013
1012
1011
1010
109
108
107
106
105
104
103
102
101
100 0
b = 1.0, d = 0.5
1014
1013
1012
1011
1010
109
108
107
106
105
104
103
102
101
1000.2
numerical
theoretical
τ
τ
numerical
theoretical
20
40
60
80
N
100
120
140
160
0.3
0.4
0.5
0.6
0.7
d
0.8
0.9
1.0
図 3.3: 左：b = 1.0, d = 0.5 での τ の N 依存性 右：N = 50, b = 1.0 での τ の d 依存性
全ての bounce 解からの寄与をまとめると（activation trajectory に対する作用を Sa とおいて）、
t3
T /2
∞ t2
$
%n
p(ns , T /2|ns , −T /2) dt1
dt2 · · ·
dtn Ke−N Sa
n=0
∞
1
n!
n=0
=
−T /2
−T /2
T /2
−T /2
dt1 · · ·
(3.121)
−T /2
T /2
−T /2
dtn
$
Ke−N Sa
%n
∞
1 $ −N Sa %n
Ke
T
n!
n=0
%
$
= exp Ke−N Sa T
=
(3.122)
(3.123)
(3.124)
ただし、1 行目では、bounce 解からの寄与が、それぞれの activation trajectory からの寄与の積で書けると仮定
した17 （(0, ps ) → (0, 0) → (qs , 0) の軌道は H = 0 と pq̇ = 0 より作用に寄与しないことに注意）。
更に、K を |K| ∼ 1 と近似する。ただし、確率は 1 以下であることから K の符号は正であってはならないの
で、K ∼ −1 とする。以上の近似の下に、
$
%
p(ns , T /2|ns , −T /2) ∼ exp −e−N Sa T
(3.125)
τ ∼ eN S a
(3.126)
が得られる。
従って、系の寿命を τ とおくと
と表され、Sa は
Sa =
0
1− db
1− db
dq pa (q) =
dq log
0
b
(1 − q)
d
= log
b
d
−1+
d
b
(3.127)
となる。
以上の計算の結果を master 方程式から数値計算で求めた結果と比べたのが図 3.3 である18 。N 依存性では N が
大きくなるにつれて、d 依存性では d が小さくなるにつれて、近似の精度が上がり、共にそれらの領域では半古典
近似の結果と数値計算結果がよく一致していることが見て取れる。
なお、本論文では経路積分とそれに対する半古典近似を用いて確率過程の解析を行っているが、[35] など多くの
研究では Master 方程式に直接的に WKB 法 [36, 37] を適用して同様の現象の解析を行っている。この 2 つの手法
17 これは量子トンネル効果の解析に対して通常行われる近似である。
18 master
方程式を直接解くのではなく、master 方程式の右辺に現れる行列の固有値を求めた。
35
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
は言わば量子力学における経路積分法と演算子法のような関係にあり、(量子力学の場合と同様に) 経路積分法に
は dynamics をより直感的に捉えられるという利点がある。
36
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
第 4 章 モデルと解析
2.5.5 項でも述べたように、集団の構造が進化ゲームの dynamics に対してどのような影響を及ぼすかについて
の研究は、まだ緒に就いたばかりである。動的な振る舞いとしては様々な問題が考えられるが、本研究では特に
「準安定状態間の遷移」という現象に着目し、この現象、特に準安定状態の寿命に対する空間自由度 (集団の構造
として最も単純なもの) の影響を探る。
「準安定状態」とは、系が極めて長い時間留まるような状態のことである。確率過程を考えると、準安定状態
は、与えられた確率過程に対応する決定論的力学系の安定状態に対応する1 。決定論的力学系では系が安定状態に
落ち着けば永遠にそこに留まるが、確率過程では極めて低い確率ではあるものの、それらの「安定状態」の間を
遷移することが可能である。この現象を、本論文では「準安定状態間の遷移」と呼ぶ。
準安定状態間の遷移は本質的に確率的なものであり、かつ非常に稀な現象である。しかし多くの場合、異なる準
安定状態間の遷移は、考える系で状態が大きく変化することに対応するため、この遷移を考えることは応用上の
観点からも重要である: このような現象の例としては、生態系における集団の絶滅や regime shift [38]、疫病の流
行と収束 [39]、細胞における遺伝子スイッチの切り替え [40] などが挙げられる。
この章では、進化ゲーム理論のモデルにおいて、空間自由度が準安定状態間の遷移率にどのような影響を及ぼ
すかを探る。まず、well-mixed モデルと空間構造のあるモデル双方の設定を述べる (4.1 節)。その後、well-mixed
モデルについての決定論的解析 (4.2 節) と準安定状態間の遷移率の評価 (4.3 節) を行う。これを踏まえて最後に、
空間構造のあるモデルについても同様の解析を行う (4.4 節、4.5 節)。
4.1
モデル
A,B の 2 種の生物がなす系を考える。これらの生物は進化ゲームのルールに従って増殖するとする。ゲームの
利得行列は
a b
c d
(4.1)
と書かれ、
a > c, b < d
(4.2)
を満たすとする。すなわち、ゲームは双安定な coordination game とする (c.f. 2.2.2 項)。
空間自由度が系の dynamics に与える影響を探るために、全ての個体が他の全ての個体と等しく相互作用するモ
デル (well-mixed モデル) と、周りの個体としか相互作用できないモデル (spatial モデル) の 2 種類のモデルを考
える。それぞれのモデルの詳細な設定を以下に述べる。
4.1.1
well-mixed モデル
A,B の 2 種の生物が 1 つのパッチに合計 N ( 1) 個体生息し、全ての個体が全ての個体と等しくゲームを行う
とする。
繁殖ルールとしては pairwise comparison process (2.3.1 項) を採用し、また、繁殖の際には突然変異があるとす
る。突然変異率は十分小さく、この系に対応する決定論的な dynamics は双安定性を持つとする（詳細は後述）。ま
た、pairwise comparison process を考えることから、payoff の範囲は 0 < a − c < 1, 0 < d − b < 1 に制限する。
1 安定状態は複数あるとする。
37
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
図 4.1: well-mixed モデル
A の個体数を n とし、n の dynamics を考える。パラメーターを以下のように定義する:
λ
:
繁殖の起こるレート
(4.3)
w
:
自然選択の強さ
(4.4)
μA
:
A の突然変異率 (A が B を産む確率)
(4.5)
μB
:
B の突然変異率 (B が A を産む確率)
(4.6)
ゲームと繁殖の過程は pairwise comparison process に従うとし、focal 個体が role 個体の戦略を採用する確率
p : R → [0, 1] は p(Δπ) =
ける:
+
W (n)
:=
W − (n)
:=
1+wΔπ
2
で与えられるとする (2.3.1 項参照)。このとき個体数 n の遷移率は次のように書
n
W (n → n + 1) = λ(1 − μA )
N
n
W (n → n − 1) = λ(1 − μB )
N
n + 1 + w ΠA (n) − ΠB (n)
μB *
n +2
1−
+λ
1−
N
2
2
N
*
n + 1 + w ΠB (n) − ΠA (n)
μA * n +2
1−
+λ
N
2
2 N
*
(4.7)
(4.8)
ここで、ΠA (n), ΠB (n) は A の個体数が n であるときの個体 A,B の利得をそれぞれ表し、次のように与えられる:
*
n+
n
a+ 1−
b
(4.9)
ΠA (n) :=
N
N+
*
n
n
ΠB (n) :=
c+ 1−
d
(4.10)
N
N
これらの遷移率を用いると、Master 方程式は（時刻 t で A の個体数が n である確率を P (n, t) とおく）:
∂t P (n, t) = W + (n − 1)P (n − 1, t) + W − (n + 1)P (n + 1, t) − [W + (n) + W − (n)]P (n, t)
(4.11)
と書ける。ただし、n < 0, n > N に対しては P (n, t) ≡ 0 とした。
4.1.2
spatial モデル
次に、空間自由度を持つモデルを考える。M 個のパッチが 1 次元状につらなっており、周期的境界条件をみた
すとする。各パッチ上には Np ( 1) 個の個体がいるとし、パッチ間の距離を l、L̃ := M l とおく。
個体はパッチ内でゲームと繁殖を行い、また、隣接するパッチの間で移動することができるとする。移動の際
には隣接するパッチの間で個体が「交換」されることを仮定し、各パッチ上の個体数は一定に保たれるとする。
ni を i 番目のパッチ上の A の個体数 (i ∈ {1, 2, · · · , M }) とし、ni の dynamics を考える。パラメーターを以下
38
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
図 4.2: spatial モデルの概念図。パッチは左右に続いていることに注意（周期的境界条件を課す）。
のように定義する（σ 以外は well-mixed モデルと同様）:
λ
:
繁殖の起こるレート
(4.12)
w
:
自然選択の強さ
(4.13)
μA
:
A の突然変異率 (A が B を産む確率)
(4.14)
μB
:
B の突然変異率 (B が A を産む確率)
(4.15)
σ
:
パッチ間の移動の起こるレート
(4.16)
(4.17)
ゲームと繁殖の過程は前項の well-mixed モデルと同様の pairwise comparison process に従うとする。このと
き、各パッチ内での繁殖による個体数 ni の遷移率は、次のように書ける:
2
ni
ni 1 + w ΠA (ni ) − ΠB (ni )
ni
μB
1−
W + (ni ) := W (ni → ni + 1) = λ(1 − μA )
1−
(4.18)
+λ
Np
Np
2
2
Np
2
ni 1 + w ΠB (ni ) − ΠA (ni )
ni
μA ni
−
1−
W (ni ) := W (ni → ni − 1) = λ(1 − μB )
(4.19)
+λ
Np
Np
2
2
Np
ここで、ΠA (n), ΠB (n) は A の個体数が n であるのときの個体 A,B の利得をそれぞれ表し、次のように与えられ
る2 :
A
Π (ni )
:=
ΠB (ni )
:=
ni
ni
b
a+ 1−
Np
Np
ni
ni
d
c+ 1−
Np
Np
(4.20)
(4.21)
次に、隣接パッチ間で個体の移動 (入れ替え) を考える。入れ替えは隣接するパッチ間で rate σ で起こり、それぞ
れのパッチから 1 個体ずつ random に選んで交換する:
W ((ni , nj ) → (ni − 1, nj + 1)) = σ
ni Np − nj
Np Np
(4.22)
これらの遷移率により、系の確率的 dynamics が決まる。
4.2
well-mixed モデルの決定論的解析
まず、モデルの基本的な性質を理解するために、well-mixed モデルの決定論的な振る舞い (n の期待値の振る舞
い) を考える。ここでの解析から、系が双安定性を持つことが示される。
2 well-mixed
モデルの場合と同一。
39
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
q
q
1
q
2
3
0.02
+
W (q)−W−(q)
0.00
0.04
0.06
0.0
0.2
0.4
q
0.6
0.8
1.0
図 4.3: 期待値の dynamics を表す相図。緑の丸が固定点、赤の破線が dq/dτ の値を表す。 a = 1.0, b = 0.3, c =
0.6, d = 1.3, w = 0.5, μA = 0.01, μB = 0.01
n(t) の期待値を n(t) とおき、q(t) :=
n(t)
N
とする。q の従う方程式（“rate 方程式”）は master 方程式を van
Kampen 展開して最低次の項を取れば
d
1 +
q(t) =
W (N q(t)) − W − (N q(t))
dt
N
(4.23)
3
となる。 λt
N を τ と置くと
d
q
dτ
W± (q)
:=
ΠA (q)
:=
W+ (q) − W− (q)
(
' μA + μB
μB − μ A
(ΠA (q) − ΠB (q)) +
q(1 − q) w 1 −
2
2
μB
μ
A
+
(1 − q)2 −
q2
2
2
1 ±
W (N q)
λ
aq + b(1 − q)
ΠB (q)
:=
cq + d(1 − q)
=
=
(4.24)
(4.28)
(4.25)
(4.26)
(4.27)
この力学系の相図は図 4.3 のようになる。μA , μB が十分小さいとき、系は 3 つの固定点を持つ。図からも分
かるように、q1 , q3 は安定固定点、q2 は不安定固定点になっている。つまり、本質的には coordination game の
replicator dynamics と同様な、双安定性を持つ系になっている。今後、このようなパラメーター範囲に限定して
考察を進める。
4.3
well-mixed モデルにおける準安定状態間の遷移
先ほども見たように、条件 (4.2) を満たす 2 × 2 ゲームを well-mixed モデルで決定論的に考えると、系は双安定
性を持つ。今考えているモデルでは突然変異があるので、系は吸収状態を持たない。決定論では系がどちらかの
安定状態に落ち着けばそのまま永遠にそこに留まり続けるが、確率過程を考えると、ごくごくまれに系はこの状
態からもう片方の安定状態に遷移する可能性がある (図 4.4)。この意味において、決定論での 2 つの安定状態は確
率過程を考えると安定ではなく、「準安定状態」と呼ばれる（[31] の Chap. XIII）。
3 本来であれば
d
Q(τ ) := q(N τ /λ) などとおき、 dτ
Q(τ ) = W+ (Q(τ )) − W− (Q(τ )) などと書くべきである。しかしここでは物理の慣習
に従い、同一の記号 q で異なる関数を用いている。 以下、確率過程や空間モデルの解析でも同様。
40
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
䛾ಶయᩘ
䜋䜌
䜋䜌
図 4.4: Monte Carlo simulation の結果と準安定状態間の遷移の概念図。a = d = 1.5, b = c = 1.0, w = 0.4, μA =
μB = 0.01, N = 50
これらの準安定状態間の遷移率（もしくはその逆数である寿命）を評価し、系の諸々の性質にどのように依存
するかを考える。解析にあたっては 3.1 節で述べた経路積分の手法で系を記述し、半古典近似、特に 3.3.2 項で述
べた bounce 解の手法を用いる。
4.3.1
経路積分表示
まずは 3.1 節に述べた経路積分の手法でこのモデルを記述する。
時刻 0 から時刻 T̃ までの振る舞いを考えると、path n(·) の確率的な重み P[n(·)] は、3.1 節に述べたように
P[n(·)] = Dp e−S̃[n(·),p(·)]
(4.29)
と書くことができる。ただし、
T̃
S̃[n(·), p(·)] =
H̃0 (n, p)
=
dt[p(t)∂t n(t) − H̃0 (n(t), p(t))]
0
p
(e − 1)W + (n) + (e−p − 1)W − (n)
λT̃
q := n/N とおき、 λt
N を τ 、 N を T とおくと、path q(·) の確率的な重み P[q(·)] は
P[q(·)] = Dp e−N S[q(·),p(·)]
(4.30)
(4.31)
(4.32)
と書ける。ここで、作用 S[·, ·] は (4.26) 式で定義した W± を用いて、
T
S[q(·), p(·)] =
H0 (q, p)
=
dτ [p(τ )∂τ q(τ ) − H0 (q(τ ), p(τ ))]
0
p
(e − 1)W+ (q) + (e−p − 1)W− (q)
と表すことができる。以後、時間はこの無次元化した時間 τ で考えることにする。
41
(4.33)
(4.34)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
1.0
0.168
0.144
0.5
0.120
p
0.0
q1
q3
q2
0.096
0.072
0.048
0.5
0.024
0.000
1.0
0.0
0.2
0.4
q
0.6
0.8
1.0
図 4.5: 半古典近似で得た運動方程式の相図。緑の丸が固定点、右のカラーバーは Hamiltonian の値を表す。 a =
1.0, b = 0.3, c = 0.6, d = 1.3, w = 0.5, μA = 0.01, μB = 0.01
4.3.2
半古典近似
前項で導出した経路積分表示に対して半古典近似を行う。得られる運動方程式は4
q̇
=
ṗ
=
∂H0
= ep W+ (q) − e−p W− (q)
∂p
∂H0
= −(ep − 1)W+ (q) − (e−p − 1)W− (q)
−
∂q
(4.35)
(4.36)
この相図は、μA , μB が十分小さいところでは図 4.5 のようになる。3.3.2 項でも述べた通り、各軌道上で H0 の
値は一定であることに注意（energy 保存）。特に、p = 0 の軸上では常に H0 = 0 が成り立つ。また、p = 0 の軸と
交差する軌道（(q1 , 0), (q2 , 0), (q3 , 0) を結ぶ軌道）は activation trajectory と呼ばれる [35]。activation trajectory
上でも H0 = 0 なので、その表式を pa (q) とおくと、H0 の定義 (4.34) より
pa (q) = − log
W+ (q)
W− (q)
(4.37)
と表される。
半古典近似の下、T 1 なる時間 T に対して
• 系が時刻 0 に q1 にあって、そのまま一度も q3 に達せずに時刻 T にも q1 にいる確率5
• 系が時刻 0 に q3 にあって、そのまま一度も q1 に達せずに時刻 T にも q3 にいる確率
を求めれば、これらの確率の時間 T に対する依存性から、各双安定状態の寿命が求まることになる（前者は q1 近
辺の準安定状態の寿命、後者は q3 近辺の準安定状態の寿命を与える。）。
まず、後者の確率 i.e.「系が最初 q3 にあって、そのまま一度も q1 に達せずに時刻 T にも q3 にいる確率」を考え
よう。この条件を満たす古典的な軌道は、相図を見ると、q3 に留まり続ける軌道か、q3 を通る閉軌道しかないこと
が分かる6 。3.3.2 項で述べたように、この閉軌道を何度も回る軌道も上記の条件を満たす。これらの解を bounce
4 W は W の q 微分を表す
±
±
5 厳密に q である必要はなく q の近くであればよいが、ここでは簡単のために
1
1
6 T は十分長いとしているので、その他の軌道は不適当。
42
q1 をとる。q3 についても同様。
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
解と呼ぶのも、3.3.2 項で既に述べたとおりである。
従って、3.3.2 項と同様にして、求める確率 P3 (T ) は
P3 (T )
Sq3
∼
$
%
exp −e−N Sq3 T
(4.38)
[pa (q(τ ))q̇(τ ) − H0 (q(τ ), pa (q(τ )))] dτ =
:=
closed orbit through
q2
q3
pa (q)dq
(4.39)
q3
ここで、activation trajectory 上では H0 = 0 であることを用いた。(4.39) 式より、Sq3 は図 4.5 の斜線部分の面積
に等しいことが分かる。また、(4.38) 式より、q3 近辺の準安定状態の寿命を τq3 とおくと、
τq3 ∼ exp(N Sq3 )
(4.40)
となることが分かる7 。q1 近辺の準安定状態の寿命も全く同様に計算できる。それぞれの計算の結果は 5 章に述
べる。
最後に、bounce 解を構成する activation trajectory の構造について注意しておく。図 4.5 より分かるように、q3
周りの bounce 解を構成する activation trajectory は決定論的解析で得られた安定固定点 q3 と不安定固定点 q2 を
結ぶものであった。不安定固定点 q2 は q1 と q3 の “分水嶺” の役割を果たしているので、bounce 解は「系が q3 か
らスタートして、q1 に引き込まれる境界ぎりぎりのところまで進んでから、再び q3 に戻ってくる」軌道であると
も言える。q1 周りの bounce 解についても同様である。このことは、spatial モデルの解析を行う際に利用する。
4.4
spatial モデルの決定論的解析
次に、空間的に広がったモデル（4.1.2 項）の決定論的解析について述べる。
まず、個体数の期待値が従う微分方程式を導出する。これを空間が連続とみなして連続化し、定常解とその安定
性を求める。
rate 方程式
4.4.1
ここでは ni (i ∈ {1, 2, · · · , M }) の期待値をそのまま ni と書くことにする8 。ni の dynamics を表す方程式は:
dni
ni+1 N − ni
ni−1 Np − ni
ni Np − ni+1
ni Np − ni−1
+
−
−
= W + (ni ) − W − (ni ) + σ
dt
Np
Np
Np
Np
Np
Np
Np
Np
(i ∈ {1, 2, · · · , M })
(4.41)
と書ける。
qi :=
ni
Np
λt
と書き、 N
を τ と置くと、
p
dqi
dτ
=
W+ (qi ) − W− (qi ) +
=
F (qi ) +
σ
(qi+1 + qi−1 − 2qi )
λ
σ
(qi+1 + qi−1 − 2qi )
λ
(4.42)
(4.43)
ただし、
F (q)
:=
W± (q)
:=
W+ (q) − W− (q)
1 ±
W (Np q)
λ
とおいた。
7 ここでは λ で無次元化した時間を用いているが、通常は
N
8 添字 i はパッチの場所を表す。
λ = O(N )
43
(4.44)
(4.45)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
4.4.2
連続化
σ λ を仮定すると、局所的な反応に比べて移動がはるかに速く起こることになる。従って、qi は i の関数と
して滑らかに変化することが期待される。そこで、この条件下で連続的な記述を用いることにする:
q(x, τ ) := q xl (τ )
(x ∈ [0, L̃])
(4.46)
とおくと（l: パッチ間の “距離”）、
qi+1 (τ ) + qi−1 (τ ) − 2qi (τ ) = l2
となる。よって、
D :=
∂2q
(x, τ ) + O(l4 )
∂x2
σ 2
l
λ
(4.47)
(4.48)
とおいて l の高次の項を無視すれば、
∂q(x, τ )
∂ 2 q(x, τ )
+ F (q(x, τ ))
=D
∂τ
∂x2
が得られる。
√
x/ D を ξ とおいて無次元化すれば、最終的に以下の式を得る:
∂ 2 q(ξ, τ )
L̃
∂q(ξ, τ )
=
+ F (q(ξ, τ ))
ξ ∈ [0, L], L := √
∂τ
∂ξ 2
D
4.4.3
(4.49)
(4.50)
定常解
well-mixed モデルの決定論的解析では、dynamics の特徴を捉えるにあたって固定点とその (線型) 安定性が大
きな役割を果たした。また、半古典近似を用いた解析においても、bounce 解を構成する activation trajectory は、
決定論的解析で得られた安定固定点と不安定固定点を結ぶものとして与えられた (4.3.2 項の最後)。spatial モデル
の解析でも定常解とその安定性が同様の役割を果たすことが期待されるので、ここで定常解とその安定性を調べ
ておく。
q の方程式 (4.50) の定常解を qs (·) とおくと、qs の満たす方程式は
d2 qs (ξ)
+ F (qs (ξ)) = 0
dξ 2
となる。これを積分すると、
2
dqs (ξ)
+ V (qs (ξ)) = Const.
dξ
q
V (q) :=
F (q)dq
1
2
(4.51)
(4.52)
(4.53)
0
が得られる。
ここから、ξ をあたかも時間のように見なすことによって、qs を 1 次元空間上で potential V を受けて運動する
粒子の座標と解釈することができる。従って、qs は (4.52) 式を満たす (q, dq/dξ) 平面上での軌道と捉えることが
できる。(4.52) 式を満たす軌道を (q, dq/dξ) 平面上にプロットすると、図 4.6 のようになる。
ここでは周期的境界条件を課しているので、qs として実現できるのは、(q, dq/dξ) 平面上での閉軌道、または 1
点に限られる。ゆえに、定常解としてありうるのは、
• qs (ξ) ≡ q1 or q2 or q3 の一様解
• 図の q2 を囲む、非一様な閉軌道のうちの 1 つ
のいずれかに限られる。後者の非一様な閉軌道解を「臨界核」と呼ぶことにする (具体的な系においてどの閉軌道
が定常解になるかは、系のサイズ L に依存する (後述)。)。一様定常解と臨界核のそれぞれについて、その性質を
調べる。
44
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
0.20
0.0150
0.15
0.0075
0.10
0.0000
dq/d
0.05
0.00
q2
q1
q3
0.0075
0.0150
0.05
0.0225
0.10
0.0300
0.15
0.0375
0.20
0.0
0.2
0.4
q
0.6
0.8
1.0
図 4.6: 定常解の “相図”。ここでは V (q1 ) > V (q3 ) である。 （a = 1.0, b = 0.3, c = 0.6, d = 1.3, w = 0.8, μA =
μB = 0.005）
一様定常解
L の値に関わらず、一様定常解 q1 , q2 , q3 は存在する。これらの (線型) 安定性を考える。
q(ξ, τ ) = qi + δq(ξ, τ ) とおいて (4.50) 式に代入して線形化すると、
∂τ δq(ξ, τ ) = ∂ξ2 + F (qi ) δq(ξ, τ )
(4.54)
となる。この δq の振る舞いを考えたい。
δq(ξ, τ ) の ξ に対する周期性より Fourier 級数展開
δq(ξ, τ ) =
cn (τ )ei
2πn
L ξ
(4.55)
n∈Z
ができる。これを (4.54) 式に代入すると、
2
dcn (τ )
2πn
+ F (qi ) cn (τ ) =: λin cn (τ ) (∀n ∈ Z, ∀i ∈ {1, 2, 3})
= −
dτ
L
(4.56)
well-mixed モデルの相図からも分かるように
F (q1 ) < 0, F (q2 ) > 0, F (q3 ) < 0
(4.57)
なので、i = 1, 3 に対しては ∀n ∈ Z (λin < 0)、i = 2 に対しては ∃n ∈ Z (λin > 0) が成り立つ。従って、
qs (ξ) ≡ q1 or q3 は線形安定、qs (ξ) ≡ q2 は線形不安定であることが分かる。
臨界核
次に、臨界核の性質を考える。図 4.6 を見て分かる通り、閉軌道は無数にあった。しかし、今は系に周期的境界
条件を課しているので、どの閉軌道が定常解になるかは、系のサイズ L に依存する: 具体的には、ちょうど “時
間”L で一周するような軌道が解になる。
45
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
q3
−
q ξ plot)
critical nucleus (
1.0
critical nuclei (phase space plot)
0.20
0.0150
0.15
0.0075
0.8
q2
0.10
0.0000
0.05
q
dq/dξ
0.6
Lc
q1
q3
q2
0.00
0.0075
0.0150
0.4
0.05
0.0225
0.10
0.2
0.0300
0.15
L =13.8
L =14.5
0.0
0
5
10
15
20
ξ
0.0375
L =40.0
q1
25
30
35
40
0.20
0.0
0.2
0.4
q
0.6
0.8
1.0
図 4.7: 臨界核の例 。Lc は臨界核の典型的なスケール（a = 1.0, b = 0.3, c = 0.6, d = 1.3, w = 0.8, μA = μB = 0.005,
このパラメーターでは Lc 13.5）。
臨界核を解析的に求めるのは困難なので、数値的に求めた: (4.51) 式を離散化し、これを高次元の代数方程式と
見なして Newton 法を適用した。この手法によって求めた臨界核の例を図 4.7 に示す。臨界核は、L が小さくなる
ほど q2 に近づき、L が大きくなるほど q3 を通る閉軌道 (homoclinic 軌道) に近づく。また、L がある値 Lc より小
さいと臨界核は存在しないことが予想され、
2π
Lc = )
F (q2 )
(4.58)
と求められる9 。Lc は臨界核の典型的な大きさでもある。
次に、臨界核の安定性を検討する。臨界核を qs (·) と表す。一様解と同様に q(ξ, τ ) = qs (ξ) + δq(ξ, τ ) とおいて
(4.50) 式に代入して線形化すると、
∂τ δq(ξ, τ ) = ∂ξ2 + F (qs (ξ)) δq(ξ, τ )
(4.59)
L := ∂ξ2 + F (qs (ξ))
(4.60)
が得られる。
の固有値問題を考える。この作用素はいわゆる Sturm-Liouville 型作用素で、周期的境界条件の下で自己共役であ
り、その固有関数は完全系をなすことが知られている [41]。従って、この作用素の固有値が全て負なら臨界核は線
型安定になり、正の固有値を持つなら線形不安定になることが分かる。
まず、L が固有値 0 を持つことが次のように示せる: qs (ξ) は (4.51) 式の解なので、(4.51) 式の空間並進対称性
より、任意の ε ∈ R に対して qs (ξ + ε) も (4.51) 式の解である10 。よって、辺々引いて ε で割り、ε → 0 の極限を
取れば
∂ξ2 qs (ξ) + F (qs (ξ))qs (ξ) = 0
(4.61)
9 L → L + 0 では軌道は q に近づく。そこで、q(x) = q + δq(x) とおいて (4.51) 式を線形化すると、δq(x) = −F (q )δq(x) とな
c
2
2
2
る。F (q2 ) > 0 よりこの解は三角関数で書け、周期的境界条件を満たすので Lc F (q2 ) = 2π が成り立つ。
10 ξ + ε ∈
/ [0, L] である ξ に対しては周期的に拡張すればよい
46
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
となる。従って、qs (·) は L の固有値 0 の固有関数である11 。相図から分かるように、qs (·) は 2 つの零点を持つこ
とに注意。
次に、これを用いて L が正の固有値を持つことを示す。そのために、次の定理を用いる ([41]Chapter 8, section
3, Theorem 3.1)
12
定理
x ∈ [a, b] ⊂ R とし、p, q, r, y : [a, b] → R は連続かつ p > 0, r > 0 とする。更に、p は微分可能で p
は連続とする。また、y(·) は周期的境界条件
y(a) = y(b), y (a) = y (b)
を満たすとする。以上の条件の下で、λ を固有値とした Sturm-Liouville 問題
'
(
d
d
p(x) y(x) + q(x)y(x) = λr(x)y(x)
dx
dx
(4.62)
(4.63)
を考える。
このとき、(4.63) 式の固有値 λ の集合は無限に続く数列 (λn )n∈N で表され、
∞ > λ0 > λ1 ≥ λ2 · · ·
(4.64)
を満たす。また、固有値 λ0 に対応する固有関数 φ0 は [a, b] で 0 にならず、λ2i+1 , λ2i+2 に対応する固
有関数 φ2i+1 , φ2i+2 は半開区間 [a, b) でちょうど 2i + 2 個の零点を持つ。
この定理は p(ξ) = 1, r(ξ) = 1, q(ξ) = F (qs (ξ)) (∀ξ ∈ [0, L]) とおけば、今考えている問題に適用できる。上
でも指摘したように、L の固有値 0 の固有関数 qs (·) は [0, L) で 2 つの零点を持つので、λ1 = 0 または λ2 = 0 が
成り立つ。従って、上の定理の主張より、λ0 > 0 である。つまり、L は正の固有値を持ち、qs は線型不安定であ
ることが示された。
4.5
spatial モデルにおける準安定状態間の遷移
spatial モデルにおける準安定状態の遷移を考える。まず経路積分の手法で定式化を行い (4.5.1 項)、これに半古
典近似を適用する (4.5.2 項)。最後に、半古典近似で得られた運動方程式を解くための数値計算手法について述べ
る (4.5.3 項)。
4.5.1
経路積分表示
時刻 0 から十分長い時刻 T̃ までの系の振る舞いを考える13 。
(n1 (·), n2 (·), · · · , nM (·)) =: n(·)
11 q (ξ) は q (ξ) の ξ 微分を表す
s
s
12 引用元とは固有値 λ と q(·) の符号が異なる。
13 再び、時間・空間は無次元化前のスケールで考える
47
(4.65)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と定義すると、n(·) なる path の確率的重みは
P[n(·)] =
Dp exp(−S̃[n(·), p(·)])
S̃[n(·), p(·)] :=
H̃(n, p)
:=
,
dt p(t) · ∂t n(t) − H̃(n(t), p(t))
T̃
0
M
i=1
+σ
(4.66)
(4.67)
(epi − 1)W + (ni ) + (e−pi − 1)W − (ni )
M '
(e
pi −pi+1
i=1
(4.68)
ni+1 Np − ni
ni−1 Np − ni
− 1)
+ (epi −pi−1 − 1)
Np
Np
Np
Np
(
(4.69)
と書ける。
ここで、
q=
n
W ± (Np q)
, W± (q) =
Np
λ
(4.70)
λt
とおき、 N
を τ 、 λNT̃p を T と変換すると、q(·) なる path の確率的重みは
p
P[q(·)] =
Dp exp(−Np Sdis [q(·), p(·)])
Sdis [q(·), p(·)]
Hdis (q, p)
T
:=
0
M
:=
(4.71)
dτ [p(τ ) · ∂τ q(τ ) − Hdis (q(τ ), p(τ ))]
(epi − 1)W+ (qi ) + (e−pi − 1)W− (qi )
(4.72)
i=1
σ pi −pi+1
(e
− 1)qi+1 (1 − qi ) + (epi −pi−1 − 1)qi−1 (1 − qi )
λ i=1
M
+
(4.73)
で与えられる14 。
更に、σ λ の状況を考え、空間自由度を連続と見なす記述に移る。併せて
√il
D
= ξ (i ∈ {1, 2, · · · , M }) と空
間を無次元化すると（計算の詳細は付録 B.1 を参照のこと）
P[q(·, ·)] =
S[q(·, ·), p(·, ·)] :=
T
.
σ
S[q(·, ·), p(·, ·)]
Dp exp −Np
λ
(4.74)
L
h[q(·, ·), p(·, ·)](ξ, τ )
:=
dξ {p(ξ, τ )∂τ q(ξ, τ ) − h[q(·, ·), p(·, ·)](ξ, τ )}
(4.75)
0
/
h0 (q(ξ, τ ), p(ξ, τ )) − [∂ξ q(ξ, τ )][∂ξ p(ξ, τ )] − q(ξ, τ )[1 − q(ξ, τ )][∂ξ p(ξ, τ )]2 (4.76)
h0 (q, p)
:=
(ep − 1)W+ (q) + (e−p − 1)W− (q)
dτ
0
0
(4.77)
と書ける。
4.5.2
半古典近似
上記の経路積分表示に対して半古典近似を適用すると、得られる運動方程式は
∂τ q(ξ, τ )
=
ep(ξ,τ ) W+ (q(ξ, τ )) − e−p(ξ,τ ) W− (q(ξ, τ ))
0
/
(4.78)
+ ∂ξ2 q(ξ, τ ) − 2q(ξ, τ )[1 − q(ξ, τ )]∂ξ2 p(ξ, τ ) − 2[1 − 2q(ξ, τ )][∂ξ p(ξ, τ )][∂ξ q(ξ, τ )]
∂τ p(ξ, τ )
=
−(ep(ξ,τ ) − 1)W+ (q(ξ, τ )) − (e−p(ξ,τ ) − 1)W− (q(ξ, τ ))
0
/
− ∂ξ2 p(ξ, τ ) + [1 − 2q(ξ, τ )][∂ξ p(ξ, τ )]2
14 添え字の
dis は discerete の意味。
48
(4.79)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
1.0
uniform bounce solution
q3
1.0
0.9
0.9
0.8
0.8
q
q2
0.7
nonuniform bounce solution
q3
0.7
0.6
0.6
0.5
0.5
critical nucleus
0.4
0
2
4
L
6
8
10
0.4
0
5
10
L
15
20
図 4.8: bounce 解の概念図。左: 一様な bounce 解 A、右: 非一様な bounce 解 B
となる（導出は付録 B.2 参照）。この運動方程式に従う解は、Hamiltonian H :=
&L
0
dξh を保存する (導出は付録
B.3 を参照のこと)。
well-mixed モデルの解析と同様にして、T 1 なる時間 T に対して
• 時刻 0 で一様に q1 で、そのまま一度も「ほぼ一様に q3 」という状態に達せずに、時刻 T でもほぼ一様に q1
である確率
• 時刻 0 で一様に q3 で、そのまま一度も「ほぼ一様に q1 」という状態に達せずに、時刻 T でもほぼ一様に q3
である確率
を求めれば、これらの確率の T 依存性から各準安定状態の寿命が求まることになる (前者は「ほぼ一様に q1 」の
準安定状態の寿命、後者は「ほぼ一様に q3 」の準安定状態の寿命を与える。)。
まず、後者の確率 i.e.「系が最初ほぼ一様に q3 にあって、そのまま一度も『ほぼ一様に q1 』という状態に達せ
ずに時刻 T にもほぼ一様に q3 にある確率」を考えよう。
well-mixed モデルの場合はこの種の軌道の形は相図（図 4.5）を見れば一目瞭然であったのに対し、空間モデル
ではこのように容易に軌道を見つけることはできない。しかし、well-mixed モデルにおいて、考えるべき軌道が
安定固定点 q3 と不安定固定点 q2 を結ぶものであったことから類推すると、空間モデルで考えるべき軌道は、
• bounce 解 α: 「ほぼ一様に q3 」 → 「ほぼ一様に q2 」 → 「ほぼ一様に q3 」
• bounce 解 β: 「ほぼ一様に q3 」 → 臨界核 → 「ほぼ一様に q3 」
のいずれかであると考えられる15 。前者のような bounce 解を一様な bounce 解、後者のような bounce 解を非一
様な bounce 解と呼ぶことにする (図 4.8)。
このような軌道 (q(·, ·), p(·, ·)) を求めて (具体的な求め方は次項に述べる)、それに対応する作用 Sq3 を計算すれ
ば、求める確率 Pq3 (T ) は
Pq3 (T )
∼
Sq3
:=
*
+
√σ
exp −e−Np λ Sq3 T
T
L
dτ
dξ [p(ξ, τ )∂τ q(ξ, τ ) − h [q(·, ·), q(·, ·)] (ξ, τ )]
0
(4.80)
(4.81)
0
と表すことができる
従って、「ほぼ一様に q3 の状態」の寿命 τq3 は
τq3 ∼ exp Np
15 後者はもちろん、その系に臨界核が存在する場合に限られる。
49
.
σ
Sq
λ 3
(4.82)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と表せる。「ほぼ一様に q1 の状態」の寿命についても全く同様に計算できる。
ただし、bounce 解 α, β 両方が存在する場合は次のように考える: q3 の寿命を例にとり、bounce 解 α, β に対応
する作用をそれぞれ Sα , Sβ とし、Sα > Sβ とする (逆でも全く同様。)。(3.121) 式と同様にして (Kα , Kβ をそれ
ぞれ bounce 解 α, β の prefactor として)
Pq3 (T )
T
T
∞
n
*
+k *
+n−k
√
√
1 −Np σ
−Np σ
λ Sα
λ Sβ
K
K
C
dt
·
·
·
dt
e
e
1
n
α
β
n k
n!
0
0
n=0
(4.83)
k=0
=
=
∞
∞
+k +l
√
√σ
1 *
1*
−Np σ
Sα
λ
Kα e
Kβ e−Np λ Sβ T
T
k!
l!
k=0
l=0
,*
√σ
√σ + exp Kα e−Np λ Sα + Kβ e−Np λ Sβ T
(4.84)
(4.85)
となる。Kα = Kβ = −1 と近似すると、準安定状態 q3 の寿命 τ は、
√σ
√σ
1
∼ e−Np λ Sα + e−Np λ Sβ
τ
(4.86)
となる。今、Np 1 としているので、Sα > Sβ より、
τ ∼ eNp
√σ
λ Sβ
(4.87)
と近似できる。従って、作用の値が小さい方の bounce 解が効くことが分かる。
4.5.3
bounce 解の数値計算
さて、well-mixed モデルでは所望の軌道（bounce 解）を見つけるのは相図を用いれば容易であったが、空間モ
デルではこの手法は適用できない。ここでは、bounce 解を数値的に求める手法について述べる。
well-mixed モデルの軌道の類推から、求める軌道は bounce 解 α, β のいずれかと考えられることは既に述べた。
どちらも後半の「→」での遷移は p = 0（元のモデルでいうと期待値の dynamics）の軌道として与えられるので、
この部分は作用に寄与しないと考えられる。従って、
• bounce 解 α : 時刻 0 で一様に q3 → 時刻 T ( 1) で一様に q2
• bounce 解 β : 時刻 0 で一様に q3 → 時刻 T ( 1) で臨界核
のいずれかの境界条件を満たす軌道を求めれば十分であり、偏微分方程式 (4.78)(4.79) の 2 時刻境界値問題に帰着
する。この境界値問題を、代数方程式の Newton 法 [42] による解法を用いて数値的に解く（なお、このように微
分方程式の境界値問題を Newton 法を用いて解くことを緩和法と呼ぶ [42]）。
解くべき偏微分方程式は
∂τ q
=
∂τ p =
ep W+ (q) − e−p W− (q)
+ ∂ξ2 q − 2q(1 − q)∂ξ2 p − 2(1 − 2q)(∂ξ p)(∂ξ q)
1)W+ (q)
−p
−(e −
− (e −
2
− ∂ξ p + (1 − 2q)(∂ξ p)2
p
である。これを次のように離散化する:
• NT : 時間 [0, T ] の離散化ステップ数, NL : 空間 [0, L] の離散化ステップ数
• Δt := T /NT , Δx := L/NL
• τi := iΔt (i ∈ {0, 1, · · · , NT − 1}), ξj := jΔx (j ∈ {0, 1, · · · , NL − 1})
50
(4.88)
1)W− (q)
(4.89)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
• q(ξj , τi ) ↔ qi,j ,
p(ξj , τi ) ↔ pi,j
（なお、q0,j , qNT −1,j (j ∈ {0, 1, · · · , NL − 1}) は境界条件として与えられているとする）。この操作により、偏微
分方程式（4.78）
（4.79）は 2NL (NT − 1) 本の代数方程式に帰着する。更に変形すれば、差分方程式の解を求める
問題は、ある函数 f : R2NL (NT −1) → R2NL (NT −1) の零点を求める問題に帰着する（計算の詳細については付録 B.4
を参照のこと）。この高次元の代数方程式の問題に Newton 法を適用することにより、解を数値的に得ることがで
きる。以下に、数値的に得た一様/非一様な bounce 解の振る舞いを示す。bounce 解の模式図 4.8 も参照のこと。
一様な bounce 解
図 4.9 は q1 ,q3 の周りの一様な bounce 解である。境界条件は前者が「時刻 0 で一様に q1 、時刻 T ( 1) で一様
に q2 」で、後者が「時刻 0 で一様に q3 、時刻 T ( 1) で一様に q2 」である。
この軌道が well-mixed モデルの場合と全く同じ軌道の上を動いていることが、(4.78),(4.79) 式から容易に分か
る: q, p が ξ に依存しないことを仮定すると、well-mixed モデルの半古典近似の運動方程式 (4.35),(4.36) が得られ
るからである。また、このことは次のようにしても見て取れる: bounce 解 q, p の時間方向の断面図をとり、それ
らを (q, p) 平面にプロットする。得られた軌道を well-mixed モデルの相図と重ねると図 4.10 のようになる。
非一様な bounce 解
図 4.11 は q1 ,q3 の周りの非一様な bounce 解を表す。境界条件は前者が「時刻 0 で一様に q1 、時刻 T ( 1) で臨
界核」で、後者が「時刻 0 で一様に q3 、時刻 T ( 1) で臨界核」である。臨界核については 4.4.3 項の図 4.7 も参
照のこと。
51
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
図 4.9: 上:q1 の周りの一様な bounce 解。下:q3 の周りの一様な bounce 解。a = 1.0, b = 0.3, c = 0.6, d = 1.3, w =
0.8, μA = μB = 0.005, L = 10.5, NL = 21, T = 40.125, NT = 321 (q1 = 0.0031, q2 = 0.72, q3 = 0.99)
52
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
uniform bounce solutions
2.0
bounce solution around
q1
bounce solution around
q3
1.5
p
1.0
0.5
0.0
q2
q1
q3
0.5
1.0
0.0
0.2
0.4
q
0.6
0.8
1.0
図 4.10: well-mixed モデルと一様な bounce 解の軌道。a = 1.0, b = 0.3, c = 0.6, d = 1.3, w = 0.8, μA = μB =
0.005, L = 10.5, NL = 21, T = 40.125, NT = 321
53
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
図 4.11: 上: q1 の周りの非一様な bounce 解。下: q3 の周りの非一様な bounce 解。a = 1.0, b = 0.3, c = 0.6, d =
1.3, w = 0.8, μA = μB = 0.005, L = 20.5, NL = 41, T = 40.125, NT = 321 (q1 = 0.0031, q2 = 0.72, q3 = 0.99)
54
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
第 5 章 結果と考察
4 章に述べた解析の結果を示す。まず well-mixed モデルについての結果を述べた後、空間自由度を持つモデル
(spatial モデル) についての結果をまとめ、空間自由度の影響を述べる。特に、準安定状態の寿命に空間自由度が
どのような影響を及ぼすかについて、詳しく議論する。
5.1
well-mixed モデル
4.3 節で行った半古典近似を用いた解析の結果 (4.40) 式と、master 方程式からの数値計算で得られた結果を比較
する。具体的には、各準安定状態の寿命の、(1) 集団サイズ N 、 (2) 自然選択の強さ w、(3) 利得行列の値 a, b, c, d
への依存性を調べる。いずれの場合においても、半古典近似による結果が数値計算の結果と半定量的に一致する
ことを見る。
5.1.1
集団サイズ依存性
まず、各準安定状態の寿命の集団サイズ N への依存性を図 5.1 に示す。数値計算の結果と半古典近似の結果は、
定数倍だけ異なるが、N 依存性はほぼ同じである。特に、N が十分大きくなると「各準安定状態の寿命が集団サ
イズ N に対して指数的に増加する」という点においては一致している。(4.40) 式の寿命 τqi は、prefactor を省略
しなければ、
τqi |Ki |eN Sqi (i ∈ {1, 3})
(5.1)
で与えられる (ただし、Ki は qi を通る bounce 解の prefactor)。図 5.1 では、Ki ∼ −1 と近似した値を「半古典近
似の結果」として用いており、寿命の値自体が合わないことは prefactor の寄与を無視したことが原因である。ま
た、「寿命が N に対して指数的に増加する」という振る舞いは (5.1) 式の後半からも明らかである。なお、τq1 と
τq3 とでグラフの傾きが異なるのは、作用 Sq1 , Sq3 の値が異なるからである。
直感的に考えて、集団サイズ N の増加に伴って各準安定状態の寿命が長くなることは自然である: 一般に確率
過程では、N が大きくなると系の確率的なゆらぎの効果が相対的に小さくなり、より決定論的な描像に近づく (例
√
えば、van Kampen 展開では期待値は O(N ) であったのに対し、その周りの揺らぎは O( N ) であった。)。また、
決定論では各準安定状態の寿命は無限であった。このことより、N の増加に伴って準安定状態の寿命が長くなる
ことが理解される。
5.1.2
w 依存性
次に、自然選択の強さ w への依存性を見る（図 5.2）。数値計算の結果と半古典近似の結果は、値自体は一致し
ない。これは N 依存性の場合と同様に prefactor を無視したことによるものと考えられる。しかし、半古典近似・
数値計算、いずれの場合でも、各準安定状態の寿命は w の増加に伴って大きく増加することが見て取れる。
このことは直感的には次のように解釈できる: w が大きくなると、自然選択が確率的浮動に比べてより強く働
き、結果として決定論的描像に近くなる。
55
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
lifetime of q1
lifetime of q3
104
1017
1015
103
1013
1011
102
109
107
101
105
103
101 0
numerical
theoretical
10
20
30
N
40
50
60
70
100 0
numerical
theoretical
10
20
30
N
40
50
60
70
図 5.1: 各準安定状態の寿命の N 依存性（a = 1.0, b = 0.3, c = 0.6, d = 1.3, w = 0.8, μA = μB = 0.005）
lifetime of q1
10
10
10
numerical
theoretical
19
10
4
15
13
10
10
numerical
theoretical
17
10
10
lifetime of q3
5
3
11
10
10
10
10
10
9
10
2
7
5
10
1
3
1
0.1
0.2
0.3
0.4
0.5
w
0.6
0.7
0.8
0.9
1.0
10
0
0.1
0.2
0.3
0.4
0.5
w
0.6
0.7
0.8
0.9
1.0
図 5.2: 各準安定状態の寿命の w 依存性（N = 50, a = 1.0, b = 0.3, c = 0.6, d = 1.3, μA = μB = 0.005）
5.1.3
利得行列依存性
最後に、ゲームの利得行列の値への依存性を見る。4.1 節の W ± と ΠA , ΠB の定義を見ると分かるように、今考
えているゲームの dynamics は利得行列の要素 a, b, c, d にそのまま依存するのではなくそれらの差 a − c, d − b に
のみ依存する。従って、これら 2 つのパラメーターに対する依存性を見ればよい (図 5.3 1 )。図より、数値計算の
結果と半古典近似の結果で寿命の値自体は一致していないものの (これも prefactor を無視したことによるものと
考えられる)、payoff に対する依存性は再現していることが見て取れる。なお、この計算では μA = μB としている
ので、対称性より、q1 の寿命のグラフを直線 a − c = d − b について反転させれば q3 の寿命のグラフが得られる。
図 5.3 で示された寿命の利得行列の値への依存性は、定性的には次のように解釈できる：
• a − c が大きくなる。
⇐⇒ 戦略 A の個体を相手にしたとき、戦略 A の個体は戦略 B の個体よりもより有利になる。
=⇒ q3 の寿命は長くなり、q1 の寿命は短くなる。
1 （a − c, d − b < 0.2, > 0.8 の範囲では系が双安定でなくなる利得行列の組合せがあったので、準安定状態の寿命を考えることは意味をな
さず、従ってここには載せていない。）
56
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
0.8
lifetime of q1 (numerical)
27.5
lifetime of q1 (theoretical)
22.5
0
.0
5
2
0.8
0
00
.0
25.0
0.7
0
00
.
20
20
0
00
.
15
22.5
0.6
20.0
0.7
17.5
0.6
15.0
0.5
17.5
1
d−b
d−b
20.0
00
0
5.
0
00
.
10
0.5
12.5
10.0
15.0
0.4
0.4
12.5
00
.0
10
0.3
7.5
00
5.0
0.3
5.0
10.0
0.8
0.3
0.4
a−c
0.5
7.5
0.6
0.7
0.8
lifetime of q3 (numerical)
27.5
25.0
0.7
2.5
0.2
0.2
0.3
0.8
lifetime of q3 (theoretical)
0.4
a−c
0.5
0.6
0.7
0.8
5.0
0.2
0.2
24
0.7
21
22.5
0.6
18
0.6
15
00
.0
9
20
00
12.5
.0
12
00
15
10
0.4
20
0.3
25
10.0
6
0.2
0.2
0.3
0.4
a−c
0.5
7.5
0.6
0.7
0.8
25
.0
00
0.3
.00
0
0.5
15.0
.0
.0
15
10
0.4
00
00
.00
0
17.5
.0
0.5
d−b
d−b
20.0
0.2
0.2
0.3
0.4
a−c
0.5
3
0.6
0.7
0.8
図 5.3: 上:準安定状態 q1 の寿命の対数。下:準安定状態 q3 の寿命の対数。N = 50, w = 0.8, μA = μB = 0.005
• d − b が大きくなる。
⇐⇒ 戦略 B の個体を相手にしたとき、戦略 B の個体は戦略 A の個体よりもより有利になる。
=⇒ q1 の寿命は長くなり、q3 の寿命は短くなる。
5.2
spatial モデル
4.5 節の解析と、それに基づく数値計算から得られた結果を示す。
5.2.1
作用の集団サイズ依存性
4.5.3 項で得た bounce 解から求めた作用 S の L 依存性を図 5.4 に示す。数値的に求めた bounce 解から作用を
計算し、一様な bounce 解と非一様な bounce 解が両方存在する領域では、より作用の値が小さい方を用いた。結
果、5 ≤ L ≤ 12 の範囲では一様な bounce 解を用いて計算した値、L ≥ 14 の範囲では非一様な bounce 解を用い
て計算した値を採用した。その結果を図 5.4 に示す。
57
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
q1
12
q3
1.0
0.9
10
action
S
0.8
8
0.7
0.6
6
0.5
4
0.4
2
4
6
8
10
12
L
14
16
18
20
22
0.3
4
6
8
10
12
L
14
16
18
20
22
図 5.4: 作用の L 依存性。左: q1 の周りの bounce 解 右: q3 の周りの bounce 解。直線は well-mixed モデルから
予想される値。a = 1.0, b = 0.3, c = 0.6, d = 1.3, w = 0.8, μA = μB = 0.005
一様な bounce 解を用いて計算した作用は、どちらも L について線型に増加している。このことは、q, p が位置
ξ に依らないことと、(4.81) 式からも分かる。図 4.10 にも示したように、一様な bounce 解の軌道は well-mixed
モデルの軌道と同一のものとなるので、この直線の傾きは well-mixed な場合の作用で与えられる。図中の直線は
well-mixed モデルから予想される直線を表し、数値計算の結果とよく一致している。
また、非一様な bounce 解を用いて計算した作用は、q3 の周りの bounce 解については L に対してほぼ一定値を
取るが、q1 の周りの bounce 解については L に対して線型に増加する。このことは、直感的には次のように説明
できる: 図 4.7 でも見たように、臨界核は特徴的なスケールを持つ2 。従って、q3 から臨界核に遷移するにはその
スケールの範囲内でだけ q が変化すれば良いため S は L にほぼ依存しない。対して、q1 から臨界核に遷移するに
は全系で q が変化する必要があるため、S は L に依存して増える。
5.3
考察
図 5.1 より分かるように、well-mixed モデルにおいては、準安定状態 q1 , q3 の寿命はいずれも集団サイズ N が
増加すると指数的に増加する。このように「well-mixed な集団において系の何らかの状態の『寿命』が集団サイ
ズ N に指数的に依存する」という現象は、生態学や進化ゲーム理論のモデルにおいていくつもの例が知られてい
る: 生態系などの population のモデルについては集団が絶滅するまでの平均時間の研究が数多く行われ、[43] は
あるクラスのモデルについての厳密解を、[35] はより広いクラスのモデルについて WKB 法を用いた計算を行っ
た。また進化ゲームの文脈では、[28] は突然変異のない anti-coordination ゲームを考え、2 戦略が共存する状態の
寿命 (= 戦略が固定するまでの待ち時間) を求めた。この他にも様々な研究があるが、いずれの場合も系の準安定
な状態の寿命 τ は、集団サイズを N とすると
τ C exp(bN ) (b > 0)
(5.2)
と書けることが知られている [44, 45]。ここで、b > 0 は定数、C > 0 は N の函数だが N に指数的には依存しな
いものである。
2 ただし、(4.53)
式の V について V (q1 ) = V (q3 ) であれば、特徴的なスケールは存在しない。
58
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
䡚Lc
(a)L < Lc
䡚Lc
(b)L > Lc での q3 から q1 への遷移
(c)L > Lc での q1 から q3 への遷移
図 5.5: spatial モデルで準安定状態間の遷移が起きる様子。
一方、空間自由度があるモデル (spatial モデル) においては、準安定状態の寿命はこれとは異なる振る舞いを見
せる。(4.82) 式より、spatial モデルにおける準安定状態 q1 , q3 の寿命 τq1 ,τq3 はそれぞれ作用 Sq1 , Sq3 を用いて
τqi ∼ exp Np
.
σ
Sq (L)
λ i
(i ∈ {1, 3})
(5.3)
と表せる。ただし、L は系のサイズ (パッチの個数に比例)、Np はパッチ 1 つあたりの個体数であり、Np は L に
は依存しない。このことと図 5.4 より、準安定状態の寿命は L の増加に伴って次のように変化することが分かる:
L が増えると寿命は最初は指数的に増加する。臨界核が存在するスケール Lc を超えると q3 の寿命は一定値をと
るようになるが、q1 の寿命は (傾きを変えて) 直線的に増加する。つまり、L が小さいときは well-mixed モデルと
同様に振る舞うが、L が臨界核のスケールを超えると q1 と q3 とで全く異なる振る舞いが見られる。
Meerson ら [46] は 1 次元格子上の生態系モデルでの絶滅現象を考えて同様の解析を行ったが、そこでは準安定
状態が 1 つしかなく、集団が絶滅した状態が安定状態として存在した。そのため、本研究に見られるような 2 つ
の準安定状態の寿命の非対称性は、具体的には計算されていない。
このような非対称性の起源は、次のように理解することができる: 異なる準安定状態の間には「障壁」があり、
これを超えれば容易にもう片方の準安定状態に遷移することができる。具体的には、spatial モデルでは不安定定
常解が「障壁」の役割を果たす3 。というのも、これらを超えれば後は決定論的に (つまりは非常に高い確率で) も
う一方の準安定状態に遷移するからである。つまり、準安定状態から障壁を超えるのがどれほど困難か (稀か) が
準安定状態間の遷移の難しさ (遷移がどれほど稀か) を決めると言える。
• L < Lc では、不安定解 q2 が障壁の役目を果たし、片方の準安定状態からもう片方の準安定状態に遷移する
には系全体が変化して q2 を超える必要があった (図 5.5(a))。従って、遷移は系のサイズ L が増すほど困難
になる。つまり、各準安定状態の寿命は L が増すと大きく増える。
• L > Lc では臨界核も障壁の役割を果たし、臨界核を超えさえすればもう一方の準安定状態に移ることがで
きる。
– q3 から q1 の遷移のためには系のうち Lc 程度の部分だけが変化すればよく (図 5.5(b))、従って系のサ
イズ L が増しても遷移の難しさは変わらない。つまり、L が増しても準安定状態 q3 の寿命は増えない。
– q1 から q3 への遷移のためには系のうち L − Lc 程度の部分の変化を必要とし (図 5.5(c))、従って、遷
移は系のサイズ L 増すほど困難になる。つまり、準安定状態 q1 の寿命は L が増すと大きく増える。
ここでは特定のパラメーターの組を用いて計算したために上述のような「q3 の寿命が L に対して増えなくなる」
という結果が得られたが、q1 , q3 のうちどちらの準安定状態がこのような性質を示すかはパラメーターに依存する。
3 well-mixed
モデルでは不安定固定点が障壁となる。
59
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
しかしいずれの場合でも、「どちらかの準安定状態の寿命が L に対して増えなくなる」という結論は変わらない。
これは、一定のスケール Lc を持った臨界核が系に存在することによるものであり、更に言えば図 4.3 や図 4.6 の
ような相図が書けることによるものである。
また、このモデルでは p(Δπ) =
1+wΔπ
2
の pairwise comparison process を採用したが、p の函数形を変えた
pairwise comparison process や Moran process を採用したとしても (利得行列が条件 (4.2) 式を満たし、突然変異
率が十分小さければ) 図 4.3 や図 4.6 のような相図が描ける。従って系には特徴的なスケールを持った臨界核が存
在し、準安定状態の寿命に関してもこのモデルと同様な結論が得られると考えられる。
以上より、双安定な進化ゲーム一般に対して、空間自由度は
• 空間自由度がないと準安定状態の寿命は集団サイズ N に対して指数的に増える。
• 空間自由度が存在すると、一方の準安定状態の寿命は集団サイズ L に対して指数的に増えるが、もう一方の
準安定状態の寿命は (L が一定値を超えれば)L が大きくなっても増えず一定値を取る。
という形で、準安定状態の寿命の集団サイズ依存性を質的に変えることが分かった。
60
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
第 6 章 まとめと展望
本研究のまとめ
6.1
我々は双安定な進化ゲームの確率的 dynamics、具体的には準安定状態間の遷移を考え、準安定状態の寿命に空
間自由度の有無が及ぼす影響を調べた。解析にあたっては、確率過程として定式化されたモデルを経路積分によっ
て記述し、ここに半古典近似を適用して各準安定状態の寿命を計算した。
空間自由度がない場合については集団サイズ N や種々のパラメーター依存性を求め、master 方程式からの数値
計算結果と比較した。空間自由度がある場合については、系のサイズ L への寿命の依存性を求めた。また、空間
自由度がある場合とない場合の結果の比較から、次のことが分かった:
• 空間自由度がない場合、各準安定状態の寿命 τ は集団サイズ N に指数的に依存し、τ ∼ eaN (a > 0) となる。
• 空間自由度がある場合、系のサイズ L がある一定の値 Lc を超えるまでは τ は L の増加に伴って指数的に増
加する。L > Lc では一方の準安定状態では L に対して寿命が指数的に増えるが、もう一方の準安定状態の
寿命は L に依らない一定値を取る。
まとめると、本研究では
1. 確率過程における準安定状態間の遷移という稀な現象を取り扱うための、一般的な方法論を整理した。
2. この方法論を適用することにより、進化ゲームの確率的 dynamics に空間自由度が及ぼす影響の 1 つの側面
を解明することができた。
6.2
6.2.1
展望
他の系への適用
我々は準安定状態を持つ進化ゲームのモデルを考えたが、4 章の解析からも分かる通り、ここで用いた解析の手
法はモデルの詳細には依存しない: well-mixed モデルの相図が図 4.5 のような構造を持つ確率過程であれば、全て
同様に扱うことができる。
従って、本研究で用いた手法は、生態学や細胞生物学、化学反応など多種多様な系についての応用的研究にも
有用であると考えられる。
6.2.2
より一般の構造上での進化ゲーム
本研究では 1 次元格子上での進化ゲームを考えたが、より一般に 2,3 次元の格子も考えられる。また、一般の
network 上での進化ゲームを考えることもできる [21, 23, 25]。
特に、格子から一般の network に拡張することにより、通常の物理のモデルでは見られないような振る舞いが
発見されることが期待される: というのも、格子上では物理的に近い範囲で相互作用が起きるが、network 上では
非常に遠くの個体との相互作用も起きるからである。これは、通常の物理系には見られない性質である。
61
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
6.2.3
より精緻な近似手法による解析
本研究では解析にあたって半古典近似を用いたが、経路積分で記述されたモデルに対しては他の近似手法も適
用可能である: 系統的な摂動やくりこみなどの、物性物理学で発達してきた手法が適用できると考えられる [32]。
これらの手法を利用することにより、準安定状態の寿命だけでなく、他の動的な側面に空間自由度が与える影響
を探ることができると期待される。
62
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
付 録A
A.1
確率過程の近似手法
master 方程式の van Kampen 展開
ここでは、3 章で利用した、van Kampen の Ω 展開による linear noise 近似（[31] の X 章）の手法を解説する。
大雑把にいうと、van Kampen 展開は master 方程式をある大きな parameter に対して展開し、系統的な近似を行
う手法である。多くの場合 2 次までの展開を取り、求める確率分布を期待値の周りの Gauss 分布として評価する。
まず簡単な例として崩壊過程の展開について述べ、その後に一般の one-step process（n が一度の反応で ±1 し
か変化しない過程）について述べる。
A.1.1
崩壊過程
ここでは崩壊過程を例とする:
A → ∅ (with rate λ per particle)
(A.1)
時刻 t での A 粒子の数を n(t) とし、n(0) = N ( 1) とする1 。
この確率過程の master 方程式と初期条件は、時刻 t で粒子数が n である確率を P (n, t) とすると、
∂t P (n, t) = λ(n + 1)P (n + 1, t) − λnP (n, t), P (n, 0) = δn,N
(A.2)
となる。
van Kampen 展開では、粒子数を表す確率変数 n が
n(t) = N φ(t) +
√
N ξ(t)
(A.3)
と分解できると仮定する。N に比例した項 N φ は n の期待値を表し、ξ はその周りでの揺らぎを表す。n につい
ての初期条件より、
φ(0) = 1, ξ(0) = 0
(A.4)
とする。更に、ξ の確率密度関数を Π(ξ, t) とおく:
√
1
Π(ξ, t) = √ P (N φ(t) + N ξ, t)
N
2
このとき、
∂t P (n, t) =
√ ,
√
N ∂t Π(ξ, t) − N φ̇(t)∂ξ Π(ξ, t)
(A.5)
(A.6)
が成り立つ。
1 上述したとおり、van Kampen 展開はこのように何らかの大きな parameter の存在を利用する（e.g. 体積、粒子数など。van Kampen
はこの「大きな parameter
。
√ 」を Ω と表したので、この近似法は「van Kampen の Ω 展開」とも呼ばれる）
√
2 右辺全体にかかる
N は規格化条件から。van Kampen は Π の規格化は考慮せずこの N を省略しているが、後の計算は全く同様で
ある。
63
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
これらを元の master 方程式に代入し、N について展開する:
√
√
∂t Π(ξ, t) − N φ̇(t)∂ξ Π(ξ, t)
√
√
1
= λ(N φ + N + 1)Π ξ + √ , t − λ(N φ + N ξ)Π(ξ, t)
N
'
(
√
√
1
1 2
− 32
∂ξ Π + O(N ) − λ(N φ + N ξ)Π
= λ(N φ + N ξ + 1) Π + √ ∂ξ Π +
2N
N
√
φ 2
− 12
= λΠ + λ
N φ∂ξ Π + ξ∂ξ Π + ∂ξ Π + O(N )
2
√
1
λ 2
=
N λφ∂ξ Π + λ∂ξ (ξΠ) + φ∂ξ Π + O(N − 2 )
2
(A.7)
(A.8)
(A.9)
(A.10)
(A.11)
N に比例する項が相殺することを要請すると
dφ
= −λφ
dt
(A.12)
が得られる。これは、期待値の dynamics を定める常微分方程式である。初期条件より φ(0) = 1 なので、
φ(t) = e−λt
(A.13)
が得られる。
従って、Π についての方程式は
∂t Π = λ∂ξ (ξΠ) +
となり、この解は
λ −λt 2
e ∂ξ Π
2
(A.14)
(
'
ξ2
1
, s2 (t) := e−λt − e−2λt
exp − 2
Π(ξ, t) = )
2s (t)
2πs2 (t)
(A.15)
と求められる3 。
よって、元の分布 P の近似として
P (n, t) )
(
'
(n − N e−λt )2
, σ 2 (t) := N (e−λt − e−2λt )
exp −
2σ 2 (t)
2πσ 2 (t)
1
(A.16)
が得られる。
A.1.2
one-step process
一般の one-step process の master 方程式は次のように書ける。
∂t P (n, t) = W + (n − 1)P (n − 1, t) + W − (n + 1)P (n + 1, t) − [W + (n) + W − (n)]P (n, t)
(A.17)
系に特徴的なスケール N があり、N 1 を満たすと仮定する。更に、
W ± = O(N ), j N に対して W ± (n ∓ 1) = W ± (n) ∓ V ± (n) + O(
1
), V ± = O(1)
N
(A.18)
と書けると仮定する。このとき、
W± (x) :=
1 ±
W (N x), V± (x) := V ± (N x) (x = O(1))
N
(A.19)
と定義する。
ここで、
n(t) = N φ(t) +
3 一般論から、(A.14)
√
Nξ
式の解が Gaussian になることが言える（[31]）ので、平均と分散を求めればよい。
64
(A.20)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
と書けると仮定する。φ(t) は n の期待値を表し、ξ はその周りのゆらぎを表す。加えて、ξ の確率密度関数を Π と
おく：
√
1
Π(ξ, t) := √ P (N φ(t) + N ξ, t)
N
このとき、
∂t P (n, t) =
√ ,
√
N ∂t Π(ξ, t) − N φ̇(t)∂ξ Π(ξ, t)
(A.21)
(A.22)
と、
√
N ξ − 1) =
√
−
W (N φ + N ξ + 1) =
W + (N φ +
√
W + (N φ) + ( N ξ − 1)V + (N φ) + O(N −1/2 )
√
W − (N φ) + ( N ξ + 1)V − (N φ) + O(N −1/2 )
(A.23)
(A.24)
が成り立つ。
これらを元の master 方程式に代入する4 ：
√
∂t Π(ξ, t) − N φ̇(t)∂ξ Π(ξ, t)
(A.25)
-
,
√
1
1 2
∂ξ Π + O(N −3/2 )
= W + (N φ) + N ξV + (N φ) − V + (N φ) + O(N −1/2 ) Π − √ ∂ξ Π +
2N
N
-
,
√
1
1 2
−
+
−
−1/2
−3/2
∂ Π + O(N
) Π + √ ∂ξ Π +
)
+ W (N φ) + N ξV (N φ) + V (N φ) + O(N
2N ξ
N
,
√
√
(A.26)
− W + (N φ) + N ξV + (N φ) + W − (N φ) + N ξV + (N φ) + O(N −1/2 ) Π
,
√
1
1
W + (N φ) + N ξV + (N φ) ∂ξ Π +
= −√
W + (N φ)∂ξ2 Π − V + (N φ)Π
2N
N
√
1
1 , −
W − (N φ)∂ξ2 Π + V − (N φ)Π + O(N −1/2 )
(A.27)
W (N φ) + N ξV − (N φ) ∂ξ Π +
+√
2N
N
'
(
√
1
1
= − N W+ (φ) + √ V+ (φ) ∂ξ Π + W+ (φ)∂ξ2 Π − V+ (φ)Π
2
N
'
(
√
1
1
+ N W− (φ) + √ V− (φ) ∂ξ Π + W− (φ)∂ξ2 Π + V− (φ)Π + O(N −1/2 )
(A.28)
2
N
√
1
= − N [W+ (φ) − W− (φ)] ∂ξ Π − [V+ (φ) − V− (φ)] ∂ξ (ξΠ) + [W+ (φ) + W− (φ)] ∂ξ2 Π + O(N −1/2()A.29)
2
√
N に比例する項が相殺することを要請すると、
φ̇(t) = W+ (φ(t)) − W− (φ(t))
(A.30)
Π の満たす方程式は、O(1) まで取って、
∂t Π(ξ, t) = − [V+ (φ) − V− (φ)] ∂ξ (ξΠ) +
1
[W+ (φ) + W− (φ)] ∂ξ2 Π
2
(A.31)
が得られる。
A.2
半古典近似の詳細
3 章で半古典近似について述べたときは、
p(nT , T |n0 , 0) K e−N S[qcl (·),pcl (·)]
(A.32)
と書き、K の具体的な評価は行わなかった。この appendix では簡単な例で K を具体的に評価し、その寄与が小
さいことを見る。
4 引数は適宜省略する。
65
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
確率過程として、Poisson 過程を考える：
∅ → A (with rate λ)
5
(A.33)
初期条件は n(0) = N ( 1) とする。λ = O(N ) と仮定し、λ := λ/N とおく。
これを 3 章で述べたように経路積分で記述すると、
p(nT , T |N, 0) =
dt[p(t)∂t n(t) − H̃(n(t), p(t))]
S̃[n(·), p(·)] :=
H̃(n, p)
DnDp e−S̃[n(·),p(·)]
λ(ep − 1)
:=
DqDp e−N S[q(·),p(·)]
T
dt[p(t)∂t q(t) − H(q(t), p(t))]
0
λ (ep − 1)
:=
(A.38)
q(0)=1,n(T )=qT
:=
H(q, p)
(A.36)
(A.37)
従ってこれは、q := n/N 、qT := nT /N とおけば、
p(nT , T |N, 0) =
S[q(·), p(·)]
(A.35)
n(0)=N,n(T )=nT
(A.39)
(A.40)
と書けて半古典近似を適用できる。
古典軌道の満たす運動方程式は
∂t q = λ ep , ∂t p = 0, q(0) = 1, q(T ) = nT /N =: qT
で、この解は
qcl (t) = (qT − 1)
t
+ 1, pcl (t) = log
T
qT − 1
λ T
(A.41)
(A.42)
と書ける (ただし、qT ≥ 1 に限る。)。これらを代入すると、
qT − 1
qT − 1
qT − 1
log
−
λ
−
1
T
λ T
λ T
0
qT − 1
(qT − 1) log
− (qT − 1) + λ T
λ T
S[qcl (·), pcl (·)] =
=
T
dt
(A.43)
(A.44)
なので、
e
−N S[qcl (·),pcl (·)]
=e
−λT
(λT )
nT −N
e
nT − N
nT −N
(A.45)
となる。ここまでの寄与だけを取れば、
p(nT , T |N, 0) ∼ e
−λT
(λT )
nT −N
e
nT − N
nT −N
(nT ≥ N )
(A.46)
次に、古典解周りの揺らぎを（2 次まで）求める。
q(t) = qcl (t) + r(t), p(t) = pcl (t) + is(t)
5 master
(A.47)
方程式は、時刻 t で A の個数が n である確率を P (n, t) とすると、
∂t P (n, t) = λP (n − 1, t) − λP (n, t)
となる。
66
(A.34)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
とおいて S に代入し (r, s ∈ R、r(0) = r(T ) = 0 を仮定)、r, s の 2 次まで取ると、
S[q, p]
=
=
=:
dt (pcl + is) ∂t (qcl + r) − λ (epcl +is − 1)
0
(
T '
1
dt pcl ∂t qcl − λ (epcl − 1) + pcl ∂t r + is (∂t qcl − λ epcl ) + is∂t r + λ epcl · s2
2
0
qcl , pcl は運動方程式の解なので、
(
T
T '
λ pcl 2
pcl
dt [pcl ∂t qcl − λ (e − 1)] +
dt is∂t r + e s
2
0
0
S[qcl , pcl ] + S [r, is]
T
これより、
p(nT , T |N, 0) e−N S[qcl ,pcl ]
DrD(is)e−N S
[r,is]
(A.48)
(A.49)
(A.50)
(A.51)
(A.52)
と表される。この式の後半の経路積分を評価する (経路積分の離散の仕方については (3.44) 式を参照のこと)。
DrD(is)e−N S [r,is]
(A.53)
M
−1
M
M
−1 idsk
λ
rl+1 − rl
exp N Δt
sl+1
=
lim
− epcl,l+1 s2l+1 − i
(A.54)
N drj
M →∞
2πi
2
Δt
j=1
=
=
lim
M →∞
lim
=
=
k=1
N drj
j=1
M
−1
M →∞
M
k=1
N drj
j=1
M
1
.
k=1
l=0
'
λΔt pcl,k 2
dsk
exp −
e
sk − iN (rk − rk−1 )sk
2π
2
'
N 2 (rk − rk−1 )2
exp
−
2πλepcl,k Δt
2λepcl,k Δt
1
1
、 xj
2λepcl,k Δt
M
−1
M
(2
(A.55)
(
:= N rj とおくと、
.
Ak
lim
dxj
exp −Ak (xk − xk−1 )2
M →∞
π
j=1
k=1
M
M .
Ak
exp −Ak (xk − xk−1 )2
dxj δ(xM )
lim
M →∞
π
j=1
k=1
M
M .
Ak
1
iwxM
dwe
dxj
lim
exp −Ak (xk − xk−1 )2
M →∞
2π
π
j=1
Ak :=
=
M
−1
(A.56)
(A.57)
(A.58)
(A.59)
k=1
=
=
=
=
yk := xk − xk−1 (k ∈ {1, 2, · · · , M }) とおくと、
M
M .
M
%
$
1
Ak
dweiw l=1 yl
lim
dyj
exp −Ak yk2
M →∞
2π
π
j=1
k=1
.
M
$
%
1
Aj
dw
dyj exp −Aj yj2 + iwyj
lim
M →∞ 2π
π
j=1
w2
dw
exp −
4Aj
j=1
⎛
⎞
M
1
λΔt 2
lim
dw exp ⎝−
epcl,j ⎠
w
M →∞ 2π
2
j=1
1
lim
M →∞ 2π
⎡
=
lim ⎣2πλΔt
M →∞
M
M
(A.60)
(A.61)
(A.62)
(A.63)
⎤− 12
epcl,j ⎦
(A.64)
j=1
67
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
M → ∞ では
Δt
M
e
pcl,j
T
→
dte
pcl
0
j=1
n−N
=
(∵ pcl = log
λ
n−N
λT
(A.65)
なので、
p(nT , T |N, 0)
=
1
e−N S[qcl ,pcl ] )
2π(nT − N )
nT −N
1
e
)
e−λT (λT )nT −N
nT − N
2π(nT − N )
(A.66)
(A.67)
nT を n に、T を t に書き換えれば
p(n, t|N, 0) e
−λt
(λt)
n−N
e
n−N
n−N
1
)
2π(n − N )
(A.68)
となる。
一方、Poisson 過程の厳密解は
p(n, t|N, 0) =
⎧
⎨e−λt (λt)n−N
(n ≥ N )
⎩0
(otherwise)
(n−N )!
(A.69)
で与えられる6 。
従って、n > N において
厳密解:
prefactor なしの半古典近似:
prefactor ありの半古典近似:
(λt)n−N
(n − N )!
p(n, t|N, 0) ∼ e−λt (λt)n−N
p(n, t|N, 0) = e−λt
p(n, t|N, 0) e
−λt
(λt)
n−N
(A.70)
e
n−N
e
n−N
n−N
(A.71)
n−N
1
)
2π(n − N )
となる。prefactor なし/ありの半古典近似の結果は、厳密解をそれぞれ Stirling の近似式
* n +n
* n +n
√
n! ∼
, n! 2πn
e
e
(A.72)
(A.73)
で近似したものであることが分かる。
これらの結果 (と van Kampen 展開による結果7 ) を比較したグラフを図 A.1 に示す。prefactor があると近似の
精度は上がるが、精度が大きく上がるわけではないことが見て取れる。
n
n
:=
n∈Z x P (n, t) とおき、master 方程式の両辺に x をかけて n についての和を取ると、F についての偏微分方程式
∂t F (x, t) = λ(x − 1)F (x, t) が得られる。この方程式と初期条件 F (x, 0) = xN より、F (x, t) = xN eλ(x−1)t と書ける。この F を x につ
いて展開すれば、解 P (n, t) が得られる。
7 van Kampen 展開による結果は、A.1.2 項の結果より
1
(n − N − λt)2
p(n, t|N, 0) √
exp −
2λt
2πλt
6 F (x, t)
となる。
68
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
Poisson process
0
exact
without prefactor
with prefactor
van Kampen
expectation value
10
log p
20
30
40
50
6020
25
30
35
n
40
図 A.1: λ = 1.0, N = 20, t = 2.0
69
45
50
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
付 録B
空間モデルの Hamiltonian と半古典近似
Hamiltonian の連続化
B.1
(4.72),(4.73) 式の作用と Hamiltonian は
T
Sdis [q(·), p(·)] :=
Hdis (q, p)
:=
0
M
dτ [p(τ ) · ∂τ q(τ ) − Hdis (q(τ ), p(τ ))]
(B.1)
(epi − 1)W+ (qi ) + (e−pi − 1)W− (qi )
(B.2)
i=1
σ pi −pi+1
(e
− 1)qi+1 (1 − qi ) + (epi −pi−1 − 1)qi−1 (1 − qi )
λ i=1
M
+
(B.3)
であった。これらの表式に対して空間を「連続化」したとき、(4.75),(4.76) 式の作用と Hamiltonian が得られるこ
とを示す。
σ λ の状況を考えると、qi (·) は空間座標 i の関数として滑らかに変化すると考えられる。従って、
q(x, τ ) := q xl (τ ) (x ∈ {l, 2l, 3l, · · · , M l})
(B.4)
を x ∈ [0, L] に拡張して連続的な記述に移ることができる。以下、この連続的な記述での作用と Hamiltonian の表
式を導出する。
このとき、
p(τ ) · ∂τ q(τ )
M
(epi − 1)W+ (qi ) + (e−pi − 1)W− (qi )
i=1
−→
−→
1
l
1
l
L
dx p(x, τ )∂τ q(x, τ )
(B.5)
0
L
,
dx (ep(x) − 1)W+ (q(x)) + (e−p(x) − 1)W− (q(x))
(B.6)
0
Hamiltonian Hdis の後半部分は以下のように求められる1 ：
まず、
pi±1 − pi
=
∴ epi −pi±1 − 1
=
1
p(x ± l) − p(x) = ±(∂x p)l + (∂x2 p)l2 + O(l3 )
2
2
1
2
2
(∂x p) − ∂x p l + O(l3 )
∓(∂x p)l +
2
(B.7)
(B.8)
また、
qi±1 (1 − qi ) = q(x ± l)[1 − q(x)] = q(1 − q) ± (1 − q)(∂x q)l + O(l2 )
従って、
D :=
1 以下の計算では時間変数を省略。また、q, p
σ 2
l
λ
の引数は特に明示されなければ x であるとする。
70
(B.9)
(B.10)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
とおくと、
σ pi −pi+1
− 1)qi+1 (1 − qi ) + (epi −pi−1 − 1)qi−1 (1 − qi )
(e
λ i=1
M
M
0
σ 2/
−2(1 − q)(∂x q)(∂x p) + q(1 − q) (∂x p)2 − ∂x2 p + O(l3 )
l
λ i=1
0
D L /
−→
dx q(1 − q) (∂x p)2 − ∂x2 p − 2(1 − q)(∂x q)(∂x p)
l 0
D L dx q(1 − q)(∂x p)2 − (∂x q)(∂x p)
=
l
0 L
L
2
∵−
dx q(1 − q)(∂x p) =
dx (1 − 2q)(∂x q)(∂x p)
=
0
(B.11)
(B.12)
(B.13)
(B.14)
0
を得る。
従って、全体としては
P[q(·, ·)]
Stmp [q(·, ·), p(·, ·)] :=
htmp [q(·, ·), p(·, ·)](x, τ ) :=
h0 (q, p)
Dp exp (−Np Stmp [q(·, ·), p(·, ·)])
=
:=
L̃
1 T
dτ
dx {p(x, τ )∂τ q(x, τ ) − htmp [q(·, ·), p(·, ·)](x, τ )}
l 0
0
h0 (q(x, τ ), p(x, τ ))
0
/
−D [∂x q(x, τ )][∂x p(x, τ )] − q(x, τ )[1 − q(x, τ )][∂x p(x, τ )]2
(e − 1)W+ (q) + (e
p
−p
− 1)W− (q)
0
h[q(·, ·), p(·, ·)](ξ, τ ) :=
h0 (q, p)
B.2
(B.16)
(B.17)
(B.18)
√
√
D で無次元化する: ξ := x/ D
.
σ
S[q(·, ·), p(·, ·)]
P[q(·, ·)] =
Dp exp −Np
λ
T
L
S[q(·, ·), p(·, ·)] :=
dτ
dξ {p(ξ, τ )∂τ q(ξ, τ ) − h[q(·, ·), p(·, ·)](ξ, τ )}
更に、空間次元を
(B.15)
:=
(B.19)
(B.20)
0
h0 (q(ξ, τ ), p(ξ, τ ))
0
/
− [∂ξ q(ξ, τ )][∂ξ p(ξ, τ )] − q(ξ, τ )[1 − q(ξ, τ )][∂ξ p(ξ, τ )]2
(e − 1)W+ (q) + (e
p
−p
− 1)W− (q)
(B.21)
(B.22)
運動方程式の導出
作用 S は (B.20),(B.21) 式 (または (4.75),(4.76) 式) のように
dτ
0
h[q(·, ·), p(·, ·)](ξ, τ ) :=
h0 (q, p)
:=
T
S[q(·, ·), p(·, ·)] :=
L
dξ {p(ξ, τ )∂τ q(ξ, τ ) − h[q(·, ·), p(·, ·)](ξ, τ )}
(B.23)
0
h0 (q(ξ, τ ), p(ξ, τ ))
0
/
− [∂ξ q(ξ, τ )][∂ξ p(ξ, τ )] − q(ξ, τ )[1 − q(ξ, τ )][∂ξ p(ξ, τ )]2
(e − 1)W+ (q) + (e
p
−p
− 1)W− (q)
(B.24)
(B.25)
と書ける。ここに半古典近似を適用した時に、q, p についての偏微分方程式 (4.78),(4.79) 式が得られることを示す。
q を q + δq に、p を p + δp に変化させたとする。δq, δp は ξ について周期的境界条件を満たし、
∀ξ ∈ [0, L] [δq(ξ, 0) = δq(ξ, T ) = 0, δp(ξ, 0) = δp(ξ, T ) = 0]
71
(B.26)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
を仮定する。δq, δp はこれらの条件を満たし、かつ十分良い解析的性質を持つ範囲で任意とする。このときの S
の変化は O(δq), O(δp) まで取って（ ˙ は τ による微分を、 は ξ による微分を表す）
=
S[q + δq, p + δp] − S[q, p]
T
L
dτ
dξ [(p + δp)(q̇ + δ q̇) − pq − h0 (q + δq, p + δp) + h0 (q, p)]
0
0
T
0
dξ [pδ q̇ + q̇δp − (∂q h0 )δq − (∂p h0 )δp]
dτ
0
0
T
0
T
0
T
T
0
L
dξ −q δp − p δq + 2q(1 − q)p δp + 2(1 − 2q)q p δp − (1 − 2q)(p )2 δq
L
dτ
0
(B.28)
dξ [−ṗδq + q̇δp − (∂q h0 )δq − (∂p h0 )δp]
0
dτ
0
dξ q δp + p δq − 2q(1 − q)p δp − (1 − 2q)(p )2 δq
0
L
dτ
=
L
dτ
+
=
dξ (q + δq )(p + δp ) − (q + δq)(1 − q − δq)(p + δp )2 − q p + q(1 − q)(p )2 (B.27)
0
L
T
=
L
dτ
+
dξ
/
(B.29)
−ṗ − ∂q h0 − p − (1 − 2q)(p )2 δq
0
[q̇ − ∂p h0 − q + 2q(1 − q)p + 2(1 − 2q)q p ] δp}
(B.30)
ここに停留条件 S[q + δq, p + δp] − S[q, p] = 0 を課すと、[· · · ] の部分がそれぞれ 0 に等しいことが言え、古典
軌道が満たすべき Hamilton の運動方程式
∂τ q
=
∂p h0 (q, p) + ∂ξ2 q − 2q(1 − q)∂ξ2 p − 2(1 − 2q)](∂ξ p)(∂ξ q)
ep W+ (q) − e−p W− (q)
+ ∂ξ2 q − 2q(1 − q)∂ξ2 p − 2(1 − 2q)](∂ξ p)(∂ξ q)
0
/
= −∂q h0 (q, p) − ∂ξ2 p + (1 − 2q)(∂ξ p)2
=
∂τ p
=
−(ep − 1)W+ (q) − (e−p − 1)W− (q))
0
/
− ∂ξ2 p + (1 − 2q)(∂ξ p)2
(B.31)
(B.32)
が従う。
B.3
Hamiltonian の保存
H [q(·, ·), p(·, ·), τ ] =
L
/
0
dξ h0 (q(ξ, τ ), p(ξ, τ )) − (∂ξ q(ξ, τ ))(∂ξ p(ξ, τ )) − q(ξ, τ )(1 − q(ξ, τ ))(∂ξ p(ξ, τ ))2
0
(B.33)
とおく。
Hamilton の運動方程式 (B.31),(B.32) 式に従う q(·, ·), p(·, ·) を 1 つ固定したとき、
Eq,p (τ ) := H [q(·, ·), p(·, ·), τ ]
が時間 τ に依らず一定値を取ることを示す。
72
(B.34)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
周期的境界条件に注意して部分積分を行うと、
d
Eq,p (τ )
dτ
L
=
dξ[(∂q h0 )(∂τ q) + (∂p h0 )(∂τ p) − (∂τ ∂ξ q)(∂ξ p) − (∂ξ q)(∂τ ∂ξ p)
0
+(1 − 2q)(∂τ q)(∂ξ p)2 + 2q(1 − q)(∂ξ p)(∂τ ∂ξ p)]
L
=
0
dξ[(∂q h0 )(∂τ q) + (∂p h0 )(∂τ p) + (∂τ q)(∂ξ2 p) + (∂ξ2 q)(∂τ p)
+(1 − 2q)(∂τ q)(∂ξ p)2
−2(1 − 2q)(∂ξ q)(∂ξ p)(∂τ p) − 2q(1 − q)(∂ξ2 p)(∂τ p)
L
=
0
(∵ 部分積分)
/
dξ (∂τ q) (∂q h0 ) + (∂ξ2 p)(1 − 2q)(∂ξ p)2
0
+(∂τ p) (∂p h0 ) + (∂ξ2 q) − 2(1 − 2q)(∂ξ q)(∂ξ p) − 2q(1 − q)(∂ξ2 p)
L
[(∂τ q)(−∂τ p) + (∂τ p)(∂τ q)]
=
(∵ 運動方程式)
0
= 0
(B.35)
従って、古典軌道 q(·, ·), p(·, ·) に対して、Eq,p (τ ) は τ ∈ [0, T ] で一定値を取る。
B.4
緩和法による数値解法
ここでは、4.5.3 項で用いた偏微分方程式の数値解法、特に偏微分方程式の境界値問題を Newton 法で解ける形
に帰着させるプロセスについて述べる。
解くべき偏微分方程式は
∂τ q
=
∂τ p
=
F (q, p) + ∂ξ2 q − 2q(1 − q)∂ξ2 p − 2(1 − 2q)(∂ξ p)(∂ξ q)
G(q, p) − ∂ξ2 p + (1 − 2q)(∂ξ p)2
F (q, p)
:=
e W+ (q) − e
G(q, p)
:=
−(ep − 1)W+ (q) − (e−p − 1)W− (q)
p
−p
W− (q),
であった。
これを次のように離散化する:
• NT : 時間 [0, T ] の離散化ステップ数, NL : 空間 [0, L] の離散化ステップ数
• Δt := T /NT , Δx := L/NL
• τi := iΔt (i ∈ {0, 1, · · · , NT − 1}), ξj := jΔx (j ∈ {0, 1, · · · , NL − 1})
• q(ξj , τi ) ↔ qi,j ,
p(ξj , τi ) ↔ pi,j
なお、q0,j , qNT −1,j (j ∈ {0, 1, · · · , NL − 1}) は境界条件として与えられているとする。
このとき、偏微分方程式 (B.36),(B.37) は次のように離散化される2 :
2 どのように離散化するかは任意性があるが、ここでは対称性を重視した離散化を行った。
73
(B.36)
(B.37)
(B.38)
(B.39)
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
i ∈ {1, 2, · · · , NT − 2} , j ∈ {0, 1, · · · , NL − 1} に対しては、
qi+1,j − qi−1,j
2Δt
=
+
pi+1,j − pi−1,j
2Δt
=
−
F (qi,j , pi,j )
'
qi,j+1 − 2qi,j + qi,j−1
qi,j+1 − qi,j−1 pi,j+1 − pi,j−1
− 2(1 − 2qi,j )
·
(Δx)2
2Δx
2Δx
(
pi,j+1 − 2pi,j + pi,j−1
−2qi,j (1 − qi,j )
(Δx)2
G(qi,j , pi,j )
'
(
pi,j+1 − 2pi,j + pi,j−1
(pi,j+1 − pi,j−1 )2
+
(1
−
2q
)
i,j
(Δx)2
4(Δx)2
(B.40)
(B.41)
i = 0, j ∈ {0, 1, · · · , NL − 1} に対しては、
p1,j − p0,j
Δt
=
−
G(q0,j , p0,j )
'
(
p0,j+1 − 2p0,j + p0,j−1
(p0,j+1 − p0,j−1 )2
+
(1
−
2q
)
0,j
(Δx)2
4(Δx)2
(B.42)
i = NT − 1, j ∈ {0, 1, · · · , NL − 1} に対しては、
pNT −1,j − pNT −2,j
Δt
=
−
G(qNT −1,j , pNT −1,j )
'
pNT −1,j+1 − 2pNT −1,j + pNT −1,j−1
(Δx)2
(
(pNT −1,j+1 − pNT −1,j−1 )2
+(1 − 2qNT −1,j )
4(Δx)2
(B.43)
この離散化した差分方程式を満たす q1,0 , q1,1 , · · · , qNT −2,NL −1 , p0,0 , p0,1 , · · · pNT −1,NL −1 が求まれば、元の偏微分
方程式の 2 時刻境界値問題が（近似的に）解けたことになる。
(B.40)-(B.43) 式を変形したものを用いて
Ai,j
:=
−
Bi,j
:=
+
B0,j
:=
+
BNT −1,j
:=
+
qi+1,j − qi−1,j − 2Δt F (qi,j , pi,j )
'
1
Δt
− qi,j (qi,j+1 − qi,j−1 )(pi,j+1 − pi,j−1 )
(qi,j+1 − 2qi,j + qi,j−1 ) −
2
(Δx)2
2
− 2qi,j (1 − qi,j )(pi,j+1 − 2pi,j + pi,j−1 )
(i ∈ {1, 2, · · · , NT − 2} , j ∈ {0, 1, · · · , NL − 1})
(B.44)
pi+1,j − pi−1,j − 2Δt G(qi,j , pi,j )
'
(
1
Δt
2
(1
−
2q
(p
−
2p
+
p
)
+
)(p
−
p
)
2
i,j+1
i,j
i,j−1
i,j
i,j+1
i,j−1
(Δx)2
4
(i ∈ {1, 2, · · · , NT − 2} , j ∈ {0, 1, · · · , NL − 1})
(B.45)
p1,j − p0,j − Δt G(q0,j , p0,j )
'
(
1
Δt
2
(1
−
2q
(p
−
2p
+
p
)
+
)(p
−
p
)
0,j+1
0,j
0,j−1
0,j
0,j+1
0,j−1
(Δx)2
4
(j ∈ {0, 1, · · · , NL − 1})
(B.46)
pNT −1,j − pNT −2,j − 2Δt G(qNT −1,j , pNT −1,j )
'
(
1
Δt
2
(1
−
2q
(p
−
2p
+
p
)
+
)(p
−
p
)
NT −1,j+1
NT −1,j
NT −1,j−1
NT −1,j
NT −1,j+1
NT −1,j−1
(Δx)2
4
(B.47)
(j ∈ {0, 1, · · · , NL − 1})
とおき、
f : (q1,0 , q1,1 , · · · , qNT −2,NL −1 , p0,0 , p0,1 , · · · pNT −1,NL −1 ) → (A1,0 , A1,1 , · · · , ANT −2,NL −1 , B0,0 , B0,1 , · · · BNT −1,NL −1 )
(B.48)
74
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
とおくと、差分方程式 (B.40)-(B.43) の解を求める問題は f の零点を求める問題に帰着される。
75
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
参考文献
[1] John von Neumann, Oskar Morgenstern. Theory of Games and Economic Behavior. Princeton University
Press (1944).
[2] Gyorgy Szabo, Gabor Fath. Evolutionary games on graphs. Physics Reports, 446, 97-216 (2007).
[3] John Maynard Smith, George Robert Price. The Logic of Animal Conflict. Nature 246, 15-18 (1973).
[4] Martin A. Nowak, Karl Sigmund. Evolutionary Dynamics of Biological Games. Science 303,793-799 (2004).
[5] Arne Traulsen, Christoph Hauert. Stochastic evolutionary game dynamics (in Reviews of Nonlinear Dynamics and Complexity Vol. II edited by H.G.Schuster. Wiley-VCH.(2009))
[6] Martin A. Nowak. Evolutionary Dynamics: Exploring the Equations of Life. Belknap Press (2006).
[7] William H. Sandholm. Population Games and Evolutionary Dynamics. The MIT Press (2010). ISBN:
978-0-262-19587-4
[8] Peter D. Taylor, Leo B. Jonker. Evolutionarily stable strategies and Game Dynamics. Mathematical Biosciences 40, 145-156 (1978).
[9] Josef Hofbauer, Karl Sigmund. Evolutionary Games and Population Dynamics. Cambridge University Press
(1998). ISBN: 0-521-62365-0 (hardback), 0-521-62570-X (paperback).
[10] Patrick Alfred Pierce Moran. The statistical processes of evolutionary theory. Clarendon Press (1962).
[11] Martin A. Nowak, Akira Sasaki, Christine Taylor, Drew Fudenberg. Emergence of cooperation and evolutionary stability in finite populataions. Nature 428,646-650 (2004).
[12] Christine Taylor, Drew Fudenberg, Akira Sasaki, Martin A. Nowak. Evolutionary game dynamics in finite
populations. Bulletin of Mathematical Biology 66,1621-1644 (2004).
[13] Arne Traulsen, Jens Christian Claussen, Christoph Hauert. Coevolutionary Dynamics: From Finite to
Infinite Populations. Physical Review Letters 95, 238701 (2005).
[14] Arne Traulsen, Jens Christian Claussen, Christoph Hauert. Coevolutionary dynamics in large, but finite
populations. Physical Review E 74, 011901 (2006).
[15] Tibor Antal, MartinA.Nowak, ArneTraulsen. Strategy abundance in 2 × 2 games for arbitrary mutation
rates. Journal of Theoretical Biology 257, 340-344 (2009).
[16] Tibor Antal, Arne Traulsen, Hisashi Ohtsuki, Corina E. Tarnita, Martin A. Nowak. Mutation-selection
equilibrium in games with multiple strategies. Journal of Theoretical Biology 258, 614-622 (2009).
[17] Hisashi Ohtsuki, Martin A. Nowak, Jorge M. Pacheco. Breaking the symmetry between interaction and
replacement in evolutionary dynamics on graphs. Physical Review Letters 98, 108106 (2007).
76
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
[18] Gyorgy Szabo, Csaba Toke. Evolutionary prisoner’s dilemma game on a square lattice. Physical Review E
58, 69-73 (1998).
[19] Christoph Hauert, Michael Doebeli. Spatial structure often inhibits the evolution of cooperation in the
snowdrift game. Nature 428, 643-646 (2004).
[20] Benjamin Kerr, Margaret A. Riley, Marcus W. Feldman, Brendan J. M. Bohannan. Local dispersal promotes
biodiversity in a real-life game of rock-paper-scissors. Nature 418, 171-174 (2002).
[21] Corina E. Tarnita, Hisashi Ohtsuki, Tibor Antal, Feng Fu, Martin A. Nowak. Strategy selection in structured populations. Journal of Theoretical Biology 259,570-581 (2009).
[22] Martin A. Nowak, Robert M. May. Evolutionary games and spatial chaos. Nature 359, 826-829 (1992).
[23] Erez Lieberman, Christoph Hauert, Martin A. Nowak. Evolutionary dynamics on graphs. Nature 433,
312-316 (2005).
[24] Hisashi Ohtsuki, Christoph Hauert, Erez Lieberman, Martin A. Nowak. A simple rule for the evolution of
cooperation on graphs and social networks. Nature 441, 502-505 (2006).
[25] Martin A. Nowak, Corina E. Tarnita, Tibor Antal. Evolutionary dynamics in structured populations.
Philosophical Transactions of the Royal Society B: Biological Sciences 365, 19-30 (2010).
[26] Tibor Antal, Istvan Scheuring. Fixation of Strategies for an Evolutionary Game in Finite Populations.
Bulletin of Mathematical Biology 68, 1923-1944 (2006).
[27] Arne Traulsen, Jorge M. Pacheco, Martin A. Nowak. Pairwise comparison and selection temperature in
evolutionary game dynamics. Journal of Theoretical Biology 246, 522-529 (2007).
[28] Michael Assaf, Mauro Mobilia. Large fluctuations and fixation in evolutionary games. Journal of Statistical
Mechanics: Theory and Experiment P09009 (2010).
[29] Philipp M. Altrock, Chaitanya S. Gokhale, Arne Traulsen. Stochastic slowdown in evolutionary processes.
Physical Review E 82, 011925 (2010).
[30] Philipp M. Altrock, Arne Traulsen, Tobias Galla. The mechanics of stochastic slowdown in evolutionary
games. Journal of Theoretical Biology 311, 94-106 (2012).
[31] Nicolaas Godfried van Kampen. Stochastic Processes in Physics and Chemistry, Third Edition. North
Holland (2007). ISBN: 978-0-444-52965-7
[32] Alexander Altland, Ben Simons. Condensed Matter Field Theory, Second Edition. Cambridge University
Press (2010). ISBN: 978-0-521-76975-4
[33] Alexandre Lefevre, Giulio Biroli. Dynamics of interacting particle systems: stochastic process and field
theory. Journal of Statistical Mechanics: Theory and Experiment (2007) P07024.
[34] James D. Murray. Mathematical Biology: I. An Introduction. Third Edition. Springer (2002). ISBN: 0-38795223-3
[35] Michael Assaf, Baruch Meerson. Extinction of metastable stochastic populations. Physical Review E 81,
021116 (2010).
[36] Ryogo Kubo, Kazuhiro Matsuo, Kazuo Kitahara. Fluctuation and relaxation of macrovariables. Journal of
Statistical Physics 9, 51-96 (1973).
77
《修士論文》
物性研究・電子版 Vol. 3, No. 2 032601 （2014年5月号）
[37] M.I.Dykman, Eugenia Mori, John Ross, P. M. Hunt. Large fluctuations and optimal paths in chemical
kinetics. Journal of Chemical Physics 100(8), 5735-5750 (1994).
[38] Roger. M. Nisbet, William .S.C. Gurney. Modelling Fluctuating Populations. The Blackburn Press (2004).
[39] Hakan Andersson, Tom Britton. Stochastic Epidemic Models and Their Statistical Analysis (Lecture Notes
in Statistics). Springer (2000). ISBN: 978-0-387-95050-1
[40] Mark Ptashne, A Genetic Switch, Third edition: Phage Lambda Revisited. Cold Spring Harbor (2004).
ISBN: 0-87969-716-4
[41] Earl A. Coddington, Norman Levinson, Theory of Ordinary Differential Equations. McGraw-Hill (1955).
[42] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery. Numerical Recipes 3rd
Edition: The Art of Scientific Computing. Cambridge University Press (2007). ISBN: 978-0-521-88068-8.
Section 18.3.
[43] Charles R. Doering, Khachik V. Sargsyan, Leonard M. Sander. Extinctiontimes for birth-death processes:
exact results, continuum asymptotics, and the failure of the Fokker-Planck approximation. Multiscale Modeling and Simulation. 3(2), 283-299 (2005).
[44] Otso Ovaskainen, Baruch Meerson. Stochastic models of population extinction. Trends in Ecology and
Evolution 25(11), 643-652 (2010).
[45] Carlos Escudero, Alex Kamenev. Switching rates of multistep reactions. Physical Review E 79, 041149
(2009).
[46] Baruch Meerson, Pavel V. Sasorov. Extinction rates of established spatial populations. Physical Review E
83, 011129 (2011).
78