《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

細胞の
自律力学
名古屋大学理学研究科
物質理学専攻（物理系）
非線形物理研究室（R 研）
博士前期課程２年
田辺光弘
平成 26 年 3 月 24 日

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

概要「群れ」はユビキタスで、イワシ、ムクドリ、アリ、更にはバクテリアやアメーバも群れで行動する。ここで「群れ」とは単に局在しているということではなく、生存戦略の手段としての協調的な集団行動を指す。多くの場合、集団を牽引するリー
ダーが不在で、各個体は自分の周りの局所的な情報のみに従い行動する。それでもあたかも全体が一つの意志を共有しているかのごとく振る舞う [7][8][35]。ここに「群れ」の面白さがある。我々はアメーバをモデル生物に選び、興奮性と走化性のある粒子からなる多体モデル Swarm Oscillators を考え、自発的に現れる秩序的な挙動を調べた。以下の二つの成果を得た。
まず一つ目の成果について説明する。先行する巌佐氏らの研究で本モデルが走化性のある細胞集団の動きを再現することが予言された。その後、具体的に

D.discoideum が飢餓状態に起こす凝集行動に似た振る舞いが発見された。そこ
で私は具体的に解を構成するという最も明らかな形で、その見つけられた解の存在を証明した。さらにその構成した解を詳しく調べることで、凝集過程において位相面の特異点が重要であることを発見した。また D.discoideum を用いた実験も行い、実際の系との整合性も確認した。スパイラル状の波が凝集に伴うことは 70
年代からよく知られていた事実であるが、数理モデルでそのメカニズムを解き明かしたのは本研究が初である。
次に二つ目の成果について説明する。本モデルが提案された当初はパラメータ
チューニングにより現れる粒子の様々な振る舞いに興味が持たれ、走化性粒子の示
す全ての振る舞いを説明し得ると期待された。そして少数粒子の場合に限り、それぞれの現象のメカニズムが巌佐氏により説明された。しかし多体現象は未解明だった。そこで私は十分に粒子数が多い状態で何が起こるかを予言する理論を作った。先行する結果と全く矛盾せず、あらゆるパラメータ領域でのモデルの挙動を定性的に説明できた。
本モデルは走化性粒子に対する広いクラスの数理モデルを包含し、将来的には全ての走化性粒子の研究で用いられることになるだろう。また社会現象に対しての適用を唱える研究も現れ [39]、本モデルが適用できる射程は予想を上回る広がりを見せている。本研究の成果は後続する数多くの研究に影響を与えることと信じている。以上が私の研究内容である。

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

目次
第 1 章 イントロダクション

3

第 2 章 キイロタマホコリカビの集団挙動を記述する数理モデルの設計

6
7
12
12
15
15
18
20

2.1
2.2

2.3
2.4

キイロタマホコリカビについて . . . . . . . . . . . . . . . . . . . .
先行する研究について . . . . . . . . . . . . . . . . . . . . . . . . .

2.2.1 Laub Loomis の化学振動子モデル . . . . . . . . . . . . . . .
2.2.2 Goldbeter の化学振動子モデル . . . . . . . . . . . . . . . . .
2.2.3 長野の細胞運動モデル . . . . . . . . . . . . . . . . . . . . .
Chemotactic Oscillator について . . . . . . . . . . . . . . . . . . . .
Swarm Oscillators について . . . . . . . . . . . . . . . . . . . . . .

第 3 章 Swarm Oscillator の凝集フローの解析

3.1
3.2
3.3
3.4

Swarm Oscillator の連続化について . . . . . . . . . . . . . . . . . .
特異摂動によるダイナミクス分解 . . . . . . . . . . . . . . . . . . .
自己無撞着方程式による平均場の決定 . . . . . . . . . . . . . . . .
細胞のフローの計算 . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23
27
30
46

第 4 章 安定性解析による田中・巌佐チャートの解読

4.1
4.2

50
巌佐氏の２粒子の厳密フローによるチャートの定性的解読 . . . . . 54
線形安定性解析によるチャートの解読 . . . . . . . . . . . . . . . . 57

第 5 章 キイロタマホコリカビの凝集フローの観察

5.1

5.2
5.3
5.4

実験の手順

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 D.discoideum の培養について . . . . . . . . . . . . . . . . .
5.1.2 凝集の観察準備と観察 . . . . . . . . . . . . . . . . . . . . .
5.1.3 画像解析の方法 . . . . . . . . . . . . . . . . . . . . . . . . .
実験データとその解析結果 . . . . . . . . . . . . . . . . . . . . . . .
結果の考察と課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
結論と課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 6 章 Appendix

6.1

63
64
64
64
65
66
66
68

71
Appendix A:化学反応振動子モデルの Hopf 分岐の証明 . . . . . . . 71
1

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

6.1.1 W.Liu の Hopf 分岐クライテリオン . . . . . . . . . . . . . .
6.1.2 Liu のクライテリオンの適用例 1（Brus-selator） . . . . . .
6.1.3 Liu のクライテリオンの適用例 2（Lorentz 方程式） . . . . .
6.2 Appendix B:Chemotactic Oscillator の位相縮約 . . . . . . . . . . .
6.3 Appendix C:Particle Image Velocimetry の概要 . . . . . . . . . . . .

2

72
75
76
77
85

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

第1章

イントロダクション

キイロタマホコリカビという生き物をご存知だろうか？粘性の高い液状の物質を先に付けた支持棒が地面から真っすぐ上向きにに生えた形をした子実体の写真
は、誰もが一度は見たことがあるはずである図 1.1。その球は休止状態になった細胞の集まりで、この細胞たちは動物が触れたり、雨で流されたりするなど、外界からアクションで地面に落ちるのを待っている。また支える棒も同じ細胞で出来ており、この細胞たちは球を地面から浮かせる支えとしての役割を担い、次世代に受け継がれない犠牲となる細胞たちである。その先端の球体が雨で流されたり動物に運ばれたりして地面に落ちてしばらく経つとアメーバに発芽する。アメー
バは周囲を動き回り、辺りの餌を食べ尽くす。周りの餌が無くなり飢餓状態になるとアメーバは化学物質をお互いに発し、寄せ集まってあたかも１個体のナメク
ジのごとく別の場所へと移動する。誰もがこの合理的な生存戦略に驚かれるだろう。アメーバが集まる過程、ナメクジ状態での行動、子実体が形成される過程、発芽過程。どのフェーズを見てもこの生物は面白い。特にアメーバフェーズの終わりに見られる１点に集まる過程は細胞が互いにコミュニケーションを取って行動するという観点からとても興味深い。化学走性と興奮性と呼ばれる性質だけから生まれる現象で、細胞が誰かから指示を受けて動いているように見える、いわば統率性が勝手に生じているように見える非常に不思議な現象である。ここで一つ問題を提起する。
「どのようにして特別な性質や振る舞いをするリーダー細胞無しに細胞たちは１点に集まることができるのか？」
細胞が一様に広がっていても飢餓状態になると細胞たちは cAMP 放出のスパイ
ラルを作りながら勾配に沿って１点に集まってくるのである。細胞一つ一つにはまったく個性も知性も無い。夢のある文学的な解釈をあてるならば、次のように仮説を立てられる。
「何らかのメカニズム (X1)”により集団としての知性が創発され、統率的な行動が行えているのだ。言い換えるならば個が全体意志を共有するのではなく、集団効果がダイナミクスを構造的にしているのだ。」
では、ここでいう”何らかのメカニズム (X1)”とは何なのか？私の研究はこの
何らかメカニズム (X1)”の正体に迫っている。もっと冷静な視点でこの現象を考えるならば、次のような壁にまずぶつかる。無個性な細胞を一様に広げれば、特異3

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

点なんて存在せず延々と一様なままのはずだ。そこで、考える自由度を一つ足すことを試みるべきだろう。細胞は cAMP と呼ばれる化学物質を感知し放出する機構を持っており、それはサイクルをなしている。そこで、細胞がその化学反応のどこに居るか、すなわち細胞に位相を付随させようではないかと。これが新しい
自由度である。さらに位相は互いに相互作用し、飢餓状態に位相引き込みが起きているとも仮定しよう。この仮定は極めてミニマムなものであり、無相互作用な
らばそもそも考えること自体が無駄になる。人類は局所的に結合した位相振動子
系についての物理学の知識を既に持っている。誰もが知るように、局所的に引力相互作用する移動振動子系を考え、位相が一様に乱れた状態からの時間発展を考えた場合、勝手に分布が不安定化し、自発的に位相が欠陥した点（あるいは小さな領域）が発生する。そして大域的に見ればスパイラルが発生する (Y)。そしてその位相欠陥に対して 何らかのメカニズム (X2) で粒子が引きずり込まれて集まるのではないかと考えられる。またしても 何らかのメカニズム という言葉を使った。ただし、この推論で再び現れた (X2) は問題 (X1) から蔵本モデルのスパ
イラル発生機構の知識 (Y) が差し引かれ、問題の所在がより明確になっている。

(X1) = (X2) + (Y )(1.1)

問題 (X2) は次のように言い表すことができる。
「粒子の移動が cAMP 勾配に駆動されることは間違いない。そして、集まる際にはスパイラル状の位相波が発生する。抽象的な概念である位相波が cAMP 駆動する粒子を駆動しているように見えるメカニズムは一体何なのか？」
この問題に対して、本研究で明確な答を与える。本論文では、ある極めて一般的な自己駆動多体系を数理モデルとして仮定し、その集団ダイナミクスについての詳細な解析を述べる。細胞がシグナリングしながら一点に集まる現象がどのようにして創発するのかを構成的に説明する。また同時にその正当性について、細胞性粘菌の一種であるキイロタマホコリカビをモデル生物として検証した結果も述べる。
最後に本論文の構成について説明する。まず２章でキイロタマホコリカビについての簡単な説明をする。どのようにして、外部の cAMP を感じ取り、cAMP を分泌し、勾配に沿った移動をし、そして周りに溜まりすぎないように cAMP を除去しているのか。生化学的な視点から理解する。次に議論で使用する数理モデル
について説明する。そこから３章以降は凝集過程を数式を使ってどのように理解できるのか説明する。本研究のメインの目的はここで尽きる。しかし、本研究の計算を通して、モデルの先行研究を上回る内容の結果を副産物として得たので、それについて４章で説明する。５章では実際にキイロタマホコリカビを使った実験について説明する。本モデルはキイロタマホコリカビに限らず、様々な物理系、あるいは社会学的な問題に対して適用されることが知られ [39]、モデルの取り扱い方

4《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

法を深めておくことは意義のあることである。

図 1.1: D.discoideum の子実体の写真。マクロレンズを使って上から撮影した。

5

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

第2章

キイロタマホコリカビの集団挙動を記述する数理モデルの設計

この章ではキイロタマホコリカビ（以下、学術名の D.discoideum で呼ぶ）の集
団運動を記述する数理モデルの作り方について説明する。まずはその前提知識となる D.discoideum のシグナル制御の分子生物学的な説明をし、これらの知見を踏まえた化学反応方程式からモデルの説明を進める。過去にシグナル制御が振動的に起こることを確認する研究が数多くなされている。例えば Goldbeter のモデル [13]や Laub、Loomis[19] のモデルを参照したい。素直に化学反応方程式を書き、時定
数の差を利用して簡単化したものであるが、見た目はとても複雑で、数値解析による周期軌道の研究が先行している。また、このモデルの周期軌道が Hopf 分岐由
来でありパラメータに関して安定なことを、Hopf の分岐定理と Rauth Hurwitz の定理を組み合わせた Liu 氏のクライテリオン [15][33] を用いて厳密に証明した研究もある [23][21]。このアイディアは Appendix A に書くことにするが、この研究により D.discoideum の内部自由度の振動は Hopf 分岐する力学系で考えれば良いことが推論できる。さらに我々は細胞の内部自由度のタイミングに合わせて化学物質が空間に放出される現象、その化学物質が拡散、分解される現象、さらに化学物質の濃度勾配に駆動される粒子の移動を記述する方程式を考えなくてはならない。この章ではそれらの機能が全て入ったモデルをどのようにして作るかをこの章で説明する。

6

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

2.1

キイロタマホコリカビについて

学術名 Dictyostelium discoideum。略称 Dicty。世界各地の温暖な地域に生息する。国内にも生息し、例えば北海道や長野県の高度 2000m から 3000m 程度の高所帯で、日の光の当たり、やや湿り気のある山肌で見つけることができるとされる。
しかし通常、その大きさや土の表面に混ざっているという生態からして人の目に留まることは考えにくい。
多フェーズの生活環を持ち、子実体 stalk、アメーバ amoeba、塊 mold、ナメク
ジ slug と４つの形態を順に周期的にたどる図 2.1。定常状態は子実体と呼ばれる状
態で、地面から垂直に伸びた棒状の支えの先に球形の液体として休止している。支えになっている部分も同じ細胞からできているのだが、支えになった細胞は死ぬ
ことが決まっていて次世代には受け継がれない。実はどの細胞が受け継がれるか、どの細胞が支えとして犠牲になるかということは slug の段階で既に決まっている。
その子実体は小動物が触れたり、雨風などにさらされたり、外部からの何らかの外からの作用で先端が地面に落ちる。しばらくすると、発芽してアメーバ状態になる。アメーバになると例えば大腸菌などの環境にある餌を食べながら動き回る。
そして周囲の餌が少なくなって飢餓状態になると、細胞はそれまでとは違った振る舞いを始める。周囲に cAMP と呼ばれる化学物質を放出し、また周囲の cAMP
勾配を感知して移動し、互いに集まり合って細胞の塊を作る。一定の大きさまで集まると、ナメクジ状に形を変え、別の場所へと移動し再び子実体を形成する。このサイクルの周期は種や環境に寄ってはこの限りではないが、およそ４日程度と分かっている。このように D.discoideum は確実に子孫をつなぐために、形を様々に変えながら合理的な生存戦略をとっている。
以上のような現象の豊かさに加え、D.discoideum は遺伝子操作に向く遺伝子操
作に向いており、光走行性、熱走行性、カフェインに対する応答、レセプターのノックアウトによる振る舞いの変化など様々な観点からモデル生物として 70 年代
からよく研究されてきた対象である。特に本研究で注目している凝集過程については、cAMP の放出、感知、除去のプロセスについて一部の G 共役タンパクの機構を除いて分子生物学的機構はよく分かっている。また、その撮影も 70 年代から既に成功していて [2]、現在では蛍光染色を利用した１細胞レベルでのトラッキン
グ撮影行われている [6]。１細胞レベルだけでなく集団としての挙動も実験、理論
ともに調べられており、イントロでも述べたようにリーダーなく１点に集まることができるのメカニズムは、より広い枠組みで集団の行動の問題として考えても興味深い。
ここではシグナルの分子生物学的知見、すなわち cAMP の放出、感知、除去の機構について既知の事実をまとめておく。数理モデルを設計する上で対象となる生物の実際の機構が分かっていることは、モデルの妥当性を考える上で重要だか7

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

らである。しかし、過剰な労力をかけて一つ一つ分かっている酵素の特性を説明していくことはせず、ここではどんな酵素がどのように働いて制御システム構成
しているのかを要約するにとどめる。
図 2.2 に cAMP 合成の制御機構のパスウェイを示す [28]。cAMP はアデニル酸シクラーゼ A（Adenylyl cyclase A,ACA）という膜貫通型の合成酵素により行われる。ACA は２つのドメインからなり、凝集過程において活発に発現される。細胞外の cAMP 濃度に応じた活性化機構と抑制機構が ACA に対して働いていて、興奮性と順応性の機能をなしている。抑制機構についてはよく分かっていない。活性化は以下のようなシナリオで行われる。cAR1 に cAMP が結合すると 3 量体結合タンパクである G2 が α サブユニットと βγ サブユニットに分解される。βγ サ
ブユニットはフォスファチヂルイノシトール２キナーゼ (PI3K) を活性化し、フォ
スファチヂルイノシトール２リン酸 (PIP2) をフォスファチヂルイノシトール３リ
ン酸 (PIP3) にリン酸化する。PIP3 はアデニル酸シクラーゼの細胞質レギュレー
タ (CRAC) の結合ドメイン（プレクストリン相同ドメイン）に結合し、ACA の活性化に関与する。また、その他図の下方に書いたタンパク質群（ERK2：MAP キ
ナーゼ、RasC：グアニンヌクレオチド交換因子、RasGEF：グアニンヌクレオチ
ド交換、Pia：ピアニッシモ、Eip3）も ACA の活性化に関与するが詳細は分かっていない。以上が正のフィードバック制御である。また cAR1 は外部の cAMP の結合に応じて ACA の活性を抑制する負のフィードバック制御をすることも分かっている。しかしその詳しいパスウェイは未解明である。
次に cAMP の感知機構について説明する。D.discoideum の cAMP 感知には４つの G タンパク共役受容体 cAR1、cAR2、cAR3、cAR4 が関与している。しかし、それぞれ発現するタイミングが異なり、cAR1→cAR3→cAR2→cAR2 の順番に発現する。具体的にそれぞれの多く発現されるタイミングを説明する。cAR1 は cAMPの濃度勾配に駆動された凝集過程に、cAR3 は凝集過程の終盤に、cAR2 はナメク
ジ状になって移動する過程に (slug 状態)、cAR4 は子実体の状態にそれぞれ多く発現される。特に我々は細胞癒着が進んでいない凝集過程に着目しているので、ここでいう cAR1 が発現する段階が問題となる。ただ、濃度勾配を感知して細胞移動を起こすのは良いのだが、cAMP がパルス的に伝播するパルス波が伝わって来た場合、細胞は濃度の上り坂を感じてパルス波と逆方向に動き、その後、下り坂を感じて元に戻るのではないかという議論がある。パラドキシカルなこの事実を説明する仮説として、Goldstein のトイモデル [22] と石渡氏、巌佐氏のトイモデルの２つのシナリオが提案されている [25]。
最後に cAMP の除去機構について説明する。D.discoideum には 5 つのヌクレオチ
ドホスホヂエステラーゼが関わっていて、Pde3、PdeD という２つの酵素が cGMPにのみ特異性がある酵素、PdsA が cAMP と cGMP の両方に同程度の特異性がある酵素、そして RegA と PdeE が cAMP にのみ特異性がある酵素である。従って8

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

