数理解析研究所講究録
第 828 巻 1993 年 91-102

91
走化性乱歩による集合モデル
九州大学理学部生物学科

江副日出夫、

巌佐

庸

はじめに昆虫、細胞性粘菌、 バクテリアなどがフェロモンを出して寄り集まる行動は、 現在まで主に反応拡散モデルによって研究されてきたが、各個体の乱歩をそのまま個体密度の拡散におきかえて扱ってよいかどうかについてはあまり問題にされてこなかった。 そこで、 本研究では、 とくに生物個体の数が比較的少ない (数十匹) 場合について、 生物個体の個体密度のかわりに各個体の位置の変化を追うモデルが偏微分方程式にもとつく反応拡散モデルとは異なるふるまいを示すことを明らかにする。 たとえば反応拡散モデルにおいて一様平衡解が安定であっても、 対応する乱歩モデルでは、 個体が少数のクラスタに集合し、 クラスタどうしがランダムな動きをしながらしだいに融合する。 クラスタの動態を解析的に記述することにより、 走化性集合行動モデルの新しいアプローチを試みた。 このとき、 動物が自らつくったフェロモンにひきずられて、 動きが遅くなる 「摩擦」現象が重要であることが判明した。

走化性乱歩モデル

昆虫の幼虫、細胞性粘菌、 バクテリアなどの生物では、 摂食、 防衛、 移動などそれぞれの適応的理由から、 フェロモンを出しながら同種個体同士が集合する行動が一般に知られている。 このような走化性集合行動のダイナミックスは、 Keller and Segel (1971) 以来、 主に反応拡散モデルとして研究されてきた (他に引用文献 [2],[3] など)。

しかし、 反応拡散モデルは、 本来は確率的である生物個体のランダムウオークを、 決定論的な個体密度の拡散に置き換えて近似したモデルである。 したがって、 生物個体の数が決定論的に論じ得るには少ない場合には、 そのような置き換えはもとのモデルの性質を調べる上で適当ではない。

このような観点から、 走化性集合行動のダイナミックスを記述する試みとして、 次のような非常に単純な 1 次元系モデルを用いることにする。

92(なお、 以後は簡単のため、 生物個体のことを『虫』 とよぶが、 昆虫に限定したモデルを考えているわけではない。)$\frac{d\xi_{i}}{dt}=\chi\frac{\partial c}{\partial x}+\sqrt{2\mu}\frac{dB_{i}}{dt}$

(1)

ただし、

$\xi_{i}=\xi_{i(t)}$は番目の虫の座標、 $B_{i}=B_{i}(t)$ は 番目の虫のブラウン運動、
$i$

$i$

$c=c(t, X)$ は場所 x でのフェロモン濃度である。右辺第一項はフェロモンの濃度勾配に比例する虫の走化性、 第二項は虫の動きの確率的な乱歩の効果を表す。
また、 フェロモンの濃度 $c=c(t, x)$ に関しては、 ふつうの反応拡散方程式に従うことにする。 2$\frac{\partial c}{\partial t}=D\frac{\partial c}{\partial x^{2}}-kc+fa$(2)上式で、 右辺第一項はフェロモンの拡散、 第二項はフェロモン濃度に比例したフェロモンの分解、 第三項は虫によるフェロモンの生産をあらわす。 ここで、 虫の個体密度$a=a(t, X)$ は、 デルタ関数\delta $(x)$ を用いて、$a(t, x)=\Sigma\delta(x-\xi_{l}\{t))$(3)$i$となる。パラメータ$\chi,$

$\mu$

$D$

$k,f$ は全て定数であるとする。

方程式系 (1)、 (2)、 (3) を、走化性乱歩モデルとよぶことにする。

これと比較するために、 式 (1) に対応する反応拡散モデルを、$\frac{\partial a}{\partial t}=\mu\frac{\partial^{2}a}{\partial x^{2}}-\chi\frac{\partial}{\partial x}(a\frac{\partial c}{\partial x})$(4)のように考える。第一項は拡散、 第二項は走化性の項である。パラメータ$\chi,$$\mu$

は式 (1) で用いたものにそれぞれ対応している。式 (4) と、 フェロモン濃度の式 (2) を組み合わせると、反応拡散系のモデルになる。

一様平衡解の安定性についての反応拡散モデルとの差異

