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?

2. How far does the infection spread beyond its point of origin?

3. How long does the infection take to reach across the lattice?

4. How does social distancing slow the spread?

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.

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)

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