S1.Rmd

---
title: "S1"
output: 
  html_document:
    toc: true
    toc_float:
      smooth_scroll: FALSE
    number_sections: true
bibliography: references.bib
---

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

library(dplyr)
library(ggplot2)
library(kableExtra)
library(patchwork)
library(pomp)
library(purrr)
library(readr)
library(readsdr)
library(stringr)
library(tictoc)
library(tidyr)

source("./R_scripts/plots.R")
source("./R_scripts/R0_estimation.R")

pop_val <- 17475415

beta_vals <- 0.724637681 * c(1, 1.5, 2)

bio_df <- read_csv("./Data/bio_params.csv", 
                       show_col_types = FALSE)

bio_params        <- as.list(bio_df$Value)
names(bio_params) <- bio_df$Parameter
```

This electronic supplementary material supports the results presented in the 
main text regarding the **no-intervention** sub-models (*SEI3R* and *SEI3RD*).
Specifically, this HTML file is the rendered version of a dynamic document 
(R markdown) containing the *R* code that simulates the models and produces
the plots shown in the main text. Additionally, this file includes supplementary
information to complement the discussion in the main text.

# Basic reproduction number ($\Re_0$) 

In this section, we illustrate the procedure to derive an analytical expression 
for the *SEI3R*'s basic reproduction number via the next-generation matrix 
method. This method consists of rewriting the rates of infected states in the 
form of two matrices. One matrix corresponds to the rate of appearance of new 
infections ($V$ matrix), and the other one to the rate of transitions between 
compartments of infected individuals ($F$ matrix). The product between $F$ and 
$V^{-1}$ is known as the next-generation matrix, whose largest eigenvalue 
(spectral radius) corresponds to $\Re_0$. We obtain the analytical expressions
using *Mathematica* (see *R0.nb* file in the Github repository).

## $F$ matrix

```{r, message = FALSE}
F_matrix           <- read_csv("./Data/NI_F_matrix.csv", show_col_types = FALSE)
c_names            <- colnames(F_matrix)
c_names[1]         <- ""
colnames(F_matrix) <- c_names

F_matrix |>
  kbl(escape = FALSE) |> kable_styling(full_width = FALSE)
```

## $V$ matrix

```{r, message = FALSE}
F_matrix           <- read_csv("./Data/NI_V_matrix.csv", show_col_types = FALSE)
c_names            <- colnames(F_matrix)
c_names[1]         <- ""
colnames(F_matrix) <- c_names

F_matrix |>
  kbl(escape = FALSE) |> kable_styling(full_width = FALSE)
```

## Spectral radius

\begin{equation}

  \Re_0 = \beta  \left(\omega  \left(\frac{1}{\gamma }+\frac{1}{\nu
   }\right) + \frac{\eta -\eta  \omega }{\kappa }\right)
\end{equation}


# Simulations

## SEI3R

We simulate the SEI3R model ~1000 times to produce Fig 2 in the main text.

```{r}
set.seed(1550)

samples <- runif(997, beta_vals[[1]], beta_vals[[3]])
betas   <- c(samples, beta_vals) |> sort()

c_df <- data.frame(par_beta = betas)
```

```{r}
mdl <- read_xmile("./models/1A_SEI3R.stmx", 
                  const_list = bio_params)

fn <- "./saved_objects/1A_sens.rds"

if(!file.exists(fn)) {
  
  # Without mortality (SEI3R)
  wo_mort <- sd_sensitivity_run(mdl$deSolve_components, 
                              consts_df    = c_df,
                              integ_method = "rk4", start_time = 0,
                              stop_time    = 150, 
                              timestep     = 1 / 16,
                              multicore    = TRUE)
  
  SEI3R_df <- wo_mort |> select(time, C, iter, par_beta)
  
  saveRDS(SEI3R_df, fn)
   
} else SEI3R_df <- readRDS(fn)
```

```{r}
SEI3R_df <- SEI3R_df |>
  mutate(c     = C / pop_val,
         model = "SEI3R")
```

```{r}
g <- plot_fig_02(SEI3R_df)

ggsave("./plots/Fig_02_Attack_rate.pdf", plot = g, 
       height = 3, width = 5)
```


## SEI3RD 

We compare the simulations from the *SEI3R* against those produced by a similar
model that accounts for disease-induced mortality (*SEI3RD*). This comparison
indicates that the dynamics created by both models are nearly equivalent at a 
mortality level similar to that of the 1918 influenza pandemic 
[@Taubenberger_Morens_2006]. This comparison is not performed to indicate that 
the number of deaths is negligible. In fact, the *SEI3RD* estimates between 
*290 000* and *360 000* disease-induced deaths. However, omitting 
disease-induced mortality does not compromise the overall insights derived from 
this work and facilitates the incorporation of additional structures.

```{r}
mdl2 <- read_xmile("./models/1B_SEI3RD.stmx",
                  const_list = c(bio_params, 
                                 list(par_rho = 0.025)))
```

```{r}
fn <- "./saved_objects/1B_sens.rds"

if(!file.exists(fn)) {
  
  # With mortality (SEI3RD)
  mort_df <- sd_sensitivity_run(mdl2$deSolve_components, 
                              consts_df    = c_df,
                              integ_method = "rk4", start_time = 0,
                              stop_time    = 150, 
                              timestep     = 1 / 16,
                              multicore    = TRUE)
  
  SEI3RD_df <- mort_df |> select(time, C, iter, par_beta)
  
  saveRDS(SEI3RD_df, fn)
   
} else SEI3RD_df <- readRDS(fn)
```

```{r}
SEI3RD_df <- SEI3RD_df |>
  mutate(c     = C / pop_val,
         model = "SEI3RD")
```

```{r, fig.cap = "Fig 1. Comparison between the SEI3R and the SEI3RD", fig.height = 4}
df <- bind_rows(SEI3R_df, SEI3RD_df)
plot_fig_02(df)
```

### Final size

To further illustrate that incorporating disease-induced mortality does not 
substantially impact the dynamics of the within-host profile, we compare the
attack rates obtained from both variants at day 150.

```{r, fig.cap = "Fig 2. Distribution of the final size by model"}
df2 <- df |> filter(time == 150)

plot_final_size(df2)
```

# References