Abstract
Susceptible-Infected-Recovered (SIR) models have long formed the basis for exploring epidemiological dynamics in a range of contexts, including infectious disease spread in human populations. Classic SIR models take a mean-field assumption, such that a susceptible individual has an equal chance of catching the disease from any infected individual in the population. In reality, spatial and social structure will drive most instances of disease transmission. Here we explore the impacts of including spatial structure in a simple SIR model. In particular we assume individuals live on a square lattice and that contacts can be ‘local’ (neighour-to-neighbour) or ‘global’ or a mix of the two. We combine an approximate mathematical model (using a pair approximation) and stochastic simulations to consider the impact of increasingly local interactions on the epidemic. We find that there is a strongly non-linear response, with small degrees of local interaction having little impact, but epidemics with susbtantially lower and later epidemics once interactions are predominantly local. We also show how intervention strategies to impose local interactions on a population must be introduced early if significant impacts are to be seen.
1. Introduction
The classic Susceptible-Infected-Recovered (SIR) model has long been used to model the spread of infectious disease in human, animal and plant populations (Kermack and McKendrick, 1927; Anderson and May, 1979). More recently it has formed a central pillar of much of the modelling of the Covid-19 pandemic (Ferguson et al, 2020; Kucharski et al, 2020; Firth et al, 2020). In its standard form, the SIR model has a mean-field assumption, such that individuals in the population have purely random, ‘global’ interactions (Boots and Sasaki, 2000) and there is no spatial structure. In reality, individuals in a population are more likely to contract disease from infected individuals who are closer to them, both physically and socially. Incorporating this spatial structure into mathematical models is extremely challenging. In some cases, large datasets of known contact networks have been used to replicate epidemics to excellent effect (Ferguson et al, 2020; Firth et al, 2020). While such models have a high degree of realism and thus predictive power, they cannot be readily modelled by a simple set of equations and require significant computational exploration to capture possible outcomes and feedbacks.
One common approach to incorporating a degree of regular spatial structure, and particularly ‘local’ near-neighbour interactions, in to infectious disease models is to use a lattice-based probabilistic cellular automata (Sato et al, 1994; Rand et al, 1995). These stochastic individual-based models have also been combined with an analytic pair-approximation method (Matsuda et al, 1992; Sato et al, 1994), where the full spatial dynamics are approximated by a set of ordinary differential equations based on the classic SIR model. Such models have been applied to infectious disease systems both with (Keeling et al, 1997; Webb et al, 2007a,b; Best et al, 2012) and without (Keeling, 1999; Sharkey, 2008) demography. These have found that local interactions reduce the value of R0, slowing or even preventing an epidemic that would occur when interactions are global (Keeling, 1999). These approaches largely insist on a strict degree of spatial structure, where infection and/or host reproduction can only be through near-neighbour interactions. While this is useful for comparison with the mean-field case, interactions are unlikely to be entirely local or global in reality, and we may be missing important features of systems where the interaction structure lies between these two extremes.
The ability to move between local, near-neighbour interactions and global, mean-field interactions has been considered in a few spatial models of infectious disease, primarily in evolutionary (Boots and Sasaki, 1999, 2000; Kamo et al, 2007; Best et al, 2011; D’ebarre et al, 2012) and ecological (Ellner, 2001; Webb et al, 2007a) contexts. This multiscale method is commonly achieved by allowing a proportion of transmission and/or reproduction to occur locally and the rest globally. We might interpret this, for example in a human population, as an individual mostly interacting within their household or community (local interactions), but also travelling some distance for work, holidays or visiting friends or family (global interactions). These studies have shown that there is increased potential for ecological cycles and disease-driven extinction as interactions become predominantly local (Webb et al, 2007a), while evolutionary selection is generally towards lower levels of infection in both host and parasite as interactions become more local (Boots and Sasaki, 1999; Best et al, 2011), but not necessarily monotonically (Kamo et al, 2007). Most recently, this multiscale method has been applied to a human epidemiology model with equal births and deaths (Maltz and Fabricius, 2016), showing that pronounced (but damped) oscillations in infection may result after a sudden shift to local interactions. However, this simple mechanism to investigate the impacts of varying the degree spatial structure has yet to be applied to simple human epidemic models over short-term scales such that demography does not impact the dynamics.
In this study we present a combination of a stochastic individual-based model and a pair approximation of epidemics on a lattice. We explore how changing the proportion of local-to-global interactions alters the course of an epidemic and investigate whether increasing the degree of local interactions - which we may define as restrictions on movement - can lessen the impact of an epidemic.
2. Model
2.1. Mean-field model
The underlying dynamics of the model are based on the classic Susceptible-Infected-Recovered (SIR) epidemiological framework (Kermack and McKendrick, 1927), with no demographic processes (births/deaths). We first consider the model under a mean-field assumption with no local interactions. All individuals in the population are either susceptible (S), infected (I) or recovered (R). The total population size N = S + I + R is constant (assume N = 1 for consistency with what follows), meaning we only need to track the dynamcis of S and I densities, given by the following ordinary differential equations,
Transmission is assumed to be density-dependent with coefficient β, while recovery occurs at rate γ and immunity is assumed to be permanent.
2.2. Pair-approximation model
To account for spatial structure and local transmission, we use a pair-approximation model (Matsuda et al, 1992). Assume individuals live on a square lattice, where each site is always occupied by one susceptible, infected or recovered individual. We define the probability that a site is occupied by a susceptible individual as PS, an infected individual as PI and a recovered individual as PR. The dynamics of these ‘singlet’ densities mirror those of the mean-field model above, with the following ordinary differential equations, with PR = 1 − PS − PI. Here we have introduced our key parameter, L, which determines the proportion of transmsision that occurs ‘locally’ between neighbouring individuals, with the remainder of transmission (1 − L) occurring ‘globally’ between random individuals on the lattice. This corrersponds to individuals’ interactions being predomiantly local (with their near neighbours) or global (randomly across the population). The conditional probability, called the ‘environs density’, that an infected individual has a neighbour that is susceptible is denoted qS/I = PSI /PI. Therefore there are two routes to transmission:
global: (1 − L) βPSPI
local: LβqS/IPI.
This system of equations is not closed, since to calculate the conditional probabilities we need to know the ‘pair’ density, PSI, e.g. the probability that a randomly chosen pair of neighbouring sites are a susceptible and an infected. The dynamics of these pair densities are governed by an additional set of ordinary differential equations, and PRR = 1 − PSS − PII − 2PSI − 2PSR − 2PIR. Again, this system of equations is not closed as we have further conditional probabilities that depend on ‘triplets’ (e.g. qI/SI = PSII /PSI). One can appreciate that this pattern will continue and that the equations will never form a closed system. We thus apply a ‘pair approximation’ (Matsuda et al, 1992) where we assume that, for example, qI/SI = qI/S, allowing us to close the system.
2.3. Basic reproductive ratio
The basic reproducive ratio, R0, is the well-known quantity that measures the average number of secondary infections caused by an infected individual in an otherwise disease-free population (Anderson and May, 1981). For the mean-field (global) case where L = 0, this is simply given by R0 = β/γ. When interactions are fully local with L = 1, we have R0,l = βqS/I /γ. In the limit where the population is indeed entirely disease-free, the environs density qS/I = PS = 1, and the two basic reproductive ratios will be equal. However, in the early stages of an epidemic the environs density qS/I rapidly shrinks as the contact network is formed, meaning it quickly becomes that R0,l < R0, leading to a slower epidemic (Matsuda et al, 1992; Keeling, 1999). Given the total reproductive ratio will be, it is clear that the initial growth rate of an epidemic will be slower the greater the degree of local interactions.
2.4. Stochastic simulations
Alongside these mathematical models we additionally conduct stochastic individual-based simulations using a probabilistic cellular automata. Similarly to the model described above, a lattice of sites is established, now of fixed size, where each site is again occupied by one individual. A Gillespie algorithm (Gillespie, 1977) is implemented for tau-leaping between events of recovery and transmission (local or global). At each step, exactly one of these events occurs, with probabilities proportional to their rates, and a suitable host is chosen randomly from the lattice for it to occur to (e.g. recovery requires an infected host to be selected). After an event occurs, the lattice is updated and a new tau-leap calculated for the next event. This approach is fully spatially explicit, unlike the approximation present in the mathematical methods above. Code for the models are provided as electronic supplementary material.
3. Results
3.1. Epidemic curves
We begin by examining the epidemic curves predicted by the pair approximation and stochastic simulations for different values of L (0.1 and 0.9) and different mean-field basic reproductive ratios, R0 (2, 5 and 10). Recent work has highlighted the pitfalls of combining mutliple stochastic individual-based models into simple static statsitics of means and variances (Juul et al, 2020). We follow the methods of Juul et al (2020) by finding the ‘most central’ 50% of 100 simulated curves to present here (see appendix for details).
Focussing on the effect of increasing the proportion of local interactions, from figure 1 it is clear viusally that the higher value of L produces a lower and later peak of infection. Restricting global interactions may therefore, in itself (without further reductions to transmission probability), slow down and limit the spread of an epidemic. Increasing R0 not only moves the epidemics earlier and higher, but also reduces the effect of local interactions. Comparing the plots, we can see that control mechanisms that both shift interactions from predominatly global to predominatly local and reduce R0 (for example, through both movement restrictions and ‘social distancing’ and hygiene measures) are predicted to have a signfiicant effect on reducing the epidemic.
We can also compare the fit of the pair-approximation to the stochastic models. As we might expect, when L is small the pair approximation appears to present a reasonable ‘average’ of the stochastic model runs. As L becomes larger we find that, while the pair approximation often sits within the most central runs, for larger R0 at least, it tends to predict that the epidemic peak is rather earlier and higher than seen in most of the fully spatially-explicit model runs. The discrepancy between the pair approximation and stochastic simulations is most pronounced at low values of R0. In particular, in this case a number of the stochastic simulations produce ‘failed’ epidemics, as evidenced by the lower bound of the 50% central curves running close to 0. In the online appendix we show that for L = 1 and R0 = 2 around 30% of stochastic simulations do not result in an epidemic.
3.2. Descriptive statistics
We now explore the behaviour as we vary local interactions across the full range of L from 0 (fully global) to 1 (fully local). Five descriptive statistics were evaluated, with three presented here (percentage of the population infected by day 300, percentage of the population infected at the peak and the day of the peak) with two further statistics in the online appendix (days till less than 1% of the population were infected and days with greater than 15% of the population infected). In order to take results from the individual-based simulations, 100 runs of the model were created and the mean and variance of the aforementioned statistics were taken.
Two clear trends emerge from all of the results. Firstly, the impact of local interactions is significantly reduced the higher R0 is. For every statistic investigated, varying the value of L has little effect on the plots for R0 = 10. Secondly, there is an accelerating impact of local interactions, with little effect seen as L is first increased from 0, but the impact growing as L moves towards 1. Both effects are clear when plotting the percentage of the population infected by day 300 (broadly, the total infected during an epidemic in this case) in figure 2(a)-(c). When R0 = 2, there is a slow decrease in the percentage infected for 0 < L < 0.6, but once L is greater than this there is a large accelerating decrease in the number of individuals infected during the epidemic. However for larger values of R0, there is a smaller decrease which doesn’t happen until L is almost 1. A similar pattern can be seen in plots (d)-(f): increasing L initially causes little change in the proportion infected at the peak of the epidemic, but then decreases more rapidly as L approaches 1, but this is less evident the higher R0 is.
Figure 2(g)-(i) shows that the number of days until the peak increases with L, again accelerating as L increases. There is an exception to this when R0 = 2 as L approaches 1. Here, the peak moves significantly earlier because the infection fails to spread through the population meaning the peak of the epidemic is both very early on and very low, as confirmed in figure 2a. Obviously, the larger R0 is, the faster the disease will be able to spread through the population and therefore the faster it will die out (with no susceptible individuals left to infect).
In general, the pair approximation appears to be a good fit to the results from the stochastic model and is almost always within a standard deviation of the mean, but this fit appears to be least good as L approaches 1. The pair approximation is less accurate for R0 = 2 than for higher values of R0, and this is likely due to the large proportion of infections which fail to become established in the stochastic model when the disease spreads slowly, resulting in a lower mean and larger standard deviation, as described in the online appendix.
3.3. Using local interactions as a control mechanism
We now explore how enforcing movement restrictions, resulting in more localised interactions, might impact the spread of an epidemic. We assume that initially a population has predominantly global interactions (L = 0.1). We then assume that when a threshold of percentage infected (here, 5%) is reached, interactions immediately switch to being predominantly local (L = 0.9) and remain so until the infected percentage returns below the threshold. Figure 3 shows that compared to the case where interactions remain predominantly global throughout (red), if movement restrictions are imposed (blue) the peak of the epidemic is reduced, but less substantially than if interactions had always been predmoniantly local, particularly for the lower R0 (see figure 1 and table 1). Interestingly, we also see a second wave emerging for lower R0 once restrictions are lifted since the herd-immunity threshold has not been reached, suggesting further and/or longer restrictions may need to be imposed.
We further investigate by varying the threshold at which restrictions are imposed and the value of L moved to under the restrictions (figure 4). Changing the bound at which L increases seems to have relatively little effect on the course of the epidemic, with there being little change in the the number of people infected, but a modest decrease in the peak and a somewhat later peak for lower thresholds. When varying L, again there is little change in the total proportion infected, a modest decrease in the peak for higher restrictions but almost no change to the peak day. The variance in the results from the stochastic simulations is large, suggesting that it may be more difficult to predict the outcome of a disease once restrictions on global interactions are implemented, but the mean of these results is close to the pair approximation. When the higher L = 1, the PA results change dramatically, due to the emergence of the second peak. As the higher L increases towards 0.9, it can be seen that the impact of the epidemic is mitigated, with the biggest change seen in the peak infected proportion of the population almost halving from L = 0.1 to L = 0.9.
This relative lack of impact is because of the speed with which the lattice becomes correlated in the early stages of an epidemic. The correlation between S and I sites on the lattice is given by,
At the start of an epidemic with predominantly global interactions then the lattice is uncorrelated and CSI = 1. During the early stages, the correlation rapidly approaches a quasi-equilibrium as the contact network forms (Keeling, 1999), which we show in the appendix can be approximated as, Figure 5 shows that increasing L leads to much stronger early-time S-I correlation, due to increasingly spatially-localised contact networks. If an epidemic begins in a population with predominatly local interactions, the lattice quickly becomes correlated, qS/I falls and the infection slows itself down due to a lack of locally available susceptible individuals. In contrast, if an epidemic has established with predominatly global interactions, the network is already highly uncorrelated before the movement restrictions are imposed. The late implementation of local interactions therefore cannot cause high correlation of the lattice, and a large number of local epidemics can still occur.
4. Discussion
In this study we have used a pair approximation alongside stochastic simulations to investigate the impact of local interactions on an epidemic. Our results show that epidemics where interactions are predominantly local will result in fewer infections spread over a longer period of time than those where interactions are global, in line with previous studies (Keeling, 1999). Importantly, we find that the trends as we move from purely global to purely local interactions are not linear. Instead, our results consistenly show initially flat responses in various infection statistics as L is increased, with rapid changes as L approaches 1. This suggests that the course of an epidemic in a population with relatively high proportions of local interactions (even 50:50) will be roughly the same as an epidemic in a population with purely global interactions. Even at relatively low proportions of global interactions, enough long-range infections can occur in the early stages of an epidemic to seed large numbers of local epidemics, allowing the infection to spread throughout the population. For example, if R0 = 2 and L = 0.5, on average an infected individual passes the disease to one local and one global contact, allowing the disease to become established across the lattice and to then form a series of outbreaks. It is only as L becomes close to 1 and almost all interactions are local that the likelihood that an infected individual transmits the disease globally is small enough to have a significant impact. Interestingly, in the similar model by Maltz and Fabricius (2016) that includes simple demographics (and thus yields an endemic equilbrium), the infected equilibrium is initially fairly static as interactions become more local before rapidly falling as local interactions become more dominant, suggesting this non-linear trend is robust in simple epidemic models.
Our results have important implications for attempting to limit an epidemic through restricting movement. In particular, such restrictions must be considerable, with almost all global interactions removed, if significant effects are to be seen. It is important to note that in our model restricting movement does not lead to lowered per-individual contacts, as might be assumed under ‘lockdown’ scenarios (for example due to social distancing, regular hand-washing, wearing masks, etc). We found that restrictions that both make interactions more local and infectious contacts less frequent (through lowered R0) can substantially reduce the impact of an epidemic. Moreover, we found that if the population starts from a position of having predomninantly global interactions, movement restrictions must be imposed very early on in the course of an epidemic or they will have minimal effect. This is due to the fact that, if a disease has already begun to spread randomly through a population with global contacts, when restrictions are put in place there will already be large numbers of local outbreaks forming across the lattice. If an infection has a particularly high R0, and therefore rapid growth, it may be that infection is already too widespread for movement restrictions to take effect by the time public health officials realise an epidemic has begun. In this study we assumed a simple switch such that interactions returned to the default after the infected proportion fell back below the threshold. More realistic approaches might be to gradually ease restrictions or enact further restrictions in cases where a “2nd wave” emerges. In the similar study by Maltz and Fabricius (2016), they found a simple switch to a different proportion of local interactions led to pronounced (damped) oscillations and significant periodic outbreaks as the system was effectively moved such that it was no longer at its steady state. Further investigation in to the use of movement restrictions as a control mechanism is needed to explore the best strategies.
Combining mathematical analysis, using the pair approximation (Matsuda et al, 1992; Sato et al, 1994), and stochastic simulations has allowed us to explore the dynamics of our model in depth. Interestingly we found that the deterministic results from the pair approximation model provide a good ‘average’ of the dynamics from fully-spatial stochastic simulations. The weakest ‘fits’ were for low values of R0, where a proportion of simulations lead to failed epidemics, whereas the analytical model always assumes an outbreak occurs. Given the problems in accurately depicting ‘averages’ of stochastic simulations (Juul et al, 2020), such analytic approximations may provide a useful guide.
We have deliberately focussed on the simplest possible epidemic model in this study, with the only two mechanisms being transmission and recovery. This has allowed us to draw clear conclusions and insight in to the behaviour of the model, but it clearly cannot and should not be used as an accurate predictive model for a particular epidemic. In an earlier study, Maltz and Fabricius (2016) considered the same model with simple demographics, finding that the infected equilibrium reduces with more local contacts, while (Webb et al, 2007a) examined the impact of varying local interactions on a fully ecological model, noting the potential for disease-induced extinctions and endemic cycles of disease. Clearly, however, there are many further elements that could be considered to make the model appropriate for specific infections or systems. A standard extension for many disease models is to add an exposed compartment, separating out those that are infected from those that are also infectious (see Keeling and Rohani, 2008). It may also be instructive to consider the dynamics if immunity to infection wanes over time, since the non-spatial model would then yield an endemic equilibrium, unlike our model. If we wish to consider a disease persisting over the long-term, we should not only add demographics but also consider seasonal-forcing (Aron and Schartz, 1984; Schwartz, 1985; Altizer et al, 2006). Finally, more realistic spatial and social networks would be needed for any conclusions around movement/interaction restrictions in specific circumstances to be considered, such as in recent models of Covid-19 in the UK (Ferguson et al, 2020; Kucharski et al, 2020; Firth et al, 2020). As it is, our model suggests that significant movement restrictions may be a useful strategy in tackling an epidemic.
Data Availability
N/A
Acknowledgements
LW received an Undergraduate Research Internship stipend from The School of Mathematics and Statistics at Sheffield University.