📂동역학

롯카-볼테라 포식자-피식자 모델

개요

롯카-볼테라 포식자-피식자 모델은 종간의 상호작용을 시스템으로써 모델링하며, 특히 포식자-피식자 모델은 두 종의 포식관계를 나타낸다. 두 종에 대해서만 다루면 그 확장은 끝이 없기 때문에 먹이사슬을 표현하기엔 충분하다.

모델¹

$$ \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*} $$

$x$ 에 대한 $y$ 의 포식행위는 $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) $$

이 고정점에 대한 분석이 상당히 흥미로운데, 포식자와 피식자를 구분하지 않고 공평하게 죽일 경우 피식자의 개체수는 오히려 늘어난다는 것이다. 이를 볼테라 원리^{volterra’s Principle}라 하는데, 동역학적인 분석이 간단한 것에 비해 외부 개입에 대해 명확한 인사이트를 주고 있음을 확인할 수 있다.

코드

다음은 본문의 이미지들을 만드는 줄리아 코드다.

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")