📂동역학

동역학에서의 먹이사슬 시스템

모델 ¹

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

변수

• $R(t)$: $t$ 시점에서 자원^resource의 밀도다.
• $C(t)$: $t$ 시점에서 소비자^consumer의 밀도다.
• $P(t)$: $t$ 시점에서 포식자^predator의 밀도다.

파라미터

• $K = 0.94$: 자원의 환경 용량^{carrying capacity}이다.
• $x_{c} = 0.4$, $y_{c} = 1.7$, $R_{0} = 0.16129$: 소비자가 자원을 포식하는 것에 관계된 값들이다.
• $y_{p} = 5.0$, $x_{p} = 0.08$, $C_{0} = 0.5$: 포식자가 소비자를 포식하는 것에 관계된 값들이다.

설명

소개된 $3$차원 먹이사슬 시스템^{food-chain system}은 자원이 로지스틱 성장모델을 따르며 자원-소비자와 소비자-포식자라는 두 가지 롯카-볼테라 포식자-피식자 모델이 커플링된 다이내믹 시스템이라 볼 수 있다.

원래 이렇게 간단한 시스템을 조합한 유형의 시스템은 복잡할 수가 없지만, $R-C-P$ 사이에 홀링 타입 2의 반응 함수를 가미함으로써 캐어릭한 시스템이 된다. 실제로 $K, y_{c}, y_{p}$ 를 바꿔가며 그린 바이퍼케이션 다이어그램은 다음과 같다.

연구적인 측면에서 이런 시스템을 알아두는 것은 논문의 가치를 더하는 방법이 될 수도 있겠다. 똑같이 $3$차원이고 캐어릭한 시스템의 예로는 로렌츠 어트랙터가 있지만 너무 유명하고 많이 쓰인 나머지 어떤 벤치마크의 대상으로써 매력적이지 않으며, 현재에 이르러서는 $2$차 이내의 다항함수로 표현된다는 점에서 그 복잡성도 떨어지는 것으로 비칠 수 있다. 반면 먹이사슬 시스템은 그 자체로 이미 생태학^ecology의 직관적인 의미를 담고 있으며 캐어릭하면서 수식도 유리함수를 포함하고 있어 상대적으로 복잡하다².

과도적 혼돈 ^{3 4}

한편 $x_{c} = 0.4$, $y_{c} = 2.009$, $x_{p} = 0.08$, $y_{p} = 2.876$, $R_{0} = 0.16129$, $C_{0} = 0.5$ 일 때 $K \approx 0.99976$ 근방에서는 과도적 혼돈이 발생한다. 이 바이퍼케이션 다이어그램을 그리기 위해서는 과도적 혼돈을 고려하여 여러 초기조건에 대한 시뮬레이션이 반복되는 것에 주의해야 한다.

다음은 $K = 0.99976 + 2\times 10^{-4}$ 일 때 초기 조건이 바뀜에 따라 과도적 수명의 빈도가 지수분포를 따른다는 걸 확인할 수 있는 도표다.

다음은 이를 재현하기 위한 줄리아 코드다.

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 factory_foodchain2(K::Number; ic = [0., 0.4rand() + 0.6, 0.4rand() + 0.15, 0.5rand() + 0.3], tspan = 0:1e-2:10)
     xc,    yc,   xp,    yp,      R0,  C0 = (
    0.4, 2.009, 0.08, 2.876, 0.16129, 0.5)
    function sys(v::AbstractVector)
        _,R,C,P = v

        dt = 1
        dR = R*(1 - frac(R,K)) - xc*yc*frac(C*R, R + R0)
        dC = xc*C*(frac(yc*R, R + R0) - 1) - xp*yp*frac(P*C, C + C0)
        dP = xp*P*(frac(yp*C, C + C0) - 1)
        return [dt, dR, dC, dP]
    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_foodchain2(T::Type, args...; kargs...) =
DataFrame(factory_foodchain2(args...; kargs...), ["t", "R", "C", "P", "dt", "dR", "dC", "dP"])[:, Not(:dt)]

hrzn = []
vrtl = []
@time for k in 0.9:0.0001:1
    for _ in 1:500
        traj = factory_foodchain2(DataFrame, k, tspan = 0:1e-1:5000)[40000:end,:]
        bits = arglmin(traj.P)
        if maximum(traj.P) > 0.55
            push!(hrzn, fill(k, length(bits)))
            push!(vrtl, traj.P[bits])
            break
        else
            println("k: $k should be tried again")
        end            
    end
end
using LaTeXStrings
scatter([hrzn...;], [vrtl...;], ylims = [0.55, 0.8], yticks = 0.55:0.05:0.8, xticks = 0.9:0.02:1,
    ms = 1, msw = 0, color = :black, alpha = 0.5, legend = :none, xlabel = L"K", ylabel = L"\min P")
png("bifurcaiton.png")

tau_ = zeros(Int64, 100000)
@showprogress @threads for k in eachindex(tau_)
    traj = factory_foodchain2(DataFrame, 0.99976 + 2e-4, ic = rand(4), tspan = 0:1e-0:5000)
    tau = findlast(diff(cumsum(traj.P .< 0.55)) .== 0)
    tau_[k] = !isnothing(tau) ? tau : -1
end
freq = [count(x .≤ tau_ .< x + 250) for x in 250:250:4750]
scatter(250:250:4750, freq ./ sum(freq), yscale = :log10, xticks = 0:1000:5000, xlims = [0, 5000], yticks = [1e-1, 1e-2], color = :white, legend = :none, xlabel = "transient lifetime", ylabel = "frequency")