cAMP の加水分解に関わるのは PdsA、RegA、PdeE の３つの酵素である。PdsA
は飢餓状態に発現する酵素であり、凝集過程の制御機構を考える上で最も重要と考えられる。RegA も飢餓状態に発現するのだが、飢餓させて 4 時間程度経った後に発現することが知られていて、子実体の軸となる前の細胞（予定枝細胞）あるい
は胞子となる前の細胞（予定胞子細胞）において重要であるとされ（図 2.1）、この酵素は我々の問題から除外できるだろう。また、PreE は cAMP に対してアフィ
ニティが小さく、PreE を発現しないミュータントにおいても通常と同じ細胞の振るまいが見られるという。つまり PdeE は外部 cAMP を増加させまいとする穏やかな抑制制御の働きをしている。

9

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 2.1: D.discoideum のライフサイクル。この図は [24] を参考に描いた。子実体(spores) が発芽し、アメーバ (amoebae) となり分裂増殖し、飢餓状態になると集まって塊を作り (c)、その後ナメクジ状の状態 (d) を形成し移動する。しばらく移動した後、突起が現れ (e)、軸を上へと伸ばして (f) 子実体を形成する (g)。slug(d)の後部の黄色く書かれた細胞は子実体を形成する時に軸を形作り、後世には引き継がれない。

10

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 2.2: D.discoideum のシグナル機構の模式図。[28] を参考に描いた。cAMP のcAR1 への結合から ACA による cAMP 合成までのシグナルの流れを示した図。下のクエスチョンは回路が未だ決定されていない酵素である。外部 cAMP の濃度上
昇が cAMP の合成を抑制するフィードッバク機構となっている。

11

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

2.2

先行する研究について

ここでは我々のモデルと共通する所の多い長野氏のモデル [27] についてここで説明する 1 。彼のモデルは我々のモデルと共通する点が式の形だけでなく目的においても多く、ここで彼の研究をレビューすることは我々の成果を自己評価するのに必要でと思うからである。しかし、彼のモデルを理解するためにはそれに先行す
る Goldbeter モデルを理解する必要がある。しかし、Goldbeter モデルは抽象的で分かりづらい嫌いがあるため、ここではより具体的で分かりやすい Laub Loomis
モデルと比較して説明する。

2.2.1

Laub Loomis の化学振動子モデル

Laub、Loomis の提案したモデルは次の化学反応方程式で定式化される。
d[ACA] = k1 [ERK2] − k2 [ACA]dtd[P KA] = k3 [i.cAM P ] − k4 [P KA]
dtd[ERK2] = k5 [CAR1] − k6 [ERK2][P KA]
dt
d[REG − A]dtd[i.cAM P ]dt
d[e.cAM P ]dtd[CAR1]
dt= k7 − k8 [REG − A][ERK2]= k9 [ACA] − k10 [REG − A][i.cAM P ]= k11 [ERK2] − k12 [e.cAM P ]= k13 [e.cAM P ] − k14 [CAR1][P KA](2.1)

定数値は次の通り。

k1 = 1.4min−1 、k2 = 0.9min−1 、k3 = 2.5min−1 、k4 = 1.5min−1 、k5 = 0.6min−1 、
k6 = 0.8min−1 μM −1 、k7 = 2.0min−1 μM −1 、k8 = 1.3min−1 μM −1 、k9 = 0.7min−1 、
k10 = 1.0min−1 μM −1 、k11 = 0.3min−1 、k12 = 3.1min−1 、k13 = 1.8min−1 、
k14 = 1.5min−1 μM −1 。
これらのデータはさらに先行する研究の実験データに基づいている。この値で数値計算すると自発的に振動する解が得られる。
図 2.3 に示した経路を仮定し、この化学反応方程式からモデルが設計されている 2 。各反応の意味は次の通りに明確である。cAR1 が外部の cAMP との結合で活
1

最近、Loomis モデルを使った Nagano モデルの論文も発表された [20]。
Laub Loomis モデルは cAMP の振動を導くタンパク質レベルのシグナルネットワークモデル
としては成功を収めた。その後で、細胞内のリン酸ジエステラーゼ（RegA）と PKA が中心的な
2

12

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

性化し、ERK2 を活性化、引き続いて ERK2 が ACA を活性化し cAMP が合成される。すると、外部の cAMP 濃度が上昇し、正のフィードバックがかかる。同時に内部の cAMP 濃度も上昇する。内部の cAMP が上昇すると、PKA が活性化される。PKA が ERK2 を抑制制御し、cAMP 合成経路が阻害される。ERK2 は REG

A を抑制する機能があるのだが、抑制が弱まり REG A が活性化する。REG A により内部 cAMP は分解される。このようにして、内部の cAMP 濃度が上昇しすぎないように制御されている。
外部 cAMP の方程式には粒子配置などの情報は平均化されている。また ERK2周りの反応は前述の通りまだ実験的には証明されていない部分があり仮定を含んでいる。とはいえ、D.discoideum の化学反応を記述するモデルとしては、最も自然で理解しやすい化学反応モデルである。

役割をもって振動するという Laub Loomis とは異なるシグナル経路の模型も澤井氏により提案されている [29]。

13《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 2.3: Laub Loomis モデルで仮定されているシグナル経路。[19] より引用。cAR1に細胞外部の cAMP が結合すると、ERK2 を活性化し、続いて ACA が活性化され cAMP が合成される。ERK2g 活性化すると、同時に REG A が抑制され、ACAから PKA への経路が抑制される。PKA は ERK2 の活性を抑え、結果的に cAMPの合成を抑える作用をする。一方で、細胞内部の cAMP の濃度が上昇すると PKAの活性が上がって ERK2 が抑制される。冒頭で述べたシグナル機構とは少し異なり、細胞外部の cAMP 濃度上昇による抑制制御を考えていない。

14

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

2.2.2

Goldbeter の化学振動子モデル

Goldbeter 氏が提案したモデルは次の式 (2.2) で定式化される。細胞内 cAMP、細胞外 cAMP、レセプターの３変数で記述する極めてシンプルなモデルである。
Michaelis Menten kinetics をベースに作られており、反応経路で関わる酵素反応などは繰り込まれている。歴史的な順序で言えば逆だが、このモデルはいわば Laub
氏と Loomis 氏のモデルをマクロに簡約したものと見なすこともできる。
k1 + k2 [e.cAM P ]k1 L1 + k2 L2 c[e.cAM P ]
d
ρ = −
ρ+
(1 − ρ)dt
1 + [e.cAM P ]1 + c[e.cAM P ]2ρ[e.cAM P ]α(λθ + 1+[e.cAM P ] )d[i.cAM P ] = qσ
− (ki + kt )[i.cAM P ]2dtρ[e.cAM P ]1 + αθ + (λθ + 1+[e.cAM P ] )(1 + α)
kt [i.cAM P ]d[e.cAM P ] =− ke [e.cAM P ]
dt
h

(2.2)

ρT は活性化している受容体の全体での割合、[i.cAM P ] は細胞内の cAMP 濃度、[e.cAM P ] は細胞外の cAMP 濃度を表す。細胞外 cAMP の濃度変化は化学反応方
程式は内部の cAMP 濃度に比例して増加速度が上がり、一定の確率で分解されていくという方程式になっている。また、Goldbeter 氏らが数値計算に用いた各定数の値は次の通り。
k1 = 0.036min−1 、k2 = 0.666min−1 、ke = 5.4min−1 、L1 = 10、L2 = 0.005、
c = 10、q = 4×103 、σ = 0.6min−1 、λ = 10−2 (no experimental data)、θ = 10−2 (noexperimental data)、 = 1(no experimental data)。
最後の３つを除き、これらのデータはさらに先行する研究の実験データに基づいている。この値で数値計算するとやはり自発的に振動する解が得られる。それ
ぞれのモデルでは細胞外 cAMP 分解の方程式を拡散方程式で置き換えて作られて
いる。

2.2.3

長野の細胞運動モデル

長野氏は次のようなモデルを提案した。Goldbeter のモデルをベースに、細胞外部 cAMP 濃度が空間不均一であるとして拡散方程式に取り替えて拡張したモデル

15

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）である。まず、cAMP の分泌および空間拡散を記述する方程式はつぎの３つ (2.3)。


k1 
4
∂S(t, r)2
= −ke S(t, r) + DΔS(t, r) +βμ (t) exp − 2 (r − rμ (t))
∂t
h μ
σ0

2
ρμ S(t,rμ )0.0001+
1+S(t,rμ )dβμ (t)= 18002 φ(ρμ , S(t, rμ (t))) − (ki + kt )βμ (t)dtρμ S(t,rμ )0.2575 + 1+S(t,rμ)

dρμ (t)k−1 + 10k−2 S(t, rμ (t))
k1 + k2 S(t, rμ (t))=(1 − ρμ (t)) −ρμ (t) (2.3)
dt
1 + 10S(t, rμ (t))
1 + S(t, rμ (t))
ここに、運動方程式を組み合わせる。各細胞が細胞外の cAMP 濃度勾配に駆動され、減衰項と細胞癒着ポテンシャルを入れる。具体的には次の方程式 (2.4) を仮定
する。

ma

 ∂
drμ
d2 rμ
∂γ
=
(t,
r
)
−
v
(|r
−
r
|)
−
η
1
μ
m,n
μ
ν
dt2
∂rμ
∂rμ
dt
ν(=μ)  mn
1 σ0 m 1  σ0 nvm,n (r) = 0−
m−n m rn r

(2.4)

このモデルを以下のパラメータ設定で数値計算するとつぎのような現象が見られる。各細胞内で振動的な化学反応が発生し、それに応じて振動的に細胞が cAMPを外部に分泌する。それによって、空間にできる cAMP 濃度勾配に駆動されて細胞は空間を移動する。その過程はスパイラルを描くという（長野氏の論文ではス
パイラルの波が発生するメカニズムまでは明らかにされていない）。自然と細胞同士が集まり、ある程度の距離まで細胞同士が近づくと、細胞癒着ポテンシャルにより σ0 の距離に細胞同士がトラップされる。全体としてスパイラルを描きながら凝集する過程がみられるという。その様子を図 2.4 に示す 3 。

ke = 5.4min−1 、ki = 1.7min−1 、kt = 0.9min−1 、k1 = 0.036min−1 、k−1 =0.36min−1 、k2 = 0.666min−1 、k−2 = 0.00333min−1 、h = 0.104、D = 4.0 ×
10−4 cm2 /s、ma = 10−9 g 、 0 = 10−13 g mm/min2 、 1 = 10−17 g mm/min2 、m = 9、n = 3、m = 9、η = 0.5 （space unit:0.333σ0 、time step:0.0033min）

3

この結果については長野正道先生より転載許可をいただきました。

16

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 2.4: 長野氏のモデルの計算結果。[27] より引用。5 min の図には t = 0min、

t = 3min における位置がそれぞれ三角、丸で示されている。スパイラルを描きな
がら粒子は集まっている。

17

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

2.3

Chemotactic Oscillator について

スパイラル状の波を発生させて凝集することを示すだけなら、長野氏のモデルを調べればいいじゃないかと考える方もいると思う。しかし、彼のモデルはミニマルに設計されているとは言え、まだまだ複雑で解析的な手法を使うにはほど遠い。そこで、ここからは解析的に取り扱えるモデルを作ることを考える。更にミニマルに情報を落としていくと言うと、それはアーティフィシャルな操作が入るのではないかと懸念する人もいるだろう。その心配は無用である。ここでの導出
は一度普遍的なモデルを考え、不必要な情報を系統的に落としていく方針をとる。
アーティフィシャルさを一切感じない系統的な方法である [5]。
まずは、カビの内部で起こる化学反応をどのように記述するかを考える。ここまでに説明した Goldbeter や Loomis のモデルのようなものを考えればよい。そこ
で、両者のモデルを含む以下の方程式 (2.5) で内部自由度が記述できると仮定する。

dXμ(t) = f (Xμ ; λ) + kg (S (rμ ))(2.5)
dt
X は内部自由度を記述する多次元の変数であり、例えば内部の化学反応に関わる化学物質の濃度のリストを考えれば良い。もっと複雑な変数でも構わない。関数 S(t, r) は時刻 t、二次元平面上の位置 r における cAMP の濃度を表す。f の関数
形はあるパラメータ λ に依存し、λ = 0 で Hopf 分岐し λ > 0 で周期軌道があると仮定する。定数 k や関数 g に対しては特に要請はなく、g は S = 0 近傍で C 1 以上の滑らかさがあれば十分である。ミカエルス・メンテン型だろうが、量子論的か
ら導出する現代的な反応速度論を考えようが関係ない。振動的に振る舞う幅広い
化学反応方程式のモデルを含んでいる。
次に、シグナル伝達を記述する方程式について説明する。シグナルは周囲に２次元的に cAMP を拡散することで行われる。従って、cAMP が拡散方程式に細胞から放出されることを表す湧き出し項とホスホヂエステラーゼによる分解を表す
項を入れればよい。すると次の方程式 (2.6) が考えられる。
τ∂S(t, r) = −S(t, r) + dΔS(r, t) +h (Xμ ) δ (r − rμ (t))
∂t
μ

(2.6)

第一項が cAMP の分解、第二項が cAMP の拡散、第三項が cAMP の細胞からの放出を表す。τ は時間スケールを表す定数、d は拡散係数である。また、関数 h は放出のタイミングに値を持つ変数であり、X について原点周りで C 1 程度以上に滑らかであると仮定する。この式は長野氏の設計とほとんど同じである。
最後に粒子の移動を記述する方程式を説明する。cAMP 勾配によって加速度が決まる減衰のある次の２階の常微分方程式 (2.7) を考えることにする。

18

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

m

drμ2
drμ
dS
− σ (Xμ )(t) = −γ(t, rμ )
2
dt
dt
drμ
r ∈ M (M = Rdim , T dim , etc.)(2.7)

式 (2.7) の右辺第一項が減衰であり γ は減衰率を表す定数、第二項は走化性を表す項で σ は走化性の強度を表す定数で内部変数に依存しているものとした。後にSwarm Oscillator に移行する際にオーバーダンプ極限をとってしまうため、Swarm
Oscillator に縮約する場合は加速度の効果は考えない。空間 M は例えば dim 次元
ユークリッド空間 Rdim や周期境界条件を課した有限サイズの空間 T dim など何を考えても良い。また右辺に短距離の強いファンデルワールス型のポテンシャル項を足すことで、長野氏のモデルをも含めることができる。おそらく、縮約後の方程式の現象はシンプルになるだろうが、我々はそのモデルについては本研究では調べない。なぜなら、ここで考えるのは凝集とスパイラル状の波の関係であり、凝集の初期過程だけを考えても十分だからである。凝集が進んだ後の細胞癒着までを議論に入れるのは、後続する研究者のヒューチャーワークである。
さて、以上の３つの式を連立して解くことで走化性粒子の挙動を決めることができる。ここで再度振り返ってみると f や h などの関数形を具体的に決めておらず、Chemotactic Oscillators は D.discoideum に限らず広く走化性粒子の挙動を記述していることが分かる。当然ながらこのままでは数値計算はできないし、議論を進めるのは難しい。ここから先は、この普遍的なモデルから不要な情報をいかに系統的に削ぎ落としていくかを説明する。このモデルを位相縮約することで、実際に数値計算や手計算の解析が可能な具体的なモデルに落としこむことができる。

19

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

2.4

Swarm Oscillators について

前節では Chemotactic Oscillator モデルについて説明した (2.5)(2.6)(2.7)。細胞癒着効果以外の D.discoideum の行動を完全に記述していると思われる一般的な方程式であった。
ここで、式 (2.5) は極端に一般的な方程式で、振動するということだけを仮定した。振動するというのは曖昧な仮定であり、もう少し具体的な仮定をおく。章の冒頭で説明したように、Goldbeter や Loomis の化学反応モデルはロバストな周期軌
道を持っていることが Liu のクライテリオンから証明されている。見方を変えればそれと同時にパラメータ空間に Hopf 分岐点（面）を持つことも分かっている。
仮にその周期軌道が Hopf 分岐してすぐのあまり成長していない周期軌道であったとしよう。つまり、周期軌道の形があまり複雑に崩れていないということを仮定する。分かりにくい仮定だが、決して強過ぎることはなく自然な仮定である。
さて、このように内部状態の発展方程式が超臨界 Hopf 分岐直上であると仮定する。すると位相縮約と呼ばれる操作、つまり力学系の言葉で言えば中心多様体縮
約 4 に相当する操作を施すことで、扱いやすい方程式に変形することができる。さらに、周期軌道が整っているという性質を他の方程式にも使って、上手く変形することで Swarm Oscillators を得ることができる。
その導出は蔵本モデルの導出を参照すれば自然で馴染みやすいものだろう。ここ
では手順の要約を述べる。まず、式 (2.5)) が Hopf 分岐の標準形である複素 Ginzburg

Landau 方程式に座標変換する。式 (2.7) はオーバーダンプを仮定して慣性項を落としておく。その座標変換を式 (2.6) と式 (2.7) に適用し、式 (2.6) の形式解を Fourier
変換を使って求め、漸近解析的に近似する。ここで得られた化学物質濃度分布 S(t, r)
を式 (2.5) と式 (2.7) に代入し、定数を整えると求めている方程式 (2.8) を得ることができる。具体的な計算はやや長いので Appendix B で詳しく説明する。以上のよ
うな処理の結果が次の方程式である。
dθμ
=
dt



