01_sir_MPOX.Rmd

---
title: "Simulating MPOX disease SIR dynamics, considering contact mixing"
author: "Shadrach Mintah, Benjamin Tommy Bavug, Kama Mary Ofuru, Victoria Nakalanzi"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE
)

# Load the necessary libraries
library(deSolve)
library(ggplot2)
library(tidyr)
```

MPOX (Monkeypox) is a viral zoonotic disease caused by the Monkeypox virus, primarily affecting central and West African regions but recently spreading to other parts of the world. It transmits from animals to humans and between humans through close contact with lesions, body fluids, or respiratory droplets. Symptoms include fever, rash, and swollen lymph nodes.

This model simulates the spread of MPOX in two distinct age groups—children and adults—using an SIR (Susceptible, Infected, Recovered) framework. It explores disease dynamics within and across these groups, focusing on how the virus spreads, the rate of recovery, and immunity formation over time.

The model divides the population into three compartments: Susceptible (S), Infectious (I), and Recovered (R) individuals. The model equations describe the flow of individuals between these compartments based on the transmission rate ($\beta$) and the recovery rate ($\gamma$).



## ASSUMPTIONS

The SIR (Susceptible-Infected-Recovered) models for children and adults make several key assumptions about the population, transmission dynamics, and disease progression. Here is a list of the main assumptions:

1. Population Assumptions
- Homogeneous Mixing: The model assumes that individuals within each age group (children or adults) mix homogeneously, meaning every susceptible individual has an equal chance of coming into contact with an infected individual.
- Fixed Population Size: The total population size is constant over time (N = S + I + R), meaning no births, deaths, or migration in or out of the population.
- Closed Population: The population is closed to external influences such as new infections from outside the population or changes due to natural death rates not related to the disease.
- Age Segmentation: The population is segmented into two distinct groups (children and adults) with no crossover between these groups in terms of contact or infection unless explicitly modeled.

2. Disease Transmission Assumptions
- Direct Transmission: The disease is transmitted through direct contact between susceptible and infected individuals. For children, this transmission is modeled with parameter beta_c, and for adults, it's modeled with beta_a.
- Constant Transmission Rate: The transmission rate (represented by beta) is constant over time and uniform across the population within each age group.
- No Latent Period: The model assumes that once an individual is infected, they are immediately capable of transmitting the disease (i.e., no latent or exposed state in the model).

3. Recovery and Infectiousness Assumptions
- Fixed Infectious Period: Infected individuals recover at a constant rate (gamma), which is the inverse of the infectious period (assumed to be 14 days in this model).
- Complete Recovery: Once recovered, individuals are assumed to gain complete immunity and cannot be reinfected (hence the transition from I to R).
- No Disease-Induced Death: The model assumes no deaths due to the disease (i.e., infected individuals either recover or remain infected but do not die).

4. Initial Conditions Assumptions
- Initial Infected Proportion: A small portion of the population is infected at the start (I0), and the rest are susceptible (S0 = N - I0). No one starts in the recovered (R) category.
- Equal Infectiousness Across Age Groups: While the population is split into children and adults, it is assumed that individuals within each group have the same transmission and recovery rates.

5. No Behavior Changes Over Time
- No Intervention: The model does not account for behavioral changes, medical interventions, vaccination, or other non-pharmaceutical interventions (e.g., quarantine, social distancing) that could alter the course of the disease.

6. Time Scale
- Fixed Time Horizon: The model runs over a fixed period (365 days in this case) with no changes to model parameters or structure during this time.


## The model function

The SIR model equations is given by:

\begin{align}
\frac{dS}{dt} & = \color{orange}{\frac{-\beta S I}{N}} \\
\frac{dI}{dt} & = \color{orange}{\frac{\beta S I}{N}} - \color{blue}{\gamma I} \\
\frac{dR}{dt} & = \color{blue}{\frac{\gamma I}{N}}
\end{align}

where:

- $N$ is the population
- $\frac{dS}{dt}$ is the change in number of $Susceptible$ over time,
- $\frac{dI}{dt}$ is the change in number of $Infected$ over time,
- $\frac{dR}{dt}$ is the change in number of $Recovered$ over time,
- $\beta$ is the transmission rate,
- $\gamma$ is the recovered rate.


but

\begin{align}
\beta_c = (b1 + b2) ---(1) \\
\beta_a = b3 ------(2) 
\end{align}

- $b1$ is the transmission by direct contact
- $b2$ is the transmission by airborne
- $b3$ is the transmission by sexual contact

- $\beta_1$ is the transmission rate via direct contact/airbone
- $\beta_2$ is the transmission rate via sexual contact


which translates into:

Transmission from children

\begin{align}
\frac{dS_c}{dt} & = \color{orange}{-\beta_c S_c I_c} \\
\frac{dI_c}{dt} & = \color{orange}{\beta_c S_c I_c} - \color{blue}{\gamma I_c} \\
\frac{dR_c}{dt} & = \color{blue}{\gamma I_c}
\end{align}

Transmission from adults

\begin{align}
\frac{dS_a}{dt} & = \color{orange}{-\beta_a S_a I_a} \\
\frac{dI_a}{dt} & = \color{orange}{\beta_a S_a I_a} - \color{blue}{\gamma I_a} \\
\frac{dR_a}{dt} & = \color{blue}{\gamma I_a}
\end{align}


where:

- $S_a$ is the number of susceptible individuals(adults),
- $I_a$ is the number of infectious individuals (adults),
- $R_a$ is the number of removal individuals(adults),
- $S_c$ is the number of susceptible individuals(children),
- $I_c$ is the number of infectious individuals (children),
- $R_c$ is the number of removal individuals(children)

$However$

To account for CONTACT MIXING among children and adults after being infected:

for children

\begin{align}
\frac{dS_c}{dt} & = \color{orange}{-(\beta_cc S_c I_c + \beta_ac S_c I_a)} \\
\frac{dI_c}{dt} & = \color{orange}{(\beta_cc S_c I_c + \beta_ac S_c I_a)} - \color{blue}{\gamma I_c} \\
\frac{dR_c}{dt} & = \color{blue}{\gamma I_c}
\end{align}

for adults

\begin{align}
\frac{dS_a}{dt} & = \color{orange}{-(\beta_aa S_a I_a + \beta_ca S_c I_a)} \\
\frac{dI_a}{dt} & = \color{orange}{(\beta_aa S_a I_a + \beta_ca S_c I_a)} - \color{blue}{\gamma I_a} \\
\frac{dR_a}{dt} & = \color{blue}{\gamma I_a}
\end{align}

where:

- $\beta_cc$ is Child to child transmission rate ($\beta_cc$ = $R_0$ * gamma / $N_c$)
- $\beta_aa$ is Adult to adult transmission rate ($\beta_aa$ = $R_0$ * gamma / $N_a$)
- $\beta_ca$ is Child to adult transmission rate ($\beta_ca$ = $R_0$ * gamma / $N_a$)
- $\beta_ac$ is Adult to child transmission rate ($\beta_ac$ = $R_0$ * gamma / $N_c$)

and

- $R_0$ is the Reproduction number of MPOX ($R_0$ = 1.8 for our discussion)
- $N_c$ is population of children below 18 years (<18yrs) in the country (Nigeria)
- $N_a$ is population of children above 18 years (>18yrs) in the country (Nigeria)


## Model for children below 18yrs
```{r SIR_below_18yrs}
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE
)


