---
title: Epidemic Models (A Requested Special Topic...)
author: 36-462/662, Spring 2020
date: 16 April 2020 (small revisions, 20 April 2020)
output: slidy_presentation
bibliography: locusts.bib
---

```{r, include=FALSE}
library(knitr)
set.seed(2020-04-16)
```

\newcommand{\Prob}[1]{\mathrm{Pr}\left( #1 \right)}
\newcommand{\Expect}[1]{\mathbf{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}



## Simple Epidemic Models

- **Compartments**: Every person can be in one of three states:
  + Susceptible ($S$): not sick, could get sick
  + Infectious ($I$): sick, can make others sick
  + Removed ($R$): either recovered or dead; not sick, can't get sick, can't make others sick
- **Contagion**: When an $S$ meets an $I$, some probability of the $S$ also becoming $I$
- **Removal**: $I$s spontaneously change into $R$s
  + Maybe after some random time
  + Maybe at some steady rate per unit time
- **Mixing** or **mass action**: probability of an $S$ meeting an $I$ just depends on total numbers of $S$s, $I$s and $R$s

## Stochastic Form of the SIR Model (1)

- $S(t) =$ number of $S$s at time t, similarly for $I(t)$ and $R(t)$
- Total (initial) population is $n$
- Small time increment $h$, during which each $S$ person encounters $k h$ people
- The probability of contagion from an encounter between an $S$ and an $I$ is $c$
- Probability that a random person is infected is $I(t)/n$
- So the probability of an $S$ becoming $I$ is $kch I(t)/n$
  + Strictly it's $1-(1-c I(t)/n)^{hk} \approx kch I(t)/n$ for small $h$
- Number of _new_ infections in time $h$ is $\mathrm{Binom}(S(t), kch I(t)/n)$
  + If $h$ is very small but $S(t)$ is very large, this is approximately $\mathrm{Poisson}(\frac{kch}{n} S(t) I(t))$
- Every $I$ gets removed with probability $\gamma h$ so the number of removals in time $h$ is $\mathrm{Binom}(I(t), \gamma h)$
  + If $h$ is very small but $I(t)$ is very large, this is approximately $\mathrm{Poisson}(\gamma h I(t))$

## Stochastic Form of the SIR Model (2)

- Over-all equations:
\begin{eqnarray}
C(t+h) & = & \mathrm{Binom}(S(t), \frac{\beta}{n} h I(t))\\
D(t+h) & = & \mathrm{Binom}(I(t), \gamma h)\\
S(t+h) & = & S(t) - C(t+h)\\
I(t+h) & = & I(t) + C(t+h) - D(t+h)\\
R(t+h) & = & R(t) + D(t+h)
\end{eqnarray}

- This bundles up multiple parameters we don't get to measure separately
into some "effective" parameters $\beta$ and $\gamma$
  + I left in the factors of  $h$ and $n$ for reasons which will become clear later

## SIR is easy to simulate

```{r}
sim.SIR <- function(n, beta, gamma, s.initial=n-1, i.initial=1, r.initial=0, T) {
    stopifnot(s.initial+i.initial+r.initial == n)
    states <- matrix(NA, nrow=3, ncol=T)
    rownames(states) <- c("S", "I", "R")
    states[,1] <- c(s.initial, i.initial, r.initial)
    for (t in 2:T) {
        contagions <- rbinom(n=1, size=states["S",t-1],
                             prob=beta*states["I",t-1]/n)
        removals <- rbinom(n=1, size=states["I",t-1],
                           prob=gamma)
        states["S",t] <- states["S",t-1] - contagions
        states["I",t] <- states["I",t-1] + contagions - removals
        states["R",t] <- states["R",t-1] + removals
    }
    return(states)
}    
```

- The code doesn't need $h$ here because it's (implicitly) taking $h=1$

## What One Simulation Looks Like

```{r, echo=FALSE}
one.SIR.run <- sim.SIR(n=1000, beta=1e-1, gamma=1e-2, T=800)
plot(one.SIR.run["S",], xlab="t", ylab="People (count)",
     type="l", lty="solid",
     main="n=1000, beta=0.1, gamma=0.01")
lines(one.SIR.run["I",], lty="dashed")
lines(one.SIR.run["R",], lty="dotted")
legend("right", legend=c(expression(S(t)),
                        expression(I(t)),
                        expression(R(t))),
       lty=c("solid", "dashed", "dotted"),
       cex=0.5)
```

## Multiple Simulations with the Same Settings

```{r, echo=FALSE}
# Run the simulation b times, with the given settings, and plot them all
# together
# Inputs: number of replicates (b)
  # population size (n)
  # infection parameter (beta)
  # removal parameter (gamma)
  # duration of simulation (T)
  # optional further arguments to sim.SIR (...)
# Output: the 3D array of simulation runs, invisibly
# Side-effect: the plot is made
# Calls: sim.SIR()
plot.many.SIR.runs <- function(b, n, beta, gamma, T, ...) {
    # Run the simulation b times, with the same parameters
    many.runs <- replicate(b, sim.SIR(n=n, beta=beta, gamma=gamma, T=T))
    # Create an empty plot with the right limits, labels, etc.
    plot(0, type="n",
         xlab="t", ylab="People (count)",
         main=paste("n=", n, ", beta=", beta, ", gamma=", gamma, sep=""),
         ylim=c(0,n),
         xlim=c(1,T))
    # Now plot all 3 state variables for each simulation run, distinguishing
      # variables by line-type and the runs by color
    for (run in 1:b) {
        lines(many.runs["S",,run], lty="solid", col=run)
        lines(many.runs["I",,run], lty="dashed", col=run)
        lines(many.runs["R",,run], lty="dotted", col=run)
    }
    legend("right", legend=c(expression(S(t)),
                        expression(I(t)),
                        expression(R(t))),
       lty=c("solid", "dashed", "dotted"),
       cex=0.5)
    invisible(many.runs)
}
```

```{r, echo=FALSE}
plot.many.SIR.runs(b=10, n=1000, beta=1e-1, gamma=1e-2, T=800)
```

- Most of these show _similar_ behavior, but not _quite_ the same
  + What's up with the horizontal black line at the top?

## What If We Make Contagion Harder?

```{r, echo=FALSE}
plot.many.SIR.runs(b=10, n=1000, beta=1e-2, gamma=1e-2, T=800)
```
```

- The one initial $I$ gets removed before it can start an epidemic

## Maybe Not Quite That Much Harder?

```{r, echo=FALSE}
plot.many.SIR.runs(b=10, n=1000, beta=2e-2, gamma=1e-2, T=1600)
```

- Note the change on the horizontal axis
- Sometimes the epidemic never gets started, sometimes it takes off but doesn't get everyone before dying out...
  + Notice that when it does break out, it stops when about 80% of the population has been infected, not 100% like in the first set of simulations

## Go Back to the Original Ease of Contagion, But Make Removal Faster

```{r, echo=FALSE}
plot.many.SIR.runs(b=10, n=1000, beta=1e-1, gamma=5e-2, T=1600)
```

- Same long-run number of infected as the last simulations, though we get there faster

## Some Suggestions from Those Simulations

- Epidemic outbreak when $\beta$ is sufficiently large, compared to $\gamma$
  + At least, outbreaks with high probability
- Epidemics don't have to affect everyone before they die out
- It may be the ratio $\beta/\gamma$ is what matters for the probability of
  an outbreak and for the ultimate fraction ever infected
  + Rather than, say, the difference $\beta-\gamma$ or even $\beta^2-\gamma^2$


## Deterministic limit

- Imagine making the population size really big ($n \rightarrow \infty$), and the time step really small ($h \rightarrow 0$)
- The fluctuations get (relatively) negligible, and we're left with _deterministic_ equations
  + See backup

\begin{eqnarray}
S(0) + I(0) + R(0) & = & n\\
\frac{dS}{dt} & = & -\frac{\beta}{n} S(t)I(t)\\
\frac{dI}{dt} & = & \frac{\beta}{n} I(t)S(t) - \gamma I(t)\\
\frac{dR}{dt} & = & \gamma I(t)
\end{eqnarray}

- When people talk about "_the_ SIR model", it's often ambiguous between the stochastic model we started with and these ordinary differential equations (ODEs)
- See backup for manipulations to these equations which show how to re-write
  them so the _only_ parameter is $\beta/\gamma$

## Simulating the Deterministic Limit

Go back to the original parameter values:

```{r, echo=FALSE}
library(deSolve) # For solving differential equations

# Function to calculate the rate of change in the deterministic SIR model
  # the "t" argument, not used here, is required by deSolve
SIR <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dS <- -beta*S*I/(S+I+R)
    dI <- beta*S*I/(S+I+R) - gamma*I
    dR <- gamma*I
    list(c(dS, dI, dR))
  })
}

