📂動力学

ロトカ＝ヴォルテラ 捕食者-被食者モデル

概要

ロトカ-ボルテラの捕食者-被食者モデルは種間の相互作用をシステムとしてモデリングし、特に捕食者-被食者モデルは二種の捕食関係を示す。二種についてのみ扱うとその拡張は無限にあるから、食物連鎖を表現するには十分だ。

モデル¹

$$\begin{align*} \dot{x} =& a x - b y \cdot x \\ \dot{y} =& c x \cdot y - d y \end{align*}$$

変数

• $x(t)$: $t$の時点での被食者集団$x$の個体数を表す。
• $y(t)$: $t$の時点での捕食者集団$y$の個体数を表す。

パラメータ

• $a>0$: $x$の繁殖率だ。
• $b>0$: $x$に対する$y$の捕食率だ。
• $c>0$: $y$が$x$を捕食することで得られる繁殖率だ。
• $d>0$: $y$の死亡率だ。

導出

$$\dot{N} = rN$$

まず、マルサス成長モデルから始めよう。被食者$x$は、天敵がいない場合、繁殖率$a>0$で成長すると仮定し、捕食者$y$は、餌がない場合、死亡率$d>0$で減少すると仮定しよう。これは厳密に連立方程式ではないが、次のように式で表すことができる。

$$\begin{align*} \dot{x} =& a x \\ \dot{y} =& - d y \end{align*}$$

$y$による$x$の捕食行為は、$x$、$y$両方に比例すると仮定し、どんな定数$h>0$に対して項$hxy$として表すことを試みる。$x$は捕食行為が起こるほど個体数が減少し、$y$は捕食行為をすることで生き延び、繁殖することができるので、次のように式を修正する。

$$\begin{align*} \dot{x} =& a x - hxy \\ \dot{y} =& hxy - d y \end{align*}$$

ここで、捕食行為が必ずしも100%の効率を出すとは限らないことに注意が必要だ。つまり、食べた分がそのまま繁殖につながるわけではなく、その前に$x$と$y$は単に個体数を表しているので、正確に「摂取量」を反映していないこともある。したがって、$x$、$y$それぞれに適した$h$が別途適用されるように式を修正する場合、ある$b,c >0$について次のように表せばいい。

$$\begin{align*} \dot{x} =& a x - b y \cdot x \\ \dot{y} =& c x \cdot y - d y \end{align*}$$

ここで、$by \cdot x$は$x$の瞬間死亡率が天敵の数に比例した$by$であり、$cx \cdot y$は$y$の瞬間成長率が餌の数に比例した$cx$という意味で、理解しやすいようにわざと別に書いたにすぎない。

■

カスタム

被食者$x$がプランクトンや雑草のような種の場合、天敵がいないときにマルサス成長すると言うのも意味があるかもしれないが、ウサギや羊のように草をめぐって同種間でも競争が起こると考えると、マルサス成長が非現実的な仮定である可能性がある。

$$\dot{N} = {{ r } \over { K }} N ( K - N)$$

それならば、単にマルサス成長をロジスティック成長モデルに変えることでモデルを改善できる。$x$が種内競争をするならば、例えばどんな$K > 0$に対して次のように修正すればいい。

$$\begin{align*} \dot{x} =& {{ a } \over { K }} x (K - x) - b y \cdot x \\ \dot{y} =& c x \cdot y - d y \end{align*}$$

定点

$$\overline{(x,y)} = \left({{ d } \over { c }} , {{ a } \over { b }}\right)$$

ロトカ-ボルテラの捕食者-被食者システムにおける非自明な定点は上記の通りだ。

存在性

$$\begin{align*} 0 =& a x - b y \cdot x \\ 0 =& c x \cdot y - d y \end{align*}$$

を満たす定点を見つけよう。$(x,y) = (0,0)$は自明解なので$x \ne 0, y \ne 0$を仮定して、一つ目の式から$x$を、二つ目の式から$y$を消去して

$$\begin{align*} 0 =& a - b y \\ 0 =& c x - d \end{align*}$$

が得られる。これを$x,y$について整理すると

$$\overline{(x,y)} = \left({{ d } \over { c }} , {{ a } \over { b }}\right)$$

■

安定性

非自明な定点についてだけ確認しよう。

$$\begin{align*} & \begin{bmatrix} a - by & -bx \\ cy & cx - d \end{bmatrix}_{\overline{(x,y)} = \left({{ d } \over { c }} , {{ a } \over { b }}\right)} \\ =& \begin{bmatrix} 0 & -b {{ d } \over { c }} \\ c{{ a } \over { b }} & 0 \end{bmatrix} \end{align*}$$

$\left({{ d } \over { c }} , {{ a } \over { b }}\right)$でのヤコビアンは上記のように計算される。その固有値は $$\begin{align*} \det \left[ \begin{bmatrix} 0 & -{{ bd } \over { c }} \\ {{ ac } \over { b }} & 0 \end{bmatrix} - \lambda I \right] =& \det \begin{bmatrix} -\lambda & -{{ bd } \over { c }} \\ {{ ac } \over { b }} & -\lambda \end{bmatrix} \\ =& \lambda^{2} - \left( - {{ bd } \over { c }} \right)\left( {{ ac } \over { b }} \right) \\ =& \lambda^{2} + ad \\ =& 0 \end{align*}$$ の解である$\lambda_{1,2} = \pm i\sqrt{ad}$だ。全ての固有値が純虚数なので、この非自明な定点はセンターだ。