# Define the SIR model for children
sir_model_children <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    dS_c <- -beta_c * S_c * I_c
    dI_c <- beta_c * S_c * I_c - gamma * I_c
    dR_c <- gamma * I_c
    return(list(c(dS_c, dI_c, dR_c)))
  })
}

# Population and parameters for children
N_c <- 95000000  # Number of children under 18 years (45% of total population in Nigeria)
I0_c <- 20       # Initial number of infected children
S0_c <- N_c - I0_c  # Initial number of susceptible children

# Define R0, infectious period, and gamma
R0_c <- 1.8
infectious_period_c <- 14
gamma_c <- 1 / infectious_period_c  # Recovery rate for children

# Calculate beta using the formula: R0 = beta / gamma
beta_c <- R0_c * gamma_c / N_c  # Transmission rate for children

# Parameters for the children's model
params_c <- c(beta_c = beta_c, gamma = gamma_c)

# Initial conditions for children
inits_c <- c(S_c = S0_c, I_c = I0_c, R_c = 0)

# Time points
time_c <- 1:365

# Solve the model for children
results_children <- deSolve::lsoda(y = inits_c, times = time_c, func = sir_model_children, parms = params_c)

# Convert the results to a data frame
results_children <- as.data.frame(results_children)
head(results_children)
tail(results_children)

# Plotting the results for children
results_long_children <- results_children |> 
  pivot_longer(cols = c(S_c, I_c, R_c), names_to = "compartment", values_to = "value")

ggplot(results_long_children, aes(x = time, y = value, color = compartment)) +
  geom_line(linewidth = 1) +
  labs(title = "SIR Model for Children (Under 18 years)", x = "Time", y = "Population Size")

