在集合种群模型中纳入时间步长相关速率/参数

Question

我想使用 B. Raynor 的代码将随时间变化的速率合并到元种群模型中。该代码最初发布于 https://rpubs.com/bhraynor/MetapopulationModel

在下面的代码中，我希望 mu 取向量中的值

mu = seq(from = 0, to = 1, length = 100)

下面是运行良好的代码，但如果我取消注释

mu = seq(from = 0, to = 1, length = 100)

向量并注释向量

mu = c(1/1000, 1/1000, 1/1000)

，则不会运行。

我收到错误消息： eval(substitute(expr), data, enclos = Parent.frame()) 中的错误： dims [产品 3] 与对象 [100] 的长度不匹配

#clean environment
rm(list = ls()) 

#load libraries
library(dplyr)
library(ggplot2)
library(deSolve)

#initial state conditions
rfun.FormatInit <- function(init.df){
  
  #create flexible naming structure for vectors
  l <- nrow(init.df)
  state.names <- as.vector(colnames(init.df))
  varnames <- NULL
  for(i in 1:length(state.names)){
    varnames <- cbind(varnames, paste(state.names[i], seq(1:l), sep=""))
  }
  as.vector(varnames)
  
  #vectorize matrix input
  init.vector <- as.vector(as.matrix(init.df))
  
  
  names(init.vector) = varnames
  return(init.vector)
}

#parameters
rfun.FormatParams <- function(parm.df){
  
  #create flexible naming structure for vectors
  l <- nrow(parm.df)
  parm.names <- as.vector(colnames(parm.df))
  varnames <- NULL
  for(i in 1:length(parm.names)){
    varnames <- cbind(varnames, paste(parm.names[i], seq(1:l), sep=""))
  }
  as.vector(varnames)
  
  #vectorize matrix input
  parm.vector <- as.vector(as.matrix(parm.df))
  
  
  names(parm.vector) = varnames
  return(parm.vector)
}

#Contact matrix
rfun.FormatContact = function(df.contact, beta.contact_constant){
  l = length(df.contact)
  matrix.contact  = as.matrix(df.contact, ncol=4)
  return(beta.contact_constant*matrix.contact)
}



MODEL <- function(time, state, parameters, beta.contact) {
  with(as.list(c(state, parameters)), {
    
    #number of cells
    l1 = length(state)
    l2 = length(unique(substr(names(state),1,1)))
    l = l1/l2
    
    #Define initial conditions
    S = matrix(state[1:l], ncol=1)
    I = matrix(state[(l+1):(2*l)], ncol=1)
    R = matrix(state[(2*l+1):(3*l)], ncol=1)
    
    #Define params
    beta <- matrix(parameters[paste0("beta", 1:l)], ncol=1)
    gamma <- matrix(parameters[paste0("gamma", 1:l)], ncol=1)
    alpha <- matrix(parameters[paste0("alpha", 1:l)], ncol=1)
    mu <- matrix(parameters[paste0("mu", 1:l)], ncol=1)
    
    #mu = seq(from = 0.001, to = 0.0015, length = 100)
    
    #System of eqns
    #bSI = (beta) %*% t(S) %*% I
    bSI = beta*S*I
    bSI.contact = beta.contact%*%I*S
    birth = mu*(S+I+R) + alpha*I
    
    
    
    dS <- birth -bSI - bSI.contact-mu*S
    dI <-  bSI + bSI.contact - gamma*I -alpha*I - mu*I
    dR <-  gamma*I - mu*R
    
    #Output
    #return(list(c(dS, dI, dR),mu))
    
    output <- c(dS, dI, dR)
    list(output, Y = mu)
  })
}




#set initial conditions
init <- data.frame(S=c(100,100,100),
                   I= c(10,10,10),
                   R= c(0,0,0)) %>%
  rfun.FormatInit()

#parameters


parms <- data.frame(beta  = c(0.01, 0.05, 0.001),
                    gamma = c(1/5, 1/4, 1/3),
                    mu = c(1/1000, 1/1000, 1/1000),
                    alpha = c(0.2, 0.3, 0.4))%>%
  rfun.FormatParams()



#contact matrix
beta.contact_constant = 0.001

Time = 100
dt = 1 #Step size dt
times <- seq(0, Time, by = dt)



#set up contact matrix for example
df.contact <- data.frame(q1 = c(0,1,1),
                         q2 = c(1,0,1),
                         q3 = c(1,1,0))
beta.contact = rfun.FormatContact(df.contact, beta.contact_constant)

#Run simulation
out <- ode(y=init, times=times, func=MODEL, parms=parms, beta.contact=beta.contact)

Answer 1

参数通常是固定值。与时间相关的值可以通过多种方法实现：

1. 将依赖项集成到模型函数中，例如通过为参数添加一个额外的微分方程，
2. 使用强制函数并从外部“信号”进行插值，
3. 使用事件，修改状态变量。此方法可以与方法 #1 结合使用。

方法 #2 和 #3 在 deSolve 包的

?forcings

和

?events

帮助页面以及 在线教程 页面中进行了描述。

要获得更详细的答案，请避免发布长代码并将您的问题简化为最小可重现示例。