e−|rν −rμ | sin (θν − θμ + α|rν − rμ | − c1 ) ,ν(=μ)

drμr ν − rμ= c3
sin (θν − θμ + α|rν − rμ | − c2 )e−|rν −rμ |dt|rν − rμ |

(2.8)

ν(=μ)

0 ≤ θμ (t) < 2π は時刻 t における粒子の内部自由度を表す。２つ目の方程式はのダイナミクスを決める方程式である。相互作用の距離は e−|rν −rμ | で決まることを容易に見て取れるが、空間をここの係数が１になるように規格化している。規4中心多様体とは非双曲不動点から中立安定方向に伸びる不変多様体であり、この多様体の上に力学系の変数を制限する操作を中心多様体縮約と呼ぶ。安定方向の成分は指数関数的に素早く一様に減衰するので、この縮約は力学系の振る舞いをうまく抽出できるのである。位相縮約は複素
Ginzburg Landau 方程式に帰着させる系統的な変形方法である。[12][31][34][38][37] が詳しい。

20

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

格化のスケールは意外に複雑で、内部自由度の角振動数、化学物質の拡散方程式の拡散定数、そして拡散方程式の時定数の３つの変数が組合わさった定数となっている。また、sin の中身に内部自由度の位相と空間変数が入るのが特徴的で、このモデルが様々な性質を示す根源になっている。本来は無次元な変数しか入らな
さそうなここに空間変数が混ざってくるのは化学物質の拡散の波打ちに起因する。
定数は α、c1 、c2 、c3 の 4 つに減っており、基本的にはこれらの定数を動かすこと
で縮約前の一般的なモデルの全ての振る舞いを網羅的に調べることができる。うまくパラメータをチューニングすると、図 2.5 のように、スパイラル状の位相を伴いながら、中心へ粒子が集まっていく様子を見ることができる。
ただパラメータの意味解釈は難しく、縮約前の系からのそれぞれの変数変換は非常に複雑なので Appendix にその変数変換を示す。α は各粒子が見える視野のようなものであり、α が小さい程に各粒子は近くの粒子の影響しか受けない。また、

c1 、c2 は何とも表現し難い系の特徴を示す変数で、D.discoideum の場合は c1 ≈ 0、
π < c2 < 2π とおおよそ見積もられている [25]。以上が Swarm Oscillators5 の導出である。非常に広いクラスの内部振動反応のモデルを含んでおり、様々ある先行研究の内容を含む普遍的な方程式である。以下このモデルを D.discoideum を記述するモデルとして採択し議論を進める。

5

Swarm Oscillators は [5] で初めて導出された。以後、このモデルに関する研究としては [32]、[17]、[11] などがある。少数粒子の振る舞いについて非常に詳しく議論されている。

21

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 2.5: Swarm Oscillators の凝集の様子。パラメータの設定は c1 = 0.1、c2 = 4.0、

c3 = 0.01、α = 0.5。スパイラル状の位相波が発生し、中心に集まる時の初期段階のスナップショット。初期条件によって発生したり、しなかったりする。c3 を小さめに設定すると粒子が広がった状態が長く保たれるので、簡単にスパイラル状の位相波を作ることができる。

22

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

第 3 章 Swarm Oscillator の凝集フローの解析
さて、前章で導いた方程式 Swarm Oscillators を解きかたを以下で説明する。最も自然なアイディアは蔵本モデルとの類似性を利用である。蔵本モデルと同じ処
方箋で解析できると期待し、連続化から安定性解析、あるいは自己無撞着方程式
による解析に持ち込むのが最も自然なアプローチだろう。しかし、素直には蔵本
モデルの処方箋を使うことはできない。Swarm Oscillators は空間移動の自由度があるため、蔵本モデルでいう結合が時間ごとに変化するということである。そのような理論は蔵本モデルの処方箋には無い。蔵本モデルの理論に結合の時間変化を入れる新しいアイディアが必要である。空間移動が位相のダイナミクスに比べて十分に遅いことを利用すれば、漸近解析的な摂動法で近似的に解けることを私は思いついた。以下の節でその具体的な解析の方法と結果を説明する。
スパイラル波が発生するとどうして粒子が集まるのか、どのように集まるのかが具体的に分かる。その後、中心に特異点を持たない回転対称なターゲット型の波では細胞が集まらなくて、特異点を持っているスパイラル型の波の場合には中心付近で加速して効率よく細胞が集まることを示す。ここでは当たり前に認識されている「位相勾配に逆らって移動する」という定説だけでは説明できない効果が見られる。ターゲット型もスパイラル型も動径方向位相の勾配は全く同じであるが、ターゲット型の波には原点に特異点があり、今回初めて導入する平均場の絶対値にその違いが現れる。つまり、位相の勾配ではなく平均場としての勾配が粒子の移動を決めるということである。

3.1

Swarm Oscillator の連続化についていま注目している現象すなわち cAMP 放出の波伝播を伴う凝集過程は本質的に多体現象である。従って、方程式を連続化して考えることは妥当な第一ステップ
である 1 。まず次のような位相と空間の直積空間における状態密度分布を与える密度関数 (3.2) を導入する。
1

ここでの計算は [10] の 10 章を参考にしている。この本には普通の Kuramoto モデルの連続化と大偏差原理によるその正当化の議論が説明されている。

23

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

ρ(t, r, θ) =δ(r − rμ (t))δ(θ − θμ (t))(3.1)μ=1,··· ,N(t ∈ R, r ∈ M, θ ∈ S 1 = [0, 2π)).
従って、次の規格化条件 (3.2) を満たす。

1N

drdθρ(t, r, θ) = 1(3.2)

M ×S 1

また、粒子分布 S(t, r) は次のように与えられる (3.3)。

S(t, r) =dθρ(t, r, θ)

(3.3)

S1

（記号が重複したが、cAMP の濃度分布とはもちろん全く異なる関数である。）
ここで、連続化された Swarm Oscillators の導出方法を説明する。密度関数 (3.2)
を微分して次の式 (3.4) を得る。


∂δ
dθμ
∂ρ
(t, r, θ) = −(t)δ(r − rμ (t)) (θ − θμ (t))
∂t
∂θ
dt
μ
 ∂δ
drμ
−
(r − rμ (t))(t)δ(θ − θμ (t))
∂r
dt
μ

(3.4)

ここに Swarm Oscillators の運動方程式 (2.8) を代入するのだが、電子系の DFTなどと同様に、自己相互作用の問題が起こる。この場合は位相変数を

θμ(new) (t) := θμ(old) (t) − t sin c1(3.5)のように取り直して、更に単位ベクトルの定義について、対角成分を０と定義する。つまり、rν (t) − rμ (t)|rν (t) − rμ (t)|(new):= (1 − δνμ )rν (t) − rμ (t)|rν (t) − rμ (t)|

(old)

(3.6)

のようにすることで、運動方程式 2.8 の右辺の和が”ν : ν = μ”(μ を排除) から”ν”（自分を含めた全粒子の和）に置き換わり、この問題を回避することが出来る 2 。
2

注意点として位相の定義が変わっていることに気をつける必要がある。この定義で計算すると後に導入する Ω が sin c1 だけシフトするので、もとの Swarm Oscillators の振動数と対応させるには、シフト分を考慮して Ω(ref ined) = Ω(bear) + sin c1 と補正する必要がある。本論文での振動数はsin c1 シフトを補正していない。

24

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）√
以下のように変形する 3 。途中で c̃1 = c1 , c̃2 = c2 +tan−1 (−α)+ π2 , c̃3 = c3 / 1 + α2と置き換える。

∂ρ∂t(t, r, θ)∂ = −δ(r − rν (t))δ(θ − θν (t))e−|rν −rμ | sin (θν − θμ + α|rν − rμ | − c1 )∂θ ν∂ rν − rμ−c3
sin (θν − θμ + α|rν − rμ | − c2 )δ(r − rν (t))δ(θ − θν (t))e−|rν −rμ |∂r ν|rν − rμ |∂= − ρ(t, r, θ) dθ drδ(θ − θν )δ(r − rν )e−|r −r| sin (θ − θ + α|r − r| − c1 )∂θν∂
−c3 ρ(t, r, θ) dθ drδ(r − rν (t))δ(θ − θν (t))∂rνr − rsin (θ − θ + α|r − r| − c2 )|r − r|∂∂δ(θ − θν )δ(r − rν ) e−|r −r| cos (θ − θ + α|r − r| − c̃1 )= − ρ(t, r, θ) dθ dr
∂θ
∂θ
ν

∂
∂

−c̃3 ρ(t, r, θ) dθ drδ(θ − θν )δ(r − rν ) e−|r −r| cos (θ − θ + α|r − r| − c̃2 )
∂r
∂r
ν


e−|r −r|∂∂ρ(t, r, θ)dθ dr ρ(t, r , θ )e−|r −r| cos (θ − θ + α|r − r| − c̃1 )
∂θ
∂θ
∂
∂

dθ dr ρ(t, r , θ )e−|r −r| cos (θ − θ + α|r − r| − c̃2 )−c̃3 ρ(t, r, θ)∂r∂r −ic̃1∂∂∂e−ic2 ∂eρ(t, r, θ) + c̃3
ρ(t, r, θ)= −dθ dr ρ(t, r , θ )
2 ∂θ
∂θ
2 ∂r
∂r
exp [i(θ − θ) + (−1 + iα)|r − r|] + c.c.

= −

(3.7)

演算子、関数を以下に定義する (3.10)、(3.9)、(3.11) に従って置き換えると、連続化された Swarm Oscillator(3.8) が得られる。

∂ρ= −Dζ(t, r)e−iθ + c.c.,
∂t
3

(3.8)

本論文では 2 次元の変数”r = (x, y)”がベクトルであることを特に明示しない。微分に関しても同様に、∂/∂r := ∇ = (∂/∂x, ∂/∂y) と省略する。また、ベクトル量が二つ以上現れているときは内積をとるものとし、内積記号”·”を明示しない。例えば、




∂
∂
∂
∂
∂
∂
a(x, y) b(x, y) +a(x, y) b(x, y)a(r) b(r) :=
∂r
∂r
∂x
∂x
∂y
∂y

25

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

with
D :=
ζ[ρ](t, r) := ζ(t, r) :=∂e−ic̃1 ∂∂c˜3 e−ic˜2 ∂ρ(t, r, θ) +ρ(t, r, θ) ,
2 ∂r
∂r
2 ∂θ
∂θ


dt dr dθ ρ(t, r , θ )g(r − r)eiθ ,g(r) := exp ((−1 + iα)|r|) .

(3.9)(3.10)(3.11)

ζ(t, r) には点 r から半径オーダー 1 程度にある粒子の情報が含まれ、点 r にある粒子が感じる外場を表しており、それに対して演算子 D を作用させることで密度関数の時間発展を知ることができる。以下の章で、この ζ が凝集過程においてどのような形になるかを計算する。これを計算することで、凝集過程における位相の分布や速度勾配を知ることができる。
以下の議論では ζ が重要な役割を果たす。位相と絶対値因子を分けて考えるために、次のように表記を定義しておく。
ζ(t, r) =: R(t, r) exp[iΘ(t, r)](3.12)

特に位相がクラスター的に振動数 Ω で定常的に回っているときは、その回転を差し引くと議論しやすい。そこで、

ζ(t, r) =: R(t, r) exp[iΘ(t, r)] = exp[i(Ω(r)t + Θ0 (r))]のように定常回転を差し引いた関数 Θ0 (r) を定義する。

26(3.13)

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

3.2

特異摂動によるダイナミクス分解

ここで特異摂動論を使う理由は Swarm Oscillators のモデルに対して自己無撞着方程式を作ろうと試みればすぐに分かると思う。そのままでは特殊なパラメータを選ばない限り、閉じた方程式まで変形できないのである。だから別の方法を考えることにした。時定数のズレを利用した特異摂動論で位相と空間移動のダイナ
ミクスを分解し、近似的に計算すれば良いことを思いついた。複スケール法が特に今回の問題にマッチする [30][36]。まず時定数がどれくらい説明する必要があるが、この研究は巌佐氏と石渡氏により行われた [25]。彼らのロジックはシンプル
で、D.discoideum の実験論文に書かれていた凝集のスピードとスパイラルの伝播速度の比率から c3 ≈ 0.1 ∼ 0.01 の値をとっていると評価できるという。この結果
から空間移動のダイナミクスは位相のダイナミクスに比べて 10 倍から 100 倍程度遅いということになる。このタイムスケールの乖離を利用するのであり、思想としては Born Oppenheimer 近似のようなものである。
実は偏微分方程式系に対して複スケール解析を適用した例は少なく、ここでの計算自体も多くの人が目新しく感じるのではないかと思う。アイディアを要約す
れば次の通り。関数を c3 で展開して方程式を階層別に分解する時に、摂動解をベ
キ的に発散させる項（永年項）が現れるのだが、遅い時間変数 c3 t と速い時間変数t の２つの独立な時間変数があるとして、余計に現れる時間微分項をその永年項と打ち消させて発散を回避するのである。この操作により得られる方程式は繰り込
み群方程式と呼ばれるが、非摂動方程式の解をこの繰り込み群方程式に代入する
ことで、摂動の効果を無限に長い時間まで繰り込むことができる。粗い言い方をすれば、遅い時間変数の時計でみた各時刻で速い自由度の計算を行うということである。
それでは具体的な計算方法を以下で説明する。まず２つのスケールの異なる時
間変数を用意する。
t0 = t, t1 = c3 t

(3.14)

それに応じて密度関数も展開する。時間依存性が (t0 , t1 ) となっていることに注意したい。

ρ = ρ0 (t0 , t1 ) + c3 ρ1 (t0 , t1 )(3.15)

この展開式を運動方程式 (3.8) に代入する。

∂ρ0l.h.s =
+ c3
∂t0

27



∂ρ1 ∂ρ0
+
∂t0
∂t1


(3.16)

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

 −ic1∂e−ic1 ∂∂∂e
ρ0 g · ρ0 + c.c. − c3
ρ1 g · ρ0r.h.s = −2 ∂θ ∂θ
2 ∂θ ∂θ

−ic1
−ic˜2
∂
∂
e
e
∂
∂
ρ0 g · ρ1 +
ρ0 g · ρ0 + c.c. + O(c23 ) (3.17)+
2 ∂θ ∂θ2 ∂r ∂rwith g(r − r , θ − θ ) = exp [(−1 + iα)|r − r| + i(θ − θ)](3.18)

g を定義し直している 4 。· は空間と位相の両方に関するコンボリューション積分を意味する 5 。各次数で項を比較することで、摂動方程式の系列を得ることができる。0 次の方程式は次の通り。
e−ic1 ∂
∂ρ0∂
ρ0 g · ρ0 + c.c.
=−∂t02 ∂θ ∂θ(3.20)

この方程式は蔵本方程式と一致している。１次方程式は次の通り。

e−ic1 ∂
∂ρ1 ∂ρ0∂
ρ1 g · ρ0
+
=−
∂t0
∂t1
2 ∂θ ∂θ
−ic1
−ic˜2
∂
∂
e
e
∂
∂
ρ0 g · ρ1 −
ρ0 g · ρ0 + c.c.
−2 ∂θ ∂θ2 ∂r ∂r(3.21)

右辺の３つ目の項は解を１次オーダーで発散させてしまう。この発散は物理的ではなく、永年項と呼ばれるそこで次の繰り込み群方程式 (3.22) を仮定してこの永年項を取り除く。

e−ic˜2 ∂∂ρ0∂ρ0 g · ρ0 + c.c.
=−∂t1
2 ∂r ∂r(3.22)

空間変数のみの方程式が現れた。この方程式 (3.22) が空間移動を記述する。この方程式に S(t, r)ρ0 (t, r, θ) = S(t1 , r)ρ0 (t0 , r, θ) := ρ0 (t0 , t1 ) という変数変換を行うと、この方程式が粒子の移動を表していることがより理解しやすくなる。まず、変数変換により繰り込み群方程式は次のようになる。

∂∂S
e−ic˜2 ∂(t, r) = −c3
S(t, r) ζ[Sρ0 ](t, r)e−iθ + c.c.
∂t
2 ∂r
∂r
ここで、次の規格化 (3.24) が保存されることが分かる。
45(3.23)

位相因子を分けて書くと記号が煩雑になるためである。
ここで言うコンボリューション積分は次の意味で定義する。

g · ρ(t, r, θ) :=drdθg(r, r − r , θ − θ )ρ(t, r , θ )(3.19)

その他のコンボリューション積分もこれに準じて計算する。この積分を積と見なせば、分配律や結合律を満たす可換な積であることが分かる。

28

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

S1

dθρ0 (t, r, θ) = 1(3.24)

初期条件で右辺が 1 になるように選べば、いかなる時刻でもこの式が成り立つ。
また応じて、S についても次の規格化条件 (3.25) が保存される。

1
N

drS(t, r) = 1(3.25)

M

従って、S は粒子分布を表していることがわかる。ここで紛らわしく密度関数の繰り込み部分を S(t, r) と書いたのはその理由からである。
これで、位相変化と空間移動のダイナミクスの分解が完成した。２つの微分方程式が得られたが、まずは速い時間変数の方の方程式から解く。解き方は後に議論するとして、仮に ρ0 (t, r, θ) と解けたとする。この解を用いて記述される繰り込み群方程式 (3.22) を考える。この式の解 S(t, r) を用いれば、近似解 ρ(t, r) ≈ S(t, r)ρ0 (t, r)が得られる。これで方程式が位相部分と空間移動の部分に分解された。以下ではこの２つの方程式 (3.20)(3.23) の解き方を説明する。

29

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

3.3

自己無撞着方程式による平均場の決定