```


## Model for adults above 18yrs
```{r SIR_above_18yrs}
# Define the SIR model for adults
sir_model_adults <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    dS_a <- -beta_a * S_a * I_a
    dI_a <- beta_a * S_a * I_a - gamma * I_a
    dR_a <- gamma * I_a
    return(list(c(dS_a, dI_a, dR_a)))
  })
}

# Population and parameters for adults
N_a <- 125000000  # Number of adults above 18 years (55% of total population in Nigeria)
I0_a <- 20        # Initial number of infected adults
S0_a <- N_a - I0_a  # Initial number of susceptible adults

# Define R0, infectious period, and gamma for adults
R0_a <- 1.8
infectious_period_a <- 14
gamma_a <- 1 / infectious_period_a  # Recovery rate for adults

# Calculate beta using the formula: R0 = beta / gamma
beta_a <- R0_a * gamma_a / N_a  # Transmission rate for adults

# Parameters for the adult model
params_a <- c(beta_a = beta_a, gamma = gamma_a)

# Initial conditions for adults
inits_a <- c(S_a = S0_a, I_a = I0_a, R_a = 0)

# Time points
time_a <- 1:365

# Solve the model for adults
results_adults <- deSolve::lsoda(y = inits_a, times = time_a, func = sir_model_adults, parms = params_a)

# Convert the results to a data frame
results_adults <- as.data.frame(results_adults)
head(results_adults)
tail(results_adults)

# Plotting the results for adults
results_long_adults <- results_adults |> 
  pivot_longer(cols = c(S_a, I_a, R_a), names_to = "compartment", values_to = "value")

ggplot(results_long_adults, aes(x = time, y = value, color = compartment)) +
  geom_line(linewidth = 1) +
  labs(title = "SIR Model for Adults (Above 18 years)", x = "Time", y = "Population Size")

```



```{r contact_mixing_model}
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE
)

# Define the new SIR model with contact mixing between adults and children
sir_model_mixed <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    # For children
    dS_c <- -(beta_cc * S_c * I_c + beta_ac * S_c * I_a)
    dI_c <- (beta_cc * S_c * I_c + beta_ac * S_c * I_a) - gamma * I_c
    dR_c <- gamma * I_c
    
    # For adults
    dS_a <- -(beta_aa * S_a * I_a + beta_ca * S_a * I_c)
    dI_a <- (beta_aa * S_a * I_a + beta_ca * S_a * I_c) - gamma * I_a
    dR_a <- gamma * I_a
    
    return(list(c(dS_c, dI_c, dR_c, dS_a, dI_a, dR_a)))
  })
}

# Population size and parameters
N_c <- 95000000   # Number of children under 18 years
N_a <- 125000000  # Number of adults above 18 years
I0_c <- 10        # Initial infected in children
I0_a <- 10        # Initial infected in adults

# Initial conditions for both groups
inits_mixed <- c(S_c = N_c - I0_c, I_c = I0_c, R_c = 0, S_a = N_a - I0_a, I_a = I0_a, R_a = 0)

R0 <- 1.8
infectious_period <- 14
gamma <- 1 / infectious_period

# Transmission rates
beta_cc <- R0 * gamma / N_c  # Child to child
beta_aa <- R0 * gamma / N_a  # Adult to adult
beta_ca <- R0 * gamma / N_a  # Child to adult
beta_ac <- R0 * gamma / N_c  # Adult to child

# Parameters for the mixed model
params_mixed <- c(
  beta_cc = beta_cc,
  beta_aa = beta_aa,
  beta_ca = beta_ca,
  beta_ac = beta_ac,
  gamma = gamma
)

# Time points
dt <- 1:365

# Solve the model
library(deSolve)
results_mixed <- deSolve::lsoda(
  y = inits_mixed,
  times = dt,
  func = sir_model_mixed,
  parms = params_mixed
)

# Manipulate and interpret the results
results_mixed <- as.data.frame(results_mixed)
head(results_mixed)
tail(results_mixed)

# Plot the results
library(ggplot2)
library(tidyr)

# Reshape the data for plotting
results_long_mixed <- results_mixed |>
  pivot_longer(
    cols = c(2:7),
    names_to = "compartment",
    values_to = "value"
  )

sir_plot_mixed <- ggplot(
  data = results_long_mixed,
  aes(
    x = time,
    y = value,
    color = compartment
  )
  ) +
  geom_line(linewidth = 1) +
  labs(
    title = "SIR model with contact mixing between children and adults",
    x = "Time",
    y = "Number of individuals"
  )
plot(sir_plot_mixed)

```