This is a simple implementation of a SARS model (Gumel et al. 2004) adapted from code given by <http://biom300.weebly.com/simulating-disease-models-examples.html>. The model is not realistic. It would be very difficult to estimate values for all the parameters. Running sensitivity analysis on a model with so many assumptions would be challenging. However the results provided in the paper appear to be reproducible using the model structure provided.
library(tidyverse)
library(deSolve)
library(ggplot2)
library(dygraphs)
library(reshape)
library(xts)
sars_model <- function(t, state, parameters) {
P <- unname(parameters['P']) # growth of susceptible population
B <- unname(parameters['B']) # transmission rate
e_E <- unname(parameters['e_E']) # relative transmission from asymptomatic individuals
e_Q <- unname(parameters['e_Q']) # relative transmission from quarantined individuals
e_J <- unname(parameters['e_J']) # relative transmission from isolated individuals
m <- unname(parameters['m']) # background mortality rate
p <- unname(parameters['p']) # migration rate of asymptomatic individuals
g_1 <- unname(parameters['g_1']) # rate of quarantine
K_1 <- unname(parameters['K_1']) # rate of symptom development for asymptomatic individuals
K_2 <- unname(parameters['K_2']) # rate of symptom development for quarantined individuals
g_2 <- unname(parameters['g_2']) # rate of isolation
s_1 <- unname(parameters['s_1']) # rate of recovery for symptomatic individuals
d_1 <- unname(parameters['s_1']) # rate of SARS-death for symptomatic individuals
s_2 <- unname(parameters['s_2']) # rate of recovery for isolated individuals
d_2 <- unname(parameters['d_2']) # rate of SARS-death for isolated individuals
S <- unname(state['S'])
E <- unname(state['E'])
Q <- unname(state['Q'])
I <- unname(state['I'])
J <- unname(state['J'])
R <- unname(state['R'])
D <- unname(state['D'])
# total population size
N <- S + E + Q + I + J + R
dSdt <- P - S * B * (I + e_E * E + e_Q * Q + e_J * J) / N - m * S
dEdt <- p + S * B * (I + e_E * E + e_Q * Q + e_J * J) / N - (g_1 + K_1 + m) * E
dQdt <- g_1 * E - (K_2 + m) * Q
dIdt <- K_1 * E - (g_2 + s_1 + d_1 + m) * I
dJdt <- g_2 * I + K_2 * Q - (s_2 + d_2 + m) * J
dRdt <- s_1 * I + s_2 * J - m * R
dDdt <- d_1 * I + d_2 * J
list(c(dSdt,dEdt,dQdt,dIdt,dJdt,dRdt,dDdt))
}initial_sars <- c(S=4000000, E=6, Q=0, I=1, J=0, R=0, D=0)
times <- seq(0, 1000, 0.1)
## parameters with no quarantine or isolation
parameters_sars_1 <- c(P=0, p=0.06, m=0.000034, K_1=0.1, K_2=0.125, s_1=0.0337, s_2=0.0386, d_1=0.0079, d_2=0.0068, g_1=0, g_2=0, B=0.2, e_J=0.36, e_E=0.1, e_Q=0)
## simulate the model
out_1 <- ode(y=initial_sars, times=times, func=sars_model, parms=parameters_sars_1)
## parameters with quarantine but no isolation
parameters_sars_2 <- c(P=0, p=0.06, m=0.000034, K_1=0.1, K_2=0.125, s_1=0.0337, s_2=0.0386, d_1=0.0079, d_2=0.0068, g_1=0.2, g_2=0, B=0.2, e_J=0.36, e_E=0.1, e_Q=0)
## simulate the model
out_2 <- ode(y=initial_sars, times=times, func=sars_model, parms=parameters_sars_2)
## parameters with isolation but no quarantine
parameters_sars_3 <- c(P=0, p=0.06, m=0.000034, K_1=0.1, K_2=0.125, s_1=0.0337, s_2=0.0386, d_1=0.0079, d_2=0.0068, g_1=0, g_2=0.2, B=0.2, e_J=0.36, e_E=0.1, e_Q=0)
## simulate the model
out_3 <- ode(y=initial_sars, times=times, func=sars_model, parms=parameters_sars_3)d <- as.data.frame(out_1)
dygraph(d)d <- as.data.frame(out_2)
dygraph(d)d <- as.data.frame(out_3)
dygraph(d)