📂동역학

동역학계로써의 식물-수분자 모델

모델 ¹

$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}}$ 로 변화하는 동안 멸종에 관한 베이신의 면적이 어떻게 달라지는가가 달라지는 것이다. 이 차이가 크다는 것은 그만큼 멸종의 가능성이 높다는 의미가 된다.

코드

다음은 식물-수분자 모델을 시뮬레이션 할 수 있는 줄리아 코드다.

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