---
title: 'VBD: building a simple model for Zika'
author: "Zulma Cucunuba & Pierre Nouvellet"
authors: ["Zulma Cucunuba", "Pierre Nouvellet"]
date: 2019-06-10
image: img/highres/mosquito.jpg
topics: ["zika", "compartmental models"]
categories: practicals
licenses: CC-BY
---

```{r options, include = FALSE, message = FALSE, warning = FALSE, error = FALSE}
library(knitr)

opts_chunk$set(collapse = TRUE)
full_version <- FALSE
```



This practical aims to illustrate the basics of **vector borne disease (VBD) modelling** using R, with an emphasis on how the methods work. We will use a basic model for an arboviral infection as an example. In this practical we will begin by gaining some understanding of the components which contribute to `$R_0$` and how potential interventions influence transmission. Later in the practical you will construct a model of Zika transmission to look at the effects of several parameters.



## Core Concepts 

From the previous lecture, we will futher develop these concepts:

- Herd effect 
- Biology of the mosquito
- Natural history of the infection in humans
- Contact rate
- Density dependence 
- Immigration-death and age-structured models
- Infection and morbidity control / elimination of infection
- Control strategies (on vectors and on humans)

## Required packages


```{r eval = FALSE}
# install.packages("deSolve", dep = TRUE)
# install.packages("ggplot2")
# install.packages("gridExtra", dep = TRUE)
```

Then load the packages using:

```{r}
library(deSolve)
library(ggplot2)
library(gridExtra)
```

## The basic Zika model

- Sh : Susceptible Humans
- Ih : Infected/Infectious humans
- Rh : Humans recovered from infection (with lifelong immunity)
- Sv : Susceptible vectors
- Ev : Exposed vectors
- Iv : Infected vectors


## The flow diagram (part I)

In this section, please make a diagram to connect the different compartments of the model



## The parameters 

We will need several parameters to connect the different compartments of our model. 

Please, look at the suplementary material of the paper http://science.sciencemag.org/content/early/2016/07/13/science.aag0219
and look at the parameter table of this model. 

Let's find the parameter values for the model. Note that we are using all the parameters in the same time unit (days)


```{r eval = FALSE}
Lv       <-        # life span of mosquitos (in days)
Lh       <-        # life span of humans (in days)
Iph      <-        # Infectious period in humans (in days)
IP       <-        # Infectious period in vectors (in days)
EIP      <-        # Extrinsic incubation period in adult mosquitos
muv      <-        # mortality of mosquitos
muh      <-        # mortality of humans
gamma    <-        # recovery rate in humans
delta    <-        # extrinsic incubation rate
betah    <-        # Probability of transmission from vector to host
betav    <-        # Probability of transmission from host to vector
Nh       <-        # Number of humans (Population of Cali 2.4 million)
m        <-        # Vector to human ratio
Nv       <-        # Number of vectors
R0       <-        # Reproductive number
b        <-        sqrt((R0 ^2 * muv*(muv+delta) * (muh+gamma)) /
                   (m * betah * betav * delta)) # biting rate

TIME     <-        # Number of years to run the simulation for 
```


## The model (Equations)

###  Humans

`$$\ \frac{dSh}{dt}  = \mu_h N_h - \frac {\beta_h b}{N_h} S_h  I_v - \mu_h  S_h $$`    
`$$\ \frac{dIh}{dt}  = \frac {\beta_h b}{N_h}S_h I_v - (\gamma_h + \mu_h) I_h $$`    
`$$\ \frac{dRh}{dt}  = \gamma_h I_h  - \mu_h R_h$$`    

###  Vectors      

`$$\ \frac{dSv}{dt}  = \mu_v N_v  - \frac{\beta_v b} {N_h} I_h S_v - \mu_v Sv$$`    
`$$\ \frac{dE_v}{dt}  = \frac{\beta_v b} {N_h} I_h S_v - (\delta + \mu_v) Ev$$`    
`$$\ \frac{dI_v}{dt}  = \delta Ev - \mu_v I_v$$`    


## Estimating `$R_0$` (Reproductive number)

We need a formula to estimate `$R_0$`: 

`$$ R_0^2 = \frac{mb^2 \beta_h \beta_v \delta}{\mu_v (\mu_v+\delta)(\mu_h+\gamma_h)} $$`

## The flow diagram (part II)

Now that you know the equations, complete the flow diagram with the parameters and the correct connection between the different compartments.


```{r, include = full_version}
# Parameters
Lv       <- 10        # life span of mosquitos (in days)
Lh       <- 50 * 365  # life span of humans (in days)
Iph      <- 7         # Infectious period in humans (in days)
IP       <- 6         # Infectious period in vectors (in days)
EIP      <- 8.4       # Extrinsic incubation period in adult mosquitos

muv      <- 1/Lv      # mortality of mosquitos
muh      <- 1/Lh      # mortality of humans
gamma    <- 1/Iph     # recovery rate in humans
delta    <- 1/EIP     # extrinsic incubation rate


# Population size
Nh       <- 2.4 * 10^6# Number of humans (Population of Cali, Colombia)
m        <- 2         # Vector to human ratio
Nv       <- m * Nh    # Number of vectors

betah    <- 0.7       # Probability of transmission from vector to host
betav    <- 0.7       # Probability of transmission from host to vector
R0       <- 3         # Reproductive number
b        <- sqrt((R0 ^2 * muv*(muv+delta) * (muh+gamma)) /
                   (m * betah * betav * delta)) # bitting rate

TIME     <- 100       # Number of years to run the simulation 


```

