Introduction

Spatial models are critical to the study of infectious disease. Epidemics spread across geographic environments and understanding the role that geography has to play provides important information for protecting population from disease. There are a number of approaches to spatial modeling (Keeling and Rohani 2008). Some researchers use lattice environments (Kao 2003) with discrete sub-populations connected by migration along cardinal directions. Others use probabilistic point-pattern dynamics, which assign infection status based on pair-wise distances to other infected (Hamada and Takasu 2019). Some even ignore explicit geographic distances and build network models of disease spread (Kiss, Green, and Kao 2005). Each of these types of models are useful for representing the spread of infection, though complexity often limits generalizability.

One of the most common lattice models in mathematical biology is the Stepping Stone model by Motoo Kimura (Kimura and Weiss 1964). The model, designed for population genetics, represents gene flow over a set of discrete “demes” or sub-populations connected by a migration rate \(m\). Migration is only allowed between nearest neighbor demes \(r+1\) and \(r-1\) on the lattice. Each deme contains a constant \(N\) individuals, with evolutionary forces such as mutation, natural selection and genetic drift acting on allele frequencies. The model captures how these forces change the distribution of allele frequencies over geography.

A schematic of the one-dimensional Stepping Stone model by Kimura.
A schematic of the one-dimensional Stepping Stone model.

Project aims

In this project, we examine SIS model dynamics in a Stepping Stone migration framework. By incorporating lattice migration, we observe how an epidemic spreads across space. Our goal is to study the relationships between our parameters and the spread of infection. First, we replicate the Stepping Stone model dynamics in one spatial dimension and incorporate SIS dynamics using \(\beta\) and \(\gamma\) for COVID-19. Then, we expand the model to two spatial dimensions and compute important quantities from our simulations.

Finally, we introduce a “social distancing” parameter \(d\) that decreases the typical Stepping Stone migration rate. We examine pre-equilibrium dynamics to measure how distancing slows the spread of infection. We argue that SIS dynamics are sufficient to observe this process, since before equilibrium, SIS essentially represents the early stages of an outbreak, before enough removals can occur to slow the infection. These early stages of an epidemic are critical for the geographic spread of a disease like COVID-19, where everyone is initially susceptible.

Main questions

We present four main questions of interest that we will answer throughout this project.

  1. What fraction of the population is expected to become infected?
  1. How far does the infection spread beyond its point of origin?
  1. How long does the infection take to reach across the lattice?
  1. How does social distancing slow the spread?

Theory

The standard SIS model

Recall the differential equations of the standard SIS model. \[ \frac{dS}{dt} = -\beta S I \\ \frac{dI}{dt} = \beta S I - \gamma I \] Where \(S\) and \(I\) represent the fraction of the population of size \(N\) that are susceptible or infected, respectively. \(\beta\) is the transmission rate and \(\gamma\) is the recovery rate.

We can write this with just one equation if we only want to simulate change in the fraction of infected.

\[ S+I = 1 \\ \implies \frac{dI}{dt} = \beta I (1-I) - \gamma I \]

One-dimensional habitat

As mentioned, we incorporate SIS dynamics in a Stepping Stone model framework. In the Stepping Stone model, we have \(L\) many discrete demes indexed by position \(r\). A symmetric migration rate \(m\) allows for spread of infection between nearest neighbor demes.

We can track changes in infection frequency at each deme \(r\) between time generations as \(I_r(t+1) - I_r(t)\). The Stepping Stone model incorporates changes in infection across the lattice as follows:

\[ I_r(t+1) - I_r(t) = m(I_{r-1}+I_{r+1}-2I_r) \] Where each generation \(I\) fraction of infected individuals from nearest neighbors \(r+1\) and \(r-1\) move into deme \(r\). There are also two directions for individuals in \(r\) to leave from, so we subtract \(2I_r\). This discrete form is the main contribution of the Stepping Stone model. We will incorporate this into simulations, but will focus on our addition of SIS dynamics.

The infection spread over space and time is represented by the following difference equation:

\[ I_r(t+1) - I_r(t) = \beta I_r(1-I_r) - \gamma I_r + (m-d)(I_{r-1}+I_{r+1}-2I_r) \] Where \(d\) is the “social distancing” parameter that decreases the migration rate.

We also want to consider stochastic changes to infection frequency. We model this by considering the conditional probability of having \(Y\) number of infected individuals in deme \(r\) at time \(t+1\). This only depends on the current number of infected in that location. Therefore:

\[ Y_r(t+1)|\vec{I}(t) \sim Binom(N,I_r(t+1)) \\ \, \\ I_r(t+1)|\vec{I}(t) \sim \frac{1}{N}Binom(N,I_r(t+1)) \] Where \(\vec{I}(t)\) is a vector of infection frequencies across the spatial lattice and \(I_r(t+1)\) is given by: \[ I_r(t+1) = I(t) + \beta I_r (1-I_r) - \gamma I_r + (m-d)(I_{r-1}+I_{r+1}-2I_r) \] This allows for random changes due to finite population size. Also known as “drift.”

All together, this equation captures the four main forces of interest:

  1. Infectivity with rate \(\beta\)

  2. Recovery with rate \(\gamma\)

  3. Migration with rate \(m\)

  4. Drift with rate \(\frac{1}{N}\)

We consider each of these forces in our discrete simulations.

Simulations

Here we employ our previous equations to simulate changes in the infection frequency over space and time.

One-dimensional habitat

spatial_sis_1d <- function(num_gens, L, pop_size, beta, gamma, m,d){
  f0 = 1/pop_size # set initial fraction of infected
  mid<-floor(L/2) # set initial location 
  I <- array(data=0,dim = c(num_gens,L+1)) # initialize matrix of I values
  for(i in 1:(num_gens-1)){ # loop of all time
    if(sum(I[i,] > 0)==0){ # check if infection is extinct, if so reset
        I[i,mid] <- f0
      }
    for(j in 2:L){ # loop for all space
      # difference equation
      dI = beta*I[i,j]*(1-I[i,j])  - gamma*I[i,j] + (m-d)*(I[i,j-1] + I[i,j+1] - 2*I[i,j])
      p = I[i,j] + dI # update frequency
      p[which(p>1)] = 1 # bound to be less than 1
      p[which(p<0)] = 0 # bound to be greater than 0
      I[i+1,j] = rbinom(n = length(p),size = pop_size, prob = p)/pop_size # sample using binomial
    }
  }
return(I)
}
num_gens = 1e3 # number of time generations to simulate
L = 50 # number of demes
pop_size = 1e4 # population size N
beta = 0.2875 # transmission rate
gamma = 1/8 # recovery rate
m = 1e-1 # migration rate
d=0 # distancing

I_1d <- spatial_sis_1d(num_gens,L,pop_size,beta,gamma,m,d)

Time series in one deme

We plot changed in infection frequency at one location over time.

mid <- floor(L/2) # identify starting location of epidemic
plot(I_1d[,mid+10],type="l", main ="I vs. time in 1 deme",
     xlab="time (generations)", 
     ylab ="frac