📂力学系

ヘイスティングス・パウエル系

モデル ¹

$$ \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}$: ホーリング型IIの機能応答に従う飽和^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]$を分岐パラメータとして分岐図を描いたものである².

コード

以下は上の分岐図を再現する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 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])