设定R解旋中导数的下限

Question

我正在使用R的软件包deSolve创建数学模型。我有三个派生词：dx（感染宿主的比例），dY（感染载体的比例）和dm（这是总体中载体与宿主的比例）。我的模型的重点是显示某种杀虫剂处理对种群的影响（效果表示为参数“ z”。为了将此时间相关协变量纳入模型中，使用了近似函数功能。该模型正常工作，但是我想为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)

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

这种方法是否存在根本无法使用的基本原则？还是更多的编码错误？谢谢！

Answer 1

您可以为此尝试Monod型防护装置。 MM是您的阈值，km是Monod常数。它的值不是很严格，但是它的值不能太小，以避免对求解器造成数值困难。

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)