前述の走化性乱歩モデルが反応拡散モデルとどのように異なるふるまいをするのかをみ93るために、 一様平衡解の安定性を反応拡散モデルと走化性乱歩モデルで比較した。

反応拡散モデルについては、 ノイマン条件のもとでの安定性条件が、 Keller andSege1(1971)\d ag er )$t$$\mu\{D(\begin{array}{l}mL\end{array})+k\}>\chi\overline{a}f_{(n=}1,$2,)(5)である。 ただし、 L は領域の長さ、 a は領域で一様な虫密度を表す。

この条件が走化性乱歩モデルでも成立するかどうかをみるために、 いくつかのパラメータについてシミュレーションをおこない、 また、 比較のため反応拡散モデルでも同じパラメータを用いて数値解析をおこなった (図 l-la\tilde 1-4) 。ただし、 走化性乱歩モデルでは虫の座標を虫密度に直したグラフに示した。
まず、 条件 (5) が成り立つような領域では、 反応拡散モデル (図 l-la) でも走化性乱歩モデル (図 l-lb)

でも共に虫が集合し、 高密度のクラスタをつくる。 ただし、後で論じるが、図 l-lb に示したものはまだ定常状態には至っていない。次に、 条件 (5) が成り立たないようなパラメータでは、 反応拡散モデルでは一様平衡状態からの微小なずれは減衰する$2a)$$2b)$(図 1-。一方、走化性乱歩モデルでは虫密度は一様にはなっていないようにみえるが (図 1、この場合の時間変化をみると、 この図にみえている小さなクラスタは一時的なものであり、 持続的なクラスタは形成されることがないので、 反応拡散モデルの結果と較べて顕著な差異はないと思われる。
しかし、 これらの場合に対して、 実は、 条件 (5) が成立しなくても走化性乱歩モデルでは虫が強く集合するパラメータ領域があることが, シミュレーションによって確かめられた (例えば図 1-3b) 。また、 この場合、 反応拡散モデルの数値解析では、 かなり非一様な初期状態 (例えば図 1-4) からでも一様平衡状態 (図 1-3a) に到達した。

このように走化性乱歩モデルが反応拡散モデルと質的に違った性質を示すパラメータ領域は、 かなり広い範囲にわたっている。 そこで、走化性乱歩モデルの性質をさらに調べるために、 (このモデルで虫が集合するようなパラメータにおいて) 空間を縦軸、 時間を横軸にとり、 各虫の座標変化をおってみた (図 2) 。これより、 走化性乱歩モデルによる虫の集合は次のような性質をもっていることがみてとれる。

94

$1$

、比較的早い段階で、 互いに近いところにいる虫が寄り集まって局所的なクラスタをつくる。 このクラスタは安定であり、 いったん属した虫が離れることはほとんどない。
$2$、互いに充分近くにあるクラスタは、 フェロモンによって引き合い、 やがて融合して、より大きなクラスタをつくる。
$3$、近くに他のクラスタがない場合、 クラスタは全体としてまとまりを保ちつつランダムな動きをする。

ランダムな動きによってたまたま近づいたクラスタ同士は融合するので、 最終的な状態、(安定状態) では、 クラスタが融合して 1 個ないし 2 個の大きいクラスタをつくり、 それが片方または両方の境界に貼りつくと考えられる。 しかし、 クラスタの融合がすすむにつれて、残ったクラスタ同士の距離は大きくなり、 また後で述べるように、 クラスタの規模が大きくなるとランダムな動きは急に小さくなるので、 それ以上はクラスタの融合がなかなかすすまない状態、 いわば『準安定』 な状態が長く保たれることになる。
この準安定状態に比べると、 各個体の独立した動きは非常に初期の段階にしかみられない。 したがって、 このモデルを記述するためには、 クラスタを単位としてそれらのふるまいを調べることが有効であると考え、、クラスタが他のクラスタの出すフェロモンにひかれて動く (定向移動)、他のクラスタの影響が無視できる場合にクラスタがランダムな動きをする (ランダム$1$$2$ドリフト)というクラスタの特徴的なふるまいについて、 さらに詳しく調べることにした。

定向移動
いま、 $N$ 個体からなるクラスタが、 そのクラスタの外部に由来するフェロモンの勾配から影響を受けている場合を考える。 クラスタに属する虫が生産したフェロモンの濃度を$c(x, t)$
、クラスタの外に由来するフェロモンの濃度を $c_{out}(x, t)$ とする。 このとき、 虫の走化性による速度成分は$\chi\frac{\partial}{\partial x}(c+c_{out})=\chi\frac{\partial c}{\partial x}+\chi\frac{\partial c_{out}}{\partial x}$(6)

もし、$\chi(\partial c/\partial x)=0$

なら、 クラスタを構成している虫は平均速度$v_{0}=\chi(\partial c_{out}/\partial x)$

95で移動し、 結果としてクラスタ全体が速度$v_{0}$で移動するはずである。 しかし、 シミュレーションによれば、 クラスタの移動速度はクラスタが大きくなるにしたがって急に小さくなった。 このことから、 クラスタが外からのフェロモンに引かれて移動するとき、 自分自身の出したフェロモンの影響をうけて動きが抑制されていることを示す。 つまり、 自分自身がしばらく前に生産したフェロモンは、 その時点の位置で濃度が高いためにクラスタの動きをとどめる方向にかかり、 あたかも動摩擦のような作用をする。以下ではこの現象を『摩擦』 とよぶことにする。
この摩擦を定量的に理解するために、 次のような仮定をおいて、 反応拡散方程式をもとにした解析を試みた。 すなわち、、虫はほとんど一点と見なせるくらいに狭い範囲に集中している。

、クラスタ外部由来のフェロモンの濃度勾配は一定である。

$1$

$2$

仮定 2 は、外部由来のフェロモンの生産地点が十分離れていると考えることによって正当化される。 この仮定により、 クラスタが動きはじめてからしばらくすると、 クラスタの移動速度は一定になり、 このときクラスタの出したフェロモンの濃度分布は形を保ちながらクラスタと同じ速度で進む進行波になると思われる。 つまり、 このときの速度をくとフェロモンの分布は $C(X, t)=c$ (x- $vt,$$0$$v$とお) となる。

そこで、 $s=x- vt$ とおくと、 仮定 1 より、 クラスタにはたらく摩擦は、 クラスタの中心の座標におけるフェロモンの勾配のみを考慮すればよいから、$v=v_{0}- \chi(\frac{dc}{ds})_{s=0}$(7)

また、 式 (2) は に関する常微分方程式になって、$s$-v

何 $=D$2$d$$ds$$- kc+ffl\delta(s)$$ds^{2}$と書ける (ただし、(8)$\alpha_{S}$) はデルタ関数) 。これを解くと、 定常状態におけるクラスタ由来のフェロモンの分布は、$c(s)=\{M_{exp\{-(v- p)s\}}^{\rho}L_{\rho^{-exp\{-\frac{1}{\frac 2^{1}D2D}(v+\rho)s\}}}^{N}$$(s>0)(s<0)$

(9)

96
ただし\mbox{\boldmath $\rho$}$=\sqrt{v^{2}+4kD}$

この分布が図 3 である。 この分布はクラスタ後方に偏った非。

対称な形であり、 クラスタは自分の出したフェロモンに引き戻されている。
しかし、 フェロモンの分布は原点 $(s=0)$ において不連続であり、 式 (7) は計算できない。 そこで、 まず、虫が充分狭い有限の範囲 [-e, ] に分布し、 そこでは滑らかな分布関$\epsilon$数をもつと考える。 すなわち、 式 (7) のかわりに、$- v=- \frac{1}{N}\{\int_{\epsilon}^{0}\chi(\frac{dc}{ds})a(s)ds+\int_{0^{\epsilon}}\chi(\frac{dc}{\ })a(s)d_{S} \}$$v_{0}$(10)を用いる。式 (9) と (10) を用い、 $a(s)/N$ がデルタ関数に近づく極限をとると、$v_{0}=v+\underline{\chi vffl}$$v=v_{0}-\underline{x\Psi^{N}}$$D\rho$$Dp$または

この式では実現された速度$v$$\text{し^{}\rho=}\sqrt{v^{2}+4kD}$

ただが外から引かれる速度

ラスタの動きが遅いことを考慮して、

$v$

(11)。

$v_{0}$の陰関数になっている (図 4) 。クが十分小さい近似により$v=(1+ \frac{\chi ffl}{2k^{\iota a}D^{312}})^{-\iota_{v_{0}}}$(12)とかける。 この式から、 摩擦による速度の減少の程度は、 個体数 N、 フェロモンの感度\chiおよび生産速度 f が大きいほど、 またフェロモンの分解速度 k や拡散係数 D が小さいほど、強いことが読みとれる。

以上の結果を、 走化性乱歩モデルのシミュレーションと比較してみた (図ミュレーション結果は$v_{0}$$5- 1^{-}3$) 。シの 3 つの値についていずれも (12) の理論値がよい近似であることを示しているが、 理論値の速度の方がやや小さい。 これはクラスタをなす虫の分布にわずかの広がりがあるため、 もしくはシミュレーションのクラスタの動きが確率的要素を含むため、 などが考えられるが、 さらなる研究が必要である。

ランダムドリフト

次に、 クラスタが他からの影響を受けない場合を考える。前に述べたように、 構成する虫の数が有限であるために、 クラスタは全体としてまとまりを保ちつつもランダムな動き

97をする。
ここで、 このランダムドリフトの大きさを示す量として、 一定時間 $T$ の問にクラスタ(13)各個体のブラウン運動の係数か村 7\mbox{\boldmath $\varpi$} であったことを考えると、個体間の相互作用が無視できれば、 $N$個体からなるクラスタの中心位置のランダムドリフトは、 $V_{T}=2\mu T/N$である。 しかし実際の走化性乱歩モデルでは、$V_{T}$は $T$ には比例しているものの、 比例定数は $2\mu/N$ よりはるかに小さい (図 6) 。よってランダムドリフトもまた、 前節の定向移動と同じく、 フェロモンによるクラスタの動きの抑制、 つまり摩擦の影響を受けていることがわかる。
そこで、 まず、各パラメータ ($\chi,$$\mu,$$D,$$k,$$N_{)}$のさまざまな値について $V_{T}/T$ を計算し、 両対数グラフで表示した (図 7) 。すべての回帰直線がほぼ平行であることは、$V_{T}/T$ が各パラメータの巾に比例していることを示す。 そこで重回帰により指数および定数要素を推定すると、

20

1.1

$V_{T}=4.8 \frac{D\mu}{(\chi f)^{1.7}N^{2.7}}T$(14)となった。

式 (14) で注目すべきは、 $N$ に関する依存性である。 $V_{T}/T$ は、 フェロモンがなかった場合の $N^{-1}$ にではなくて、 $N^{-2.7}$ に比例する。 したがって、 クラスタ同士の融合などにより $N$ が大きくなると、 $V_{T}/T$ は急激に減少する。 このことは図 2 でもみることができる。

考察
この研究では、 生物個体の走化性集合行動について、 反応拡散モデルの適用が不適当である場合を示した。 したがって、 反応拡散モデルがどの範囲までどれだけの信頼性をもって適用できるかを調べることが重要になるが、 確率論にもとつく解析が有効であろう。 また次のようなテーマがとくに重要な課題である。

98
一方、 走化性乱歩モデルにおいては、 局所的なクラスタの動態が重要である。 そのうち今回はクラスタの定向的な移動とランダムドリフトについて解析したが、 さらに、 クラスタの形成と融合の過程の解析も興味深い。 しかし、一様に虫が分布した状態から局所的クラスタの位置や規模を初期状感から予測することは、各個体の動きのランダム性により困難である。一方、複数のクラスタが融合していく過程は、集団遺伝学における coalescentprocess に類似している。
今回シミュレーションを行なった範囲では、 いったん形成されたクラスタはばらばらになることなく非常に安定に存在していたが、 虫の乱歩のために充分長い時間がたてばクラスタが崩壊する可能性がある。形成されたクラスタの安定性についても確率論による解析が必要である。
さらに、 ランダムドリフトの単位時間あたり移動距離の分散$V_{T}/T$ は、 フェロモンがない極限では $2\mu/N$ になるはずだが、 式 (14) はそうなっていない点で不十分である。

謝辞
本研究は、 江副の京都大学理学部における卒業研究課題から発展したものである。

ご指導いただいた重定南奈子教授、 梅田民樹博士、 および理論生物物理学研究室の皆様に感謝いたします。 また以下の方々からは、 有益な助言をいただきました。松田博嗣、 三村昌泰、 飯塚勝、 佐々木顕、 および九大数理生物学研究室の皆様。

引用文献[1] Keller, E. F., Segel, L. A. (1970) Initiation of slime mold aggregation viewed asinstability. J. Theor. Biol. 26, $399A15$[2] Deneubourg, J. -L., Gr\’egorie, J. -C., Le Fort, E. (1990) Kinetics oflaval gregarious

behavior in the bark beetle Dendroctonus micans (Coloptera:Scolitidae).
J.Insect Behav. 3, 169-182[31 Myerscough, M. R., Murray, J. D. (1992) Analysis ofpropagating pattern in achemotaxis system. Bull. math. Biol. 54, 77-94[41 Okubo,A.,(1980) Diffusion and Ecological Problems: Mathematical Models.

Springer-Verlag, Berlin Heidelberg New York

99

$0$

図 l-la40$20121648 \frac{\ovalbox{\tt\small REJECT}_{\#}[\Uparrow_{[1}|\Vert\iota\Lambda_{\Lambda}\iota\Vert\Lambda_{111_{\Lambda}}|\ovalbox{\tt\small REJECT}}{0\backslash \text{図^{}\backslash }1- 2b40}$$0$40

図 1-4. 図 1-3a に到達する初期状態の例.
$0$

図 1-3a40

図 1; 反応拡散モデルと走化性乱歩モデルの比較.
; 反応拡散モデル,b; 走化性乱歩モデル.
$a$$D=4.8,fi– 6.0,$ $k=6.0,$ $N=40$l-lab; $\chi=2.4,$ $\mu=1.2$ ,1-2ab; $\chi=0.6,$ $\mu=4.8$ ,1-3ab および 1-4; $\chi=1.2,$ $\mu=2.4$

100

図 2. 虫座標の時間変化. 横軸が時間、 縦軸が空間座標.
$N=20,$ $\chi=2.4,$ $\mu=1.2,$ $D=4.8,$ f–6.0, $k=3.0$$D=9.6,f-6.0,$ $k=1.5,N=40,$ $\Gamma-3.0$$\chi=2.4$, その他のパラメータは図 3 と同じ.

101

図 6.

$T$と $Vr$ との関係.

$\chi=2.4,$ $\mu=1.2,$$D=4.8$$k=6.0,$ $f=6.0,$ $N=40$

図 5-2.

図 5. 解析的モデルとシミュレーシ
ョン結果の比較. 時間に対するクラスタの移動距離を表す. 太実線が理論値、細実線が各 3 回のシミュレーション結果、点線は摩擦がな

図 7. パラメータと $V/Vr$ との関係.
図では異なる個体数 $N$ に対してかったとした場合 $(x=vot)$$k$.

図 5-1. $V\sigma 0.06$
図 5-2.

$\mathcal{V}\sigma 0.3$

$N=40$.

についても同様なグラフを描き、 それぞれほぼ平行な回帰直線を得た.
$(\chi, \mu, D, N)$

図 5-3. $Vt=0.6$$\chi=2.4,$ $\mu=1.2,$を変化させた回帰直線をとっている. これらの直線はほぼ平行になっている. その他のパラメータ$D=9.6,f– 6.0,$ $k=1.5$ ,

102

$Ag$

gregation with Chemotaxis-Random Walk

Hideo Ezoe and Yoh Iwasa
Department of Biology, Faculty of Science, Kyushu University, Fukuoka 812, Japan

Group formation of insect larvae by secreting aggregation pheromone is studied by achemotaxis-random walk model (CR-model), in which the location of eachindividual is followed. If the number of individuals is mildly large (several tens),
the model’s behavior is different from the one predicted by the corresponding
reaction diffusion model (RD-model). For example, even if uniform distribution is astable steady state for the RD-model, in the corresponding CR-model with uniform
initial distribution, animals may quickly form many small clusters; then theseclusters drift randomly, gradually merge with each other, and finally form one or a
few large clusters. Next, the directed and random movement of a single cluster ofanimals is analyzed in detail. Pheromone produced by a cluster tends to slow down
the moment of the cluster. A formula for the maginitude of slowing down effect ofattracting pheromone, we call ‘Triction“, is derived using an approximation andconfirmed by computer simulation.