times <- seq(from=0, to=800, by=1)
state <- c(S=999,
           I=1,
           R=0)
params <- c(beta=1e-1,
            gamma=1e-2)
ode.sir <- ode(y=state, times=times, func=SIR, parms=params)
plot(0, type="n", xlab="t", ylab="People (count)",
     xlim=c(0,800), ylim=c(0,1000))
lines(S ~ time, data=ode.sir, lty="solid")
lines(I ~ time, data=ode.sir, lty="dashed")
lines(R ~ time, data=ode.sir, lty="dotted")
legend("right", legend=c(expression(S(t)),
                        expression(I(t)),
                        expression(R(t))),
       lty=c("solid", "dashed", "dotted"),
       cex=0.5)
```

## Is There Going to be an Epidemic?

- **Epidemic transition**: Typically, the contagion either dies out rapidly _or_ it takes over a large fraction of the population
  + We saw evidence for this in the simulation
  + More mathematically: As $n$ grows, does the total number who are ever infected grow $n$ or not?
  + I.e., is $R(\infty) = O(n)$ (an epidemic) or $R(\infty) = o(n)$ (no epidemic)?
- One way to approach this is through the **basic reproductive number** $R_0$:
\[
R_0 \equiv \Expect{\text{number of new infections directly produced by adding one infected into an all-susceptible population}}
\]
  + Don't confuse this with $R(0) =$ number removed at time $0$!
  + Some people write $\mathcal{R}_0$ to avoid this confusion
- Each of those $R_0$ new infections _also_ faces an (almost) entirely susceptible population, so they'll cause $R_0$ new infections on average as well
  + Obviously can't be true for ever, but is true initially
- After $d$ rounds of contagion, we expect $R_0^d$ infected
  + Again, initially
- _The early stage of an epidemic has exponential growth_
- There is an **epidemic threshold** at $R_0=1$
  + $R_0 < 1$: contagion dies out before getting very far (**sub-critical**)
  + $R_0 > 1$: contagion will grow exponentially, epidemic outbreak  (**super-criticial**)
  + $R_0 = 1$: Will die out, but can live a long time before doing so (**critical**)

## The Epidemic Threshold, Illustrated

```{r, include=FALSE}
library(network)
# Simulate a branching process with a geometric distribution for the
  # number of offspring, and a given number of generations,
  # plotting the result as a network
