Realistic age structured model using ODE from the deSolve package

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!

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:

1. 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.
2. 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.
3. 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!