r - 在 R desolve 中设置导数的下限

Question

我正在使用 R 的包 deSolve 创建一个数学模型。我有三个导数：dx（受感染的宿主的比例）、dY（受感染的载体的比例）和 dm（载体与宿主的比例）。我的模型的重点是显示某种杀虫剂处理对种群的影响（效果表示为参数“z”。为了将这个时间相关的协变量纳入模型，使用了 approxfun 函数。模型工作正常, 但是我想为 dm 设置一个下限(假设不是人口中的所有向量都会受到影响). 如果不设置下限, 我的代码和图表如下所示：

Initial vectors for days post treatment, % killed
x <- c(4, 30, 60, 90, 120, 210, 360)
z <- c(1.0, 0.99, 0.99, 0.79, 0.7, 0.02, 0) 

plot(z ~ x)

#=============================================
#  fit data with logistic curve
#       -extract fit values using equation y = Asym / (1 + exp((xmid - input) / scal))
#=============================================
fit2 <- nls(z ~ SSlogis(x, Asym, xmid, scal), data = data.frame(x, z))
summary(fit2)

lines(seq(0, 400, length.out = 400),
      predict(fit2, newdata = data.frame(x = seq(0.5, 400, length.out = 400))))

Asym<-summary(fit2)$parameters[1,1]
xmid<-summary(fit2)$parameters[2,1]
scal<-summary(fit2)$parameters[3,1]

times <- seq(0, 1000, by = 1)
signal <- data.frame(times = times, import = rep(0, length(times)))
signal$import=  Asym / (1 + exp((xmid - times) / scal))

#Force time dependent covariate into the model
input <- approxfun(signal, rule = 2)


RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y)) 
      dm <- ((R*(1-m/K)*m )+(-m*a*z))
      return(list(c(dX, dY, dm)))
    }
  )
}


initTx2 <- c(X = 0.01, Y= 0, m=40) 
parametersTx2 <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09)
outTx2 <- as.data.frame(ode(y = initTx2, times = times, func = RMTx2, parms = parametersTx2))
RESULTS2<-data.frame(outTx2$X,outTx2$Y)
RESULTS2m <-data.frame(outTx2$m, outTx2$Y*outTx2$Y)
timesTx2 <- seq(0, 10000, by = 1)

这是情节。我想做的是限制 dm 随着时间的推移而下降，以便在治疗期间它不会下降到某个值以下

我正在尝试设置参数 dm，例如，该值不能低于 15。我尝试了一些代码来执行此操作，包括：

RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y))
      dm <- if (isTRUE (((R*(1-m/K)*m )+(-m*a*z)) > MM)) ((R*(1-m/K)*m )+(-m*a*z)) else 15 
      return(list(c(dX, dY, dm)))
    }
  )
}

initTx2 <- c(X = 0.01, Y= 0, m=40) 
parametersTx2 <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM= 15)
outTx2 <- as.data.frame(ode(y = initTx2, times = times, func = RMTx2, parms = parametersTx2))
RESULTS2<-data.frame(outTx2$X,outTx2$Y)
RESULTS2m <-data.frame(outTx2$m, outTx2$Y*outTx2$Y)

不幸的是，由于某种原因，这导致人口无限期地增加：

这种方法是否有一些根本不起作用的东西？或者这更像是一个编码错误？谢谢！

标签： rodederivative

Answer 1

您可以为此尝试使用 Monod 类型的保护措施。MM是你的阈值和kmMonod 常数。它的值是相当不重要的，但它不应该太小以避免求解器的数值困难。

RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y)) 

      #dm <- (R*(1-m/K)*m )+(-m*a*z)
      dm <- (R*(1-m/K)*m )+(-m*a*z) * (m - MM) / (km + (m - MM))  

      list(c(dX, dY, dm), dm=dm) ## return derivative dm as additional output
    }
  )
}


initTx2 <- c(X = 0.01, Y= 0, m=40) 

parameters_old <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM = 0, km=0)
parameters_new <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM = 15, km=0.1)

out_old <- ode(y = initTx2, times = times, func = RMTx2, parms = parameters_old)
out_new <- ode(y = initTx2, times = times, func = RMTx2, parms = parameters_new)

plot(out_old, out_new)