# Inputs: mean number of descendants (mu)
  # maximum number of generations to run (generations, default 3)
  # further optional graphical parameters (...)
# Output: the corresponding network, invisibly
# Side-effects: The network is plotted, color-coded by generations
# ATTN: It'd be better coding style to split this into two functions, one to
  # do the simulation and the other to do the plotting
branching.process.genealogy <- function(mu, generations=3, ...) {
    # We'll need the "network" package for plotting
    stopifnot(require(network))
    # Convert from average number of descendants to the parameter of a geometric distribution
    p <- 1/(1+mu)
    # Create a list to store the number of offspring in each generation
    gens <- vector("list", length=generations+1)
    # Create a vector to store the number of individuals produced in each
    # generation
    n.gens <- rep(0, times=generations+1)
    # The first generation begins with one individual
    gens[[1]] <- 1
    n.gens[1] <- 1
    for (k in 2:length(gens)) {
        # If the branching process has not died out,
        if (n.gens[k-1] > 0) {
            # Each individual in the previous generation gets to leave a
            # random number of descendants
            gens[[k]] <- rgeom(n=n.gens[k-1], prob=p)
        }
        # Record how many individuals we made in this generation
        n.gens[k] <- sum(gens[[k]])
    }
    # If we died out before the maximum allowed number of generations,
    # truncate at the maximum generation
    max.gen <- max(which(n.gens > 0))
    n.gens <- n.gens[1:max.gen]
    gens <- gens[1:max.gen]
    # How many individuals did we make in total?
    n <- sum(n.gens)
    # Create a big vector of labels for individuals by generations
    gen.labels <- rep(1:max.gen, times=n.gens)
    # Make a big, initially empty, matrix for displaying parent/child relations
    adj <- matrix(0,nrow=n,ncol=n)
    # Did the process branch at all?
    if (max.gen > 1) {
        # For all but the last generation
        for (k in 1:(max.gen-1)) {
            # Which of the n nodes belong to generation k?
              # They are potential parents of the next generation
            parents <- which(gen.labels==k)
            # Which nodes are their children?
            children <- which(gen.labels==(k+1))
            # Now make a vector to track  where each parents' children begin
            # in the vector of children
              # e.g., if there are 3 individuals in the generation, who have
              # 2, 0 and 4 children, then children 1 and 2 belong to the first
              # individual, and children 3, 4, 5 and 6 belong to the third
              # individual
            cumgen <- c(0, cumsum(gens[[k+1]])) + 1
            # For each potential parent
            for (l in 1:length(parents)) {
                # Create edges between the parent and the right number of
                # children
                # First check that this node has children
                if (gens[[k+1]][l] > 0) {
                    # What's the index number of this individual's first child?
                    first.child <- cumgen[l]
                    # How many children did they have?
                    num.children <- gens[[k+1]][l]
                    # What's the index number of the last child?
                    last.child <- first.child + num.children - 1
                    # Create edges pointing from this parent to all their kids
                    adj[parents[l], children[first.child:last.child]] <- 1
                 }
            }
        }
    }
    the.genealogy <- as.network(adj)
    set.vertex.attribute(the.genealogy, "generation", gen.labels)
    plot(the.genealogy, vertex.col="generation", ...)
    invisible(the.genealogy)
}
```

```{r, echo=FALSE}
par(mfrow=c(1,3))
branching.process.genealogy(mu=0.9, generations=6, main=expression(R[0]==0.9))
branching.process.genealogy(mu=1.1, generations=6, main=expression(R[0]==1.1))
branching.process.genealogy(mu=2.6, generations=6, main=expression(R[0]==2.6))
```

(Details: tracing out to a maximum of 6 rounds of growth; assuming a
geometric distribution for the number of new infections, with the mean
given by $R_0$; color-coded by generation)

## Relating $R_0$ to SIR Parameters (1)

- Start by guestimating
- $I(t)=1$, $S(t) \approx n$, so expected number of infections per time step (of length $h$) is $\approx \beta$
- Expected number of time steps before the initial $I$ gets removed is $1/\gamma$
- So (roughly!)
\[
R_0 = \frac{\beta}{\gamma}
\]
   + Our initial simulations had $\beta=10^{-1}$, $\gamma=10^{-2}$ so $R_0 = 10$
   + For the second set of simulation, $R_0 = 1$
   + For the third and fourth sets, $R_0 = 2$

## Relating $R_0$ to SIR Parameters (2)

- There are more rigorous ways to come to the same conclusion...
  + For the stochastic model: approximate it as a branching process
  + For the deterministic model: look at the stability of the all-susceptible state
- Upshot: In the SIR model,
\[
R_0 = \frac{\beta}{\gamma}
\]
  + Holds in both the stochastic and deterministic versions

## Qualitative Results for the Deterministic Model

- $dS/dt \leq 0$ so the number of susceptibles never goes up
- $dR/dt \geq 0$ so the number removed never goes down
- If we begin with $I(0) \ll n$, $R(0)=0$, $S(0) \approx n$, then at small
times $t$,
\[
I(t) \approx I(0) e^{(\beta-\gamma) t}
\]
  + See backup for derivation
- That is, early in an outbreak we expect to see exponential growth at the rate $\beta-\gamma$
  + Notice: positive growth rate iff $R_0 > 1$
  + Notice: at this early stage, $I(t)$ is _doesn't care_ about total population size $n$, so _don't_ compare outbreaks by their per-capita infection rates
- Once $S$ gets small enough, $dI/dt < 0$
  + Specifically, $dI/dt < 0$ once $\frac{S}{n} \leq \frac{\gamma}{\beta} = \frac{1}{R_0}$
- So $I(t)$ will grow exponentially and then eventually decline
  + Implication: stop an outbreak by removing enough susceptibles
- $I(t)$ is going to be a one-humped curve, usually skewed to the right
- $R(t) = \gamma \int_{0}^{t}{I(s)ds}$, it's $I(t)$ integrated over time, so it's going to look like the logistic function $\frac{e^t}{1+e^t}$
  + Or some other similar curve --- _not_ necessarily that functional form
- There are ways of calculating $R(\infty)$ but they involve solving a transcendental equation numerically [@Miller-epidemic-final-sizes]
  + The upshot is that $R(\infty)/n$ is an increasing function of $R_0$


## Back-of-the-envelope SIR Modeling for Covid-19 in the USA

- Did this back in early March (file last updated 11 March)
- Initial state:
  + $S(0) = 3\times {10}^8$ (approx. US population)
  + $I(0) = 3\times {10}^3$ (number of known US cases at the time, $\times$ a fudge factor of about 30)
  + $R(0) = 1$ just because I wanted to plot everything on a log scale
- Parameters:
  + $1/\gamma =$ average length of time someone is infectious, say 12 days
    * An early estimate; we could adjust that now
  + $R_0=\beta/\gamma$, and one plausible estimate back then was $R_0 \approx 2.6$ [@Imai-et-al-transmissibility-of-covid-19]
  + I backed out $\beta$ from $R_0$ and $\gamma$

```{r, echo=FALSE}
SIR <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dS <- -tau*S*I/(S+I+R)
    dI <- tau*S*I/(S+I+R) - gamma*I
    dR <- gamma*I
    list(c(dS, dI, dR))
  })
}
state <- c(S=3e8,
           I=3e3,
           R=0)
