library(ecoevoapps)
library(patchwork)

There is an important category of models from epidemiology called Compartment models. These models are designed to model the spread of diseases through populations that are made of individuals in different compartments (groups). For example, a population might be made up of individuals who are susceptible to the disease, individuals who have been exposed to a disease but are not yet infected/infectious, infected individuals who can spread disease, recovered individuals, etc. Different models can include different compartments, based on the population being studied and the dynamics of the infection. For more details on compartment models, please refer to the Wikipedia page or refer to this paper.

The ecoevoapps package includes functions to simulate various specific implementations of such compartment models.

### Suscptible-Infectious-Recovered (SIR) model

One of the most common implementations of compartment models tracks the dynamics of a population comprised of individuals who are susceptible to a disease, those who are infectious, and those who have recovered. The dynamics of a such a model are governed by the following set of equations:

\begin{align} \frac{dS}{dt} &= m(S + I + R)(1 - v) - mS - \beta SI\\ \frac{dI}{dt} &= \beta SI - mI - \gamma I\\ \frac{dR}{dt} &= \gamma I - mR + mv(S + I + R) \end{align}

This model can be simulated with the run_infectiousdisease_model() function:

sir_params_vec <- c(m = .1, beta = .01, v = .2, gamma = 0)
sir_init_vec <- c(S = 100, I = 1, R = 0)
sir_time_vec <- seq(0, 25, 0.1)
sir_out <- run_infectiousdisease_model(time = sir_time_vec, init = sir_init_vec, params = sir_params_vec, model_type = "SIR")
#>   time         S        I         R
#> 1  0.0 100.00000 1.000000 0.0000000
#> 2  0.1  99.70499 1.094015 0.2009934
#> 3  0.2  99.40350 1.196509 0.3999868
#> 4  0.3  99.09479 1.308208 0.5970002
#> 5  0.4  98.77806 1.429887 0.7920533
#> 6  0.5  98.45245 1.562382 0.9851656

We can then use plot_infectiousdisease_time() to plot the trajectory of these compartments over time:

plot_infectiousdisease_time(sir_out, model_type = "SIR")

We can also use plot_infectiousdisease_portrait() to track how the compartments change relative to one another:

plot_infectiousdisease_portrait(sir_out, "S", "I") +
plot_infectiousdisease_portrait(sir_out, "S", "R") +
plot_infectiousdisease_portrait(sir_out, "R", "I")

### Susceptible-Infectious-Susceptible

In the SIR model, once an individual recovers from a disease, they are fully immune to it in the future. But in reality, we may encounter a disease from which individuals never actually develops immunity. In this case, an SIS model may be more appropriate:

\begin{align} \frac{dS}{dt} &= m(S + I) - mS - \beta SI + \gamma I\\ \frac{dI}{dt} &= \beta SI - mI - \gamma I\\ \end{align}

We can similarly simulate the dynamics of an SIS model using the same functions as above, just with different inputs:

sis_params_vec <- c(m = 0.1, beta = .01, gamma = 0)
sis_init_vec <- c(S = 100, I = 1)
sis_time_vec <- seq(0, 25, 0.1)
sis_out <- run_infectiousdisease_model(time = sis_time_vec, init = sis_init_vec, params = sis_params_vec, model_type = "SIS")
#>   time         S        I
#> 1  0.0 100.00000 1.000000
#> 2  0.1  99.90588 1.094125
#> 3  0.2  99.80301 1.196989
#> 4  0.3  99.69062 1.309385
#> 5  0.4  99.56783 1.432166
#> 6  0.5  99.43374 1.566258

plot_infectiousdisease_time(sis_out, model_type = "SIS")

### Susceptible-Exposed-Infectious-Recovered

Another structure we can explore is cases in which populations are comprised of susceptible individuals, those who are exposed but not (yet) infectious, infectious individuals, and those who have recovered. This is captured by an SEIR model:

\begin{align} \frac{dS}{dt} &= m(S + E + I + R)(1 - v) - mS - \beta SI\\ \frac{dE}{dt} &= \beta SI - aE - mE\\ \frac{dI}{dt} &= aE - mI - \gamma I\\ \frac{dR}{dt} &= \gamma I - mR + mv(S + E + I + R) \end{align}

Again, we can use the functions above with different input parameters to simulate this model:

seir_params_vec <- c(m = 0.1, beta = .1, gamma = 0.2, a = .025, v = 0)
seir_init_vec <- c(S = 100, E = 0, I = 2, R = 0)
seir_time_vec <- seq(0, 25, 0.1)
seir_out <- run_infectiousdisease_model(time = seir_time_vec, init = seir_init_vec, params = seir_params_vec, model_type = "SEIR")
#>   time         S        E        I          R
#> 1  0.0 100.00000 0.000000 2.000000 0.00000000
#> 2  0.1  98.07764 1.939821 1.943316 0.03922482
#> 3  0.2  96.26275 3.767312 1.892944 0.07699503
#> 4  0.3  94.54461 5.493511 1.848433 0.11344615
#> 5  0.4  92.91365 7.128267 1.809374 0.14870377
#> 6  0.5  91.36133 8.680395 1.775392 0.18288448

plot_infectiousdisease_time(seir_out, model_type = "SEIR")

In all the models above, we assume that individuals who are infected can eventually recover, and that there is no excess mortality associated with the disease. But we know that in some cases, disease condition is associated with increased mortality rates, which we can model by introducing a new compartment into these models:

\begin{align} \frac{dS}{dt} &= m(S + I + R)(1 - v) - mS - \beta SI\\ \frac{dI}{dt} &= \beta SI - mI - \gamma I - \mu I\\ \frac{dR}{dt} &= \gamma I - mR + mv(S + I + R)\\ \frac{dD}{dt} &= \mu I\\ \end{align}

We can simulate this model as follows:

sird_params_vec <- c(m = 0.1, beta = .1, gamma = 0.2, a = .025, v = 0, mu = 0.05)
sird_init_vec <- c(S = 100, I = 2, R = 0, D = 0)
sird_time_vec <- seq(0, 25, 0.1)
sird_out <- run_infectiousdisease_model(time = sird_time_vec, init = sird_init_vec, params = sird_params_vec, model_type = "SIRD")
plot_infectiousdisease_time(sird_out, model_type = "SIRD")