ここでは、位相部分の解き方について説明する。方針としては、定常状態での平均場 ζ(t, r) を仮定し、自己無撞着に決定する理論を作り解析する。位相部分は局所的に結合した蔵本方程式と一致しており、導出の方法は蔵本モデルでの導出とほとんど同じなので、この計算はスタンダードな蔵本モデルの計算を参照すれば理解しやすいだろう [34][37][38]。計算手順を要約すれば、各振動素子の振動状
態が位相ロックと位相ドリフトの２種類にわかれると仮定し、それに応じて分布関数をロックされた位相を指すδ関数と、ドリフト状態の不変密度分布に分ける。
そのように分けた分布関数を、平均場の定義式に代入すると平均場で閉じた方程式が得られる。以下、その導出を説明する。
次の形の方程式 (3.26) を目標の形とする。

ζ(r, t) = eiΩtdr e(−1+iα)|r−r | H[ζ](r , t)

(3.26)

以下、この形の自己無撞着方程式にどのように変形するのかを説明する。同期状態（位相ロック状態）(3.27) と非同期状態（位相ドリフト状態）(3.28) で定常分
布関数 P (r, ψ) with ψ = θ − Ωt が次のようになっていると仮定する。


P (r, ψ) = S(r)δ ψ + c1 − Θ0 (r) + sin
P (r, ψ) = S(r)

−1

Ω
R(r)



Ω2 − R(r)2 /2π
Ω + R(r) sin(ψ + c1 − Θ0 (r))(3.27)(3.28)

位相ロック状態の式 (3.27) は簡単で、ロックされた位相を指すδ関数に粒子分布の重みをつけた式を仮定する。位相ドリフト状態の式 (3.28) は連続化方程式の不変分布である。式の意味は明確で、位相の運動方程式の逆数となっていて位相が遅くなる所ほどに分布が偏っているという形をしている。この場合にも空間部分
は粒子分布関数を用いる。すると次のように H の形を決めることができる (3.29)。

H(r) =dψP (r, ψ)eiψ

(3.29)

位相ロック状態では次の式 (3.30) になる。δ関数が入っているので簡単に積分を実行できる。




Ω
−1
H(r) =
dψe S(r)δ ψ + c1 − Θ0 (r) + sin
R(r)
⎧
⎫

2 ⎬

⎨
Ω
Ω
− 1−
= S(r)ei(Θ0 (r)−c1 ) −i
⎩ R(r)
R(r) ⎭iψ

30(3.30)

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

ドリフト状態では次の式になる (3.31)。三角関数を含む積分になるが、留数定理
により積分が実行できる。ただし、ここの積分は Ω の符号によって極の選び方が変わるので注意が必要である。

Ω2 − R(r)2 /2π
Ω + R(r) sin(ψ + c1 − Θ0 (r))
⎫
⎧

2 ⎬

⎨
Ω
R(r)
1− 1−
= S(r)iei(Θ0 (r)−c1 )⎭R(r) ⎩
ΩH(r) =dψeiψ S(r)(3.31)

状況を２つに分けて、それぞれ異なる分布関数を仮定して積分したのだが、結局は全く同じ式の形をしているので、まとめて書くことができる。

⎧
⎫

2
⎨ Ω
⎬
Ω
−
H(r) = S(r )iei(Θ0 (r)−c1 )−1⎩ R(r)⎭
R(r)(3.32)

この H を平均場の定義式に代入することで、自己無撞着方程式を得ることができる。結果、次のように求まる。

R(r)eiΘ0 (r) = ie−ic1

⎧
⎫

2
⎨ Ω
⎬
Ω
 (−1+iα)|r−r | iΘ0 (r  )−S(r )e−1dr e⎩ R(r )⎭R(r )(3.33)
本来の目的ならばこの方程式の解析解が欲しいところである。なぜなら、特異摂動の無摂動解がこの自己無撞着方程式の解だからである。しかし、残念ながら
この超越方程式 (3.33) の解を手計算で求めるのは一般には無謀である。数値計算
によりこの方程式の解を計算する。しかし、この方程式は解が一意ではないため
このままでは求まらない。そこで解の形をある程度制限して解を計算することになる。
まず、定常状態を仮定する。さらに、我々の目標はスパイラル状の波を計算したいため、具体的には「スパイラルの中心を軸に 360 度回転した時に位相が 2π 変化する」という仮定をする (3.35)6 。この仮定は Cohen 氏の論文 [4] で初めて提案
され、Martens 氏らによりガウス型結合の Kuramoto モデルで適用され、キメラ 7の自己無撞着方程式による計算がなされた [16]。
6

この仮定の下では index 1 のスパイラルを計算することができる。逆無気にすれば index 2 のスパイラルも計算できる。初期条件を上手く指定すれば index の大きいスパイラルも存在し得る
し、本論文で用いた方法をわずかに変更するだけで計算できる。しかしながら凝集に対する効果は定性的には同じだろうから本論文では計算しない。
7
スパイラル状の位相波が発生している時に、中心付近に位相が特異的になっている領域が現れる現象をここではキメラと呼んでいる。

31

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

R(t, r) = R(|r|)(3.34)

Θ(t, r) = Ωt + Θ0 (|r|) + arg(r)(3.35)

更に、粒子分布を具体的に設定しなければ数値計算できない。最も簡単な粒子が一様に広がった一様分布を仮定する。この仮定 (3.35) の下、適当な関数の初期形を決めて繰り返し代入して収束させる逐次代入法により計算する。方程式 (3.33)
は Ω を適切に決めないと収束しないので、ここでの計算から Ω も同時に決めることができる。スパイラルを特徴付ける量は Ω と波長 λ であるが、波長は Θ0 (r) の関数形から計算することができる。自己無撞着方程式の収束結果 R(r)、Θ0 (|r|) を
図 3.1、図 3.2、図 3.3、図 3.4 に示す。図 3.6 に、Ω と波長 λ の c1 依存性を示す。さらに、分布 (3.27)、(3.28) を用いることで位相の分布を作ることもできる。それを図 3.5 に示す。
比較のため、等位相面に特異性の無いターゲット型の平均場も計算する。平均場の関数形を次のように仮定する。

R(t, r) = R(|r|)(3.36)

Θ(t, r) = Ωt + Θ0 (|r|)(3.37)

Θ が回転対称になっている点が特徴的である。この場合の収束結果 R(r)、Θ0 (|r|)
を図 3.7、図 3.8、3.9、図 3.10 に示す。スパイラル状位相波の場合と同様に、振動数 Ω と波長 λ の c1 依存性を図 3.12 に、位相分布を図 3.11 に示す。
方程式 (3.33) に c1 が含まれるが、結果を見ると R と Θ は c1 依存性がない。一方で、Ω は c1 については見たところ線形に依存している。実は、この Ω は c1 の sin関数であることが以下の計算から分かる。まず、振動数はスパイラルの中心から離れたエリアで決まると考え、r をスパイラルの中心から十分離れた位置に取る。
すると、位相波は平面波で近似できる。この事情はターゲット型の位相波でも同じことであり、結果は全く同じになる。具体的には次のような関数形を仮定する。

R(t, r) = R
Θ(t, r) = Ωt +(3.38)
2π
r · ewave + arg(r)λ(3.39)（ここで ewave ∈ R2 は波の進行方向を表す単位ベクトルである。どの方向を選んでも結果は同じ。）計算の便宜だが、cos χ := R/Ω とおくと、自己無撞着方程式
は以下のように変形できる。

32

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

R = ie

−ic1

∞

ρ0

2π





r dr dφe0= −e−i(c1 +χ) ρ0 

0(−1+iα)r 2πi(1 − iα))2(1 − iα)2 + ( 2π
λ




2π 
exp i − r cos φ − χ
λ3/2(3.40)

方程式 (3.40) の両辺を比較すると次の等式が得られる。

R = 2πρ0 (1 − iα) 2 3/2(1 − iα)2 + 2π
λ 2 
π
32π
0 = c1 + χ − + tan−1 α + tan−1 (1 − iα)2 +22
λ(3.41)(3.42)

この式では残念ながら情報不足で R、Ω、λ を決定することはできない。しかし、少し変形すると Ω が c1 について sin 関数であることを見いだすことができる。

Ω =(1 − iα) 2 3/2(1 − iα)2 + 2π
λ 2 
32π
× sin c1 + tan−1 α + tan−1 (1 − iα)2 +2λ

(3.43)

図 3.6 と図 3.12 にこの理論曲線を重ねて表示した。ただし、ここで作った Ω は
波長 λ の情報が欠落しているので、λ = 22 を代入して理論曲線を描いた。すると、理論とよく一致していて、この計算が妥当なことが分かった。

33

《修士論文》

R(r)

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

14
12
10
8
6
4
2
0

c1 = 0.2
c1 = 0.4
c1 = 0.6
0

20

40

60

80

100

r:distance from the center
図 3.1: スパイラル状位相波の平均場:R(r)。原点付近に位相特異点があり、原点を中心に位相スパイラルが回っている。この図はスパイラル位相波の絶対値を中心を通からある方向に切った断面図である。原点には S 1 一周分の位相面が寄せ集まっており、原点の非常に狭い近傍では位相分布がブロードした特異的な領域がある。
従って、同期の具合を表す秩序変数でもある R は 0 をとり、そこから離れるに連れて値が大きくなり次第には一定値に近づく。原点付近では、中心から離れる向きに正の強い勾配が生じている。本文で述べるように、この強い勾配が粒子の凝集を加速する役割をする。また、グラフの右端の中心から離れたところで波打ちは境界で起こる数値計算上の問題から生じた乱れであり物理的な意味は無く、無視すべき乱れである。実際には一定値に漸近する。

34

《修士論文》

Theta(r)物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

4
3
2
1
0
-1
-2
-3
-4

c1 = 0.2
c1 = 0.4
c1 = 0.6
0

20

40

60

80

100

r:distance from the center
図 3.2: スパイラル状位相波の解 Θ0 (|r|)。図 3.1 が平均場 ζ(t, r) の絶対値であるが、このグラフはそれと対になる平均場の位相部分である。この図も同じくある ζ の位相部分は各振動数 Ω(r) で振動するのだが、その時間変化分を差し引いたものがここで示している Θ0 (|r|) である。スパイラル中心から外に離れるに従って一様に減少する関数である。つまり波長は位置によらず一定である。中心の位相の値が 0になっているが、ここの値は時間の原点の定義の仕方によってどんな値にもできるため全く意味がない。見やすくするために 0 に合わせた。図 3.1 と同じように、
グラフ右端で振る舞いが変わっているのは数値計算上の問題で無意味である。

35

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 3.3: スパイラル状位相波の解 R(r) の２次元プロット。平均場 ζ(t, r) の断面の絶対値成分を図 3.1 で描いたが、この図はその断面をとる前の２次元的な R(r) である。回転対象な形をしている。中心で 0 をとり、中心から離れたエリアでは一定値をとっており、図 3.1 のグラフが中心から離れたところで平らになっているのと対応している。ここでは中心を見やすくするために図 3.1 よりも狭い領域を描いている。そのため、この図には遠方で起こってしまった数値計算の無意味な乱れは描いていない。

36《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 3.4: スパイラル状位相波の解 Θ0 (|r|) の２次元プロット。図 3.3 と対になる平均場 ζ(t, r) の位相部分である。原点を中心にスパイラルを巻いている。原点周りを一周回ると位相が 2π 変化する。実際には原点に位相特異領域があるのだが、平均場ではその領域は見えない。この図も図 3.3 と同じく、図 3.1 と図 3.2 よりも狭い領域の平均場を描いている。

37

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 3.5: スパイラル状位相波の解 Θ0 (|r|) から構成した、粒子の位相分布の２次元プロット。ほとんど Θ0 (|r|) と同じに見えるのは、粒子は外場に c1 に応じて決まるわ
ずかな遅れをとって同期しているからである。平均場との特徴的な違いとしては中心の近傍に周りと不連続な狭い円形の領域が見えることである。このエリアでは位相が乱れていて分布関数 P (φ) が式 (3.28) となりブロードをもつ。このプロットで示したこの位相特異領域での値は P (φ) のピーク値をとっている。

38

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

25

4

wave length

20

3
2.5

15

2
10

1.5
wave length
phase velocity
p.v. (theory)

5

1
0.5

0

0
0

0.2

0.4

0.6

0.8

1

c1
図 3.6: スパイラル状位相波における Ω と波長の c1 依存性。Ω を正しい値にチュー
ニングしないと平均場を決める自己無撞着方程式が収束しないので、これを解く
時に同時に手に入る。λ は Θ0 を描けば、その傾きから計算することができるのでこれも手に入る。Ω は c1 に関して線形的に増加し、Ω は c1 には依存していない
ことが分かった。Ω の c1 依存性については式 3.43 により解析的に説明することができる。一方 λ は空間スケールに依存する量なので、モデルの設定で相互作用がexp(−|r|) でダンプするようにしている限りこの量は変わらないと考えられる。また振動数については理論曲線を重ねてプロットした。よく一致している。
（なお図中の p.v. は phase velocity の略。）

39

phase velocity

3.5

《修士論文》

R(r)

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

13
12
11
10
9
8
7
6

c1 = 0.2
c1 = 0.4
c1 = 0.6
0

20

40

60

80

100

r:distance from the center
図 3.7: ターゲット状位相波の平均場 R(r)。ターゲット状位相波の場合はスパイラルの場合と違って、中心に位相特異点が無い。むしろ、中心付近では位相が同じエ
リアが円形に形成されるので、同期の程度は強くなり秩序変数 R(r) は遠方より大きい値をとる。中心から離れると、スパイラルの場合とほとんど状況は同じになるので一定の同じ値に漸近していく。ただし、ターゲット型の場合は漸近が上からであって傾きが負になっている点が注目に値する。本文中で述べるように、この傾きは位相勾配の凝集駆動を妨げる働きをし、ターゲット型の位相波で粒子が集まりにくい原因となる。このグラフの右端の乱れは、以前に示したグラフと同じように非物理的な計算上の問題でおこる乱れであり、無視して良い。

40

《修士論文》

Theta(r)物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

4
3
2
1
0
-1
-2
-3
-4

c1 = 0.2
c1 = 0.4
c1 = 0.6
0

20

40

60

80

100

r:distance from the center
図 3.8: ターゲット状位相波の解 Θ0 (|r|)。ターゲット状位相波もスパイラル状位
相波（図 3.2）と特に大きな違いは見られない。単調に直線的に減少する関数形をとっている。この場合もターゲットの中心の位置が 0 ∈ [−π, π] となるように位相を合わせている。この場合も同様に c1 の依存性は見られない。

41

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 3.9: ターゲット状位相波の解 R(|r|) の２次元プロット。ターゲット状位相波の絶対値成分（図 3.7）を２次元的にプロットしたもの。当然、回転対象な形をしている。ターゲットの場合はどこにも特異性が無いため、中心が少し大きな値をとって、中心から離れるに従って単調になだらかに一定値に上から漸近する形をしている。図 3.3 の場合と同様に、中心付近を拡大した図である。

42

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 3.10: ターゲット状位相波の解 Θ0 (|r|) の２次元プロット。ターゲット状位相波の位相成分（図 3.8）を２次元的にプロットしたもの。回転対象な形をしている。
図 3.4 の場合と同様に、中心付近を拡大した図である。スパイラルの場合（図 3.4）と比較してなだらかな形をしているように見えるが、それはスパイラルが回転方向の勾配が急なためであって、凝集に効く動径方向の勾配は同じである。

43

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 3.11: ターゲット状位相波が生じているときの位相分布の２次元プロット。スパ
イラル位相波が発生しているときの粒子の位相分布（図 3.5）と同様の方法で構成した。図 3.4 や図 3.5 と同様に中心付近を拡大して描いている。スパイラルの場合（図 3.5）の場合に見られた中心の特異的な領域はもちろん見られない。

44

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

25

4

wave length

20

3
2.5

15

2
10

1.5
wave length
phase velocity
p.v. (theory)

5

1
0.5

0

0
0

0.2

0.4

0.6

0.8

1

c1
図 3.12: ターゲット状位相波における Ω と波長の c1 依存性。この結果はスパイラルの場合（図 3.6）とほとんど同じに見える。この場合にも理論曲線は、自己無撞着方程式の計算結果とよく一致している。（なお図中の p.v. は phase velocity の略。）

45

phase velocity

3.5

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

3.4

細胞のフローの計算

前章で計算した平均場 (Θ(r), R(r)) を繰り込み群方程式に代入することで、細胞のフローを表すベクトル場を計算することができる。繰り込み群方程式 (3.23) を連続の式と見ることで、速度ベクトル場 v(t, r, θ) が次の式でかけることが分かる。

∂
ζ[Sρ0 ](t, r)e−iθ + c.c.
∂r∂
= −c3 R(t, r)ei(Θ(t,r)−θ) + c.c.
∂r∂R∂Θ(t, r) + R(t, r)iei(Θ(t,r)−θ)(t, r) + c.c. (3.44)= −c3 ei(Θ(t,r)−θ)∂r∂rv(t, r, θ) = −c3

変数 θ については、自己無撞着方程式を導出する過程で作った分布関数で平均をとることで、実際に起こるであろうフローを計算することができる。

v(t, r) =dψP (r, ψ)v(t, r, Ωt + ψ)

(3.45)v(t, r) の計算結果をスパイラル状の波、ターゲット状の波、それぞれの場合について図 3.13、図 3.14 に示す。弱くスパイラルを巻きながら中心に向かって移動
し、中心付近で強く引き込まれているのが分かる。
このメカニズムは次の通り。式 (3.44) を見ると分かるように、速度ベクトル場
は平均場の絶対値部分 R の勾配と位相部分 Θ0 (r) の勾配に比例する項の和で出来ていることが分かる。スパイラル状の位相波では中心から離れたエリアでは R の
勾配はほとんど無視でき、位相から定常振動を差し引いた場 Θ0 (r) の勾配からのみ決まる。この事情はターゲット状の位相波でも同じことであり、実際２つの図
3.13、図 3.14 を比較してみても、中心から離れたエリアのベクトル場はほとんど同じである。ところが、中心付近では事情が異なる。スパイラル状位相波の中心には位相特異点がありこの近くでは R(r) が外向きに正の傾きで急激に変動する。
符号を注意深く吟味する必要があるが、考えているパラメータの条件では R(r) の勾配とは逆向きにベクトル場が発生するので、スパイラルの中心付近では強く引き込みが起きている。一方、ターゲット状位相波の中心には特異点は無く、むしろ周辺よりも中心が強く同期する。従って、R(r) は中心付近では外向きに緩く負の傾きに変動しており、R(r) の勾配は中心から湧き出す向きの粒子フローを作る。
これが Θ0 (r) 勾配の粒子フローと打ち消し合って、中心付近での引き込みは非常に弱くなる。
ただし、多くの先行論文では cAMP 勾配と位相勾配がほとんど同一視され、位相勾配が粒子のフローを作ると考えられている。この考え方は平均場として位相と絶対値を考える私のアイディアとは異なる。位相勾配のみに着目する立場ではスパイラル状位相波もターゲット状位相波も状況はほとんど同じであり、粒子の46

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

引き込みが同様に起こると考えるのが定説である。しかしターゲット状の波で集まるという報告は未だ出ておらず 8 、本研究の結果と矛盾する報告は無い。本研究の結果は定説を破り、かつ未だ実験との矛盾が無い、理論計算から得られた全く新しい成果である。ただ、D.discoideum に対してこのことを確認した実験の論文は未だ存在しない。
今後、なんらかの実験でターゲット波で凝集が起こらないことが証明されることを期待する。

8

凝集を観察した実験は膨大にあるが、例えば [1][26][14][2] などを読めば分かるように、撮影されているのはスパイラルを巻いた凝集のみであり、ターゲット型の波を伴う凝集を報告する論文は見当たらない。

47

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

y

15
10
5
0
-5
-10
-15
-15-10 -5 0 5 10 15x

図 3.13: スパイラル状の位相波が発生しているときの粒子フロー。中心から離れた所では弱く渦を巻きながら引き込まれる方向の流れが生じていて、中心付近で強く引き込まれている。

48

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

y

15
10
5
0
-5
-10
-15
-15-10 -5 0 5 10 15x

図 3.14: ターゲット状の位相波が発生しているときの粒子フロー。中心から離れた所では弱く引き込まれる向きに流れていて、中心に近づくほどに、流れが弱くなっている。

49

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

第4章

安定性解析による田中・巌佐
チャートの解読

前章ではパラメータ設定をある領域に限定し、その設定で具体的にスパイラルの解が存在について、解の構成による直接的な証明を与えた。さらに解がどのような性質をもっているかを詳しく調べた。これがここまでに説明した成果であり、既に凝集過程については完全に理解できたように思える。しかし、「なぜこのパラ
メータ設定でスパイラルが起こるのか」という疑問には一切答えていない。この章ではその疑問に対する答を与える。
まずはパラメータ設定に関する先行研究の説明をする。Swarm Oscillators が提唱されてすぐ、[11] において、パラメータ空間 (c1 , c2 ) ∈ S 1 × S 1 の大域的な調査が行われ、パラメータ空間を１３種類の領域に分類した相図（以後チャートと呼ぶ）が作成された（図 4.2）。その研究により、位相、空間分布に関して非常に豊かな現象がこのモデルに含まれることが明らかになった。その後、２粒子の厳密解から導かれたフローダイアグラムを使って、チャートの様子が定性的に説明された。以上が先行研究の内容である。

50

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.1: 田中・巌佐チャートでの分類。粒子の動き、空間配置、位相同期から分類している。A は多数の粒子が１点に同期して集まって静止。B は少数の粒子が位相が離れた状態で集まって静止。CD は規則的な配置をとって同期して静止。E は規則的な配置で一方向に流れながら、位相も同じ方向に粒子のフローとは異なる速度で伝播。FG は特徴的な配置で同期して静止。H は線状の配置をとってうねりながら動く。I は３角形の配置でキャッチボールをするように粒子が流れる。J は
小さな円を作って回る。K は位相が近い粒子が近づいては離れるを繰り返す。L は
51
位相が遠い粒子が近づいては離れるを繰り返す。
M は局在しながらも、不規則にクラスター内部で動き回る。同期も起こらない。HIJKLM は同期しない。

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.2: 田中・巌佐チャート。

52

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

しかし、たくさんの（103 オーダーの）粒子の数値計算で相図の境界を詳しく調べてみると、その境界は不明瞭で２粒子の議論を見直す必要があることが示唆された。多体系の境界は２粒子系よりも複雑な物になっているのではないかと私は直感し、前章で導いた連続 Swarm Oscillators(3.8) を使って、チャートの再現を目指すことにした。方程式 (3.8) を使って、密度関数を粒子が平面状を一様に分布し、位相が全く揃っていない位相も一様な分布をした状態の周りで展開し、線形安定性を見た。波数依存性を調べることで、長距離秩序発生による一様状態不安定化
が位相引き込みに対応し、短距離秩序の発生による一様状態の不安定化が粒子間引力に対応することが分かった。またたとえ連続化された系であっても、１粒子の視野半径に対応する α を小さくした時に巌佐氏の結果と一致することが分かった。
これは我々の解析がある粒子を取り囲む周りの粒子全ての効果を足し合わせた物であり、巌佐氏の解析が１つ隣までの効果を取り入れたものであったため、視野を小さくして結果が一致するのはとても自然なことである。我々の結果は巌佐氏の２粒子の解析結果をサポートし、さらに強い主張を与えている。
また、「スパイラル波を発生させる凝集が起こると考えられるパラメータ設定は何か」、すなわち「位相引き込みと凝集が起こるパラメータ領域がどこか」という問題にも答えることができる。以下ではまず、巌佐氏の結果をレビューし、その後で解析方法とそこから得られた結果を説明する。

53

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

4.1

巌佐氏の２粒子の厳密フローによるチャートの定性的解読

巌佐氏は２粒子の方程式のフローが厳密にかけることを利用し、粒子同士が接近した時にどのような振る舞いが起こるかを８種類に分類した。なお、ここでの内容は彼とのディスカッションで教えて頂いたアイディアである [18]。彼は、相対座標 δ = θ2 − θ1 、ρ = |r2 − r1 | を用いることで２粒子系を１つの方程式に書き直した。

δ = −2e−ρ cos(αρ − c1 ) sin δ(4.1)

ρ = −2ec3 e−ρ sin(αρ − c2 ) cos δ(4.2)

このように書くことで、δ-ρ 平面にフローダイアグラムを描くことができる。
（図4.3) cos(c1 − c2 ) の符号で不動点の安定性が変わるので、２つの状況に分けて図を描いている。特にそこに示しているのは左側の図がが c2 = c1 の時のフロー、右側の図が c2 = c1 − π の時のフローである。