params <- c(tau=2.6*(1/12),
            gamma=1/12)
```

## What Did This Predict?

```{r, echo=FALSE}
# Solve the ODE to get the time series of the state variables
trajectory <- ode(y=state, times=seq(from=0, to=365, by=0.01), func=SIR,
                  parms=params)
# Add a column showing, over time, how many are or were infected
trajectory <- cbind(trajectory, cum.inf=sum(trajectory[1,]) - trajectory[,"S"])
# Add a column showing, over time, the cumulative proportion ever infected
trajectory <- cbind(trajectory, cum.prop=trajectory[,"cum.inf"]/sum(trajectory[1,c("S","I","R")]))


par(mfrow=c(1,1))
plot(cum.prop ~ time, data=trajectory,
     ylim=c(0,1), xlab="days",
     ylab="Cumulative fraction ever infected",
     type="l",
     sub="Assuming 3000 infectious on day 0, 12 day infection period, R0=2.6")
abline(h=0.5, lty="dashed")
abline(h=0.90, lty="dotted")
legend("bottomright", legend=c("Model projection", "50%", "90%"),
       lty=c("solid", "dashed", "dotted"))
```

- About 90 days to reach 50% of the population, or early-mid June
  + More exactly, `r min(trajectory[trajectory[,"cum.prop"] >= 0.5, "time"])` days
- About 180 days to reach 90% of the population, or early-mid September
  + More exactly, `r min(trajectory[trajectory[,"cum.prop"] >= 0.90, "time"])` days
- Level off at a bit over 90% of the population
  + More exactly, `r max(trajectory[,"cum.inf"])/(3e8+3001)`
- Assuming a death rate of 1% among those who contract it, $\approx `r signif(max(trajectory[,"R"])*0.01,2)`$ deaths

## How Well Did This Do?

- Use [Kieran Healy's `covdata` package](https://kieranhealy.org/blog/archives/2020/04/10/covdata-package/) to get data on _reported_ cases and deaths for the US

```{r, include=FALSE}
library(tidyverse)
library(covdata)
usa <- covnat[covnat$iso3=="USA",]
day.of.100th.case <- with(usa, min(date[cu_cases>99]))
usa$days.elapsed <- with(usa, date-day.of.100th.case)
```

- The 100th confirmed case in the US was reported on `r day.of.100th.case`, about when I did the simple model

```{r, echo=FALSE}
plot(cu_cases ~ days.elapsed, data=usa[usa$days.elapsed>=0,], log="y",
     xlab="Days elapsed since 100th confirmed case",
     ylab="Cumulative number count",
     ylim=c(1,1e6))
points(cu_deaths ~ days.elapsed, data=usa[usa$days.elapsed>=0,], pch=4)
lines(trajectory[,"time"], trajectory[,"cum.inf"], lty="dotted")
legend("topleft", legend=c("Confirmed cases",
                           "Confirmed deaths",
                           "n-S(t) from SIR model"),
       pch=c(1,4,NA),
       lty=c(NA,NA,"dotted"))
