使用 deSolve 包中的 ODE 的真实年龄结构模型

如何解决使用 deSolve 包中的 ODE 的真实年龄结构模型

我正在尝试使用 deSolve 包中的 ODE 模拟一个真实的年龄结构模型，其中所有个体都可以在时间步长结束时转变为以下年龄组（而不是以给定的速率连续老化）。

以一个模型为例，它具有易感 (S) 和传染 (I) 两种状态，每个状态被分为 4 个年龄组（S1、S2、S3、S4 和 I1、I2、I3、I4），所有个体S1中的应该在时间步长结束时进入S2，S2中的应该进入S3，以此类推。

我尝试分两步进行，第一步是求解 ODE，第二步是在时间步长结束时将个体转移到以下年龄组，但没有成功。

以下是我的尝试之一：

library(deSolve)

times <- seq(from = 0,to = 100,by = 1) 
n_agecat <- 4

#Initial number of individuals in each state
S_0 = c(999,rep(0,n_agecat-1))
I_0 = c(1,n_agecat-1))

si_initial_state_values <- c(S = S_0,I = I_0)
# Parameter values 
si_parameters <- c(beta = 0.01) #contact rate assuming random mixing


si_model <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
    n_agegroups <- 4
    
    S <- state[1:n_agegroups]
    I <- state[(n_agegroups+1):(2*n_agegroups)]
    
    # Total population
    N <- S+I
    
    # Force of infection
    lambda <- beta * I/N
    
    # Solving the differential equations 
    dS <- -lambda * S 
    dI <- lambda * S 

    # Trying to shift all individuals into the following age group
    S <- c(0,S[-n_agecat]) 
    I <- c(0,I[-n_agecat]) 
    
    return(list(c(dS,dI)))
  })
}

output <- as.data.frame(ode(y = si_initial_state_values,times = times,func = si_model,parms = si_parameters))

任何指导将不胜感激，在此先感谢您！

解决方法

我看过你的模型。在一个事件函数中实现 shift 原理上是可行的，但是主模型仍然存在几个问题：

消亡：如果年龄组每时间步移动并且第一个元素刚填充为零，则所有内容都将在 4 个时间步内移到末尾，人口就会消亡。立>
感染：在您的情况下，感染者只能感染同一年龄组，因此您需要在计算 lambda 之前对“年龄”组进行汇总。
最后，什么是“年龄”组？你想要感染后的时间吗？

总而言之，有几种选择：我个人更喜欢离散模型进行此类模拟，即差分方程、年龄结构矩阵模型或基于个体的模型。

如果你想保持它是一个 ODE，我建议让易受感染的人一起成为一个状态，并且只将受感染的人作为阶段结构化。

这里有一个简单的例子，请检查：

library(deSolve)

times <- seq(from = 0,to = 100,by = 1) 
n_agegroups <- 14
n_agecat    <- 14

# Initial number of individuals in each state
S_0 = c(999)                     # only one state
I_0 = c(1,rep(0,n_agecat-1))    # several stages

si_initial_state_values <- c(S = S_0,I = I_0)
# Parameter values 
si_parameters <- c(beta = 0.1) # set contact parameter to a higher value


si_model <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
    
    S <- state[1]
    I <- state[2:(n_agegroups + 1)]
    
    # Total population
    N <- S + sum(I)
    
    # Force of infection
    #lambda <- beta * I/N         # old
    lambda <- beta * sum(I) / N   # NEW
    
    # Solving the differential equations 
    dS <- -lambda * S 
    dI <-  lambda * S 
     
    list(c(dS,c(dI,n_agegroups-1))))
  })
}

shift <- function(t,p) {
  S <- state[1]
  I <- state[2:(n_agegroups + 1)]
  I <- c(0,I[-n_agecat]) 
  c(S,I)
  
}

# output time steps (note: ode uses automatic simulation steps!)
times     <- 1:200

# time step of events (i.e. shifting),not necessarily same as times
evt_times <- 1:200   

output <- ode(y = si_initial_state_values,times = times,func = si_model,parms = si_parameters,events=list(func=shift,time=evt_times))

## default plot function
plot(output,ask=FALSE)

## plot totals
S <- output[,2]
I <- rowSums(output[,-(1:2)])

par(mfrow=c(1,2))
plot(times,S,type="l",ylim=c(0,max(S)))
lines(times,I,col="red",lwd=1)

## plot stage groups
matplot(times,output[,-(1:2)],col=rainbow(n=14),lty=1,ylab="S")

注意：这只是一个技术演示，不是一个有效的阶段结构SIR模型！