■

視覚的理解

$$\begin{align*} \dot{x} =& {{ 2 } \over { 3 }} x - {{ 4 } \over { 3 }} y \cdot x \\ \dot{y} =& x \cdot y - y \end{align*}$$

上記のシステムを視覚化してみよう。400種類の初期値を無作為に抽出してその軌跡を描くと、次のようになる。

一目で非自明な定点$\displaystyle \left( 1, {{ 1 } \over { 2 }} \right)$があり、全てのフローが周期を持ちながら接近も離れもしないことが確認できる。

軸ごとにタイムエボリューションを見ると、上記のようになる。捕食者が減少すると被食者が増加し、被食者が増加すると捕食者が増加し、再び被食者が減少するサイクルを繰り返す。

応用

ボルテラの原理²

捕食者であれ被食者であれ、外部の介入によって$\varepsilon > 0$だけ収穫されたり除去されたりすると考えてみよう。例えば、捕食者と被食者が魚であり、同じ網の影響を受けるとか、両方とも害虫であり、同じ殺虫剤の影響を受けると仮定するのだ。

式を見る前に考えてみると、被食者は外部の介入によって死ぬことを考慮しても、それによって捕食者も減少し、抑制が弱くなり、捕食者はすでに外部によって死んでいくのに餌も減少するだろう。しかし、その逆が再び被食者の増加と捕食者の減少につながる可能性があるので、このシステムがどうなるか直感的に予想するのは簡単ではない。

今度は式をもう一度見てみよう。外部の介入は、単に個体数に掛けて引くだけで反映することができる。

$$\begin{align*} \dot{x} =& a x - b y \cdot x - \varepsilon x \\ \dot{y} =& c x \cdot y - d y - \varepsilon y \end{align*}$$

そして、変更されたシステムで新しく得られる定点は次のようになる。

$$\overline{(x,y)} = \left({{ d + \varepsilon } \over { c }} , {{ a - \varepsilon } \over { b }}\right)$$

この定点に対する分析は非常に興味深いが、捕食者と被食者を区別せずに公平に殺す場合、被食者の個体数は実際に増加するということだ。これをボルテラの原理と呼び、動力学的な分析が単純であるにもかかわらず、外部介入に対して明確なインサイトを提供していることが確認できる。

例示コード

• ソースコード

コード

以下は、本文の画像を作成したJuliaのコードだ。

cd(@__DIR__)

@time using Plots
@time using DifferentialEquations

function lotka_volterra_prey_pradator!(du,u,p,t)
    x, y = u
    a, b, c, d = p
    du[1] = dx = a*x - b*y*x
    du[2] = dy = c*x*y - d*y
end
u0 = [0.9, 0.9]
p = [2/3, 4/3, 1., 1.]
end_time = 20.

tspan = (0.,end_time)
prob = ODEProblem(lotka_volterra_prey_pradator!,u0,tspan,p)
sol0 = solve(prob; saveat=(0.0:0.01:end_time))

time_evolution = plot(sol0, vars=(0,1),
  linestyle = :solid,  color = :red,
  size = (400,300), legend=:bottomright, label = "prey");
plot!(time_evolution, sol0,vars=(0,2),
  linestyle = :solid,  color = :blue,
  size = (400,300), legend=:bottomright, label = "predator")
png(time_evolution, "lotka_volterra_prey_pradator.png")

Δt = 50
end_time = 50.
tspan = (0.,end_time)

senario = 400
J = 0.1:0.1:0.5
sol_ = []
for k in 1:senario
    u0 = [1.0rand()*cos(2π*rand())+1.,0.5rand()*sin(2π*rand())+0.5]
    prob = ODEProblem(lotka_volterra_prey_pradator!,u0,tspan,p)
    push!(sol_, solve(prob; saveat=(0.0:0.01:end_time)))
end

anim = @animate for t = 100:10:4950
    print('|')
    frame = plot();
    for j in 1:5
        u0 = [1.,J[j]]
        prob = ODEProblem(lotka_volterra_prey_pradator!,u0,tspan,p)
        frame = plot!(frame, solve(prob; saveat=(0.0:0.01:9)), vars=(1,2),
                        label = :none, color = colorant"#AAAAAA", linestyle = :solid)
    end
    for k in 1:senario
        frame = plot!(frame, sol_[k][t:(t+Δt)], vars=(1,2),
                    label = :none, color = :black, linealpha = (1:(1+Δt))/Δt)
    end
    frame = plot!(frame,
            size = (450,300), xlims = (0,3), ylims = (0,2), aspect_ratio=:equal,
            xlab = "prey", ylab = "predator")
end
gif(anim, "lotka_volterra_preypredator.gif")