図 4.3: ２粒子の方程式から得られるフローダイアグラム。左側の図が cos(c1 −c2 ) >0 の時のフローで特に c2 = c1 の時のフローを描いている。右側の図が cos(c1 −c2 ) <0 の時のフローで特に c2 = c1 − π の時のフローを描いている。
まず既に、(case a) cos(c1 − c2 ) > 0、(case b)cos(c1 − c2 ) < 0 で様相が全く異なるフローが現れることがこの図から分かっている。では、フローダイアグラムが

c1 、c2 の設定でさらに詳しくどのように変化するか考える。Swarm Oscillator では54

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

相互作用が指数関数的にダンプするので粒子が空間的に十分に接近した初期条件で考える。従って、このダイアグラムの δ 軸の近くを注目すれば良い。２つの図は c1 と c2 を同じ値だけ動かすと、フローが αρ 軸方向にシフトする。見方を変えれば、粒子のゼロ距離状態を表す δ 軸が αρ 軸方向に上下に動くと考えることもできる。このような考え方で巌佐氏はフローの分類を行った。
まずは、(case a) を考える。図は c1 = c2 でそれぞれ正の小さい値をとっている
設定でのフローを示している。同期状態を表す点 (αρ, δ) = (0, 0) は空間方向に不安定で、空間方向に少し離れた点 (αρ, δ) = (c2 , 0) へとフローは流れ込んでいる。
従ってこのパラメータ設定の時は２粒子は空間的に少し離れて位相同期する。このままフローを上向きにシフトさせていく、すなわち c1 = c2 の値を大きくしてい
くことを考える。すると、c1 = c2 = π/2 のところで、δ 軸の下にあった中立不動
点が顔を出し、フローが質的に変わる。c1 = c2 を破って c1 、c2 を独立に少しずつ動かしてもフローの質的な性質は変わらない。どこまで変わらないかを考えると、次の領域で 2 粒子系は質的に同じ振る舞いをするだろうと考えられる。

{−π/2 < c1 − c2 < π/2} ∩ {c2 < π/2}(4.3)

この領域はチャートの H、I 領域に当たる。このようにして、分類していくと

case-a,b がそれぞれ４種類に分類される。従って全部で８種類に振る舞いが分類される。これ境界を c1 × c2 図に示したのが図 4.4 である。
十分に接近した２粒子ならどんな α を選んでも振る舞いが変わらないため、２粒子解析では α の意味はこの議論では分かりにくい。ただ、多体系では α の選び方により振る舞いが変わる。ちなみに、A、B、· · · 、M の分類がなされた時の値
は α = 1.0 であった。おおかた境界の位置を説明できているが、例えば C と B の境界（c1 = π/2 に見える境界）が検出できないなど微妙な問題がある。次の節で紹介する多体効果を入れた理論でこの問題を解決できる。

55

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.4: ２粒子解析から導かれた分類。

56

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

4.2

線形安定性解析によるチャートの解読

蔵本モデルにおいて、一様に非同期な状態は線形安定であることが知られている。それを踏まえ、連続 Swarm Oscillators を一様非同期な状態の周りで展開して、線形解析を試みる。ρ(t, r, θ) = ρ0 /2π + δρ(t, r, θ) ように分布関数を展開する。これを式 (3.8) に代入し、２次以上の項を捨てる。さらにフーリエ変換をすると次の式を得る。

δρ(t, k, 1) = L(k)δρ(t, k, 1)(4.4)

ρ0 Δc (k, 1)(1 + iα)(4.5)

L(k) =

3(|k|2 + (1 + iα)2 ) 2c˜3 e−ic˜2 2 e−ic˜1 2
|k| +
n
2
2
係数 L の実部が不安定性を決める。図の見方は図 4.5 を参照せよ。不安定化する所をその波長に応じた色で表示する。不安定モードが存在しないエリアではプ
ロットしない。白い領域がそれである。短波長程に黄色、長波長程に紫や黒になるグラデーションで表示した。
結果は α に依存していて、２粒子解析の結果と、田中・巌佐チャートの食い違いを上手く埋める結果を得た。というのは α = 1.0 の時はチャートと整合性のとれた結果となり、α がそれより小さい値のときに２粒子解析の結果に近づくことが分かった。そのメカニズムについては後で説明するとし、先に結果を述べる。まずは先行研究で用いられていた値 α = 1.0 の場合について説明する。この場合の線形安定性解析の結果は図 4.6、図 4.7 に示した通り。図を見ると、領域 {0 < c1 < π/2}において短距離モードが不安定化することが分かる。これはこの範囲の c1 を選ぶと近くの粒子同士が引き込まれて一様分布でなくなる。特徴的なのは図の右下に見える三角形状の白い領域で、この領域では同期も集合もしないのだが、境界の位置をよく検出できている。
一方、α = 0.1 の図 4.8 を見ると、領域 {0 < c1 < π/2, 3π/2 < c1 < 2π} で長距離秩序の不安定化が、{π < c2 < 2π} で短距離秩序の不安定かが起こることが分かる。実際２粒子解析では、前者の領域は粒子が引き合う領域、後者は位相引き込みが発生する領域で巌佐氏の結果とも整合する。そして領域 {π/2 < c1 <3π/2} × {0 < c2 < π} は同期も粒子間引力もない。ここの領域もよく検出できている。
ただ、巌佐氏は一般の α でそうなると主張しているが、説明したように我々の結果はそうではなかった。すなわち α が小さい値の時に 2 粒子のフローの結果が正しいことを線形安定性解析の結果は示している。本論文には示さないが、我々の確認した限りではここでの主張に反する結果は数値計算では確認されていない。
Δc (k, n) :=57

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

では、このような α に依存した結果になったメカニズムを説明する。２粒子解
析の結果である図 4.3 で説明する。２粒子の振る舞いを調べるには横軸（δ 軸）の近くの振る舞いが重要であることを前に説明した。多粒子系では各粒子の見える範囲（相互作用が及ぶ範囲）が exp(−|r|) で特徴付けられることから（すなわちr < 1 程度まで見えることから）、巌佐氏は多粒子においても、２粒子での結果がそのまま適用できるだろうと考えたのである。ところが、このグラフ（図 4.3）の
縦軸は αρ であり、このグラフで考えるべき範囲は αρ < α である。α が小さいときは確かに δ 軸近傍だけを見れば良いのだが、α が大きい時（α = 1 は決して小さい値ではない）δ 軸から α 程度離れた所まで見えてしまうのである。従って、この場合には２粒子解析の結果が成り立つことは期待されず、実際２粒子解析の結果
は α = 1.0 の時に作ったチャートとはあまり一致していない。
さて、線形安定性解析の結果はどうだったか、α = 0.1（α が小さい時）は２粒
子解析の結果と整合し、α = 1.0 の時はチャートと整合した。これは今述べた説明とよく整合している。線形安定性解析は注目している点から見える全ての粒子の効果を繰り込んだ効果を見ることができ、先の議論から言えば距離 α 程度の中に入った粒子の効果を見ていることになる。

58

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.5: 安定性解析の数値計算の見方の説明。最も不安定なモードの値を色にしてプロットしている。色に応じて波長が示されている。黒や紫色は長波長モードの不安定化、黄色やオレンジは短波長モードの不安定化を表している。白い所は不安定モードが存在せず、一様に非同期な無秩序状態が安定なパラメータエリア。このエリアでは粒子の位相波乱雑なまま同期せず、そして粒子の凝集も起こらない。
なおここに示したのは α = 1.0 での安定性解析の結果。

59

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.6: α = 1.0 での安定性解析の結果に粗くした田中・巌佐チャートの分類を重ね合わせて描いたもの。多粒子にして区別がつかないと判断される線を削除した。
線形安定性解析は無秩序状態の安定エリアをよく検出している。また図で言う左半分 (0 < c1 < 2π) のエリアでは短距離モードが不安定化している。このエリアでは粒子が互いに引き寄せ合う現象が見られる。また、上半分では長距離モードの不安定かが見られるが、このエリアでは粒子の位相が引き込み合う現象が見られる。

60

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.7: α = 1.0 での安定性解析の結果に２粒子解析による分類を重ね合わせて描いたもの。線形安定性解析の結果で c1 = π のところに境界線が現れているが、２粒子解析の場合にこの境界は無い。両者の結果はあまり一致していないことが見
て取れる。

61

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

図 4.8: α = 0.1 での安定性解析の結果に２粒子解析による分類を重ね合わせて描いたもの。図の下半分の中央にある、無秩序状態の安定エリアを見事に検出している。確かにこのエリアでは粒子が集まったり、位相が揃ったりすることは無い。
なお、線形安定性解析の結果の右側の境界が曲がっているが、α の値をもっと 0 に近づけることでシャープになる。

62

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

第5章

キイロタマホコリカビの凝集
フローの観察速度ベクトル場が予想通りに弱いスパイラルを巻きながら中心を向き、中心付
近で急に加速していると理論では予想している。しかし、百聞は一見にしかずである。アメーバが加速する所を実際に見てみたい。それが実際の系ではどうなっているのかを確かめるために、D.discoideum を使って実際に凝集過程のフローを観察した。結果、我々はスパイラルを描きながら凝集をする様子を撮影することに成功した。さらに我々は動画を連続的な画像データに分解してそれに対して、複数の画像の間の相関をとって速度を解析する PIV というアルゴリズムを適用して速度ベクトル場を計算した。技術的にはとりわけ新しい技術は無いのだが、過去に何か具体的な目的を持って絶対値を含めて速度ベクトル場を解析した研究は無く、実験においても本論文が第一報となるだろう。この章ではまず実験の手順を説明し、その後で結果の解析方法を説明、さらに結果の考察を行う。

63

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

5.1

実験の手順

5.1.1

D.discoideum の培養について

D.discoideum はバクテリアを餌に成長する。まず、KK2 バッファーに D.discoideumを撹拌し、少量のバクテリア（我々は K.aerogenes を用いた）を混ぜて寒天栄養培
地に広げる。我々は栄養培地として SM 培地を用いた。我々は先行研究のプロトコ
ル [3] を参考にして SM 培地を用いたが、最近よく用いられる LB 培地を用いても問題ない。２日程度の間、温度 20◦ C で湿度が適切に保たれた環境で放置すれば、分裂が進み細胞を培養することができる。我々が用いたインキュベータは湿度管理が出来ないため、水を入れたコップをインキュベータ内に置いた。湿度に関してはこの程度の調整で培養できる。
なお、K.aerogenes はグルコースを餌に増殖するため、グルコースが入っている栄養培地を使えば増殖し続けるため餌が不足することは無い。ちなみに K.aerogenes
は液体 SM 培地に混ぜておくだけで培養し続けることができる。
SM 培地
Glucose:5.0g, Bacto Peptone (BD REF No.211677):5.0g, K2 HP O4 :0.5g, KH2 P O4 :0.75g,
M gSO4 :0.25g, Agar:10.0g, Distilled H2 O 0.5l
KK2 バッファー
K2 HP O4 :20mM, KH2 P O4 :20mM
K.aerogenes、D.discoideum(AX2) は NBRP 細胞性粘菌中核機関提供の株を使用
した。

