📂동역학

헤이스팅-파웰 시스템

모델 ¹

$$ \begin{align*} \frac{dX}{dT} =& R_{0} X \left( 1 - \frac{X}{K_{0}} \right) - C_{1} \frac{A_{1} X}{B_{1} + X} Y \\ \frac{dY}{dT} =& \frac{A_{1} X}{B_{1} + X} Y - \frac{A_{2} Y}{B_{2} + Y} Z - D_{1} Y \\ \frac{dZ}{dT} =& C_{2} \frac{A_{2} Y}{B_{2} + Y} Z - D_{2} Z \end{align*} $$

변수

• $X(T)$: $T$ 시점에서 $Y$ 에 먹히는 종의 개체수다.
• $Y(T)$: $T$ 시점에서 $X$ 를 먹고 $Z$ 에 먹히는 종의 개체수다.
• $Z(T)$: $T$ 시점에서 $Y$ 를 먹는 종의 개체수다.

파라미터

• $R_{0}$: $X$ 의 내재적 성장률^{intrinsic growth rate}이다.
• $K_{0}$: $X$ 의 환경 용량^{carrying capacity}이다.
• $C_{1}^{-1}, C_{2}$: 각각 $Y$ 와 $Z$ 에 대한 전환률^{conversion rate}이다.
• $A_{1}, A_{2}, B_{1}, B_{2}$: 홀링타입 2 함수형 반응에 따른 포화^saturation에 관계된 계수들이다.

설명

먹이사슬 시스템: $$ \begin{align*} \dot{R} =& R \left( 1 - {\frac{ R }{ K }} \right) - {\frac{ x_{c} y_{c} C R }{ R + R_{0} }} \\ \dot{C} =& x_{c} C \left( {\frac{ y_{c} R }{ R + R_{0} }} - 1 \right) - {\frac{ x_{p} y_{p} P C }{ C + C_{0} }} \\ \dot{P} =& x_{p} P \left( {\frac{ y_{p} C }{ C + C_{0} }} - 1 \right) \end{align*} $$

이 포스트가 작성된 이유는 오로지 하나, 먹이사슬 시스템과 헤이스팅-파웰 시스템 간에 본질적인 차이가 없다는 걸 일러주기 위해서다. 어떤 논문에서는 이 둘을 구분하기도 하고 어떤 논문에서는 별도의 명시 없이 둘을 같은 것으로 취급하기도 하는데, 따지자면 후자에 가까운 것이다.

물론 방정식이 완전히 같냐면 그렇지는 않고, 일부 파라미터의 자유도가 다르기는 하지만 결국 둘 다 전문화된^specialized 포식자가 있음을 상정하는 생태계 모델이라는 것은 마찬가지다.

$$ \begin{align*} \dot{v_{1}} =& a_{0} v_{1} - b_{0} v_{1}^2 - \frac{w_{0} v_{1} v_{2}}{d_{0} + v_{1}} \\ \dot{v_{2}} =& \frac{w_{1} v_{1} v_{2}}{d_{1} + v_{1}} - \frac{w_{2} v_{2} v_{3}}{d_{2} + v_{2}} - a_{1} v_{2} \\ \dot{v_{3}} =& \frac{w_{3} v_{2} v_{3}}{d_{3} + v_{2}} - c_{3} v_{3} \end{align*} $$ 다음은 $X = v_{1}$, $Y = v_{2}$, $Z = v_{3}$ 로 나타내고 계수를 $b_{0} = 0.05$, $w_{0} = 1$, $d_{0} = 10$, $w_{1} = 2$, $d_{1} = 10$, $w_{2} = 1.5$, $d_{2} = 10$, $a_{1} = 1$, $w_{3} = 2$, $d_{3} = 20$, $c_{3} = 0.7$ 로 두고 $a_{0} \in [1.5, 2.0]$ 를 바이퍼케이션 파라미터로 두고 바이퍼케이션 다이어그램을 그려낸 것이다².

코드

다음은 위의 바이퍼케이션 다이어그램을 재현하는 줄리아 코드다.

function RK4(f::Function, v::AbstractVector, h=1e-2)
    V1 = f(v)
    V2 = f(v + (h/2)*V1)
    V3 = f(v + (h/2)*V2)
    V4 = f(v + h*V3)
    return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end

function arglmax(x)
    bits = circshift(x, 1) .< x .> circshift(x, -1)
    bits[1] = false
    bits[end] = false
    return findall(bits)
end

frac(x, y) = x ./ y
function factory_hastingpowell(a0::Number; ic = [0., 1, 1, 1], tspan = 0:1e-2:10)
      b0, w0, d0, w1, d1,  w2, d2, a1, w3, d3,  c3 = (
    0.05,  1, 10,  2, 10, 1.5, 10,  1,  2, 20, 0.7 )
    function sys(v::AbstractVector)
        _,v1,v2,v3 = v

        dt = 1
        dv1 = a0*v1 - b0*v1^2 - w0*frac(v1*v2, d0 + v1)
        dv2 = w1*frac(v1*v2, d1 + v1) - w2*frac(v2*v3, d2 + v2) - a1*v2
        dv3 = w3*frac(v2*v3, d3 + v2) - c3*v3
        return [dt, dv1, dv2, dv3]
    end
        
    len_t_ = length(tspan); h = tspan.step.hi
    v = ic; DIM = length(v)
    traj = zeros(2+len_t_, 2DIM); traj[1, 1:DIM] = v

    tk = 0
    while tk ≤ len_t_
        t,_,_,_ = v
        v, dv = RK4(sys, v, h)

        if t ≥ first(tspan)
            tk += 1
            traj[tk+1,         1:DIM ] =  v
            traj[tk  , DIM .+ (1:DIM)] = dv
        end
    end
    return traj[1:(end-2), :]
end
factory_hastingpowell(T::Type, args...; kargs...) =
DataFrame(factory_hastingpowell(args...; kargs...), ["t", "v1", "v2", "v3", "dt", "dv1", "dv2", "dv3"])[:, Not(:dt)]

hrzn = []
vrcl = []
for a0 = 0.15:0.001:2.0
    data = factory_hastingpowell(DataFrame, a0, tspan = 0:1e-2:1000)[50001:end, :]
    bits = arglmax(data.v1)
    push!(hrzn, [a0 for _ in bits])
    push!(vrcl, data.v1[bits])
end
scatter([hrzn...;], [vrcl...;], ms = 1, msw = 0, color = :black, alpha = 0.5, legend = :none, xlims = [1.5, 2.0], ylims = [25, 40])