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

Question

我正在尝试使用 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,rep(0,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))

非常感谢任何指导，提前致谢！

Answer 1

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

消亡：如果年龄组每时间步移动并且第一个元素刚好填充为零，则所有内容都会在 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, rep(0, n_agegroups-1))))
  })
}

shift <- function(t, state, 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, type="l", ylab="S")

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

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

Realistic age structured model using ODE from the deSolve package

r

ode

differential-equations

infection

desolve