5.1.2

凝集の観察準備と観察培養した D.discoideum の細胞を白金耳で少量拾い、エッペンドロフチューブの中で KK2 バッファーに撹拌する。このときにボルテックスなど強い刺激がかかる
撹拌をすると細胞が弱るため、撹拌はピペッティングで行う。また、遠心器による沈殿操作も細胞に負担をかけるため、行わない方が良いデータを得ることができる。セルカウンターを使って細胞数を数え、108 (/ml) 程度の細胞密度に調整する。これを寒天培地に KK2 バッファーを固めた寒天培地（もしくは SM 培地）の表面に広げ、2 時間程度放置すると細胞が集まり始める。ただし培地に広げるときの表面の水分の具合には繊細な調整が必要で、適切な状態を作らないとうまく集まらない。多すぎては細胞が流れてしまうし、少な過ぎては細胞が動かなくなっ64

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）てしまう。バッファーのほうが透明度が高いが、簡便のため SM 培地を用いても同
じデータを得ることができる。

40 倍から 100 倍程度の拡大率の顕微鏡で観察するとうまく撮影することができる。寒天の裏からでも表からでも撮影可能である。ただし、撮影する時にはシャーレのふたは空けない方が良い。なぜなら、撮影にはうまく凝集のタイミングに始めても連続で 3 時間以上かかり、ふたが開いていると寒天が乾いて厚みが薄くな
りピントがずれてしまうからである。また、更に良い画像を撮るための工夫としては、シャーレのふたの内側の端の方にに滅菌水で浸したガーゼをはり付けると内部の湿度を保ち乾燥を抑えることができるので、先に述べたピントのズレを抑
えられる。

5.1.3

画像解析の方法

撮影した画像データに、連続する画像の相関から粒子の速度を評価する PIV（Par-ticle Image Velocimetry）というアルゴリズムを使って、粒子が動く速度を定量的に評価した。PIV の原理については Appendix で紹介する。我々は動画でも撮影したが、適切な時間間隔（我々の考えている系では 10 秒から 30 秒程度）で静止画サンプリングをすれば同様に解析ができる。
速度評価の精度を上げるためにはサンプリングは細かくとるということも考えられる。確かに PIV の相関に使う画像数も増えるので精度が増すだろうが、同時に解析時間が長くなり扱うデータの大きさも莫大になる。細胞の移動速度や実験精度を考慮し試行錯誤した結果、我々の目的には上記のサンプリング間隔で最適と考えられた。
PIV には相関をとるメッシュサイズの調整と最大・最小のしきい値設定によるノイズ除去のチューニングをし、細胞のフローをベクトル場で表現することができた。その結果は次の章で説明する。
また、各細胞の重心を各時刻で計算して速度を評価する PTV（Particle Tracking
Velocimetry）という方法もあるが、D.discoideum の凝集を見る場合は細胞が明確に分離して見えないため、我々の目的に対しては PTV は適切ではない。

65

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

5.2

実験データとその解析結果

以上のような方法で観察をした結果、スパイラルを描きながら凝集する様子が見える連続的な画像を撮影することができた。図 5.1 にそのキャプチャを示す。ス
パイラルを巻きながら集まっている様子を確認することができる。ここには２つのスパイラルが映っていて、その２つの特異点を中心に細胞が回転移動している様子が映っている。加速と減速を繰り返しながら集まる様子が見られた。

図 5.1: D.discoideum の撮影（左側）と、それに PIV 解析をかけた結果（右側）。
実際には連続的に画像を撮影していて、そこから PIV によりベクトル場を作成した。ここに示したのはある時刻でのキャプチャである。PIV の結果を見ると、画面内には２つのスパイラルが発生している。PIV のベクトル場が小さく表示されているため矢印の向きは見にくいが、これは中心方向を向いている、速度分布がスパイラルを描いている様子を確認できる。赤い矢印は黄色に比べて速度が速いことを表しているが、中心付近にいくらか確認できる。

5.3

結果の考察と課題

スパイラルを巻きながら集まるという描像は正しかった。スパイラルと凝集の明らかな時定数のズレをも確認することができた。さらにスパイラルの勾配に比べて進行方向の巻きが緩い点もモデルと整合していた。微妙ではあるが、中心付近の加速も見える。

66

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

モデルと異なる点としては、実際のカビは加速、減速してを繰り返しながら集まっていたことだろう。この振る舞いは cAMP がパルス状に伝播していて、パルスが届いたタイミングで前に進むために見える。Swarm Oscillators は内部振動の振動を単純化して考えているため、縮約操作の時にこの加減速の情報が落ちたのだろう。
我々の撮影での課題点は、中心付近の加速の検証が微妙になってしまったことである。ただし、我々の理論は中心付近で粒子が加速することを示唆しているため、それが確認されれば過去のあらゆるモデルから一歩前進したことになる。しかし、実験のセッティングが非常に難しく、中心付近の速度まで測定できたのは本論文に載せたスナップショットの動画しかないため、このデータだけから中心付
近での加速を主張するのは苦しいだろう。だが PIV の結果でスパイラルの中心付
近には速い速度ベクトルがわずかに見えているので、このデータを楽観的に見るならば、同じ解像度のデータを繰り返しとり続ければ求めている速度分布が見えるだろうと期待できる。今後、周到に計画された実験でスパイラル中心での加速が証明されることを待望する。

67

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

5.4

結論と課題

ここまでの議論で、Swarm Oscillators の研究を通して、D.discoideum がどのようにして集まるのかを説明し、そしてこのモデルがパラメータに応じてどのような振る舞いをするのかということを詳しく説明してきた。結果、cAMP を放出す
るタイミング（位相）の波が螺旋状に伝播し、結果的に生じた cAMP の波が粒子を１点に集める向きに駆動するというシナリオを数理モデルで再現した。
さらに、反応拡散系で現れる別の位相波では集まらないのかという疑問に対しても答を出した。ターゲット型の波を作って粒子の動きを見たところ、中心付近
で粒子が集まりにくくなることが分かった。数理モデルの中の抽象的な言葉で解釈するならば、平均場の位相因子の勾配と絶対値因子の勾配の２つのファクター
で粒子の駆動が駆動されるため、中心で同期が強くなるターゲット型位相波の場
合には中心付近で絶対値因子の勾配が凝集を妨げたということになる。
実際の系で何が起きているかという、より現象に近い言葉で解釈しよう。ターゲット位相波では中心付近での同期が強くなることが分かった、これは現象としては中心付近で cAMP の放出タイミングが一致するということを意味する。一斉に中心付近での cAMP パルスが放たれて、まず中心に円形に cAMP 濃度が高い領域が現れ、それを感じ取った周りの D.discoideum から cAMP が放出され、円形がドーナツ型になって広がる。それが中心から遠くに次々と伝播し cAMP パルスが発生する。そのパルスを個体が感じ取ってその伝播方向に逆らう向きすなわち中心向きに動く。cAMP 放出の波が中心から離れた所へ伝わっていく。細胞は発生した cAMP 濃度上昇のパルスを感じとって中心に向かって移動することができる。
だが、細胞が中心付近に近づくとそこには同じタイミングで cAMP 濃度が振動的に変化する円形の領域があり、そこでは cAMP の伝播を感じることができない。
細胞はどちらに動けば良いのか迷ってしまう。
一方、スパイラル状の位相波が出ているときは、特異点近くで平均場の絶対値部
分が急激に減少するため、その勾配が粒子の移動を更に加速させていた。これを現実の系でのシナリオに焼き直そう。中心から離れた所ではターゲット状位相波の場合と同様に、cAMP 濃度勾配のパルス的伝播を感じて細胞は中心へと移動する。そして中心近くにたどり着いてからスパイラルの性質が顕著になる。等位位相線が中心まで束なるように伸びていて、中心まで cAMP の勾配があるため細胞はcAMP パルスの波を感じながら中心まで集まることができる。ここで注目したいのが、等位相線は中心に近づくにつれて曲率を増していくが、理論で描いた速度ベクトル場はそれほど急なカーブを描いていないということである。これは等位相面が cAMP 放出パルスの形と一致しているという仮説が、スパイラルの中心付
近で破れていることを示唆している。つまり、等位相線は急な角度で中心に巻き込まれているが、cAMP の勾配の形は違うのではないかということ。本論文の理論

68

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

は平均場から cAMP 濃度分布を構成する逆問題までは至っていないので、cAMPの勾配の形については具体的に言及できない。ただ、見かけ上 D.discoideum はス
パイラル状の平均場の波を感じながら中心へと集まるように動くのである。
ここで本論文の冒頭で提示した問題 (X2) に対する答をまとめる。我々は D.discoideum
が cAMP の濃度と濃度勾配を感じて動いているということを知っている。しかし、それならばスパイラル状の位相波であってもターゲット状の位相波であっても細胞は集まることができるはずであるし、むしろ、回転方向の余計な動きの無いター
ゲットの方が理にかなっているように思える。しかし現実は違う。細胞は常にス
パイラルを描きながら集まるのである。両者の違いは平均場により特徴付けられた。ここに (X2) の対する答がある。
「D.discoideum は cAMP の濃度勾配から自分の周囲の個体の状態、すなわち平均場で駆動される。」
次のようなより強烈な言い方もできる。
「D.discoideum はある点での cAMP 濃度とその勾配という非常に局所的な情報から周囲の位相を積分した平均場を読みとり、その勾配に従って動いている。」
本研究において cAMP の濃度と勾配の情報だけで動く Chemotactic Oscillators
を、平均場の勾配で動く Swarm Oscillators に帰着させることで、凝集のフローが簡単に理解できることが分かった。いわば Chemotactic Oscillator で見えていた奇妙な凝集の振る舞いが、Swarm Oscillators で見ることで分かりやすいシンプルな現象に化けたのである。見かけ上、細胞は周りの化学物質の濃度や勾配の情報をもとに、周囲の細胞の状態を感知していると考えられるのである。
しかし D.discoideum に cAMP の分布と勾配からソースの情報を復元する逆問題が解ける高度な情報処理機能が備わっているということではない。ただ興奮系特有の位相引き込みの性質と走化性が備わっているだけである。平均場勾配の認識はこの２つの性質が絶妙に絡みあって実現した機能である。
今回扱った数理モデルは複数のダイナミクスが絡んでいて、力学系（反応拡散系）の問題として十分に楽しめる素材だった。まだいくらでも複雑な挙動が出てくるため、理論的な探求の可能性には際限がない。ただ、D.discoideum の凝集プ
ロセスの理論的な理解という面では満足な到達点まで来たと考えている。凝集プ
ロセスの理論は実利的な可能性もあると私は期待している。
そこで次の２つの課題を後続する研究に託したい。まずは本研究で不十分になってしまった D.discoideum の速度分布の観察を進め、理論を検証を完全にすることである。もう一つは D.discoideum で平均場の振る舞いを操作する技術の開発である。前者は当然クリアすべき課題であり、時間をかけて実験すれば必ず解決できるだろう。それよりも後者を期待したい。本研究で立てた理論が正しいならば、細胞は平均場が決まれば決まった動きしかしないはずである。となれば、問題は狙った平均場を発生させることに移行する。こうした技術が開発されれば、生きた細69

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

胞集団の操作が可能となり、アメーバが新機能性物質としての立ち位置を獲得するだろう。アメーバがマイクロスケール操作のツールとなり、新素材の開発や外科治療に活躍する日が来ることを期待する。

70

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

第 6 章 Appendix
以下では本文中で省略した議論や計算を説明する。ここにはどんな知識が念頭にあって、実際にどんな思考をして、どんな計算をしたのか、研究の泥臭い部分をあからさまに書く。以下の説明には非常に面倒な計算が含まれていたり、細かい知識の収集が含まれていたりするかもしれない。しかし、この Appendix により
本文での論理飛躍を全て埋めようと意図しているからである。ここでは計算上の技術の応用を考える方、この研究の追計算をする方に向けて丁寧に記す。

6.1

Appendix A:化学反応振動子モデルの Hopf 分岐の証明

Hopf の分岐定理と制御工学では有名な多項式の解の実部の符号を判定する Routh
Hurwitz の定理を組み合わせる Liu 氏の方法を用いると、Hopf 分岐の存在を示すことができる。本論文のストーリーとしては、ここで Hopf 分岐を証明することが、
Chemotactic 振動子に Hopf 分岐を要請すべき根拠となる。
Laub Loomis モデルに対しては R.Ghaemi 氏らの論文でホップ分岐点からの距離が議論されており、Goldbeter モデルに対しては P.E.Rapp 氏らの論文で Hurwitzの多項式が計算されている。ただ、後者については不動点の位置が数値計算によって計算されているため、Loomis モデルのような議論をするには、不動点の計算か
ら、線形係数の決定、多項式の計算、そして行列式の計算までを自動的に計算するプログラムを書く必要があるだろう。本研究では数理モデルを作る時に Hopf 分
岐近傍であることが仮定されているため、有名な代表的なモデルでその計算を行うことが重要であるが、既に過去に行われていて非常に煩雑な計算を必要とするので、その論文を紹介するにとどめる。ここでは彼らが共通して用いた Liu により開発された Hopf の分岐定理に Routh Hurwitz の定理を援用した Hopf 分岐のク
ライテリオンを紹介する。この方法はタフな計算を必要とするものの、汎用性が高く Hopf 分岐を証明するためのツールとして有用である。

71

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

6.1.1

W.Liu の Hopf 分岐クライテリオン

Hopf の分岐定理に、制御工学で有名な Routh Hurwitz の定理を組み合わせることで W.Liu 氏は力学系に現れている周期軌道が超臨界ホップ分岐により現れたものであることを証明する定理を導いた (Liu 1992)。この章ではそのアイディアを説明する。まずは定理が適用できる条件設定を説明する。
条件（W.Liu (1992)[15], R.Ghaemi et.al (2009)[23], X.Yang (2002)[33]）：p 次元パラメタ付き n 次元可微分力学系 f (x, k) を考える。
ẋ = f (x, k) with f ∈ C ∞ (Rn × Rp → Rn )(6.1)

この力学系は不動点 f (xf ix (k)) = 0 を持ち、その点の周りの線形係数を J(k) とおく。

J(k) :=∂f(xf ix (k), k) : Rp → M at(Rn × Rp )
∂x

(6.2)

係数 A(k) の振る舞いを調べることで分岐の判定をすることができるのだが、このままでは面倒である。パラメータ空間中にパラメータを動かす開曲線 Γ ⊂ Rp を設定して、パラメータ空間を実質的に１次元と見なして議論することにする。この仮定を置いても自由に Γ を選ぶことができるので議論の一般性は破れない。つまり、この章で考えるのは結局のところ 1 次元パラメタ付き n 次元可微分力学系

f (x, μ) である。

ẋ = f (x, μ) with f ∈ C ∞ (Rn × Γ → Rn )

(6.3)

∂f
(xf ix (μ), μ) : Γ → M at(Rn × Γ)
∂x

(6.4)

J(μ) :=

さて、まずは Hopf の分岐定理を述べる。
定理（Hopf の分岐定理）:
以下に述べる条件 (SH1)、(SH2) が成立する時、またその時に限り力学系 f (x, μ)は点 μ0 にて Hopf 分岐する。

(SH1)J(μ0 ) は 1 組の純虚数固有値ペア λ(μ0 )、λ(μ0 ) を持ち、その他の固有値の実部は負である。
(SH2)(SH1) の純虚数固有値に対応する J(μ0 ) の固有値 λ(μ) は μ0 周りで滑らかであり、以下の条件を満たすとする。
dλ(μ0 ) = 0
dt
72

(6.5)

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

これはいわば超臨界 Hopf 分岐の定義でもあり、直接的に固有値を計算しない限
り Hopf 分岐かどうか判定することはできない。現在では高次元の行列であっても数値的に固有値を評価する技術が非常に発達しているため機会は少ないと思われるが、今回のような何か特別な理由があって解析的な議論が必要な場合には、高次の多項式の解分布を取り扱うための工夫が必要である。先に述べたように Liu はこの問題を Rauth Hurwitz の定理を援用することで解決した。
定理（Routh Hurwitz の定理）

n 次方程式pn λn + pn−1 λn−2 + · · · + p1 λ + p0 = 0 (pn = 1)

(6.6)

の全ての根が負の実部を持つための必要十分条件は、以下の条件（RH）を満たすことである。(RH) p0 > 0 であり、かつ以下で定義される Hurwitz の行列式

D1 , D2 , · · · , Dn の値が全て正である。

Dn =

p1

p0

0

0

p3

p2

p1

p0

···
...
...

p5

p4

p3

p2

p7
..
.

p6
..
.

p5
..
.

p4
..
.

...
...

0

0

0
..
.
..
.
..
.

0
..
.
..
.
..
.

(6.7)

p2n−3 p2n−4 p2n−5 p2n−6 · · · pn−3 0
p2n−1 p2n−2 p2n−3 p2n−8 · · · pn−1 pn
ただし、pj = 0 (j < 0 or j > n) とする。

定理（W.Liu 1992）
特性方程式の係数を次のように置き、対して Hurwitz の行列式 (D1 (μ), D2 (μ), · · · , Dn (μ))
を考える。

det(λIdn − J(μ)) = pn λn + pn−1 λn−1 + · · · + p1 λ + p0

(6.8)

この時 Hopf の分岐条件 (SH1)、(SH2) は特性方程式に対して次の２つの条件：

p0 (μ0 ) > 0, D1 (μ0 ) > 0, · · · , Dn−1 (μ0 ) = 0

