#' @include generics.R
#'
NULL
#' An S4 object representing the SEIaImIsRD.
#'
#' This class defines the SEIaImIsRD model, an extension of the SEIRD model.
#' The model shows how populations of susceptible, exposed, infectious and
#' recovered individuals evolve over time, in which the infectious individuals
#' are subgrouped into compartments according to different severity of symptoms, i.e. asymptomatic, mild and severe.
#'
#' Notes:
#' 1. Total initial population size is normalised to 1.
#' 2. The current model does not include natural death or birth.
#'
#' @slot initial_condition_names list of names of initial conditions
#' (characters). Default is list("S0", "E0", "I_asymptomatic0", "I_mild0", "I_severe0", "R0", "D0").
#' @slot transmission_parameter_names list of names of transmission parameters
#' (characters). Default is list("beta", "kappa", "p_symptom", "gamma", "mu").
#' @slot initial_conditions list of values for initial conditions (double).
#' @slot transmission_parameters list of values for transmission parameters (double).
#' @slot output_names list of compartments name which are used by the model and
#' incidence.
#'
#' @import deSolve
#' @import ggplot2
#' @import reshape2
#' @importFrom methods new
#' @export SEIaImIsRD
SEIaImIsRD <- setClass(Class = "SEIaImIsRD",
slots = c(
initial_condition_names = "list",
transmission_parameter_names = "list",
initial_conditions = "list",
transmission_parameters = "list",
output_names = "list"
),
prototype = list(
initial_condition_names = list("S0", "E0", "I0_asymptomatic",
"I0_mild", "I0_severe", "R0", "D0"),
transmission_parameter_names = list("beta", "kappa",
"p_symptom", "gamma", "mu"),
initial_conditions = vector(mode = "list", length = 7),
transmission_parameters = vector(mode = "list", length = 5),
output_names = list("S", "E", "I_asymptomatic",
"I_mild", "I_severe", "R", "D",
"Incidence", "Deaths")
))
#' @describeIn SEIaImIsRD Retrieves initial conditions of SEIaImIsRD model.
#'
#' @param object An object of the class SEIaImIsRD.
#' @export
setMethod("initial_conditions", "SEIaImIsRD",
function(object) object@initial_conditions)
#' @describeIn SEIaImIsRD Retrieves transmission parameters of SEIaImIsRD model.
#'
#' @param object An object of the class SEIaImIsRD.
#' @export
setMethod("transmission_parameters", "SEIaImIsRD",
function(object) object@transmission_parameters)
#' @describeIn SEIaImIsRD Setter method for initial population sizes (in fraction)
#' of the SEIaImIsRD model.
#'
#' @param object An object of class SEIaImIsRD
#' @param value A numeric named list containing the initial fraction of population groups:
#' S - susceptible,
#' E - exposed,
#' I_asymptomatic - infected with no symptom,
#' I_mild - infected with mild symptoms,
#' I_severe - infected severe symptoms that need further hospitalization,
#' D - dead due to the infection,
#' R - recovered.
#'
#' @return An object of class SEIaImIsRD with initial population and parameters
#' @export
#'
setMethod(
"initial_conditions<-",
"SEIaImIsRD",
function(object,
value = list(S = NA_real_, E = NA_real_,
I_asymptomatic = NA_real_, I_mild = NA_real_, I_severe = NA_real_,
R = NA_real_,
D = NA_real_)) {
init_pop_list <- value
names(init_pop_list) <- object@initial_condition_names
object@initial_conditions <- init_pop_list
# check if values are valid
errors <- character()
# check whether all required initial population groups are set
is_na_pop <- is.na(init_pop_list)
if (sum(is_na_pop) != 0) {
msg <- paste("Missing initial setting for population group:",
paste(names(init_pop_list)[is_na_pop], collapse = ", ")
)
errors <- c(errors, msg)
}else{
# check whether the sum of initial population is normalized to 1
sum_init_pop <- sum(unlist(init_pop_list))
if (sum_init_pop != 1) {
msg <- "Sum of initial population is not 1, please normalize"
errors <- c(errors, msg)
}
}
if (length(errors) == 0) {
object@initial_conditions <- init_pop_list
object
}
else stop(paste(errors, ", please check and rerun initial_conditions<-.\n"))
})
#' @describeIn SEIaImIsRD Setter method for transmission parameters
#' of the SEIaImIsRD model.
#'
#' @param object An object of class SEIaImIsRD
#' @param value A numeric named list of values for transmission parameters:
#' beta - a named list of the effective contact rate from each infected group (i.e. rate at which an infected individual exposes susceptible), each element with value in \code{[0,1]}
#' kappa - rate of progression from exposed to infectious (the reciprocal is the incubation period),
#' p_symptom - a named list of the probability of exposed individuals moving into each of the different infected groups,
#' Here, only probabilities for the mild (p_symptom.mild) and severe (p_symptom.severe) groups need be specified and
#' the asymptomatic probability is the remainder (1 - p_symptom.mild - p_symptom.severe). Thus we require
#' p_symptom.mild + p_symptom.severe <= 1.
#' gamma - a list of the rate of removal of each infected group (i.e. recovery rate of an infected individual), each element with value in \code{[0,1]}
#' mu - a list of the rate of disease-caused mortality of each infected group, each element with value in \code{[0,1]}
#' @return An object of class SEIaImIsRD with initial population and parameters
#' @export
#'
setMethod(
"transmission_parameters<-",
"SEIaImIsRD",
function(
object,
value = list(beta = list(asymptomatic = NA_real_, mild = NA_real_, severe = NA_real_),
kappa = NA_real_,
p_symptom = list(mild = NA_real_, severe = NA_real_),
gamma = list(asymptomatic = NA_real_, mild = NA_real_, severe = NA_real_),
mu = list(asymptomatic = NA_real_, mild = NA_real_, severe = NA_real_))) {
param_list <- value
names(param_list) <- object@transmission_parameter_names
object@transmission_parameters <- param_list
# check if values are valid
errors <- character()
# check whether all required parameters are set
is_na_params <- is.na(param_list)
if (sum(is_na_params) != 0) {
msg <- paste("Missing parameters:",
paste(names(param_list)[is_na_params], collapse = ", ")
)
errors <- c(errors, msg)
}
# check whether the lengths of beta, gamma and mu correspond to the number of infected groups
n_beta <- length(param_list$beta)
if (n_beta != 3) {
msg <- paste0(
"Length of parameter beta,", n_beta, ", is not equal to the setting ", 3)
errors <- c(errors, msg)
}
n_gamma <- length(param_list$gamma)
if (n_gamma != 3) {
msg <- paste0(
"Length of parameter gamma,", n_gamma, ", is not equal to the setting ", 3)
errors <- c(errors, msg)
}
n_mu <- length(param_list$mu)
if (n_mu != 3) {
msg <- paste0(
"Length of parameter mu,", n_mu, ", is not equal to the setting ", 3)
errors <- c(errors, msg)
}
n_p_symptom <- length(param_list$p_symptom)
if (n_p_symptom != 2) {
msg <- paste0(
"Length of parameter p_symptom,", n_mu, ", is not equal to the setting ", 2)
errors <- c(errors, msg)
}
else if (param_list$p_symptom %>% unlist %>% sum > 1){
msg <- paste0("Sum of p_symptom, ", param_list$p_symptom %>% unlist %>% sum, ", is greater than ", 1)
errors <- c(errors, msg)
}
if (length(errors) == 0) {
object@transmission_parameters <- param_list
object
}
else stop(paste(errors, ", please check and rerun transmission_parameters<-.\n"))
})
#' @describeIn SEIaImIsRD Solves the ode system.
#'
#' \deqn{\frac{dS(t)}{dt} = -S(t)(\beta.asymptomatic I_{asymptomatic}(t) + \beta.mild I_{mild}(t) + \beta.severe I_{severe}(t))}
#' \deqn{\frac{dE(t)}{dt} = S(t)(\beta.asymptomatic I_{asymptomatic}(t) + \beta.mild I_{mild}(t) + \beta.severe I_{severe}(t)) - \kappa E(t)}
#' \deqn{\frac{dI_{asymptomatic}(t)}{dt} = p_symptom.asymptomatic \kappa E(t) - (\gamma.asymptomatic + \mu.asymptomatic) I_{asymptomatic}(t)}
#' \deqn{\frac{dI_{mild}(t)}{dt} = p_symptom.mild\kappa E(t) - (\gamma.mild + \mu.mild) I_{mild}(t)}
#' \deqn{\frac{dI_{severe}(t)}{dt} = p_symptom.severe\kappa E(t) - (\gamma.severe + \mu.severe) I_{severe}(t)}
#' \deqn{\frac{dR(t)}{dt} = \gamma.asymptomatic I_{asymptomatic}(t) + \gamma.mild I_{mild}(t) + \gamma.severe I_{severe}(t)}
#' \deqn{\frac{dD(t)}{dt} = \mu.asymptomatic I_{asymptomatic}(t) + \mu.mild I_{mild}(t) + \mu.severe I_{severe}(t)}
#'
#' @param object An object of class SEIaImIsRD
#' @param times A list of time points of the simulation period
#' @param solve_method A string indicating which ode integrator to use. Default is set to 'lsoda'
#' @return A list of two dataframes: one with the time steps, time series of S,
#' E, I_asymptomatic, I_mild, I_severe and R population fractions, and one with the time steps,
#' time series of incidences and deaths population fraction
#' @export
#'
setMethod("run",
"SEIaImIsRD",
def = function(
object,
times,
solve_method = "lsoda") {
# initial population groups
pop_groups <- c(S = object@initial_conditions$S0,
E = object@initial_conditions$E0,
I_asymptomatic = object@initial_conditions$I0_asymptomatic,
I_mild = object@initial_conditions$I0_mild,
I_severe = object@initial_conditions$I0_severe,
R = object@initial_conditions$R0,
D = object@initial_conditions$D0,
C = 0)
# parameters
params <- c(beta = object@transmission_parameters$beta,
kappa = object@transmission_parameters$kappa,
p_symptom = object@transmission_parameters$p_symptom,
gamma = object@transmission_parameters$gamma,
mu = object@transmission_parameters$mu)
# ODE system RHS
ode_symptome_rhs <- function(t, pop_groups, parameters) {
with(
as.list(c(pop_groups, params)), {
dS <- -S * (beta.asymptomatic * I_asymptomatic + beta.mild * I_mild + beta.severe * I_severe)
dE <- S * (beta.asymptomatic * I_asymptomatic + beta.mild * I_mild + beta.severe * I_severe) - kappa * E
dI_a <- (1.0 - p_symptom.mild - p_symptom.severe) * kappa * E - (gamma.asymptomatic + mu.asymptomatic) * I_asymptomatic
dI_m <- p_symptom.mild * kappa * E - (gamma.mild + mu.mild) * I_mild
dI_s <- p_symptom.severe * kappa * E - (gamma.severe + mu.severe) * I_severe
dR <- gamma.asymptomatic * I_asymptomatic + gamma.mild * I_mild + gamma.severe * I_severe
dD <- mu.asymptomatic * I_asymptomatic + mu.mild * I_mild + mu.severe * I_severe
dC <- S * (beta.asymptomatic * I_asymptomatic + beta.mild * I_mild + beta.severe * I_severe)
# return the rate of cI_severenge
list(c(dS, dE, dI_a, dI_m, dI_s, dR, dD, dC))
})
}
# solving ode
output <- ode(y = pop_groups,
times = times,
func = ode_symptome_rhs,
parms = params,
method = solve_method)
output <- as.data.frame(output)
# Compute incidences and deaths
cases <- c(0, diff(output$C))
deaths <- c(0, diff(output$D))
output$Incidence <- cases
output$Deaths <- deaths
output <- output[, c("time", unlist(object@output_names))]
# Create long format of output
output <- melt(output, id.vars = "time")
output <- output[, c("time", "value", "variable")]
names(output) <- c("time", "value", "compartment")
# Added for consistency of output format across models
output$age_range <- rep("0-150", length(output$time))
# Split output into 2 dataframes: one with S,E,I, and R and one with C and D
states <- subset(output, !output$compartment %in% c("Incidence", "Deaths"))
states <- droplevels(states)
changes <- subset(output, output$compartment %in% c("Incidence", "Deaths"))
changes <- droplevels(changes)
list("states" = states, "changes" = changes)
})
#' @describeIn SEIaImIsRD Calculate the basic reproduction number (\code{R_0}) of the system using the next generation matrix approach.
#' @seealso \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6002118/pdf/main.pdf} mathematical
#' details of the next generation matrix approach.
#'
#' The R0 parameter is given by:
#' \deqn{R_0 = \rho(FV^{-1})},
#' where \deqn{F=\frac{\partial F_i{(x_0)}}{\partial x_j}} and \deqn{V=\frac{\partial V_i{(x_0)}}{\partial x_j}}
#' and \code{rho} represents the spectral radius of a matrix (i.e. the largest absolute of eigenvalue).
#' The \code{F_i} are the new infections, while the \code{V_i} transfers of infections from one compartment to
#' another. \code{x_0} is the disease-free equilibrium state.
#'
#' @param model A model of an SEIaImIsRD class object with initial_conditions and transmission_parameters set.
#' @return An R0 value
#' @export
#'
setMethod("R0", "SEIaImIsRD", function(model) {
# get required parameters:
S <- model@initial_conditions$S
beta <- model@transmission_parameters$beta
kappa <- model@transmission_parameters$kappa
p_symptom <- model@transmission_parameters$p_symptom
gamma <- model@transmission_parameters$gamma
mu <- model@transmission_parameters$mu
# define matrices F and V:
Fmat <- matrix(0, 4, 4)
Vmat <- matrix(0, 4, 4)
Fmat[1, 2:4] <- c(S * beta$asymptomatic, S * beta$mild, S * beta$severe)
Vmat[1, 1] <- kappa
Vmat[2, 1] <- -kappa * (1 - p_symptom$mild - p_symptom$severe)
Vmat[2, 2] <- gamma$asymptomatic + mu$asymptomatic
Vmat[3, 1] <- -kappa * p_symptom$mild
Vmat[3, 3] <- gamma$mild + mu$mild
Vmat[4, 1] <- -kappa * p_symptom$severe
Vmat[4, 4] <- gamma$severe + mu$severe
# calculate R0 as the spectral radius for the matrix F x V^(-1):
eigVals <- eigen(Fmat %*% (solve(Vmat)))$values
R0 = max(abs(eigVals))
return(R0)
})
#' @describeIn SEIaImIsRD Prints a compartmental diagram for the SEIaImIsRD model
#'
#' @param model an SEIaImIsRD model
#'
#' @return An ODE-compartmental structure diagram object of class html
#'
#' @export
setMethod("ode_structure_diagram", "SEIaImIsRD", function(model) {
htmltools::HTML(comomodels:::seiaimisrd_structure)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.