需要帮助在 R 中编码具有最小值的 ODE

Question

我正在为我的论文研究浮游植物资源模型。为了熟悉该模型，我开始研究一种物种（浮游植物）和一种资源（例如磷）的常微分方程。我为此运行了解析解和数值模拟。

解析解：

q <- function(inputs) {
  g <- inputs[3]
  r <- inputs[2]
  d <- inputs[1]
  return((g*r)/(r-d))
}

s <- function(inputs) {
  g <- inputs[3]
  r <- inputs[2]
  d <- inputs[1]
  k <- inputs[5]
  v <- inputs[4]
  return(-(d*g*k*r)/(d*g*r+d*v-r*v))
}

n <- function(inputs) {
  g <- inputs[3]
  r <- inputs[2]
  d <- inputs[1]
  k <- inputs[5]
  v <- inputs[4]
  s0 <- inputs[6]
  (d - r)*(-1*(s0/(g*r)) - (d*k)/(d*g*r + d*v - r*v))
}

数值模拟：

library(deSolve)
library(tidyverse)
one <- function (t, x, params) {
  ## first extract the state variables
  N <- x[1]
  Q <- x[2]
  S <- x[3]
  ## now extract the parameters
  r <- params["r"]
  g <- params["g"]
  d <- params["d"]
  v <- params["v"]
  k <- params["k"]
  s0 <- params["s0"]
  
  
  ## now code the model equations
  dNdt <- N*(r*(1-(g/Q))-d)
  dQdt <- v*(S/(S+k))-r*(Q-g)
  dSdt <- d*(s0-S)-(N*v*(S/(S+k)))
  ## combine results into a single vector
  dxdt <- c(dNdt,dQdt,dSdt)
  ## return result as a list!
  list(dxdt)
}
parms <- c(r=0.1, g=0.001, d=0.03, v=0.1, k=0.01, s0=1) 
times <- seq(from=0,to=1000,by=0.1)
xstart <- c(N=800,Q=1,S=1)

ode(
  func=one,
  y=xstart,
  times=times,
  parms=parms,
) %>%
  as.data.frame() -> out

out %>%
  gather(variable,value,-time) %>%
  ggplot(aes(x=time,y=value,color=variable))+
  geom_line(size=2)+
  theme_classic()+ 
  labs(x='time',y='?')


### setting up dataframe with numerical simulation

## List of parameter values
l <- list(r = 0.1, 
          g = 0.001,
          d = c(0.03, 0.035, 0.04, 0.045, 0.05, 0.055, 0.06, 0.065, 0.07, 0.075,
                0.08, 0.085, 0.09),
          v = 0.1,
          k = 0.01,
          s0 = 1)

## Matrix with 'prod(lengths(l))' rows listing combinations
## of parameter values
P <- as.matrix(do.call(expand.grid, l))

## Wrapper around 'ode' extracting state at final time
getLast <- function(parms, y, times, func, ...) {
  X <- deSolve::ode(y, times, func, parms, ...)
  X[nrow(X), -1L]
}

## Matrix of final states corresponding rowwise to 'P'
X <- t(apply(P, MARGIN = 1L, FUN = getLast, 
             y = c(N = 800, Q = 1, S = 1), 
             times = seq(from = 0, to = 1000, by = 1),
             func = one))

cbind(P, X) |> as.data.frame() -> dilution
## setting up dataframe of analytical solutions

analytical <- data.frame(d = seq(0.01, 0.09, 0.005), r = rep(0.1,17), g = rep(0.001,17),
                         v = rep(0.1, 17), k = rep(0.01, 17), s0 = rep(1, 17))

analytical$qa <- apply(analytical, 1, q)
analytical$na <- apply(analytical, 1, n)
analytical$sa <- apply(analytical, 1, s)

## merge analytical and numerical simulation together

total <- merge(analytical, dilution, by = c("r", "d", "g", "k", "v", "s0"))

然后我将它们合并并绘制以比较数值模拟和解析解（以确保它们相等）。现在，我正在尝试添加另一种资源来模拟一个物种、两种资源，因此其中一种资源将受到限制。这个方程来自一篇论文：方程，确定哪个资源受到限制的方程是这里：限制资源方程

我的顾问已经在 Mathematica 中进行了模拟并求解了解析解，我正在尝试将其转换为 R，但我对 R 仍然有点陌生。我试图让它看起来像这样： mathematica 笔记本方程，但在 R 中。我认为 Q1、Q2、S1 和 S2 的方程非常简单，类似于我已经完成的一个物种一个资源，但我是努力将最小值纳入 N 方程中。另外，这会改变我运行数值模拟的方式吗？我的顾问说数值模拟会更容易编码，所以我想知道是什么使得该分析解决方案更难实现？我感谢任何和所有的帮助！非常感谢您抽出时间。

Answer 1

我建议在实现完整模型之前首先从一些“玩具示例”开始。这里的模型具有两个物种（X1，X2）和两个限制底物（S1，S2）。底物依赖性被表述为 Monod 函数（相当于 Holling II）。另请注意底物和群体之间的单位转换因子y。

library(deSolve)

model <- function (time, y, parms) {
  with(as.list(c(y, parms)), {
    f <- function(S, k) S/(k + S) # functional response

    grow1 <- r1 * min(f(S1, k11), f(S2, k12)) * X1
    grow2 <- r2 * min(f(S1, k21), f(S2, k22)) * X2

    dS1 <- -1/y11 * grow1 -1/y12 * grow2
    dS2 <- -1/y21 * grow1 -1/y22 * grow2
    dX1 <- grow1
    dX2 <- grow2

    list(c(dS1, dS2, dX1, dX2), 
         lim1 = f(S1, k11) < f(S2, k12), # indicates which resource is limiting
         lim2 = f(S1, k21) < f(S2, k22))
  })
}

y0     <- c(S1 = 1, S2 = 1,         # two substrates
            X1 = 1, X2 = 1)         # two populations
parms  <- c(r1 = 0.1, r2=0.12,      # maximum growth rates
            k11 = 0.12, k12 = 0.08, # half saturation constants
            k21 = 0.08, k22 = 0.12,
            y11 = 5, y12 = 1,       # yield ratios (=unit conversion factors)
            y21 = 5, y22 = 1)

times  <- seq(0, 100)
out <- ode(y0, times, model, parms)
plot(out)

模拟显示了两种物种的限制如何在资源之间切换：

对于更多物种和底物，应该考虑矢量化。一种非常紧凑的方法是 Gujer-Peterson 矩阵方法，例如在rodeo包中可用。