## Finally, the model in R

After having your flow diagram and equations, now please complete the model below with the right parameters (PAR)

```{r, eval = FALSE}
arbovmodel <- function(t, x, params) {
  
  Sh <- x[1]    # Susceptible humans
  Ih <- x[2]    # Infected humans
  Rh <- x[3]    # Recovered humans
  Sv <- x[4]    # Susceptible vectors
  Ev <- x[5]    # Susceptible vectors
  Iv <- x[6]    # Infected vectors
  
  with(as.list(params), # local environment to evaluate derivatives
       {
         # Humans
         dSh   <-  PAR * Nh - (PAR * PAR/Nh) * Sh * Iv - PAR * Sh   
         dIh   <-  (PAR * PAR/Nh) * Sh * Iv  - (PAR + PAR) * Ih
         dRh   <-  PAR * Ih  - PAR * Rh
         
         # Vectors
         dSv  <-  muv * Nv - (PAR* PAR/Nh) * Ih * Sv - PAR * Sv 
         dEv  <-  (PAR * PAR/Nh) * Ih * Sv - (PAR + PAR)* Ev
         dIv  <-  PAR * Ev - PAR * Iv
         
         dx   <- c(dSh, dIh, dRh, dSv, dEv, dIv)
         list(dx)
       }
  )
}
```



```{r, include = full_version}
arbovmodel <- function(t, x, params) {
  
  Sh <- x[1]    # Susceptible humans
  Ih <- x[2]    # Infected humans
  Rh <- x[3]    # Recovered humans
  Sv <- x[4]    # Susceptible vectors
  Ev <- x[5]    # Susceptible vectors
  Iv <- x[6]    # Infected vectors
  
  with(as.list(params), # local environment to evaluate derivatives
       {
         # Humans
         dSh   <-  muh * Nh - (betah * b/Nh) * Sh * Iv - muh * Sh   
         dIh   <-  (betah * b/Nh) * Sh * Iv  - (gamma + muh) * Ih
         dRh   <-  gamma * Ih  - muh * Rh
         
         # Vectors
         dSv   <-  muv * Nv - (betav * b/Nh) * Ih * Sv - muv * Sv 
         dEv   <-  (betav * b/Nh) * Ih * Sv - (delta + muv)* Ev
         dIv   <-  delta * Ev - muv * Iv
         
         dx    <- c(dSh, dIh, dRh, dSv, dEv, dIv)
         list(dx)
       }
  )
}
```



## Solve the system

In this section, complete and comment the code for:

- The VALUES for the initial conditions of the system

- The ARGUMENTS of the **ode** function in the **deSolve** package.


```{r, eval = FALSE}
# Time 
times  <- seq(1, 365 * TIME , by = 1)

# Specifying parameters
params <- c(
  muv      = muv,     
  muh      = muh,     
  gamma    = gamma,   
  delta    = delta,   
  b        = b,       
  betah    = betah,   
  betav    = betav,   
  Nh       = Nh,      
  Nv       = Nv
  
)


# Initial conditions of the system
xstart <- c(Sh = VALUE?,       # meaning?? 
            Ih = VALUE?,        # meaning?? 
            Rh = VALUE?,        # meaning?? 
            Sv = VALUE?,        # meaning?? 
            Ev = VALUE?,        # meaning?? 
            Iv = VALUE?)        # meaning?? 

# Solving the equations
out <- as.data.frame(ode(y      = ARGUMENT?,   # meaning??  
                         times  = ARGUMENT?,   # meaning??     
                         fun    = ARGUMENT?,   # meaning??  
                         parms  = ARGUMENT?))  # meaning??   

```



```{r, include = full_version}
# ----------- Solving the model
times  <- seq(1, 365 * TIME , by = 1)

# Specifying parameters
params <- c(
  muv      = muv,     
  muh      = muh,     
  gamma    = gamma,   
  delta    = delta,   
  b        = b,       
  betah    = betah,   
  betav    = betav,   
  Nh       = Nh,      
  Nv       = Nv
  
)


# Initial consitions of the system
xstart<- c(Sh = Nh ,      # Initial number of Sh at T0
           Ih = 0,        # Initial number of Ih at T0
           Rh = 0,        # Initial number of Rh at T0
           Sv = Nv-1,     # Initial number of Sv at T0
           Ev = 0,        # Initial number of Ev at T0
           Iv = 1)        # Initial number of Iv at TO

# Solving the equations
out <- as.data.frame(ode(y      = xstart,     # Initial conditions
                         times  = times,      # Times
                         fun    = arbovmodel, # Model
                         parms  = params))    # Parameters


# Creating time options to display
out$years <- out$time / 365
out$weeks <- out$time / 7
```





