Реалистичная возрастная структурированная модель с использованием ODE из пакета deSolve

Question

Я пытаюсь смоделировать реалистичную структурированную по возрасту модель, в которой все люди могут перейти в следующую возрастную группу в конце временного шага (а не стареть непрерывно с заданной скоростью), используя ODE из пакета deSolve.

Рассматривая, например, модель с двумя состояниями восприимчивый (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

Я смотрел вашу модель. Реализация сдвига в функции события в принципе работает, но у основной модели все еще есть несколько проблем:

1. вымирать: если возрастные группы сдвинуты на временной шаг, а первый элемент просто заполнен нулем, все сдвигается в конец в пределах 4 временных шагов, и популяция вымирает.
2. Инфекция: в вашем случае инфицированные могут заразить только ту же возрастную группу, поэтому вам необходимо суммировать «возрастные» группы, прежде чем вычислять лямбда.
3. Наконец, что такое «возрастная» группа? Вы хотите время, прошедшее с момента заражения?

Подводя итог, можно сказать, что есть несколько вариантов: я лично предпочел бы дискретную модель для такого моделирования, то есть уравнения разностей, матричную модель с возрастной структурой или индивидуальную модель.

Если вы хотите сохранить 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!