📂力学系

動力学における食物連鎖システム

モデル ¹

$$ \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$ 間に ホーリング型IIの反応関数 を加えることで カオス的 なシステムになる. 実際に $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}$ のとき初期条件が変化するにつれて過渡寿命の頻度が 指数分布 に従うことを確認できる図だ.

以下はこれを再現するための Julia コードだ.

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