(6.9)

dDn−1
(μ0 ) = 0dμ(6.10)

73

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

が成り立つための必要かつ十分な条件である。
従って、我々は Dm の符号変化を検査することで、Hopf 分岐を起こす面（以下、

Hopf 分岐面と呼ぶ。）がどこにあるか探すことができる。用いているパラメータからどれくらい離れているのかを評価することができる。以下の章でどのようにこの定理を使うのか簡単な２つの具体例を通して説明する。

74

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

6.1.2

Liu のクライテリオンの適用例 1（Brus-selator）

ここでは Liu のクライテリオンの適用方法を説明するために簡単な振動モデル
を取り上げる。１つめは Brus-selator。代謝系の化学反応の簡単なモデルである。
式の形が非常に単純なため Hopf 分岐を考察する際の例題としては最適である。

Brus-selator は次の方程式で与えられる。
dx= a − (1 + b)x + x2 ydtdy
= bx − x2 ydt

(6.11)

この力学系の不動点は (a, b/a)。線形化した力学系は次の通り。

dx= (b − 1)x + a2 ydtdy
= −bx − a2 ydt

(6.12)

固有値多項式は次の通り。具体的に解くのも易しいが、今回は解かない。

λ2 + (a2 − b + 1)λ + a2 = 0(6.13)

Hurwitz 行列式を計算すると次のようになる。
p0 = a2 > 0(6.14)

D1 = a2 − b + 1(6.15)(6.16)

従って、a2 −b+1 = 0 にて Hopf 分岐が起こり安定不動点が不安定化し、a2 −b+1 > 0
にてリミットサイクルが発生する。ただ、２次元では Hurwitz 行列式が１成分の行列式となってしまうため、例題としては簡単すぎた。

75

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

6.1.3

Liu のクライテリオンの適用例 2（Lorentz 方程式）

２つ目の例は Lorentz 方程式。流体のモデルが由来になっている３次元の方程式である。Lorentz 方程式は次で与えられる。

dx= −10x + 10ydtdy= ρx − y − xzdt3dz
= − z + xy
dt
8

(6.17)

方程式系 6.17 の不動点は次の３つ。
  




8
8
8
8
(ρ − 1),(ρ − 1), ρ − 1 , −(ρ − 1), −(ρ − 1), ρ − 1(0, 0, 0),3333
固有値多項式は次の通り。

λ3 +

160
41 2 8λ + (10 + ρ)λ +(ρ − 1)3
3
3(6.18)

Hurwitz 行列式は次の通り。
160(ρ − 1)
3
8(10 + ρ)=3
8
(10 + ρ)3=
1

(6.19)

p0 =
D1
D2

160
(ρ − 1)3
41
3(6.20)


3760
808
4
152
ρ+
=−
ρ−
(6.21)=−
99919

従って、ρ = 470/19 で超臨界 Hopf 分岐が起こることが分かる。このような手順で系統的に Hopf 分岐を計算することができる。

76

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

6.2

Appendix B:Chemotactic Oscillator の位相縮約

Swarm Oscillator は general な走化性粒子の運動を決める方程式を Chemotactic
Oscillator を位相縮約することで導出される。非常に長い計算が必要なため、本文ではその導出を省略した。ここでは位相縮約の説明も兼ねて、その導出過程を説明する。ここでは他分野の方にも分かりやすいように、中心多様体縮約と呼ばれる方針をとる。その他に簡便法と呼ばれる位相縮約の方法もあるのだが、実際には中心多様体縮約よりも分かりにくいように感じる。ここでは他分野の方にも理解しやすいように、より系統的で分かりやすい中心多様体縮約の方法を紹介する。
ここの計算過程を追うことで、Swarm Oscillators に現れる定数が物理量とどのように関係付けられるかを把握することができる。
dXμ(t) = f (Xμ ; λ) + kg (S (rμ ))
dt
drμ2
drμ
dS
− σ (Xμ )(t, rμ )m 2 (t) = −γ
dt
dt
drμ

∂S
(t, r) = −S(t, r) + dΔS(r, t) +h (Xμ ) δ (r − rμ )τ∂tμ
ここで、Xμ ∈ Rn 、rμ ∈ M D はそれぞれインデックス μ ∈ {1, 2, · · · , N } で指定される細胞の内部状態、空間座標を表している。（V n は n 次元実ベクトル空間であり内部状態を指定する点の空間である。M D は D 次元微分多様体で細胞が動くことができる空間である。）
実関数 S(t, r) ∈ R は時刻 t ∈ R、空間座標 r ∈ M D における化学物質の濃度（D.discoideum の場合は cAMP の濃度）を表す。ベクトル値関数 g(S) ∈ Rn (g(0) =0) は化学物質環境に対する、内部状態の応答を表す関数。f (X; λ) は内部状態の関数で、内部状態の自己応答を表す。
実パラメタ λ はこの関数の形を変える実パラメタで、μ = 0 で Hopf 分岐を起こすことを仮定する。σ(X) ∈ RD×D は内部状態の関数で、化学勾配に対する空間移動の応答を表す。h(X) ∈ R は内部状態の関数で、細胞からの化学物質放出を表す。
M D は何を選んでも構わないが、簡単のため 2 次元トーラス T 2 = S × S を選ぶ。
また実パラメタ k 、σ 、m、γ 、τ 、d はそれぞれ、外部化学物質に対する内部状態への応答強度、質量、運動粘性、化学物質拡散の時定数、化学物質の拡散定数を表している。
その他、以下の細かい設定を仮定する。
f (0; λ) = 0, g(0) = 0, h(0) = 0

77

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

中心多様体縮約を使って、標準形に変形する。Super Critical Hopf 分岐を出すた
めに、線形項の形を制限する。まずは変数変換により線形項を実標準化する。変
換行列を Q として、


L = Q−1 Df (0; λ)Q =a(λ) −ω(λ)ω(λ) a(λ)

⊕ L1 = L0 ⊕ L1

ここで、a(0) = 0、a (0) > 0、ω(0) = 0 とする。L の固有値は a(λ) ± iω(λ) であり、固有ベクトルを U ∈ C2 、U ∈ C2 とする。これは U · U = 1、U · U = 1、

U · U = 0 を満たす。また L1 の固有値は全て実部が負であるとする。以下では初めから実標準化の既定に変換されていたとして議論する。
初めの２成分 (x, y) は回転的に振る舞うので、複素変数に書き直すと式が見やすくなる。
zμ = xμ + iyμ ∈ R ⊕ iR = C, その他の変数はζμ ∈ Rn−2
kg(S) の項は後で考えることとし（i.e k = 0）、第１式を書き換える。次のようになったとする。
dzμ= (a(λ) + iω(λ)) zμ + F1 (zμ , zμ , ζμ ; λ)dtdzμ
= (a(λ) − iω(λ)) zμ + F1 (zμ , zμ , ζμ ; λ)dtdζμ
= L1 ζμ + F2 (zμ , z¯μ , ζμ , λ)dt
F1 、F2 は２次以上の項。
F1 (z, z̄, ζ; λ) = f (re(z), im(z), ζ; λ)z 成分 −∂f(0; λ)z 成分 z∂z

F2 (z, z̄, ζ; λ) = f (re(z), im(z), ζ; λ)ζ 成分 −∂f(0; λ)ζ 成分 ζ∂ζ

中心多様体 W c は次のように決定される（ここの detail は関係ない。）

W c = (zμ , α(zμ ; λ))μ=1···N ∈C × Rn−2 | |zμ | < δ, μ = 1 · · · N
μ=1···N

with δ > 0 : small constant∂αμ(zμ ; λ) {(a(λ) + iω)zμ + F1 (zμ , zμ , αμ ; λ)} = L1 αμ (zμ ; λ) + F2 (zμ , z¯μ , αμ , λ)∂zμ中心多様体上に制限した方程式は次の通り。

dzμ= (a(λ) + iω(λ)) zμ + F1 (zμ , zμ , αμ (zμ ; λ); λ)dt78《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

標準化の方法に従って非線形項を単純化する。まず 2 次の項に対する近恒等変換を行う。

zμ → zμ + βμ (zμ , zμ ; λ) (βμ は zμ , zμ について２次式。)
方程式は次のように変換される。
dzμ∂βμ dzμ ∂βμ dzμ(2)++= (a + iω) (zμ + βμ ) + F1 + O(3)dt∂zμ dt∂zμ dtdzμ∂βμ dzμ ∂βμ(2)++(a − iω)zμ = (a + iω) (zμ + βμ ) + F1 + O(3)dt∂zμ dt
∂zμ
−1 !
"

∂βμ
dzμ
∂βμ
(2)
(a + iω)zμ −+ O(3)= 1+(a − iω)zμ + (a + iω)βμ + F1
dt
∂zμ
∂zμ

= (a + iω)zμ + (a + iω)βμ −

∂βμ
∂βμ
(2)
(a + iω)zμ −(a − iω)zμ + F1 + O(3)∂zμ∂zμこれより次の写像を考える。

φ2 : βμ → (a + iω)βμ −∂βμ∂βμ(a + iω)zμ −(a − iω)zμ∂zμ∂zμ2 次の斉次多項式空間の間の写像であり、基底を zμ2 , zμ zμ , zμ 2  と選んで表現行列を書くと次のようになる。
⎞
⎛
⎞⎛
⎞
⎛
−a − iω
0
0
zμ2
zμ2
⎟
⎜
⎟⎜
⎟
⎜
φ2 : ⎝ zμ z¯μ ⎠ → ⎝0−a + iω 0 ⎠ ⎝ zμ z¯μ ⎠
zμ 2
zμ 2
0
0
2iω
kerφ2 = 0 であり、うまく βμ を選べば全ての２次の項を消すことができる。それ(2)を βμ としよう。
(2)βμ(2) = φ−1
2 (F1 ) = −1  (2)  2
1  (2) 
1  (2) 
F1
F1
F1
zμ −
zμ z μ +zμ 2
a + iω
a − iω
2iω
zz
z z̄
z̄ z̄

次に３次を消すため、zμ を次のように変換する。

zμ → zμ + βμ(2) (zμ , zμ ; λ) + βμ(3) (zμ , zμ ; λ) (βμ は zμ , zμ について３次式)
方程式の変換を書くのだが、少し２次のときと結果が変わる。
(2)(3) −1
∂βμ
∂βμ
dzμ
=
1+
+
{(a + iω)zμ
dt
∂zμ
∂zμ
(2)

−

(3)

∂βμ
∂βμ
(a − iω)zμ −(a − iω)zμ + (a + iω)βμ(3) + F1
∂zμ
∂zμ
79

+ O(4)

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

(3)

(3)

∂βμ
∂βμ
= (a + iω)zμ + (a + iω)βμ −(a + iω)zμ −(a − iω)zμ
∂zμ
∂zμ
⎫
⎧ 

(2) 2
⎬
⎨
∂βμ
(3)
+ O(4)+ −(a + iω)zμ + F1⎭⎩∂zμ
この場合も考える写像は全く同じ形で（定義域と値域が異なる）、
(3)

∂βμ
∂βμ
φ3 : βμ → (a + iω)βμ −(a + iω)zμ −(a − iω)zμ∂zμ∂zμ２次のときと同様にして、基底 zμ3 , zμ2 zμ , zμ zμ 2 , zμ 3  に対する表現行列を計算すると、次の通り。

⎞⎛⎞⎞⎛−2(a + iω) 0
0
0
zμ3
zμ3
⎜
⎜ z2 z ⎟⎟ ⎜ z2 z ⎟
0
−2a
0
0
⎜
⎜ μ μ ⎟⎟⎜ μ μ ⎟φ3 : ⎜⎟→⎜⎟⎟⎜
2⎝⎝ zμ zμ ⎠⎠ ⎝ zμ zμ 2 ⎠00 2(a − iω)
0
zμ 3
zμ 3
0
0
0
−2a + 3iω⎛

この行列を見ると、λ = 0 にて ker φ3 = zμ |zμ |2  となることが分かる。従って、分(3)
岐点付近で zμ |zμ |2 の項が残る。一意ではないが、βμ として
⎛
⎞


(2) 2
∂β

μ
⎝−F1(3) +βμ(3) ∈ φ−1(a + iω)zμ ≡ −F (3) ⎠3∂zμ（注：後で用いる仮定だが、Re F (3) < 0 とする。）となるものを一つ選ぶ。例えば、






1
1
1
(3)
(3)
3
(3)
(3)
2
F1
F
F
βμ = −zμ +zμ zμ −zμ 3
2(a + iω)2(a − iω)2a − 3iωzzz
z z̄ z̄z̄ z̄ z̄を選べばよい。説明するまでもないが、係数は次のようにして決まっている。

 2


∂F1 ∂ 2 αμ∂F1
∂ 2 αμ
∂F1
∂ F1
(2)
2 1
F1 =
(0; λ) +(0) 2 (0) zμ +(0; λ) +(0; λ)(0) |zμ |2
∂zμ2
∂ζμ
∂zμ
2 ∂zμ ∂zμ∂ζμ∂zμ ∂zμ 2
∂ F1 (0; λ) ∂F1 ∂ 2 αμ++(0) zμ 2∂ζμ ∂zμ
∂zμ 2


 3
 3
1
∂F1 ∂ 3 αμ
∂ F1
∂ F1 ∂F1 ∂ 3 αμ
(3)
3
zμ +
z 2 zμ
F1 =
+
+
∂zμ3
∂ζμ ∂zμ3
3 ∂zμ2 ∂zμ∂ζμ ∂zμ2 ∂zμ μ


 3
 3
1
∂F1 ∂ 3 αμ
∂ F1
∂ F1 ∂F1 ∂ 3 αμ
2z μ zμ ++++zμ 3
3 ∂zμ ∂zμ 2∂ζμ ∂zμ ∂zμ 2∂ζμ ∂zμ 3
∂zμ 3

80

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

４次の項について同様の議論を行う。結局、変換行列のみが問題なので変換行列を書き下す。
⎞
⎛
zμ4
⎟
⎜ 3
⎜ zμ z μ ⎟
⎟
⎜
2 2 ⎟
φ4 : ⎜⎜ zμ z μ ⎟ →⎟⎜⎝ zμ z μ 3 ⎠
zμ 4
⎞
⎛
⎞⎛
−3a − 3iω
0
0
0
0
zμ4
⎟
⎜
⎟⎜ 3
⎜
⎟ ⎜ zμ zμ ⎟0−3a − iω
0
0
0
⎟
⎜
⎟⎜
⎜
⎟ ⎜ z 2 zμ 2 ⎟
0
0
−3a
+
iω
0
0
⎟
⎜
⎟⎜ μ
⎜
⎟⎜
3 ⎟
z
z
0
0
0
−3a
+
3iω
0
⎝
⎠⎝ μ μ ⎠
zμ 4
0
0
0
0
−3a + 5iω
よって４次の項は全て消すことができる。（実は全ての偶数項を消すことができる。）結局、方程式は次の形に標準化される。


dzμ

= (a(λ) + iω(λ))zμ + F (3)z|z|2 + O(5)dtzz z̄
変数変換で定常的な回転成分を分離する。すなわち zμ = Aμ eiωt と変数変換して、


dAμ
(3)
= a(λ)Aμ + F
Aμ |Aμ |2 + O(5)dtzz z̄
以下で kg(S) の項を取り込む計算をする。この際も X 、S は小さいと仮定して、線形性分のみ取り入れる。まず、S の方程式を解いて代入する。ここに田中氏による、計算を進めるためのやや大胆な仮定が入る。

S(t, r) =



∞

μ

−∞


≈

t

ds
−∞

dt

)
*



dD k  h(Xμ ) exp i s(t − t ) + k · (rμ (t ) − r)
1 + iτ s + d|k  |2D  exp ik · (rμ (t) − r)d k+ c.c.
(2π)D 1 + iτ ω + d|k  |2

(2π)D



∂h
(0) eiωtAμ∂zμ

空間依存因子をまとめて次のように書く。

M (r) =μAμ G(rμ − r) =






exp
ik
·
(r
(t)
−
r)
μ
d k
D
(2π)
1 + iτ ω + d|k  |2
D 

Aμ

μ

これで S の計算は終わり。この項を方程式に導入するのだが、その際にも M の複素共役の寄与を無視するという仮定が入る。この項は速い振動を与え、この近似は平均化法による繰り込みと見なすことができる。


∂h
dAμ

= a(λ)Aμ + F (3)(0)
Aμ |Aμ |2 +ηM (r)+O(5) with η = k (gx (0) − igy (0))
dt
∂z
zz z̄

81

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

パラメータの所在を見やすくするため、変数変換を行う。
 (3) 
Im
F
1
η
1
zz z̄
|Re (F (3) )zzz̄ |Aμ (t), t = a(λ)t, c =Aμ (t ) =
, χ=

(3)
a(λ)
a(λ) Re (F )zzz̄
a(λ)
すると原論文での表記に一致する。以下この変数変換を施したものとし、プライ
ム記号 ( ) を書かない。

dAμ= Aμ − (1 + ic) Aμ |Aμ |2 + χM (r) + O(5)dt
係数の調整をしたいという目的もあり、ここで座標部分の計算を行う。座標の方程式の環境物質依存部分を展開して線形化しておく。

md2 r μdrμ∂σ dS− Aμ e−iωt(t) = −γ(t, rμ )2dt
dt∂ z̄ drμ

over damp を仮定して慣性項を落とす。
1 ∂σ
drμ
dS
=−
Aμ e−iωt(t, rμ ) = −dtγ ∂ z̄
drμ

dM
1 ∂σ ∂h(0) Aμ(rμ ) + c.c.
γ ∂ z̄ ∂zdrμ係数を関数 M の中に入れこむ。（G も少し書き直す。）

1 ∂σ ∂h(0) dM (r)M (r) = bdM (r) =dγ ∂ z̄ ∂z
応じて、




