---
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