まずは,先週終わらなかったところからー
個体成長とカオス
2 ロジスティック成長
指数的成長を抑制する方法は種々あります。 個体群が大きくなれば資源がそれにし たがって減少したり,栄養状態が悪化したりして,成長率も小さくなります。 成長率は r と K を正の定数として
という ロジスティック方程式 で表わされます。(2)式との違いは,rx の項と (1-x/K) の項とがあることです。 このとき x が小さい時は, x/K が無視できるくらい小さいので,増加率 r に従って指数的に増加します。 つまり (2) に従って指数関数的に増加します。 r は個体密度が小さくて環境資源が十分にあるときの増加率を示しているとも考えられるので, 内的自然増加率 intrinsic rate of natural increase と言います。 一方 x が K に近づくと増加率は 0 に近づきます。 すなわち平衡の個体数 K に達します。
K より大きな値 x(0) で出発した時は,減少して K に近づき, 逆に x(0) が小さな値で出発した時は S 字型のカーブを描いて K に収束します。 K はその環境中に維持できる個体数という意味から, 環境収容力 carrying capacity と言います。
ロジスティック方程式は簡単に解析ができて,x=0 または x=K ならば dx/dt は 0 です。 この場合個体密度は変化しません。 0 < x < K に対しては,x は増加し x > K では減少します。 この解は
として,部分分数の和にし,
項別に積分すれば
だから結局,
対数の加法公式から
対数を元に戻して
となり
x について解くと
です。t=0 のときの x の値を x(0) とすれば
となります。この式を C について解けば,C が x(0) から求められます。
これを(12)に代入すると,初期問題の解が求められます。
念のためこの式で t=0 にしてみると
です。 (15)で,t=0 とすれば,
実習
以下のように,今日の自習教材をコピーして,
cp -rp ~asakawa/20100416 . cd ./20100416
としてください。そして,
java LogisticEquation
とすれば,ロジスティック成長のグラフを描く Java のプログラムが走ります。 ここで 初期値 x(0),環境収容力 K,内的自然増加率 r を変化させて, いろいろ遊んでみてください。
ロジスティック方程式 (4) には,x = 0 と x = K という 2 つの平衡点 equilibrium が存在します。平衡点とは,ある系がある状態におかれたとき, その後もずっとそこにとどまる点を指します。 これは dx/dt = 0 から求めることができます。 x = K の平衡点は,個体数がずれたとしても,時間がたつと元の状態に戻っていきます。 この意味で x = K は 安定な stable な平衡点であるといいます。
一方,生物の不在を表す x = 0 の点では,わずかな生物が侵入してきても, 時間とともに x = K に写ってしまいます。従って x = 0 は, 不安定な unstable 平衡点と呼ばれます。
2 つのパラメータ K と r が栄養状態や個体密度などの環境条件によって, どのように変化するのかを調べることによって,次のようなことが分かってきたそうです。
すなわち,植物の平均重量は,ロジスティック方程式に従って成長します。 その最大サイズ K は密度と反比例します。 その結果,平均個体重と密度の積が最終的に一定になる,ということです。
練習問題
r=1, K=1 とすれば,式(4)は,
となります。この式を解いてください。
低密度の影響 アリー効果
個体密度が高いことは必ずしも悪影響を与えるばかりとは限りません。 なぜなら,個体密度が低過ぎると生物個体間の協力的相互作用が得られず, 死滅してしまうかも知れないからです。 たとえば,昆虫で,一匹の幼虫だと葉への食いつきができなくて死んでしまいますが, 数個体の幼虫が一緒にいるとうまく生きられたりするそうです。 このような低密度の悪影響のことを アリー効果 Allee effect と呼んでいます。 このときの個体群動態方程式として,
を考えてみましょう。ここで 0 < a < K とします。この式は, dx/dt = 0 という平衡点から x = 0, x = a, x = K と 3 点あることになります。 下の図 からわかるとおり,x = 0 と x = K とは安定な平衡点です。 一方,中間の x = a は不安定な平衡点であると言います。
生物がいない場所にわずかな個体数が侵入してきても, 低密度のため増加できない場合があることを示しています。 しかし,侵入個体数が a を越えると,環境収容力 K まで増加して定着することができます。 このような a を しきい値 と呼びます。
実習
同様にして java による実習をしてみましょう。
java AlleeEquation
としてください。 初期値 x(0), 環境収容力 K, 成長率 r, しきい値 a を様々に変化させることによって, グラフがどう変わるのかを実感してください。
離散時間のロジスティック写像
世代が重ならない個体群の場合(例えば,1 年に1 世代を過ごす昆虫)を考えます。 Xt を t 世代目の個体数とします。 連続時間のロジスティック回帰と同じく r を成長率。K を環境収容力として, t + 1 世代目の個体数は次式
で表せます。個体数が低密度の場合には, 毎世代ごとに 1 + r 倍で増加しますが, 個体密度が高くなるにつれて増加率が落ち, 環境収容力 K を超えると減少するようになります。
安定不動点,不安定不動点,カオス
さて,先にも記したとおり,個体成長が,ある法則 f に従うとすると, 現在の状態を x として,次の世代の状態は y は, y = f(x) と表せました。 このとき,f は,個体数が少ない時には,成長し(増加し),多い時には減少する と考えれば,y = f(x) という関数は,直線によって表せるものではなく, コブをもった関数になるでしょう。
一つだけコブを持った関数としては,一番簡単なのは, おそらく中学生で学習する放物線,すなわち x に関する 2 次関数 y = a x2 + b x + c でしょう。 そこで,以下のような単純な関数に関して, 次世代がどのように変動するかを見ていくことにしましょう。
0 と 4 との間の R に対して写像
は区間 [0, 1] 上で力学系を与えます。 点 0 は明らかに不動点です。 R < 1 ならば,全ての x ∈ (0, 1] に対して xt+1 < xn です。 従って,x の起動は 0 に向かって単調に減少します。 今後は,R > 1 の場合だけを考えます。 F のグラフは点 0 と 1 で x 軸と交わり, 点 1/2 で最大値 R/4 をとる放物線です。 放物線は一辺の長さ 1 の正方形の内部で, 直線 y = x と唯一の点 P で交わります。 P の横座標値は p は pn+1 = pn を満たすので, もう一つの浮動点 p = (R−1)/R が得られます。
平衡状態は xt+1 = xt とおいて計算されるので, 図上では y = x との交点として求められます。 xt に対する次世代 xt+1 の値は曲線を使って求められますが,次に縦軸にある値を横軸に戻すには y = x を用いると直感的でわかりやすいです。このようにして, xt+2, xt+3 と順次グラフを使って, 個体数を求めることができます(下図参照)。
p の安定性に関しては,平均値の定理より,x と p との間に適当な c を置くと
が成り立ちます。x (従って c ) が p に十分近ければ, のとき が成り立ちます。従って,
となります。 すなわち,x' は x より浮動点 p に近いのです。 このことは,もし x が p のある適当な近傍から選ばれれば,x の軌道 (22 から反復的に得られる数列)は p の近傍に留まり, さらに p に収束することを意味します。 このような意味で,不動点 p は 漸近安定 であると言います。 逆に, であるならば,
となるので,x の軌道は不動点から遠ざかります。 このような場合 p は 不安定 であるといいます。
今,F(x) = R x (1-x) なので, となり,p は
- 0 < R ≤ 1 では単調収束する
- 1 < R ≤ 2 では 1 - (1/R) に単調収束する
- 2 < R ≤ 3 では 1 - (1/R) に振動しながら収束する
- 3 < R ≤ 3.5699 では R は 2 のベキ乗の周期点を振動する
- 3.5699 < R ≤ 4 では変動は不規則になり,特定の周期を持たない,すなわち カオス になる
カオスでは, ほんの少しだけ違った初期値ではじめると十分時間が経ったときに, その値は予想できなくなります。 3.5699 はカオスの境界を表しており, ファイゲンバウム点 と呼ばれたりします。
R = 4 の時カオスになる
大事な点は,計算可能生は予測可能生を意味しないということです。 決定論的な運動がランダムな動きと区別できない ことが重要です。
R=2.9 のときのグラフを次に示します。
R=4 のときカオスとなります。これを下の図に示します。
実習
端末エミュレータ(Windows ではコマンドプロンプト)から
java LogisticMapping
とタイプしてください。(22) すなわち xt+1 = R xt ( 1 - xt) の計算が始まります。R と初期値を変化させて遊んでみてください。 特に,R がカオス領域のときに,わずかな初期値の変化が, ドラスティックな変化をもたらす点に注目してください。
参考文献
- カール・シグムンド/冨田勝監訳. (1996). 数学でみた生命と進化. 東京: 講談社ブルーバッ クス.
- マーティン・ノーヴァック/竹内康博・佐藤一憲・中岡慎治訳. (2008). 進化のダイナミクス -生命の謎を解き明かす方程式. 東京: 共立出版
- ジョン・メイナード=スミス/巌佐庸・原田祐子訳. (1995). 進化遺伝学. 東京: 産業図書.
- 巌佐庸. (1990). 数理生物学入門 生物社会のダイナミックスを探る. 東京: 共立出版.
- ホッフバウアー,ジグムンド/竹内康博・佐藤一憲・宮崎倫子訳. (2001). 進化ゲームと微分方程式. 東京: 現代数学社.
ロトカ・ヴォルテラ方程式
ロトカ・ヴォルテラ方程式の発見
第一次世界大戦中, アドリア海はイタリア海軍とオーストリア・ハンガリー帝国との戦いのために, その海域での大規模な漁業は全て中止させられてしまったそうです。
その結果,予期せぬ事態が起こりました。 終戦後, 数年してからイタリア人の生態学者ダンコナ D'Ancona は, 漁業市場で漁獲量の統計を取った結果, 戦争後サメなどの肉食魚の割合が, 戦争前よりもかなり高くなっていたことに気づいたのです。 ダンコナは困惑しました。なぜ戦争はサメを好むのだろうか?と。
ダンコナは,ローマの数理物理学の教授であり, その当時第一線の数学者として広く知られていたヴォルテラ Senator Vito Volterra に問い合わせました。 その結果,今日ロトカ・ヴォルテラ方程式として知られる微分方程式を使って, アドリア海における生態系の説明モデルが産まれました。
アドリア海には多数の魚が棲息していますが,単純化すると 2 種類に分けられます。 サメなどの捕食者とその餌となる魚です。 ヴォルテラはこれら 2 種類の魚が相互に依存していることを見抜きました。 すべての捕食者は,その獲物(被食者)を食べることで,獲物の増加率を減少させます。 もし,捕食者の数があまりに増え過ぎると,その獲物である被食者の個体数は減少することになります。 逆に獲物(被食者)の個体数が大きくなると,捕食者は繁栄します。 また,もし被食者の数が少なくなり過ぎると,捕食者は子孫をつくる前に餓死してしまい, その結果捕食者の個体数は減少することになります。
このようにして捕食者であるサメの個体数は増えたり減ったりします。 しかし絶滅することはありません。 なぜなら,サメの個体数が減ると必ず餌(被食者)である魚は増殖します。 そうなるとサメ(捕食者)は餌である魚(被食者)を捕まえやすくなるからです。 すなわちサメの個体数が増えるようになります。 また,もし魚(被食者)の個体数が少なければ, サメの個体数も食料不足のために減少して, 結果的に餌である魚(被食者)にとって楽な環境になるというわけです。
つまり,サメと餌である魚の個体数は増えては減り,減っては増えるという具合に循環するわけです。
さてここで,この捕食者と被食者のゲームに漁師という登場人物を加えることにしましょう。 漁師は,捕食者であるサメも被食者である魚もどちらも捕ることができるので, 両者の総数はともに減少します。 しかし,漁の後には,捕食者ー被食者の間のバランスが崩れ, 予期せぬ結果が生じます。 漁によって被食者である魚の数が減少すると,捕食者であるサメの個体数もまた減少します。 なぜなら,サメの餌を漁師に捕られて食料不足になってしまうからである。 つまり結果として,捕食者の個体数が減少するため, 獲物である魚の数がかなり驚異的に増えてしまいます。 もし戦争によって漁が中止されると,その逆のことが起こります。 つまり,捕食者の個体数は増え,獲物である被食者の数は少なくなるわけです。
これが今日「ヴォルテラの法則」として知られるようになった現象です。 すなわち,捕食者ー被食者の両者が互いに相手の個体数を決めるという前提があるなら, 両者の成長率の低下,たとえば漁など,は被食者の増殖と捕食者の減少を招くというものです。
ヴォルテラのモデルは, アメリカ人の化学者ロトカ Alfred J. Lotka によってすでに研究されていたことが分かりました。 それで,今日では,ロトカ・ヴォルテラの方程式と呼ぶようになったのです。
捕食者-被食者方程式
ヴォルテラは,捕食者が存在しない場合,被食者の個体群の成長がある定数 a で与えられ, さらに成長率は捕食者密度 y の関数として線形に減少すると仮定しました。 これから dx/dt = a - by ただし(a,b>0) が得られる。 被食者が存在しない場合,捕食者は死亡せざるをえないので捕食者の成長率は負です。 しかし被食者密度 x が正であるならば, 成長率は回復するので dy/dt = -c + dx (c,d>0) となります。 両者をまとめると次式を得ます。
この式がロトカ・ヴォルテラ方程式と呼ばれるものです。 この式の解は,次の 3 つはすぐに分かるでしょう。
すなわち,捕食者の個体密度 y も被食者の個体密度 x も, ある時刻 t において 0 であるならば,常に 0 です。 被食者が不在の場合 x(t)=0,捕食者は絶滅に向かいます (t → +∞ で y(t) は 0 に収束する。 捕食者がいなければ y(t)=0,被食者の個体数は爆発的に増えます (x(t) → +∞)
捕食者を横軸,被食者を縦軸にして解曲線を描いたものが 下図です。このように変量 x と変量 y とを縦軸,横軸にとって その変化の様子を描いたものを相図 Phase Graph と言います。
ロトカ・ヴォルテラ方程式の相図
捕食者の個体数が増えると被食者個体数が減り,その結果捕食者個体数が減り始めます。 その結果,被食者個体数が増加する,という円環状の関係になっています。 初期値が違えばこの円環の大きさが変化します。しかし,その形は変化しません。 円環はパラメータを変化させることで形状が変わります。
実習
いつものように実習してみよー。
cd cp -rp ~asakawa/20100423 . cd 20100423
とすれば準備は完了です。ロトカ・ヴォルテラ方程式のシミュレータを動かすには
java LotkaVolterra
とタイプすれば良いです。
シミュレーションを実行してみると分かりますが, ロトカ・ヴォルテラ方程式には平衡点が存在します。 すなわち (x,y)=(0,0) という捕食者と被食者ともいない場合, 加えて, 両者の個体数が均衡して変化しない点が存在します。 このような平衡点は は と を満足しなければなりません。従って
です。しかし,この平衡点から少しでも離れれば,振動を永遠に繰り返すことになります。
ロトカ・ヴォルテラ方程式の時間発展
同じロトカ・ヴォルテラ方程式の時間発展を描いたものが 下図です。
実習としては,
java LotkaVolterraTime
としてください。
ヴォルテラの原理
捕食者と被食者の個体密度は周期的に振動し,その振幅と振動数は初期値に依存して決まります。 しかし,個体数の時間平均は一定であり休止点 に等しくなります。すなわち
ここで T は解の周期とします。
ですから,両辺を積分して
すなわち
を得ます。x(T)=x(0) より
となります。x の時間平均についても同様に が成り立ちます。
以上により,ダンコナの疑問,なぜ戦争はサメを好むのか, サメなどの捕食者の個体数が, 戦争前よりもかなり高くなっていたことに対する答えの準備が整いました。
漁師による漁業活動は,被食者の増加率を減少させます。 a の代わりに a-k になるわけです。 同時に捕食者の減少率を増加させる c の代わりに少し大きな値 c+m になります。 しかし相互作用を表す定数 b, d は不変です。 従って,漁師による漁業活動が中断されている場合と比べて, 捕食者の個体数の時間平均は (a-k)/b で少し小さくなり, 被食者の平均は (c+m)/d で少し大きくなります。 戦争により漁業活動が中断されている間,この逆のことが起こり,捕食者は増加し, 被食者は減少したというわけです。
リヤプノフ関数
ロトカ・ヴォルテラ方程式が安定であるための1つの十分条件は リアプノフ関数 Lyapunov function, V(x,y) が存在することであるとされます。 逆に言えば,リアプノフ関数を見つけることができれば,その力学系は安定です。 前出のロトカ・ヴォルテラ方程式
のリアプノフ関数は,以下のとおりです。
この場合のリアプノフ関数 V(x,y) は時間とともに変化せず, 状態点は V(x,y) の等高線をたどりながら動きます。 すなわち一周すればかならず元に戻ります。 このようにシステムの状態の関数であって, 時間変化にともなって増えも減りもしないものを保存量といいます。 平衡点のまわりからずれていたとするとシステムの挙動は, 別の軌道に移り,平衡点には戻らりません。 しかし遠くには慣れてしまうわけでもありません。 このような平衡点は 中立安定 といわれます。
リアプノフ関数を見つけることは必ずしも簡単ではありません。 しかし,見つかればシステムの安定性解析に役立ちます。
ロトカ・ヴォルテラ方程式のリヤプノフ関数 a=1, b=0.002, c=1, d=0.02
種内競争を持つロトカ・ヴォルテラ方程式
もし,捕食者がいなければロトカ・ヴォルテラ方程式(2) の被食者は となり,その解は指数関数 となって個体数が爆発的に増加するという話は以前しました。 これは,被食者の個体数の増え方にロジスティック方程式 を仮定することで修正できます。 捕食者の種内競争も同じくロジスティック方程式に従うとすると, (2)の代わりにを考えます。 この場合周期解になるとは限りません。 パラメータのとり方によってはある種が絶滅する場合もありえます。 一例を下図に示します。
種内競争をもつロトカ・ヴォルテラ方程式
実習
それでは実習です。
java LotkaVolterra2
してください。 上の図を描くシミュレータが起動します。
アイソクライン
このときの の平衡点を解析するために少し変形してみましょう。 x=0, y=0 すなわち両個体がいないときには変化がないのはすぐ分かります。 それ以外は, 傾き 0 の アイソクライン isocline (iso とは「同じ」という意味であり,cline とは「傾き」という意味)
と,同様に のアイソクライン
を考えます。 この二直線の交点が平衡点となる。平衡点の座標は となります。 このとき x と y との増減は以下の表のようになります。
x, y とも減少 | x は増加,yは減少 | |
x は減少,y は増加 | x,yともに増加 |
すなわち,アイソクライン線の上下で x と y との増減を考えることができるわけです。 この方法のことをアイソクライン法と呼んでいます。
ここで,パラメータをいじっていろいろ遊んでみましょう。 たとえば,a=8, b=1, c=-6, d=-1, e=1, f=1 を考えます。 この場合, だけが安定な平衡点になります。 パラメータにマイナスを使っているので, もはや捕食者,被食者と呼ぶのは正しくないかも知れません。 そこで x 種 と y 種と表現することにすると,この場合 y 種は滅びる。 最初にどんなに y 種の個体数が多くても は不安定平衡点であるので y 種は生き残ることができないのです。
a=8, b=1, c=-6, d=-0.3, e=1, f=1 を考えます。 この場合,両アイソクライン線の交点が安定な平衡点になって, どのような初期値から出発しても (3,5) に収束します。すなわち x 種と y 種との共存が実現するわけです。
a=8, b=1, c=-10, d=-1, e=0.5, f=1 を考えてみましょう。 この場合,両アイソクライン線の交点は不安定平衡点になる。初期値によって y 種だけが生き残る か,x 種だけが生き残る に収束する双安定になります。
このようにパラメータを変化させると,2 種の共存か, 競争によって 1 種だけが勝ち残ることが生じます。 たとえば,2 種のトカゲが多数の島に棲んでいたとしましょう。 気候条件や栄養状態, 島の大きさなどがほぼ同じであっても, 両種が生息している島がないとすれば, 競争関係のために共存できなかったと考えられるわけです。 そのため,環境の変動や個体数などのランダム性のため島の個体群が滅びることが起こりというわけです。
もう一つの安定なロトカ・ヴォルテラ方程式,リミットサイクル
ロトカ・ヴォルテラ方程式をもう少し現実的にしてみましょう。
というモデルを考えます。 オリジナルのロトカ・ヴォルテラ方程式からの変更点は, 以下のとおりです。
- 捕食者のいないときに被食者は無限に増殖することはせず, ロジスティック方程式に従って環境収容力 K に収束する。 これは y=0 とおけば確かめられます。
- 捕食者の摂食速度を としたこと。 これは,餌が多いときには捕食速度は餌の量とともに増加します, あまりに多いと食べきれなくなるため飽和し,上限値 a/h を持ちます。
K が有限で h=0 であれば被食者の個体群を安定させる働きによってロトカ・ヴォルテラ方程式 の振動は止まります。
アトラクター
一方,捕食速度の飽和傾向だけを考慮すると K = ∞, h>0 システムは大域的に不安定になり,振動の振幅が時間とともに大きくなる。
この両方の効果が加わると 極限周期道(リミットサイクル) limit cycle が現れます。
リミットサイクル r=0.18, K=30.0, a=0.02, h=0.1, b=0.0504, c=0.25
実習
java LotkaVolterra3
してみてください。アトラクターやリミットサイクルが現れるパラメータを探してください。
殺虫剤に害虫の駆除は効果があるか
害虫 x と天敵 y を考え,
に殺虫剤の効果 -mx が加わったモデルを考えてみましょう。
殺虫剤を散布すると害虫 x は減るのでしょうか。 もちろん殺虫剤を散布すれば害虫の数は減ります。 ところが(28)の平衡状態を調べると m があっても平衡状態は変化しないことが分かるのです。 つまり害虫の増殖率を下げる効果が天敵 y のレベルを下げ, その結果,害虫の増殖率を上げるように働く効果があるのです。 これらの効果が打ち消しあってしまいます。
実際には,殺虫剤の効果 m は,天敵 y の生存率も下げるでしょう。 すると(28)の右辺にも -my という項を付け加えるモデルが考えられます。
この場合,殺虫剤の散布はかえって害虫 x の個体数を増加させてしまいます。 このように, 殺虫剤の散布がかえって害虫の個体数を増加させてしまうことはしばしば見られることなのです。 害虫が突然変異によって殺虫剤への耐性を進化させたと考えるよりも, 案外こういうところに原因があるのかも知れません。
ロトカ・ヴォルテラモデルの拡張
捕食者と被食者が多様な生態系では,ロトカ・ヴォルテラ方程式の多変数拡張が 用いられます。
種間の相互作用 aij は,i 種が j 種の捕食者であるなら
a
たとえば,光合成によって有機物を作る植物は生産者と呼ばれ,それを食べる 一時消費者(草食動物),さらにそれを食べる2次消費者(肉食動物)というように 多様な生態系が形作られます。 この(31)の安定性解析を行うことで,現実の生態系の実体に近づくものと思われます。
参考文献
- カール・シグムンド/冨田勝監訳. (1996). 数学でみた生命と進化. 東京: 講談社ブルーバックス.
- バージェス, & ボリー著/垣田,大町訳. (1990). 微分方程式で数学モデルを作ろう. 東京: 日本評論社.
- ホッフバウアー, & ジグムント/竹内康宏,佐藤一憲,宮崎倫子訳. (1998). 進化ゲームと微分方程式. 東京: 現代数学社.
- 佐藤總夫. (1984). 自然の数理と社会の数理 I,II. 東京: 日本評論社.