χ =

M (r) =μAμ G (rμ −r) =

1 ∂σ ∂h(0)dγ ∂ z̄ ∂z

Aμ dG(rμ −r) =

μ

−1
χ


μ



Aμ



dD k  exp ik · (rμ (t) − r)(2π)D
ρ2 + d|k  |2

ここで、ρ = (1 + iτ ω)/d とおいた。以下、G は G と書く。さて、このように調整すると、次の式を得る。ここでもプライム ( ) は省く。

drμ
dAμ
dM
= Aμ − (1 + ic) Aμ |Aμ |2 + χM (rμ ) + O(5),= −Aμ(rμ ) + c.c.
dt
dt
drμ
リミットサイクルに捕まっているあるいは、近くに留まっていると仮定して位相縮約を行う。χ の絶対値が十分小さいとすれば成り立つ仮定である。極座標表示
(r)

(r)

Aμ = Aμ e−icφμ で、動径部分 Aμ ∈ R の方程式は安定不動点を持つ。摂動による解のずれを１次のオーダーまで取り込んで位相方程式に代入する。この操作は中82



《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

心多様体縮約と呼ばれる。具体的には次のように計算する。まずは極座標表示で方程式を書き直す。
(r) −icφμdAμdφμ− icA(r)e−icφμ = A(r)+ χM (rμ )1 − (1 + ic)A(r)2
e
μ
μ
μ
dt
dt
実部と虚部を分離して、（O(χ) のオーダーまで計算する。）
(r)



dAμ
= A(r)
1 − A(r)2+ χ Re Pμ
μ
μ
dt
dφμ
= −cA(r)3+ χ Im Pμμdtwith Pμ =e−ic(φν −φμ ) G(rν − rμ )−cA(r)
μ

ν
(r)

(r)

非摂動では動径部分は Aμ → 1 と収束する。摂動解を Aμ → 1 + χaμ として代入すると。
daμ= −2aμ + Re Pμdt
従って、摂動解は aμ → (1/2) Re Pμ となる 1 。これを位相部分に代入する。


1
dφμ
= 1 + χ Re Pμ − Im Pμ
dt
c
すなわち、

dφμ
=1+χ
dt

!

i
1
+
2 2c



"
Pμ + c.c.

結果、方程式は次のようになる。

dPμ
drμ
dφμ
= 1 + (κPμ + c.c.) ,=−+ c.c.
dt
dt
drμ



i
1
−ic(φν −φμ )+eG(rν − rμ ) , κ = χwith Pμ =22c
ν
G を少し計算すると Swarm Oscillator が得られる。

dD k  eik ·rbb
G(r) = b≈
exp (− Re ρー i Im ρ) |r| + arg(2π)D ρ2 + |k  |22ρ
ρ
すると、方程式は次のようになる。


bκ 
dφμ
bκ  − Re ρ|rν −rμ |=1+ecos c (φν − φμ ) + Im ρ|rν − rμ | − argdtρρν(=μ)

1

リミットサイクルの仮定から、軌道は収束する。

83

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）


d b  − Re ρ|rν −rμ |
b
drμ
=−
e
cos c (φν − φμ ) + Im ρ|rν − rμ | − arg
dt
drμ ρ
ρ
ν(=μ) r ν − rμbb
− Re ρ|rν −rμ |=−esin c (φν − φμ ) + Im ρ|rν − rμ | − argIm ρρ|rν − rμ |ρν(=μ)+ r ν − rμ
b
b
Re ρ
e− Re ρ|rν −rμ | cos c (φν − φμ ) + Im ρ|rν − rμ | − argρ|rν − rμ |ρν(=μ)= −|b|

 r ν − rμe− Re ρ|rν −rμ | cos (c (φν − φμ ) + Im ρ|rν − rμ | − arg b)|rν − rμ |ν(=μ)

ここで最後の変数変換を行う。

θμ = c (φμ − t) , t =cbκt, r = reρ r,ρ

パラメタも整理する。

α=

bκ π
π
Im ρ
ρ
, c1 = arg− , c2 = arg b + , c3 =
Re ρ
Re ρ
ρ
2
2
cκ

これで求めていた Swarm Oscillator が得られる。

dθμ=e−|rν −rμ | sin (θν − θμ + α|rν − rμ | − c1 )dt
ν(=μ)

drμrν − rμ= c3
sin (θν − θμ + α|rν − rμ | − c2 )e−|rν −rμ |dt|rν − rμ |ν(=μ)

84

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

6.3

Appendix C:Particle Image Velocimetry の概
要PIV(Particle Image Velocimetry) とは中程度の濃度の粒子流の動画データ（連続的な静止画データ）から、その画像の相関関数を使って速度ベクトル場を構成する方法である。もともとは、流体の速度分布を解析するために開発されたアルゴリズムで、流体中に小さて光学系でとらえやすいマーカー粒子を流して、粒子ベースで追いかけて速度を見積もることを目的とされていた [9]。
今回の実験データに対して用いたのは、連続する２次元の平面データに対しての適用である。その原理をここで説明する。まずは空間を細胞スケールの数倍程度のメッシュで切って空間をセルに分割する。時刻も適切な間隔で離散化されているものとする。
（画像の順番を指定する変数なので、時間指定は普通この状況になっている。）画像データからのパターン強度 I(t, x, y) を抽出し、セル (i,j)（代表
点 (xi , yj )）の時間相関関数 R(t, i, j; Δx, Δy) を次のように構成する。
dx dy  I(t, x , y  )I(t + Δt, x + Δx, y  + Δy) (6.22)

R(t, i, j; Δx, Δy) =cell(i,j)

セル (i, j) の領域を cell(x, y) と書いた。Δt は時間の離散間隔である。なお、この積分はパッケージの中では高速フーリエ変換を利用して効率的に行われる。相関関数 R のピークを調べることで (Δx, Δy) が決まる。これを全てのセルに対してそれぞれ計算することで (Δx(i, j), Δy(i, j)) を得る。これからセル (i, j) に含まれる粒子の速度ベクトル (vx (i, j), vy (i, j)) を次のように評価することができる。

Δx(i, j) Δy(i, j),(6.23)(vx (i, j), vy (i, j)) =ΔtΔt
ただ、粒子密度に応じて手作業でのメッシュの取り方のチューニングが必要である。粒子流の密度に対して荒すぎる空間メッシュを使うと速度が小さく見積もられてしまう。反対にこまか過ぎるメッシュを使うと細胞の変形等の細かな変化を抽出してしまい正しいベクトル場を得ることができない。またメッシュのチュー
ニングを最適化しても、画像データの荒れなどから現れる無意味な大きいベクトル場が計算結果に現れることがある。またこれも適切にチューニングされたフィルターをかけて取り除くしか無い。

85

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

謝辞
本研究を遂行するにあたり、非常に多くの方にお世話になりました。R 研の小西哲郎先生、宮崎州正先生、野崎一洋先生には私の勝手気ままな研究を暖かく見守って頂き、日々アドバイスを頂きました。生命科学科の塚田祐基先生には画像解析の方法をはじめ実験の相談にのって下さいました。情報科学科の青木摂之先生、手塚修文先生にはノウハウを提供頂き、さらには実験の器具を貸して下さいました。情報科学科の杉山雄規先生には実験道具を揃えるのに協力いただいたり、
また研究会に招待いただいたり、時には人生相談にまで載って頂き、大変お世話になりました。東大総合文化の澤井哲先生には最先端の D.discoideum のお話を伺い、さらには実験の相談にものって頂きました。東京医科歯科大学田中研様からは、実験で使う器具を一部提供して頂きました。東大大学院医学系循環器系内科
様からは、寒天培地の作り方、バクテリアの増やし方、器具の滅菌の方法など実験手法を基本的なところから教えて頂きました。
その他にも多くの先生方に暖かく見守って頂きました。先生方の多大な応援、支援に感謝しています。
また、情報科学の石渡龍輔さんとは実験を一緒に行い、実験の初段階から共同研究をさせて頂きました。愛工大の巌佐正智さんには理論の組み立てや論文の書き方などについて相談にのって頂きました。実験の経験が豊富な J 研の久保宏樹
さん、岡田薫さんにも研究を進める過程で相談にのって頂きました。G 研の三野広幸先生には学部生のときから長らくお世話になり、すれ違うたびに「おぉ君か、
ヒマそやな」とウィットに富んだ叱咤激励を頂きました。Tb 研のお隣さんの木村明洋先生のフルートの音色は素敵でした。QG 研の久木田真吾さんや私と同じ R
研の伊賀裕展さんは良いライバルとなって切磋琢磨し合い、研究の推進力となりました。また、R 研・Tb 研・S 研の仲間たち、情報科学の仲間たち、時々R 研に遊びにいらっしゃる工学研究科の伊藤夢奈さん、OB の鬼塚健太さん、石田つむぎ
さん、数多くの仲間のおかげで日々楽しく研究生活を送ることができました。物
理事務室の皆様、理学事務室の皆様、職員食堂のおばちゃん、Cafe Phonon のマ
スターにはいつもにこやかに声をかけて頂き、晴れ渡る気分で充実した学生生活を送ることが出来ました。私は人に恵まれ、多くの人々から学んだことは私にとって大きな財産となりました。
さらに、理学図書室の皆様とは３年間にわたり仕事を一緒にさせて頂き、図書86

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）のシステマティックな管理や利用者を増やす企画作りなど多くのことを学ぶことができました。また、中部大学第一高校の先生方と１年間教員の仕事を一緒にさせて頂きました。教育の難しさプロの教員として働くことの厳しさを学びました。
縁あって私の生徒となった同高校の学生の皆様からも多くのことを学びました。E
研のフォン先生には 2 年間 TA の仕事をさせて頂き、留学生たちの物理を学ぶ気迫を感じることが出来ました。障害学生支援室の瀬戸今日子先生、支援室スタッフ、
学生サポーターの仲間たちと共にボランティア活動に２年半にわたり取り組めた
ことはとても貴重な経験となり、忙しくも充実した学生生活を送ることができました。また、上海での短期留学では同済大学の皆様、名大上海事務所のリュウ・ラ
イさんとヒョウ・カニさん、GAOTIME の皆様、豊田紡織上海の皆様に大変お世話になり、大きく自分を成長させるきっかけになりました。
ここで挙げた方々を始め数多くの人々の支えがあったおかげで、私は研究を最後までやり抜くことができたのだと思います。この場を借りて皆様に感謝申し上げます。

最後に D.discoideum の株を分けて下さった NBRP 事業に関わる皆様に感謝申し上げます。

We would like to acknowledge National BioResource Project (NBRP) by MEXT.
今後とも NBRP 細胞性粘菌にご理解とご協力をよろしくお願いします。

87

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）

関連図書[1] A.J.Durston. Pacemaker activity during aggregation in dictyostelium discoideum. Developmental Biology, Vol. 37, pp. 225–235, 1973.
[2] E Alcantara and M Monk. Signal propagation during aggregation in the slimemould dictyostelium discoideum. Journal of general microbiology, Vol. 85, pp.
321–334, 1974.
[3] John Tyler Bonner and L. J. Savage. Evidence for the formation of cellaggregates by chemotaxis in the development of the slime mold dictyostelium
discoideum. Journal of Experimental Zoology, Vol. 106, No. 1, pp. 1–26, 1947.
[4] D. Cohen, J. Neu, and R. Rosales. Rotating spiral wave solutions of reactiondiﬀusion equations. SIAM Journal on Applied Mathematics, Vol. 35, No. 3,
pp. 536–547, 1978.
[5] Tanaka Dan. General chemotactic model of oscillators. Phys. Rev. Lett.,Vol. 99, p. 134103, Sep 2007.
[6] D.Dormann and C.J.Weijer. Review: Imaging of cell migration. The EMBO
Journal, Vol. 25, pp. 3480–3493, 2006.
[7] D.J.T.Sumpter. The principles of collective animal behaviour. Phil. Trans.
R. Soc. B 29, Vol. 361, , 2006.
[8] D.J.T.Sumpter. Collective Animal Behavior. Princeton Univ Pr, 2010.
[9] ヴィラート CE, コンペンハンス J. PIV の基礎と応用: 粒子画像流速測定法.
シュプリンガーフェアラーク東京, 2000.
[10] F. D. Hollander. Large Deviations (Fields Institute Monographs, 14). Amer
Mathematical Society, 2000.
[11] M. Iwasa, K. Iida, and D. Tanaka. Various collective behavior in swarm
oscillator model. Physics Letters A, Vol. 376, No. 3031, pp. 2117 – 2121,
2012.
88

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）[12] P.Holmes J.Guckenheimer. Nonlinear Oscillations, Dynamical Systems, andBifurcations of Vector Fields. Springer, 2002.
[13] J.L.Martiel and A.Goldbeter. A model based on receptor desensitization for
cyclic amp signaling in dyctyostelium cells. BIOPHYS.J., Vol. 52, pp. 807–
828, Nov 1987.
[14] K A Kellerman and J G McNally. Mound-cell movement and morphogenesis
in dictyostelium. Developmental Biology, Vol. 208, No. 2, pp. 416–429, April
1999.
[15] Wei-Min Liu. Criterion of hopf bifurcation without using eigenvalues. Journalof Mathematical Analysis and Apprications, Vol. 182, pp. 250–256, 1994.
[16] Erik A. Martens, Carlo R. Laing, and Steven H. Strogatz. Solvable model ofspiral wave chimeras. Phys. Rev. Lett., Vol. 104, p. 044101, Jan 2010.
[17] Iwasa Masatomo and Tanaka Dan. Dimensionality of clusters in a swarm
oscillator model. Phys. Rev. E, Vol. 81, p. 066214, Jun 2010.
[18] M.Iwasa. private discussion.
[19] M.T.Laub and W.Loomis. A molecular network that produces spontaneous
oscillations in excitable cells of dyctyostelium. Mol.Biol.Cell, Vol. 52, pp.
3521–3532, Dec 1998.
[20] Seido Nagano and Shunsuke Sakurai. Cell-to-cell coordination for the spontaneous camp oscillation in dictyostelium. Phys. Rev. E, Vol. 88, p. 062710,
Dec 2013.
[21] P.E.Rapp, R.A.Latta, and A.I.Mess. Parameter-dependent transitions andthe optimal control of dynamical diseases. Buletin of Mathematical Biology,
Vol. 50, No. 3, pp. 227–253, 1988.
[22] R.E.Goldstein. Traveling-wave chemotaxis. Physical Review Letters, Vol. 77,
No. 4, 1995.
[23] R.Ghaemi, J.Sun, P.A.Iglesias, and D.D.Vecchio. A method for determining
the robustness of bio-molecular oscillator models. BMC Systems Biology,
Vol. 3, , 2009.
[24] R.H.Kessin. Evolutionary biology: Cooperation can be dangerous. Nature,
Vol. 408, pp. 917–919, 2000. 図を引用。
89

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）[25] M.Iwasa R.Ishiwata. (To be published。)D.Discoideum の動きからパラメー
タ α, c1 , c2 , c3 の値を見積もった。
[26] Satoshi Sawai, Xiao-Juan Guan, Adam Kuspa, and Edward C Cox. Highthroughput analysis of spatio-temporal dynamics in Dictyostelium. Genome
biology, Vol. 8, No. 7, p. R144, 2007.
[27] Nagano Seido. Diﬀusion-assisted aggregation and synchronization in dictyostelium discoideum. Phys. Rev. Lett., Vol. 80, pp. 4826–4829, May 1998.
[28] S.Saran, M.E.Meima, K.E.Weening E.Alvarez-Curto, D.E.Rozen, andP.Schaap. camp signaling in dictyostelium : Complexity of camp synthesis, degradation and detection. Journal of Muscle Research and Cell Motility,
Vol. 23, pp. 793–802, 2002.
[29] S.Sawai, P.A.Thomason, and E.C.Cox. An autoregulatory circuit for longrange self-organization in dictyostelium cell populations. Nature, Vol. 433,
pp. 323–326, 2005.
[30] S. H. Strogatz. Nonlinear Dynamics And Chaos: With Applications To
Physics, Biology, Chemistry, And Engineering (Studies in Nonlinearity).
Westview Press, 2001.
[31] S.Wiggins. Introduction to Applied Nonlinear Dynamical Systems and Chaos(Texts in Applied Mathematics). Springer, 2003(2nd).
[32] Dan Tanaka. Swarm oscillators : Normal form of motile particles with internal
dynamics. Progress of theoretical physics. Supplement, Vol. 178, , 2009.
[33] Xiaojing Yang. Generalized form of hurwitz-routh criterion and hopf bifurcation of higher order. Applied Mathematics Letters, Vol. 15, pp. 615–621,
2002.
[34] Y.Kuramoto. Chemical Oscillations, Waves, and Turbulence. Springer, 2007.
[35] 郡司ペギオ幸夫. 群れは意識をもつ 個の自由と集団の秩序. PHP 研究所, 2013.
[36] 柴田和正. 漸近級数と特異摂動法 - 微分方程式の体系的近似解法. 森北出版,
2009.
[37] 蔵本由紀. リズム現象の世界 (非線形・非平衡現象の数理). 東京大学出版会,
2005.

90

《修士論文》

物性研究・電子版 Vol. 3, No. 2 032602 （2014年5月号）[38] 蔵本由紀, 河村洋史. 同期現象の数理―位相記述によるアプローチ (非線形科学シリーズ). 培風館, 2010.
[39] 加納剛史 (東北大), 大須賀公一 (大阪大), 小林亮 (広島大, JST CREST), 青沼
仁志 (北海道大), 石川将人 (大阪大), 杉本靖博 (大阪大), 大脇大 (東北大), 石黒章夫 (東北大, JST CREST). 交友関係の自己組織化過程の数理モデル. 第 19回交通流のシミュレーションシンポジウム, 2013.

91