SIMC_SIR_Exercise_4.Rmd

---
title: "SIMC_SIR_Exercise_4"
author: "Daniel T Citron"
date: '2023-07-09'
output: html_document
---

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

# Exercise 4: Interventions

Here we will practice altering the starting conditions and structure of our SIR model in order to represent what will happen when we add interventions. Each version of our model represents a different modeling scenario. By comparing different scenarios, we can imagine different ways that the outbreak will occur under different conditions.

# Part 1: Vaccination with a perfect vaccine
Suppose we are able to vaccinate a certain fraction of susceptibles in the population before the outbreak starts, such that they never become infected.

We represent this by removing a fraction of the susceptible population at the start, and see whether the outbreak takes off. Let p.vaccinated = fraction of susceptibles who are vaccinated.

Start by letting p.vaccinated = .1:
```{r}
# Using the optimized parameter values from before
# (if you found something different, feel free to use those values here)
parameters = c(beta = 2.559449, gamma =  2.109450)

# Define times - start at week 10
times = seq(from = 10, to = 75, by = 1)

# Adjust vaccinated proportion p.vaccinated here:
p.vaccinated = .1
# Initial state values
initial_state_values <- c(S = 55000*(1 - p.vaccinated) - 1,
                         I = 1,
                         R = 55000*(p.vaccinated)
)
sir.vaccinated <- deSolve::ode(initial_state_values, 
                        seq(from = 10, to = 75, by = 1), 
                        sir_model, 
                        parms = fitval$par, )
sir.vaccinated <- data.table::as.data.table(sir.vaccinated)



```

Plot the new model results. The original model fit is shown in red, while the new scenario with vaccination is in blue:
```{r}
p <- ggplot() + 
  geom_line(data = sir.vaccinated, 
            mapping = aes(x = time, y = I), 
            color = 'blue', size = 2) + 
    geom_line(data = sir.fit, 
            mapping = aes(x = time, y = I), 
            color = 'red', size = 2 ) + 
  geom_point(data = plague.dat, mapping = aes(x = time, y = Deaths), size = 3) + 
  ylab("Deaths") + xlab("time (Weeks)") + 
  theme(axis.text = element_text(size = 20),
        axis.title = element_text(size = 24)) + 
  theme_bw()
p
```

## Question: Cases averted through vaccination
How many people became infected overall?
```{r}
sir.vaccinated$R %>% tail(1)
```

How does that compare to the number of people who became infected before?
```{r}
sir.fit$R %>% tail(1)
```

The difference between these two numbers is the number of cases/deaths averted by intervening with the vaccine.

## Question: stopping an outbreak with vaccinations
Try adjusting p.vaccinated: how large does it need to be in order to stop an outbreak from occurring?

## Optional exercise

In this new vaccination scenario, we have altered the structure of our model.

Find an expression for R_0 as a function of beta, gamma, and p.vaccinated. Solve for p.vaccinated such that R_0 = 1. Compare the value of p.vaccinated to the value that you found which was high enough to stop an outbreak from occurring.

## Optional exercise

Suppose instead of providing perfect protection against plague, the vaccine only protected 50% of people who received it. What structural changes would we need to make to our model in order to represent this? How do your answers to each of the above questions change in this alternate scenario?

# Part 2: Treatment of infected individuals
Suppose that treatment decreases the amount of infectious time by 25%. 
The way we represent this is by dividing the "Infected" compartment into two, where one receives the treatment and the other does not. In the treated compartment, gamma -> gamma/(.75) (recovery rate increases by 1/.75).

How many cases need to be treated in order to prevent an outbreak?

Start by defining a new model with the treated compartment. This time there are four state variables, where we've added I_t to represent the number of infected people on treatment:
```{r SIR model with treatment}

sir_model_treatment <- function(time, state, parameters) {
  # Unpack Parameters
  with(as.list(c(state,parameters)), {
      # Calculate total population size
      N <-  S + I_t + I + R
      
      # Calculate force of infection
      lambda <- beta * (I + I_t) / N
      
      # Calculate derivatives
      dS <- -lambda * S
      # Notice now that there is this new compartment of people who receive treatment
      dI_t <- lambda * S * p.treatment - gamma/.75 * I_t
      # And this compartment where people do not receive treatment
      dI <- lambda * S * (1-p.treatment) - gamma * I
      dR <- gamma * I + 2*gamma * I_t
      
      # Return derivative
      return(list(c(dS, dI, dI_t, dR)))
    }
  )
}

```

p.treated is the probability that an individual receives treatment. Start by setting that equal to 0.25:
```{r}
# Using the optimized parameter values from before
# (if you found something different, feel free to use those values here)
# Also specify p.treated, the fraction of infected people who receive treatment
parameters = c(beta = 2.559449, gamma =  2.109450, p.treatment = 0.25)

# Define times - start at week 10
times = seq(from = 10, to = 75, by = 1)

# Initial state values
initial_state_values <- c(S = 54999,
                         I = 1,
                         I_t = 0,
                         R = 0
)
sir.treated <- deSolve::ode(initial_state_values, 
                        seq(from = 10, to = 75, by = 1), 
                        sir_model_treatment, 
                        parms = parameters )
sir.treated <- data.table::as.data.table(sir.treated)
```


Plot the new model results. The original model fit is shown in red, while the new scenario with treatment is in blue:
```{r}
p <- ggplot() + 
  geom_line(data = sir.treated, 
            mapping = aes(x = time, y = I+I_t), 
            color = 'blue', size = 2) + 
    geom_line(data = sir.fit, 
            mapping = aes(x = time, y = I), 
            color = 'red', size = 2 ) + 
  geom_point(data = plague.dat, mapping = aes(x = time, y = Deaths), size = 3) + 
  ylab("Deaths") + xlab("time (Weeks)") + 
  theme(axis.text = element_text(size = 20),
        axis.title = element_text(size = 24)) + 
  theme_bw()

p
```

## Question: Cases averted through treatment
How many people became infected overall when we add treatment?
```{r}
sir.treated$R %>% tail(1)
```

How does that compare to the number of people who became infected in the original scenario?
```{r}
sir.fit$R %>% tail(1)
```

The difference between these two numbers is the number of cases/deaths averted by intervening with treatment

## Question: stopping an outbreak with treatment
Try adjusting the parameters in the treatment scenario - what coverage is needed to stop an epidemic altogether?

## Optional exercise

In this new treatment scenario, we have altered the structure of our model.

Find an expression for R_0 as a function of beta, gamma, p.treatment, and the factor which changes the rate at which infected people recover with treatment. Solve for p.treatment such that R_0 = 1. Compare the value of p.treatment to the value that you found which was high enough to stop an outbreak from occurring.

## Optional exercise

Suppose that treatment does not shorten the duration of the infectious period, and instead it reduces the per-contact rate of transmission (beta) by 50%. What structural changes would we need to make to our model in order to represent this? How do your answers to each of the above questions change in this alternate scenario?