Realistic age structured model using ODE from the deSolve package

Question

I am trying to simulate a realistic age structured model where all individuals could shift into the following age group at the end of the time step (and not age continuously at a given rate) using ODE from the deSolve package.

Considering for example a model with two states Susceptible (S) and Infectious (I), each state being divided in 4 age groups (S1, S2, S3, S4, and I1, I2, I3, I4), all individuals in S1 should go into S2 at the end of the time step, those in S2 should go into S3, and so on.

I tried to make this in two steps, the first by solving the ODE, the second by shifting individuals into the following age group at the end of the time step, but without success.

Below is one of my attempts :

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))

Any guidance will be much appreciated, thank you in advance!

A shift of states within the model function is not possible, as the model function describes the derivativea, not the state. There are two methods to get this work: (1) using a "stop and go" approach, i.e. run ode for a given time, then shift states, then run ode for the next time period and so on, or (2) using the event mechanism. see `?events` in the deSolve help. — tpetzoldt, Jul 20 '21 at 19:55

score 1 · Accepted Answer · answered Jul 21 '21 at 16:32

I had a look at your model. Implementing the shift in an event function works, in principle, but the main model has still several problems:

die out: if the age groups are shifted per time step and the first element is just filled with zero, everything is shifted to the end within 4 time steps and the population dies out.
infection: in your case, the infected can only infect the same age group, so you need to summarize over the "age" groups before calculating lambda.
Finally, what is "age" group? Do you want the time since infection?

To sum up, there are several options: I would personally prefer a discrete model for such a simulation, i.e. difference equations, a age structured matrix model or an individual-based model.

If you want to keep it an ODE, I recommend to let the susceptible together as one state and to implement only the infected as stage structured.

Here a quick example, please check:

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")

Note: This is just a technical demonstration, not a valid stage structured SIR model!

This is wonderful thank you very much! Why do you think difference equations is better than the events option in ODE? My original model is more complex with 9 different states, and will include more age groups, the die out effect is what I wanted to obtain in the model. And indeed I will consider people to be infected differently according to their age group. Here I just tried to simplify it as much as possible to fully understand the aging process in the model. When you recommend to let the susceptible together as one state, is it to allow those susceptible to be infected in the same way? — Andyga, Jul 23 '21 at 09:32
I think that a difference equation model may be more consistent, as the shift of age groups occurs also at discrete time steps. But, if you have reasons to use the "mixed" approach, it is of course your decision. Implementing different infection rates of the "age" groups is of course possible. Here, just define a weighting function instead of `sum`. It is also important to define, what an "age" or "stage group" is. From the perspective of the susceptible, it may be indeed age (in years), whereas for the infected, it can be time since infection in days (as in the above example). — tpetzoldt, Jul 23 '21 at 14:03

Realistic age structured model using ODE from the deSolve package

1 Answers1