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.
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.
We present four main questions of interest that we will answer throughout this project.
- What fraction of the population is expected to become infected?
- How far does the infection spread beyond its point of origin?
- How long does the infection take to reach across the lattice?
- How does social distancing slow the spread?
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 \]
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:
Infectivity with rate \(\beta\)
Recovery with rate \(\gamma\)
Migration with rate \(m\)
Drift with rate \(\frac{1}{N}\)
We consider each of these forces in our discrete simulations.
Here we employ our previous equations to simulate changes in the infection frequency over space and time.
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)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 ="fraction infected")
We show infection spread over time on x-axis and \(I\) on y-axis.
time <- 1:num_gens
position <- 1:(L+1)
time_interval <-0:300 # choose a time interval to plot from infection freqs
image(time[time_interval],position,I_1d[time_interval,])
We see a gradual spread of the infection until an equilibrium where all demes have a large fraction of \(I\).
We’ve successfully replicated the one-dimensional Stepping Stone model and captured typical stochastic SIS dynamics.
We expand the one-dimensional Stepping Stone to two spatial dimensions. We now allow for migration along a vertical axis. The following schematic illustrates this:
A schematic of the two-dimensional Stepping Stone model by Kimura.
Each deme has a coordinate \((r,s)\). We can describe the changes in frequency from migration in the two-dimensional model as follows:
\[ I_{r,s}(t+1) - I_{r,s}(t)= m(f_{r-1,s} + f_{r+1,s} + f_{r,s+1} + f_{r,s-1} - 4f_{r,s}) \]
However, to avoid having to deal with all of the subscripts and loops, we compute the Laplacian of \(I\). The discrete Laplacian of a function \(\psi(x,y)\) (written \(\nabla^2 \psi(x,y)\)) is defined as the following:
\[ \nabla^2 \psi(x,y) = \frac{1}{h^2}(\psi_{x+h,y} + \psi_{x-h,y} + \psi_{x,y+h} + \psi_{x,y-h} - 4\psi_{x,y}) \]
Exactly the same as the form of the two-dimensional Stepping Stone! We use the R package spatialFil to directly compute the Laplacian on our matrix. When we call the applyFilter() function on \(I\), the calculation in our model equation is performed across all pairs of \((r,s)\) neighbor demes and \(dI\) is computed accordingly. All this is done in a computationally efficient, generalized fashion.
library("spatialfil")
# With laplacian
spatial_sis_2d <- function(num_gens,L,pop_size,beta,gamma,m,d){
f0 = 1/pop_size # initial freq
mid<-floor(L/2) # initial location
I <- array(data=0,dim = c(num_gens,L,L)) # intialize infection matrix
for(i in 1:(num_gens-1)){ # loop for time
if(sum(I[i,,] > 0)==0){ # check if extinct, if so reset
I[i,mid,mid] <- f0
}
# compute discrete laplace
laplace_I = applyFilter(x = I[i,,], kernel = convKernel(sigma = 1, k = 'laplacian'))
# compute diff eq
dI = beta*I[i,,]*(1-I[i,,]) - gamma*I[i,,] + (m-d)*laplace_I
p = I[i,,] + dI # update freqs
p[which(p>1)] = 1 # bound freqs
p[which(p<0)] = 0
# generate random draws for entire 2d spatial array
rand = rbinom(n = length(p),size = pop_size, prob = p)/pop_size
dim(rand) = c(L,L) # reshape 2d rng vector back into lattice shape
I[i+1,,] = rand # update infections
}
return(I)
}num_gens = 1e4
L = 50
pop_size = 1e4
beta = 0.2875
gamma = 1/8
m = 1e-1
d = 0
I_2d <- spatial_sis_2d(num_gens,L,pop_size,beta,gamma,m,d)mid <- floor(L/2)
plot(I_2d[0:1000,mid,mid],type="l", main ="I vs. time in 1 deme",
xlab="time (generations)",
ylab ="fraction infected")
position <- 1:L
time <- 130
image(position,position,I_2d[time,,])
# dir.create("gif") # make a new directory for images
# setwd("gif/") # set as wd
#
# position = 1:L # all positions
#
# for(i in 1:num_gens){ # loops for time
# png(file=paste("gif_",i,".png")) # take images
#
# image(position,position,I_2d[i,,]) # save images
# dev.off() # end
#
# }A .gif exemplifying the spatio-temporal dynamics of our model.
Here we see a problem with the discrete Laplacian function used in R. The spatial boundaries of our model are reflecting, i.e. once the infection spreads to the “walls” of the lattice, they reverberate back towards the center and effectively nullify spatial patterns. We will explore this further in the Results section. Still, this .gif demonstrates that the infection will not spread beyond \(L\), the size of the lattice. This answers our question (2)!
We’ve expanded the model to two-spatial dimensions and incorporated SIS dynamics!
We compute the equilibrium fraction of infected from our model analytically, and later verify with simulations.
Recall that we model changes in infection frequency according to the relationship:
\[ 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)[\nabla^2I_{r}] \]
We want to solve this equation for the average, or expected, infection frequency at equilibrium. Note that equilibrium implies \(I_r(t+1) = I_r(t)\). Lastly, recall that \(E[\nabla^2I_r] = 0\).
We derive the average \(E[I]\) using the mean of the binomial distribution and the Law of Total Expectation as follows:
\[ E_{r,t}[\frac{1}{N}Binom(N,I_r(t+1))] = \frac{1}{N} N \,E_{r,t}[I_r(t+1)] = E_{r,t}[I_r(t+1)]\\\ \,\\ E_{r,t}[I(t+1)] = E_{r,t}[I(t) + \beta I_r (1-I_r) - \gamma I_r + (m-d)[\nabla^2I_{r}]] \\ \, \\ E_t[I(t+1)] - E_t[I(t)] = \beta E_t[I](1-E_t[I]) - \gamma E_t[I] + 0 \\ \, \\ 0 = \beta E[I](1-E[I]) - \gamma E[I] \\ \, \\ \gamma E[I] = \beta E[I] - \beta E[I]^2 \\ \, \\ \gamma = \beta - \beta E[I]\\ \, \\ \gamma = \beta (1-E[I]) \\ \, \\ \frac{\gamma}{\beta} = 1-E[I]\\ \, \\ \therefore E[I] = 1-\frac{\gamma}{\beta} \] Which is of course equivalent to the equilibrium \(I\) in the typical SIS model covered in lecture.
We can demonstrate this relationship with our simulations. Let us run multiple replicates of the one-dimensional model and compare the calculated average \(I\) with the expected.
# Same parameters from before, but we will increase gamma
num_gens = 1e4
L = 50
pop_size = 1e4
m = 1e-1
d=0
beta = 0.2875
gamma_vec <- 1/(8:15) # increase gamma
expected_eq_I <- 1 - (gamma_vec)/beta # from eq
computed_eq_I <- rep(0,length(gamma_vec)) # hold means
i<-1
for(g in gamma_vec){ # loop for gamma
I_1d <- spatial_sis_1d(num_gens,L,pop_size,beta,gamma=g,m,d)
computed_eq_I[i] <- mean(I_1d) # store means
i <- i+1
}print("Average I from simulations:")
# [1] "Average I from simulations:"
print(computed_eq_I)
# [1] 0.5247752 0.5707780 0.6081212 0.6386061 0.6632534 0.6846534 0.7036204
# [8] 0.7196288
print("Average I from 1- gamma/beta:")
# [1] "Average I from 1- gamma/beta:"
print(expected_eq_I)
# [1] 0.5652174 0.6135266 0.6521739 0.6837945 0.7101449 0.7324415 0.7515528
# [8] 0.7681159## Plot computed vs. expected mean
plot(computed_eq_I, expected_eq_I, type="l", main ="Computed vs. Expected average I",
xlab="Average I from simulations",
ylab ="1-gamma/beta")
We confirm that the average fraction of infected from our simulations matches the expected equilibrium value of \(I\). There are some discrepancies likely from the stochasticity of our model and the burn-in time before the simulation reaches equilibrium.
We have answered our first question (1). Our results demonstrate that, on average, \(1-\frac{\beta}{\gamma}\) fraction of our population will be infected and remain infected at equilibrium. This is true regardless of spatial structure.
Here we compute the variance of infection frequencies at equilibrium.
eq_vals <- I_2d[dim(I_2d)[1],,]
comp_var <- sum((eq_vals - mean(eq_vals))^2) / (dim(I_2d)[2] * dim(I_2d)[3]) # compute variance
print(paste("Variance in infection frequency across spatial lattice is:",comp_var))
# [1] "Variance in infection frequency across spatial lattice is: 2.92553446704e-05"This quantifies our observation that after an early period of spread, the infection can essentially be found at all locations in the lattice equally. There is very little variance in infection presence across the lattice at equilibrium. This an artifact of the SIS model, where the infection persists at equilibrium. In this case, equilibrium is not of great biological relevance. We know that real epidemics burn out after enough time. Still, we can use our model to study the dynamics of the spread of the infection before we reach equilibrium. In the next section, we examine some of these non-equilibrium quantities derived from our model.
Recall the time series for the one-dimensional model in a single deme. We see that equilibrium occurs within the window plotted. Let us now quantify exactly how long the simulation takes to arrive at equilibrium.
num_gens = 1e4
L = 50
pop_size = 1e4
m = 1e-1
d=0
beta = 0.2875
gamma = 1/8
I_1d <- spatial_sis_1d(num_gens,L,pop_size,beta,gamma=g,m,d)
mid <- floor(L/2)
plot(I_1d[0:200,mid+10],type="l", main ="I vs. time in 1 deme",
xlab="time (generations)",
ylab ="fraction infected")
time_to_infection <- rep(0,L) #initialize vector
for(i in 2:L){
time_to_infection[i] <- num_gens - sum(I_1d[,i]>0)
}
plot(1:L, time_to_infection,
xlab="position", ylab="Time before first infection")
print(paste("Maximum time before first infection:",max(time_to_infection)))
# [1] "Maximum time before first infection: 104"We see that the time when the first infection occurs in the farthest deme corresponds to the equilibrium time. This describes the point at which the entire lattice is infected and spatial correlations begin to drop off. This answers our question (3)! We will use this max time to determine how the speed of infection spread changes as we limit the migration rate.
Here we explore the effect of \(d\), the “distancing” parameter that lowers the migration rate across the lattice.
num_gens = 5e2
L = 50
pop_size = 1e4
m = 1e-1
beta = 0.2875
gamma = 1/8
d_vec <- seq(0,m,by = 0.001)
max_time_to_infection <- rep(0,length(d_vec))
for(i in 1:length(d_vec)){
I_1d <- spatial_sis_1d(num_gens,L,pop_size,beta,gamma=g,m,d_vec[i])
time_to_infection <- rep(0,L)
for(j in 2:L){
time_to_infection[j] <- num_gens - sum(I_1d[,j]>0)
}
max_time_to_infection[i] <- max(time_to_infection)
}dist_effect<- d_vec/m # fraction of migration rate stopped by distancing
max_times <- max_time_to_infection/min(max_time_to_infection) # max time to infection in farthest deme scaled by case with no distancing
plot(dist_effect,max_times,
xlab= "fraction of migration rate stopped by distancing d/m",
ylab="Max time to infection / max time without distancing")
abline(h=2*min(max_times),col="red")
abline(v = 0.8, col="blue")
As we increase \(d\), we are decreasing the migration rate \(m\) and therefore slowing the spread of the infection across the lattice. To slow the spread of the infection by a factor of \(2\) we need to stop nearly \(80\%\) of all migration. However, the effect of distancing increases the time to full spread exponentially, so beyond the 0.8 threshold, we can triple to quintuple the time it takes the infection to reach the entire lattice. We’ve answered our final question (4)!
In this project, we have expanded the Stepping Stone model framework to incorporate SIS model dynamics. We expand the model to two spatial dimensions and visualize the spread of infection. Ultimately, we identified both how far and how quickly the infection spreads depending on our migration and distancing parameters. Our model is subject to several limitations. First of which is that we model discrete space and time. There are spatial SIS models that have operated with continuous space (Paeng and Lee 2017). Additionally our model considers parameters that are constant across space and ignores heterogeneity. Other work has been done to study heterogeneity in spatial disease models (Cui and Lou 2016), but ultimately our dynamics are sufficient for our scenario of interest.
Our model is most useful for representing fast-spreading epidemics in situations of no previous immunity and in locations with high enough population density to facilitate spread. Potential applications include capturing the spread of disease in urban areas. Our model also helps us consider social distancing as a quantitative force. For the COVID-19 parameters used, we illustrated the non-linear relationship between the speed of infection spread and the fraction of migration stopped by distancing. To even make a dent in the spread speed, the model recommends a near (80%) reduction in travel. However, as we increase beyond that amount of distancing, the payoff in terms of slower spread increases exponentially.
Ultimately, this model serves as a warning against the easing of distancing measures while an epidemic proceeds. Increasing migration can drastically speed up the spread of infection. Although, this model doesn’t tell us how a growing population of immune individuals impacts the spread. Other spatial models are more explicitly useful for evaluating the effects of immunity (Hayama et al. 2013).
Future directions include broadening the range of parameters simulated, expanding the model to represent time-varying or space-varying parameters, and quantifying the relationship between distancing and spread speed analytically.
Thank you to the Novembre Lab for the background in population genetics and thank you to Greg for leading a fantastic course!
Cui, Renhao, and Yuan Lou. 2016. “A Spatial SIS Model in Advective Heterogeneous Environments.” Journal of Differential Equations 261 (6): 3305–43. https://doi.org/10.1016/j.jde.2016.05.025.
Hamada, Miki, and Fugo Takasu. 2019. “Equilibrium Properties of the Spatial SIS Model as a Point Pattern Dynamics - How Is Infection Distributed over Space?” Journal of Theoretical Biology 468 (May): 12–26. https://doi.org/10.1016/j.jtbi.2019.02.005.
Hayama, Y., T. Yamamoto, S. Kobayashi, N. Muroga, and T. Tsutsui. 2013. “Mathematical Model of the 2010 Foot-and-Mouth Disease Epidemic in Japan and Evaluation of Control Measures.” Preventive Veterinary Medicine 112 (3): 183–93. https://doi.org/10.1016/j.prevetmed.2013.08.010.
Kao, Rowland R. 2003. “The Impact of Local Heterogeneity on Alternative Control Strategies for Foot-and-Mouth Disease.” Proceedings. Biological Sciences 270 (1533): 2557–64. https://doi.org/10.1098/rspb.2003.2546.
Keeling, Matt J., and Pejman Rohani. 2008. “Spatial Models.” In Modeling Infectious Diseases in Humans and Animals, 232–90. Princeton University Press. https://doi.org/10.2307/j.ctvcm4gk0.10.
Kimura, M., and G. H. Weiss. 1964. “The Stepping Stone Model of Population Structure and the Decrease of Genetic Correlation with Distance.” Genetics 49 (4): 561–76.
Kiss, Istvan Z, Darren M Green, and Rowland R Kao. 2005. “Disease Contact Tracing in Random and Clustered Networks.” Proceedings of the Royal Society B: Biological Sciences 272 (1570): 1407–14. https://doi.org/10.1098/rspb.2005.3092.
Paeng, Seong-Hun, and Jonggul Lee. 2017. “Continuous and Discrete SIR-Models with Spatial Distributions.” Journal of Mathematical Biology 74 (7): 1709–27. https://doi.org/10.1007/s00285-016-1071-8.
spatialfil was used to perform the Laplacian convolution for the expanded model. The package can be found here