📂力学系

力学系としての植物・送粉者モデル

モデル ¹

$N_{P}$種の植物^plantと$N_{A}$種の送粉者^pollinatorの間の関係を二部グラフ$\Gamma \left( N_{P}, N_{A} \right)$で表すとき、次のような結合力学系を通じて植物・送粉者の共生^mutualismをモデリングできる。 $$ \begin{align*} \dot{P}_{i} =& P_{i} \left( \alpha_{i}^{P} - \sum_{l=1}^{N_{P}} \beta_{il}^{P} P_{l} + {\frac{ \sum_{j=1}^{N_{A}} \gamma_{ij}^{P} A_{j}} { 1 + h \sum_{j=1}^{N_{A}} \gamma_{ij}^{P} A_{j} }} \right) \\ \dot{A}_{j} =& A_{j} \left( \alpha_{j}^{A} - \kappa_{j} - \sum_{l=1}^{N_{A}} \beta_{jl}^{A} A_{l} + {\frac{ \sum_{i=1}^{N_{P}} \gamma_{ij}^{A} P_{i}} { 1 + h \sum_{i=1}^{N_{P}} \gamma_{ij}^{A} P_{i}}} \right) \end{align*} $$

変数

• $P_{i}$ : $i = 1 , \ldots , N_{P}$番目の植物の豊富さ
• $A_{j}$ : $j = 1 , \ldots , N_{A}$番目の送粉者の豊富さ

パラメータ

• $\alpha_{k}$ : $k$番目の種の内部成長率
• $\beta_{kl}$ : $k$番目の種と$l$番目の種の間の競争
• $\gamma_{ij}$ : $i$番目の植物と$j$番目の送粉者の間の相互作用
• $\kappa_{j}$ : $j$番目の送粉者の減少率

説明

現実において植物がなす生態系はそれ自体が送粉者の生きる生活の場となり、送粉者は植物の繁殖に核心的な役割を果たす。植物・送粉者モデルはこのような相互作用を数学的にモデリングしたものである。

方程式を見ると一見複雑に見えるが、よく見ればロトカ・ヴォルテラ競争モデルで競争項の符号を反転させて協力項にした後、その反応をホリングタイプ1としたものにすぎない。同種内部、あるいは植物と送粉者それぞれのグループ内の競争はロジスティック成長に従う。

ここで方程式をなす係数はウェブ・オブ・ライフで提供される実際のデータに基づいている。上の図を見ると、データの種類によって植物と送粉者の間の相互作用が異なるという点、そして同時に植物同士の相互作用や送粉者の間の相互作用はなく、二部グラフの構造に従うという点がわかる。

R・ティッピング

気候変動に関連して大きな話題となっているのは、まさにミツバチなどの送粉者が消えつつあるということである。これに対し、送粉者を減らす役割をする$\kappa_{j}$が時間とともに一定に増加する現実を反映すると、このシステムは非自律的になりR・ティッピングを経験しうるようになる。 $$ \dot{\kappa_{j}} = \begin{cases} r & , \text{if } \kappa_{\text{min}} < \kappa_{j} < \kappa_{\text{max}} \\ 0 & , \text{otherwise} \end{cases} $$ 変化率^rate$r$が大きいということはそれだけ送粉者が速く消えるということを意味し、これに対して$R$・ティッピングを経験する確率を次のように定義できる。 $$ \Phi (r) := B \exp \left[ - C \left( \kappa_{\text{min}} - \kappa_{\text{max}} \right) / r \right] $$ ここで$C > 0$はある定数であり、$B > 0$は$A$が絶滅するベイスン、言い換えれば$A_{\infty} = 0$であるアトラクターのベイスンに関する値であるが、$\kappa = \kappa_{\text{min}}$のときと$\kappa = \kappa_{\text{max}}$のときの面積の差として求められ、元の高次元システムの場合は数値的に決定される。

単純化

高次元ではなく2次元の効率的なモデルとしては次のように表せる。

$$ \begin{align*} \dot{P} =& P \left( \alpha - \beta P + {\frac{ \gamma^{P} A} { 1 + h \gamma^{P} A }} \right) \\ \dot{A} =& A \left( \alpha - \kappa - \beta A + {\frac{ \gamma^{A} P} { 1 + h \gamma^{A} P}} \right) \end{align*} $$

R・ティッピングで$\Phi (r)$を説明する言葉が少し難しかったが、核心的な内容を単純化して理解してみよう。ネットワーク表現で水色(1時方向)のデータはM_PL_036というデータによって構築されており、この場合$\kappa = \kappa_{\text{min}}$のときと$\kappa = \kappa_{\text{max}}$のベイスンはそれぞれ上の図の$C$と$D$に対応する。図で灰色の部分は送粉者$A$が絶滅せず$P$と共存する固定点に対するベイスン、赤い部分は送粉者$A$が絶滅する固定点に対するベイスンである。

さて再び$\Phi (r)$を考えてみるとそれほど難しくないのだが、$\kappa_{\text{min}}$から$\kappa_{\text{max}}$へ変化する間に絶滅に関するベイスンの面積がどのように変わるかが変わるのである。この差が大きいということはそれだけ絶滅の可能性が高いという意味になる。

コード

次は植物・送粉者モデルをシミュレーションできるJuliaコードである。

using DataFrames, DifferentialEquations, OrdinaryDiffEqLowOrderRK

function factory_pollinator(κ::Number; ic = [1.,1.], saveat = 0:1e-2:10)
      α, β,  _h,   γp,   γA = (
    0.3, 1, 0.4, 1.93, 1.77 )
    function sys(du, u, p, t)
        P,A = u
        κ = p[1]

        du[1] = P*(α - β*P + frac(γp*A, 1 + _h*γp*A))
        du[2] = A*(α - κ - β*A + frac(γA*P, 1 + _h*γA*P))
        return du
    end
    sol = solve(ODEProblem(sys, ic, (0, last(saveat)), [κ]), RK4(), dt = saveat.step.hi, adaptive=false, maxiters = 1e+7)
    matrix = Matrix([sol.t'; sol[:, :]; stack([sys(zeros(2), u, [κ], 0) for u in sol.u])]')
    return matrix[sol.t .≥ first(saveat), :][1:end-1, :]
end
factory_pollinator(T::Type, args...; kargs...) =
DataFrame(factory_pollinator(args...; kargs...), ["t", "P", "A", "dP", "dA"])