## The results

In order to have more meaningful display of the results, convert time units *days* into *years* and into *weeks*


```{r, eval = FALSE}
# Creating time options to display
out$years <- 
out$weeks <- 
  
```


```{r, include = full_version}
# Creating time options to display
out$years <- out$time / 365
out$weeks <- out$time / 7
```


### General Behavior (Human Population)

```{r p1, include = TRUE}
# Check the general behavior of the model for the whole 100 years

p1h <- ggplot(data = out, aes(y = (Rh + Ih + Sh)/10000, x = years)) +
  geom_line(color = 'grey68', size = 1) +
  ggtitle('Total human population') +
  theme_bw() + ylab('number per 10,000')

p2h <- ggplot(data = out, aes(y = Sh/10000, x = years)) +
  geom_line(color = 'royalblue', size = 1) +
  ggtitle('Susceptible human population') +
  theme_bw() + ylab('number per 10,000')

p3h <- ggplot(data = out, aes(y = Ih/10000, x = years)) +
  geom_line(color = 'firebrick', size = 1) +
  ggtitle('Infected human population') +
  theme_bw() + ylab('number per 10,000')

p4h <- ggplot(data = out, aes(y = Rh/10000, x = years)) +
  geom_line(color = 'olivedrab', size = 1) +
  ggtitle('Recovered human population') +
  theme_bw() + ylab('number per 10,000')


grid.arrange(p1h, p2h, p3h, p4h, ncol = 2)
```


### General Behavior (Vector Population)

```{r p2, include = TRUE}
# Check the general behavior of the model 

p1v <- ggplot(data = out, aes(y = (Sv + Ev + Iv), x = years)) +
  geom_line(color = 'grey68', size = 1) +
  ggtitle('Total mosquito population') +
  theme_bw() + ylab('number')

p2v <- ggplot(data = out, aes(y = Sv, x = years)) +
  geom_line(color = 'royalblue', size = 1) +
  ggtitle('Susceptible mosquito population') +
  theme_bw() + ylab('number')

p3v <- ggplot(data = out, aes(y = Ev, x = years)) +
  geom_line(color = 'orchid', size = 1) +
  ggtitle('Exposed mosquito population') +
  theme_bw() + ylab('number')

p4v <- ggplot(data = out, aes(y = Iv, x = years)) +
  geom_line(color = 'firebrick', size = 1) +
  ggtitle('Infected mosquito population') +
  theme_bw() + ylab('number')

grid.arrange(p1v, p2v, p3v, p4v, ncol = 2)
```





### Proportion 

Let's take a more careful look at the proportions and discuss them

```{r p3, include = TRUE}

p1 <- ggplot(data = out, aes(y = Sh/(Sh+Ih+Rh), x = years)) +
  geom_line(color = 'royalblue', size = 1) +
  ggtitle('Susceptible human population') +
  theme_bw() + ylab('proportion')

p2 <- ggplot(data = out, aes(y = Ih/(Sh+Ih+Rh), x = years)) +
  geom_line(color = 'firebrick', size = 1) +
  ggtitle('Infected human population') +
  theme_bw() + ylab('proportion')

p3 <- ggplot(data = out, aes(y = Rh/(Sh+Ih+Rh), x = years)) +
  geom_line(color = 'olivedrab', size = 1) +
  ggtitle('Recovered human population') +
  theme_bw() + ylab('proportion')

grid.arrange(p1, p2, p3, ncol = 2)

```


### The First Epidemic

```{r p4, include = TRUE}
# Check the first epidemic

dat <- out[out$weeks < 54, ]

p1e <- ggplot(dat, aes(y = Ih/10000, x = weeks)) +
  geom_line(color = 'firebrick', size = 1) +
  ggtitle('Infected human population') +
  theme_bw() + ylab('number per 10,000')


p2e <- ggplot(dat, aes(y = Rh/10000, x = weeks)) +
  geom_line(color = 'olivedrab', size = 1) +
  ggtitle('Recovered human population') +
  theme_bw() + ylab('number per 10,000')


grid.arrange(p1e, p2e)
```


### Let's discuss some aspects

- Sensitivity of the model to changes in `$R_0$`.
- What are the reasons of the time lag between epidemics?
- How do we calculate the attack rate?

### Modelling control interventions

Now, using this basic model we are going to model the impact of three different types of interventions. 

1. Vaccination 
2. Bednets
3. Mosquito spraying

Try to find literature that explain these interventions and describe how you will parameterise the model. Are all of these interventions feasible? Are they cost effective?


# About this document

## Contributors

- Zulma Cucunuba & Pierre Nouvellet: initial version
- Kelly Charinga & Zhian N. Kamvar: edits

Contributions are welcome via [pull requests](https://github.com/reconhub/learn/pulls). The source file if this document can be found [**here**](https://raw.githubusercontent.com/reconhub/learn/master/content/post/practical-vbd.Rmd).



## Legal stuff

**License**: [CC-BY](https://creativecommons.org/licenses/by/3.0/)
**Copyright**: Zulma Cucunuba & Pierre Nouvellet, 2017


# References