```

- Remember I assumed 30 times more true infections than known infections
  + 30 was pretty arbitrary, so we should only expect to get the shape right
  + I.e., parallel curves on the log-scaled plot
  + I also show deaths because those are (arguably) less under-counted than cases; they'd lag infections but should show parallel behavior
- The model is not doing too badly up to about day 20, or late March, by which point everything was locked down
- N.B.: This is the simplest possible model and the parameter estimation was _very_ crude, hardly even worth calling "estimation"

## Extending the Basic SIR Model

- Reinfection: Change the $I\rightarrow R$ transition to be $I \rightarrow S$ (SIS model)
  + or $I \rightarrow R \rightarrow S$ (SIRS model) if there's a period of immunity after infection
- Can add an "exposed" compartment $E$ for "have the disease and will become $I$, but not yet"
  + SEIR models are very standard when the disease takes some time to develop
- Can add an "infectious but not showing symptoms" compartment that acts like $I$
  + Might be less (or more!) effective at spreading than the classic $I$s
- Can make the basic coefficients $\beta$ and $\gamma$ functions of time, to account for seasonal effects
- Can divide space into multiple regions, give each region its own count of each compartment, and include movement between regions
  + E.g., we might have two regions $A$ and $B$, with $S_A(t)$, $S_B(t)$, $I_A(t)$, $I_B(t)$, etc., and equations like
  \[
  I_A(t) = \frac{\beta}{n}S_A(t) I_A(t) - \gamma I_A(t) - \mu_{AB}I_A(t) + \mu_{BA}I_B(t)
  \]
  where the $\mu$ coefficients represent movement from region to region
  + Or use $I_A(t) = \frac{\beta_{AA}}{n}S_A(t) I_A(t) + \frac{\beta_{AB}}{n}S_A(t) I_B(t)$ if there's cross-contagion without migration, etc.

## Aside: Why diseases do _not_ always evolve to be less deadly

- Many biologists used to think diseases would tend to evolve to be more benign on their hosts
  + Killing the hosts seems bad for the pathogen's own reproduction!
  + It's explicit in @Zinsser-rats-lice-and-history (a wonderful, if eccentric, old book about epidemics, along with many other things)
    * Zinsser had a lot of unusual-for-his-time ideas but this was not one of them
  + You still find non-biologists (often economists or [legal academics influenced by economics](https://www.newyorker.com/news/q-and-a/the-contrarian-coronavirus-theory-that-informed-the-trump-administration)...) making this argument
- More mathematically:
  + if we have two variant strains of the disease, the one with the higher $R_0$ is going to dominate
    * Because: exponential growth
  + The longer the disease leaves the host alive, the smaller $\gamma$ is
  + Smaller $\gamma$ means larger $R_0$ (all else being equal)
  + So, shouldn't evolution favor less-lethal variants of the disease?
- Biologists don't believe this any more for a very simple reason:
\[
R_0 = \frac{\beta}{\gamma}
\]
  + A mutation which kills off hosts faster (increases $\gamma$) but _also_ makes the disease easier to spread (increases $\beta$) can raise $R_0$!
  + Lots of what diseases do to their hosts look like adaptations which increase spread at the cost of killing off the host faster [@Ewald-evolution-of-infectious]
  + I.e., they are adaptations which increase both $\beta$ and $\gamma$
- Moral: pathogen evolution is not our friend


## Social Network Structure

- When we derived the SIR model, we said each $S$ encounters $k$ other individuals in a time-step of length $h$
- Some people have a lot more contacts than others!
- Use a **social network** to represent this
- $a_{ij} = 1$ if person $i$ knows person $j$ (this $i \neq I$!), otherwise $a_{ij}=0$
- $k_i = \sum_{j=1}^{n}{a_{ij}} =$ number of people $i$ knows, **degree** of $i$
- $p(k) =$ fraction of people with degree $k$, the **degree distribution**

## $R_0$ on Networks

- If you get infectious, you're only infectious for so long, and you have some _limited_ probability of passing on the disease to your contacts, call this $\tau$ (for "transmission")
  + For the SIR model, you can calculate this from $\beta$ and $\gamma$ (exercise!)
- We start with _one_ transmission event, where Irene $i$ gives the disease to Joey $j$
- $R_0 =$ How many _new_ infections does Joey cause, on average?
\[
R_0 = \tau \Expect{\text{degree of Joey}-1}
\]
  + Why $-1$?

## The Degree Distribution for a Random Friend

- Might think $\Expect{\text{degree of Joey}} = \Expect{\text{degree of random node}} = \Expect{K} = \sum_{k=0}^{\infty}{k p(k)}$
- _Joey is not a random node_
- Joey is a node reached by a random _edge_
- High degree nodes are extra likely to be reached by following edges
  + If Joey has twice the degree of James, Joey is twice as likely to be reached as James
  + But there are only so many nodes of high degree
- So $\Prob{K_{\text{Joey}} = k} \propto k p(k)$, and, to normalize, $\Prob{K_{\text{Joey}} = k} = k p(k)/\sum_{k=1}^{\infty}{k p(k)}$
- The expected value of this distribution is going to be $> \Expect{K}$
- **The friendship paradox**: "On average, your friends have more friends than you do"

## Back to $R_0$ on Networks

- Plug in to the equation from two slides ago:
\begin{eqnarray}
R_0 & = & \tau \Expect{\text{degree of Joey}-1}\\
 & = & \tau \frac{\sum_{k=1}^{\infty}{(k-1)k p(k)}}{\sum_{k=1}^{\infty}{k p(k)}}\\
& = & \tau \frac{\Expect{K^2 - K}}{\Expect{K}}\\
& = & \tau \frac{\Expect{K^2}-\Expect{K}}{\Expect{K}}\\
& = & \tau \frac{\Var{K} + \Expect{K}^2 - \Expect{K}}{\Expect{K}}
\end{eqnarray}

(because $\Var{K} = \Expect{K^2} - \Expect{K}^2$ for any variable)

## Epidemic Threshold on Networks

- We're above the epidemic threshold when $R_0 > 1$, or 
\[
\tau\left(\frac{\Var{K}}{\Expect{K}} + \Expect{K} - 1\right) > 1
\]
- So 3 things help the epidemic take off:
  + High transmissibility (large $\tau$)
  + High average degree ($\Expect{K}$ big)
  + High _variance_ in degree ($\Var{K}$ big)
- Lots of social networks have _very_ high degree variance


## Implications for Disease Control

- Removing the highest-degree nodes is _especially_ effective for blocking epidemics
  + Drives down $\Expect{K}$
  + Even more effective at driving down $\Var{K}$
  + Idea goes back to @Pastor-Satorras-Vespiganani-immunization, building off @Cohen-et-al-breakdown-of-internet-under-random-attack;@Cohen-et-al-breakdown-of-internet-under-intentional-attack
- "Remove" could mean vaccinating
- "Remove" could mean isolating
  

## Complications to the Basic Network Analysis

- **Community structure**: Different sorts of people have different probabilities of knowing each other
  + Usually high probability of ties within communities compared to across communities
  + Tight communities tend to favor _localized_ outbreaks [@Smith-et-al-on-epidemic-clustering]
    * everyone in the community gets sick
	* other communities not as affected
	* still some transmission across communities so still _pan_-demic outbreaks
- **Higher-order interactions**: Everyone in the same cafeteria, airplane, classroom, office... is exposed to each other, even if they don't directly know each other
  + One approach: treat these locations as special nodes, with ties to the people who go there [@St-Onge-et-al-school-closures]
  + The locations will have very high degree
- More detailed dynamics instead of just "what's $R_0$"
  + "Agent-based" or "individual-based" models: Simulate!
  + "Mean-field" approximations: see backup
  


## Some implications for the present situation

- Do things to reduce $\tau$
  + Reduce the amount of time people spend near each other
  + Increase physical distancing
  + Wear masks (which make it harder for what comes out of you when you sneeze to spread)
- Do things to reduce $\Expect{K}$ and $\Var{K}$
  + Close down events and public places
  + Stay home
- "Contact tracing" (= isolating those who are tied to those known to have been infected) could be effective, _if_ it could be implemented well


## Connecting to Data

- The stochastic model has a likelihood, so if we could measure $S(t), I(t), R(t)$ we could estimate the parameters
  + Or any two of them, along with a knowledge of $n$
  + Initially, everything just looks like exponential growth, so it's hard to discriminate between models
  + Initially, everything just looks like exponential growth, so even assuming an SIR model, it's easier to estimate $\beta - \gamma$ than $(\beta, \gamma)$
- Without widespread _and accurate_ testing, measurements of $S, I, R$ are all very noisy
  + Need to model sheer test error (false positives, false negatives)
  + Need to model who gets tested
  + A good example for COVID-19: @Javan-Fox-Ancel-Meyers-probability-of-COVID-19-outbreaks
- Can sometimes **calibrate** the parameters from other sources
  + E.g., $1/\gamma =$ average time someone is infectious, which could be determined from clinical studies / observations


## Backup: Non-basic reproductive numbers

- If $R_0$ is the _basic_ reproductive number, what are the _other_ reproductive numbers?
- Remember that $R_0$ is the number of new infections directly caused by introducing one infected individual at time 0 of the epidemic
- $R_t =$ the number of extra infections directly caused by increasing $I(t)$ by 1
- Generally, $R_t$ declines with $t$ as there are fewer susceptibles to infect

## Backup: From the stochastic to the deterministic SIR model

- Rate of change in susceptibles is
\[
\frac{S(t+h) - S(t)}{h} = \frac{-C(t+h)}{h} \sim -\frac{1}{h}\mathrm{Binom}(S(t), \beta h I(t)/n)
\]
- Similarly
\[
\frac{I(t+h) - I(t)}{h}  = \frac{C(t+h) - D(t+h)}{h} \sim \frac{1}{h} \mathrm{Binom}(S(t), \beta h I(t)/n) - \frac{1}{h}\mathrm{Binom}(I(t), \gamma h)
\]
- So the _expected_ rates of change are
\begin{eqnarray}
\Expect{\frac{S(t+h) - S(t)}{h} | S(t), I(t)} & = & - \frac{\beta}{n} S(t) I(t)\\
\Expect{\frac{I(t+h) - I(t)}{h} | S(t), I(t)} & = & \frac{\beta}{n} S(t) I(t) - \gamma I(t)
\end{eqnarray}
- The variance of a binomial, $\mathrm{Binom}(n,p)$ is $n p(1-p)$, so as $n\rightarrow\infty$, the fluctuations around the expected values become negligible in comparison
  + This is easy to say, but the real work of a proof goes here
- Take the limit as $n\rightarrow \infty$, $h\rightarrow 0$:
\begin{eqnarray}
\frac{dS}{dt} & = & - \frac{\beta}{n} S(t) I(t)\\
\frac{dI}{dt} & = & \frac{\beta}{n} S(t) I(t) - \gamma I(t)\\
\frac{dR}{dt} & = & \gamma I(t)
\end{eqnarray}

## Backup: Exponential growth at the start of the epidemic in the deterministic SIR model

- Suppose that $I(0) \ll n$, $R(0)=0$ and so $S(0)=n-1$, and $n$ is big
- Then initially
\begin{eqnarray}
\frac{dI}{dt} & = & \frac{\beta}{n}S(t)I(t) - \gamma I(t)\\
& = & \left(\beta\frac{S(t)}{n} - \gamma\right) I(t)\\
& \approx & (\beta - \gamma) I(t)
\end{eqnarray}
and we know how to solve this differential equation:
\[
I(t) \approx I(0) e^{(\beta-\gamma) t}
\]
- So we'll get exponential _growth_ just when $\beta > \gamma$, i.e., when $R_0 > 1$
  + If $R_0 < 1$, we'll get exponential _decay_
- Notice that $n$ doesn't matter to $I(t)$!
  + Early in epidemic, there's no shortage of new susceptibles for the contagion to spread to
- If we have data on $I(t)$ over time, we can use that to estimate $\beta-\gamma$
  + ... but not $\beta$ or $\gamma$ separately
  + Works best if $I(0)$ is big enough that the deterministic approximation is holding well, so you often need to wait until there are (say) 100 cases, rather than 1

## Backup: Eliminating parameters from the deterministic SIR model

- We can change our units for time to whatever we like, i.e., define
$u=t\gamma$, so $t=u/\gamma$
- Then (chain rule of derivatives)
\begin{eqnarray}
\frac{dR}{du} & = & \left(\frac{dR}{dt}\right) \left(\frac{dt}{du} \right)\\
& = & \left( \gamma I \right) \left(\frac{1}{\gamma} \right)\\
& = & I
\end{eqnarray}
- Doing the same thing for all the variables,
\begin{eqnarray}
\frac{dS}{du} & = & -\frac{\beta}{n\gamma} SI\\
\frac{dI}{du} & = & \frac{\beta}{n\gamma} SI - I\\
\frac{dR}{du} & = & I
\end{eqnarray}
- Now re-scale our counts so they're relative to the total population, so $s=S/n$, $i=I/n$, $r=R/n$:
\begin{eqnarray}
\frac{ds}{du} & = &  \left(\frac{ds}{dS}\right)\left(\frac{dS}{du}\right)\\
& = & \left(\frac{1}{n}\right)\left(-\frac{\beta}{\gamma}\frac{S}{n}I\right)\\
& = & - \frac{\beta}{\gamma} si
\end{eqnarray}
- Doing the same thing for all the variables,
\begin{eqnarray}
\frac{ds}{du} & = & -\frac{\beta}{\gamma} si\\
\frac{di}{du} & = & \frac{\beta}{\gamma} si - i\\
\frac{dr}{du} & = & i
\end{eqnarray}
- So the _only_ parameter left is $\beta/\gamma=R_0$
  + Units: $dR/dt$ has units of $[\text{people}]/[\text{time}]$, so $\gamma$ must have units of $1/[\text{time}]$
  + $dS/dt$ also has units of $[\text{people}]/[\text{time}]$, but $SI/n$ has units of $[\text{people}]^2/[\text{people}] = [\text{people}]$, so $\beta$ must also have units of $1/[\text{time}]$
  + $\beta/\gamma$ is a unitless ratio, which makes sense for interpreting it as the _number_ of new infections directly caused by an initial infection
- The trick of manipulating variables so that parameters are combined into unitless ratios is a common way of simplifying dynamical systems
  + Often called "dimensional analysis"
  + @Boccara-complex-systems has many good examples

## Backup: Slightly more rigorous calculation of the epidemic threshold on networks (1)

(Adapted from @MEJN-epidemics-on-networks, with fewer generating functions)

- Start by infecting a single random node, Irene
- Irene is connected to a susceptible node, Joey
- Define $\alpha_d$ as the probability that this connection starts a chain of infections of length $\leq d$
- So $1-\alpha_d$ is the probability that the chain goes for _more_ than $d$ steps
- The probability that there's an epidemic = the probability that the chain goes forever = $\lim_{d\rightarrow\infty}{(1-\alpha_d})$
- We'll figure this out by recursively writing $\alpha_d$ in terms of $\alpha_{d-1}$ and taking the limit

## Backup: Slightly more rigorous calculation of the epidemic threshold on networks (2)

- So let's try to come up with the recursion for $\alpha_d$
- Irene's connection to Joey starts a chain of length $\leq d$ in one of two
ways:
  1. Irene doesn't infect Joey
  2. Irene does infect Joey, but all chains from Joey have length $\leq d-1$
- The average probability of an $S$ being infected by a neighboring $I$ is $\tau$, so
\[
\alpha_d = (1-\tau) + \Prob{\text{Joey gets infected but only starts chains of length } \leq d-1}
\]
- 2nd term: Joey has $K_{\text{Joey}}-1$ friends other than Irene (Kalpana, Karl, \ldots Kryztof), and _none_ of those potential transmissions can start a chain more than $d-1$ steps long
- Those events each have probability $\alpha_{d-1}$
- We worked out Joey's degree distribution earlier as $\Prob{\text{Joey has degree} k} = \frac{k p(k)}{\Expect{K}}$
- Plugging in,
\[
\alpha_d = (1-\tau) + \tau\sum_{k=1}^{\infty}{\frac{k p(k)}{\Expect{K}}\alpha^{k-1}_{d-1}}
\]
- Now _define_ the right-hand side as $f(\alpha_{d-1})$, i.e,,
\[
f(x) \equiv (1-\tau) + \tau\sum_{k=1}^{\infty}{\frac{k p(k)}{\Expect{K}}x^{k-1}}
\]
so we can say that
\[
\alpha_d = f(\alpha_{d-1})
\]

## Backup: Slightly more rigorous calculation of the epidemic threshold on networks (3)

#### How does any of this help???

- "Starts a chain of length $\leq d-1$" implies "Starts a chain of length $\leq d$", so $\alpha_{d-1} \leq \alpha_d$, and $\alpha_d$ is non-decreasing as $d\rightarrow\infty$
  + Alternately: use algebra to check that $f(x) \geq x$
- But $\alpha_d$ is a probability so $\alpha_d \leq 1$ for all $d$
- Since $\alpha_d$ is non-decreasing and bounded above, $\lim_{d\rightarrow\infty}{\alpha_d} \equiv \alpha_{\infty}$ must exist and be $\leq 1$
- The limit must be a fixed point of $f$:
\[
\alpha_{\infty}=f(\alpha_{\infty})
\]
- It's easy to check that $f(1)=1$ always
  + If $\alpha_{\infty}=1$, then any chain of infections _must_ stop in a finite number of steps, and the probability of an _infinite_ epidemic is 0
- If $f^{\prime}(1) < 1$, then the $f(1)=1$ solution is stable, and this is what $\alpha_{\infty}$ will converge to
- If $f^{\prime}(1) > 1$, then the fixed point $f(1)=1$ is unstable, and $\alpha_d$ will converge to the _other_ solution of $f(x)=x$
  + Proof that there _is_ another solution: what is this, a math class?
- The derivative is
\[
\frac{df}{dx} = \tau\sum_{k=1}^{\infty}{(k-1)\frac{k p(k)}{\Expect{K}} x^{k-2}}
\]
- So the epidemic criterion is
\[
f^{\prime}(1) = \tau \sum_{k=1}^{\infty}{(k-1)\frac{k p(k)}{\Expect{K}}} > 1
\]
or
\[
\tau \left(\frac{\Expect{K^2} - \Expect{K}}{\Expect{K}}\right) > 1
\]
exactly as we got before





## Backup: "Mean-field" approximations to epidemic models on networks

- Keep track of number in each compartment with each degree, so $S_2(t)$ or $I_1(t)$
  + Idea goes back to @May-Anderson-transmission-dynamics-of-HIV
- Similarly say $S_k I_l(t)$ is the number of $S$s who have degree $k$ and are connected to $I$s of degree $l$ at time $t$
- Then $\frac{dS_k}{dt} = -\sum_{l}{\beta S_kI_l(t)}$ and similarly for the other first-order variables
- But now we need to know $S_k I_l(t)$
- Option 1, the "mean-field" approximation: $S_k I_l \approx \frac{kl}{n\Expect{K}} S_k I_l$
  + The susceptibles of degree $k$ "send out" $k S_k$ edges at random
  + There are $n \Expect{K}$ incoming edges those could match
  + $l I_l$ of those incoming edges belong to infectious of degree $l$
- Option 2, the "pairwise mean-field" approximation
   + Introduce $S_k S_l I_m$ (and so forth), so
\[
\frac{d S_k I_l}{dt} = -\gamma S_k I_l + \tau\sum_{m}{S_k S_l I_m - I_m S_k I_l - S_k I_l}
\]
   + Now approximate triples like $S_k S_l I_m$ in terms of pairs and singles
     * similar, if more elaborate, counting arguments for the simple mean-field approach above
- The pairwise mean-field approximation often gives pretty good matches to individual-level simulations, with a lot less computational work
- Details: @Kiss-Miller-Simon-epidemics-on-networks

## Further Reading

- If you want to see the code for everything in these slides, it's at [http://www.stat.cmu.edu/~cshalizi/dm/20/lectures/special/epidemics.Rmd]
  + The differential equations use the `deSolve` package [@deSolve]
- SIR or "compartment" type models of epidemics go back to the 1920s, especially to work on malaria
  + See @Bacaer-on-population-dynamics, chapters 12 and 16, for some of the history
- For a modern textbook treatment of epidemic models, see @Daley-and-Gani-epidemics
  + There's also a good chapter in @Ellner-Guckenheimer, which covers models of many other biological processes
- **Branching processes** (which I used to derive $R_0$ from the SIR model) are an important part of applied probability theory; there's a good treatment in @Grimmett-Stirzaker, among many other places
  + See also @Bacaer-on-population-dynamics, chapter 9, on their origin in the 19th century
  + @Guttorp-scientific covers statistical inference for branching processes, and many other sorts of stochastic process
- The way I derived $R_0$ for the network model follows @MEJN-epidemics-on-networks, but without (explicitly) using generating functions
  + It's an example of using **percolation theory**, another important topic in applied probability [@Grimmett-on-percolation], with many surprising applications to epidemiology [@Abundance-threshold-for-plague]
  + See @Bacaer-on-population-dynamics, chapter 22 for some of the 20th century history
- @Kiss-Miller-Simon-epidemics-on-networks is an excellent textbook on epidemics on networks, but assumes more knowledge of differential equations than is usual for statistics students
  + You could always read, say, @Arnold-ODEs
- Simulating epidemic models at the individual level is a really good way to build insight and/or to handle complicated situations which defy exact solutions
  + @Bartlett-stochastic-population-models was one of the first to appreciate the power of simulations in this area
  + This is now a huge area of active research, e.g., the recent paper @Gallagher-et-al-SPEW from our department, just on the (important) topic of initializing the simulation







## References


