|
Comparison of alternative models of human movement and the spread of disease
Abstract
Predictive models for the spatial spread of infectious diseases has received much attention in recent years as tools for the management of infectious diseas outbreaks. Prominently, various versions of the so-called gravity model, borrowed from transportation theory, have been used. However, the original literature suggests that the model has some potential misspecifications inasmuch as it fails to capture higher-order interactions among population centers. The fields of economics, geography and network sciences holds alternative formulations for the spatial coupling within and among conurbations. These includes Stouffer's rank model, Fotheringham's competing destinations model and the radiation model of Simini et al. Since the spread of infectious disease reflects mobility through the filter of age-specific susceptibility and infectivity and since, moreover, disease may alter spatial behavior, it is essential to confront with epidemiological data on spread. To study their relative merit we, accordingly, fit variants of these models to the uniquely detailed dataset of prevaccination measles in the 954 cities and towns of England and Wales over the years 1944-65 and compare them using a consistent likelihood framework. We find that while the gravity model is a reasonable first approximation, both Stouffer's rank model, an extended version of the radiation model and the Fotheringham competing destinations model provide significantly better fits, Stouffer's model being the best. Through a new method of spatially disaggregated likelihoods we identify areas of relatively poorer fit, and show that it is indeed in densely-populated conurbations that higher order spatial interactions are most important. Our main conclusion is that it is premature to narrow in on a single class of models for predicting spatial spread of infectious disease. The supplemental materials contain all code for reproducing the results and applying the methods to other data sets.
The ability to predict how infectious disease will spread is of great importance in the 1 face of the numerous emergent and re-emergent pathogens that currently threatening 2 : bioRxiv preprint human well-being. We identified a variety of alternative models that predict human 3 mobility as as a function of population distribution across a landscape. These consider 4 some models that account for pair-wise interactions between population centers, as well 5 as some that allow for higher-order interactions. We trained the models using a 6 uniquely rich spatiotemporal data set on pre-vaccination measles in England and wales 7 (1944-65), which comprises more than a million records from 954 cities and towns.
Likelihood rankings of the different models reveal strong evidence for higher-order 9 interactions in the form of competition among cities as destinations for travelers and, 10 thus, dilution of spatial transmission. The currently most commonly used so-called 11 'gravity' models were far from the best in capturing spatial disease dynamics.
12 42 movement patterns [13] including a recent incarnation as the 'radiation' model [14].
While these models have proved useful in their original applications, their use in 44 describing infectious disease spread is complicated by two considerations. First, spatial 45 interaction models predict bulk migration fluxes between population centers, but the 46 impact of these movements on the dynamics of infectious disease can depend strongly 47 not only on the magnitude, but also on the composition, of the migrant pool. In 48 particular, the extent to which migrants are more or less likely to be susceptible to 49 infection, or actually infected, than the general population can be critically important, 50 as can the age profile of the migrant pool, due to pronounced age-specific patterns of 51 behavior, susceptibility, and infectiousness [5, 7, 15] . Second, infectious disease dynamics 52 and Wales using a spatial hazards approach [18] . We considered four main classes of 268 models, including the gravity, competing destinations, Stouffer rank, and radiation 269 models, in addition to pure diffusion and mean field models which are interesting special 270
Accurately predicting the geographical spread of emerging, re-emerging, and recurrent 14 epidemics in an increasingly globalized world is a matter of international urgency in the 15 wake of outbreaks of emerging pathogens such as severe acute respiratory syndrome 16 (SARS), Middle East respiratory syndrome (MERS), and ebola virus disease (EBVD) as 17 well as re-emerging and resurgent pathogens such as influenza subtype A-H1N1, 18 whooping cough, and measles. In response to this challenge, various 'spatial interaction' 19 models describing human movement as a function of population distribution have been 20 proposed. Some of these borrow from economics and human geography while others 21 adapt models from movement ecology and the physics of reaction and diffusion on 22 heterogenous landscapes. 23 In recent years, the field has seen the widespread adoption of a family of so-called 24 'gravity' models from transportation theory and human geography (see, e.g., [1, 2] ). In 25 its most common form, the gravity model posits that the migration flux between a pair 26 of cities is log-linearly dependent on their respective sizes and on their separating 27 distance. The application of this simple model to disease spread was originally proposed 28 by Murray and Cliff [3] , but over the last decade, many studies have used it to explain 29 historical, or predict future, disease spread in a range of infections including measles [4] , 30 influenza [5] [6] [7] , cholera [8] , and yellow fever [9] . While its application has yielded 31 insights, the geography literature has highlighted a prominent shortcoming of the 32 gravity model. In particular, these models ignore the potential for competitive or 33 synergistic interactions among population centers ( Fig. 1 ; [10, 11] ). Thus, for example, 34 movement between Boston and Washington DC is assumed to be unaffected by the 35 presence of the intervening city of New York. 36 As it happens, there are distinct families of models, arising from economics and 37 geography, that predict human movement patterns while allowing for higher-order 38 interactions among cities. In particular, Stouffer's [12] 'law of intervening opportunities' 39 posits that "the number of persons going a given distance is directly proportional to the 40 number of opportunities at that distance and inversely proportional to the number of 41 intervening opportunities". This idea has given rise to alternative models for human can feed back onto migration. Most obviously, disease symptoms can influence 53 movement behavior, as seen, for example, in the fact that a mild cold may induce 54 minimal changes in movement behavior but severe hemorrhagic fever or acute paralytic 55 disease will typically slow the movement of the infected hosts. These two considerations 56 do not preclude the utility of spatial interaction models in the infectious disease context; 57 they do, however, complicate their use and the importance of understanding 58 transmission relevant migration fluxes, which will be some sort of 'effective average' of 59 movement as filtered through such aforementioned complications.
Several spatial interaction models have been parameterized and tested using various 61 mobility data such as commuter flows (e.g., [5, 14] ), mobile phone geolocations 62 (e.g., [16] ), social media (e.g., [13] ), and microsimulations (e.g., [17] ). However, in view 63 of the challenges just noted, the ultimate test of the models is against data on the actual 64 spread of infection rather than bulk movement of people or cell phones. Bjørnstad and 65 Grenfell [18] proposed that for acute immunizing infections, the spatio-temporal 66 patterns of fade-outs (i.e., local disease extinction) across metapopulations provide 67 valuable information on disease spread because these patterns reflect spatial 68 transmission unclouded by local transmission. In this paper, we extract information 69 from fade-out patterns and a consistent likelihood framework to compare and contrast 70 to a suite of models including (i) the gravity model [1] , (ii) Fotheringham's competing 71 destinations model [10] , (iii) Stouffer's rank model [12] , and (iv) the radiation model [14] . 72 . We confront these models weekly data on measles incidence from all 954 cities and 73 towns in England and Wales from 1944 to 1965 [4, 19] . Comparing fits and predictions, 74 we show that while the gravity model is a reasonable first approximation, Stouffer's 75 rank model, an extended version of the radiation model, and the competing destinations 76 model all provide significantly better fits, Stouffer's model performing the best.
Data 79 Historical incidence of measles in England and Wales has been an influential testbed for 80 models and methods in disease dynamics since Bartlett's [20, 21] seminal work on its 81 recurrent epidemics. We use the spatially resolved weekly measles data across all 954 82 cities and towns of England and Wales from 1944, when notification was made 83 mandatory by the UK Registrar General (OPCS), until 1965, which saw subtle shifts in 84 political boundaries around London. Vaccination was not introduced in the UK until 85 1967, so that these data span a period where measles dynamics were unaffected by mass 86 immunization. The data set is complete except for a region-wide underreporting rate of 87 around 50% [22, 23] . Grenfell et al. [19] give a detailed description of the data; the 88 entire data set has been made available by Lau et al. [24] . 89 An important feature of the system is that, between 29% and 38% of the population 90 (c. 47M during this period) resided in a small number (13-28 depending on exact 91 definition) of communities above a critical community size (CCS) of c. 250-300k. Cities 92 larger than the CCS tend to sustain local chains of transmission. The remaining 60-70% 93 were distributed among the more than 900 communities smaller than the CCS where 94 local extinctions are more or less frequent (depending on population size and degree of 95 isolation) and, consequently, the rate of reintroduction of the pathogen via spatial 96 transmission is an important determinant of measles incidence. In ecological terms, the 97 prevaccination measles system represented an epidemic mainland-island metapopulation 98 (e.g., [25, 26] ). Our analysis exploits this fact, using the timing and spatial pattern and 99 timing of reintroductions to inform the parameters of each of the spatial interaction 100 models.
Using a spatially-extended time series susceptible-infected-recovered (TSIR) framework 103 (e.g., [4, 27] ), we can calculate, for each population center, the probability that a spatial 104 interaction happens (i.e., a contact between a resident susceptible host and a 105 non-resident infectious host) and that the contact results in a new local chain of 106 transmission [18] . In previous analyses, we showed that the signature of such events will 107 be drowned out in the presence of endemic circulation (see [23] ). Following a local 108 extinction, however, the re-colonization rate contains critical information on the spatial 109 interactions.
The spatially-extended TSIR model predicts that the expected number of new 111 infected hosts in community i in epidemic generation t + 1 will be
Here, β s(t) is the seasonal transmission rate as molded by opening and closing of schools 113 and the proportion of transmission that occurs within the school setting as the school 114 year progresses [23] ; s(t) = t mod 26 is the seasonality function; S i,t /N is the 115 probability that a local individual is susceptible; I i,t is the local number of infections; α 116 is a correction for the discrete-time approximation of the underlying continuous-time 117 process [28] ; and ι i,t is the sum of transmission-relevant spatial interactions i has with 118 the other communities in the epidemic metapopulation. Following a local extinction,
Assuming demographic stochasticity in transmission due to the underlying epidemic 121 birth-and-death process, the realized number of cases in the next generation will be [23] : 122
As a consequence, the number of susceptibles in the next generation will be
where B i,t is the recruitment of local susceptibles through births during the generation 124 interval. The 'serial interval' for measles is 10-14 days [29] , so we follow previous TSIR 125 analyses in aggregating the weekly data in 2-week increments (though [23] investigate 126 shorter intervals). Using susceptible reconstruction methods [30] , we estimate all local 127 parameters using the method of Finkenstadt et al. [31] . The online supplement contains 128 full documentation of these analyses and code for exact replication of all results 129 presented in this paper.
Conditioning on the local parameters and restricting analyses at each location to 131 biweeks in which measles incidence was zero, we use Eqs. 2 and 4 to form a negative 132 binomial likelihood to estimate all parameters of each of the candidate spatial 133 interaction models.
Spatial interaction models 135 When including special cases, such as pure diffusion and variants, we consider a total of 136 10 spatial interaction models (Fig. 1) . Each of these amounts to a specification of the ι 137 term in Eqs. 1 and 2.
Gravity model. Under the gravity model, the spatial interaction between locations i 139 and j takes the form θN τ1 i N τ2 j d −ρ ij , where θ, τ 1 , τ 2 , and ρ are non-negative parameters. 140 This leads to the following formulation for disease-relevant spatial interactions:
where the sum is across all non-self potential donor communities, and I j /N j is the 142 fraction of infected individuals in donor community j. The gravity model has two 143 important special cases: ρ = 0, τ 1 = τ 2 = 1, which is a mean-field model, and 144 τ 1 = τ 2 = 0, which is simple spatial diffusion. The φ parameter represents background 145 spatial transmission that is not predictable on the basis of distance and size [24] . 146 Xia's model. In the original analysis of the spatiotemporal dynamics of measles, Xia et al. [4] used a formulation rooted in the gravity-model literature but with a slightly different formulation:
which is the same as Eq. 5 if τ 2 = 1, but not otherwise.
Competing destinations model. Fotheringham [11] noted that gravity models may be misspecified because they only consider pairwise interactions among locations. He argued that nearby or intervening destinations may, in fact, make spatial interactions either more or less likely. A synergistic effect occurs, for example, if individuals are disproportionately inclined to do all their shopping in districts with many shops. Antagonistic effects may arise where travel between two cities, say Boston and Washington, is made less likely by the presence of an intervening city such as New York because individuals from either Boston or Washington can get their out-of-town needs fulfilled by visiting New York. Under Fotheringham's [11] competing destinations formulation the flux between i and j takes the form θ N τ1
As in the gravity model, τ 1 and τ 2 control how 'eagerness to travel' and city 'attractiveness' scales with population size; and ρ measures how the likelihood of travel decays with distance. The parameter δ quantifies how destinations k, of various sizes, at distances d jk from the donor j, modulate the spatial interaction between recipient i and and donor j. In particular, δ > 0 indicates a synergistic effect; δ < 0, an antagonistic effect. The resultant formulation for disease-relevant spatial interactions between community i and everywhere else is:
where φ again represents background, spatially unpredictable spread. 148 Stouffer's rank model. Noulas et al. [13] recently proposed that human mobility 149 patterns should be studied using Stouffer's [12] rank ("law of intervening author/funder. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available for use under
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.19.882175 doi: bioRxiv preprint population size as a proxy for "opportunities". Accordingly, we let Ω(i, j) be the 155 collection of towns located closer to town i than is town j:
156 The Stouffer model for disease-relevant spatial 157 interactions is, then
Stouffer's original model was framed in the context of continuous distribution of 159 population across a landscape. In the present metapopulation context, we therefore 164
Extended radiation model. The radiation model proposed by Simini et al. [14] was derived independently but is in spirit related to the original ideas of Stouffer [12] . Our metapopulation version of this model with respect to disease-relevant spatial interactions is
As with the Stouffer model, we consider the set Ω(i, j) to either exclude ('radiation") or 166 include ('radiation variant") community size i itself in the spatial interaction
Of the models considered, the gravity variants (including mean-field and pure 174 diffusion) include purely pairwise interactions whilst the other models allow for 175 higher-order interactions either explicitly (as in the competing destination model) or 176 implicitly (as in the Stouffer and radiation models).
Statistical inference 178 We use maximum likelihood to estimate the parameters of each spatial interaction 179 model and profile likelihoods to study correlation among parameters and possible 180 identifiability issues [27] . author/funder. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available for use under
The copyright holder for this preprint (which was not peer-reviewed) is the . Table 1 . MLE estimates. Models: SV -Stouffer variant, XR -extended radiation, CD -competing destinations, S -Stouffer, G -gravity, X -Xia, RV -radiation variant, R -radiation, MF -mean field, D -diffusion. Here, denotes the log likelihood,ĉ is the computed variance inflation factor, ∆QAIC is the QAIC relative to that of the best-fitting model, and the other parameters are as in the model equations. biweek of local disease absence [18] ) All calculations are documented in the online 189 supplement, which includes code for the detailed replication of all results.
190 Table 1 shows estimated parameters and likelihoods for each of the spatial interaction networks. To explore the drivers among model differences in fit, we therefore contrast 214 model predictions of disease import and export rates. Fig. 2 shows the mean predicted 215 export and import rates vs. population size for each of the models. The dependence of 216 mean export rate on population size contrasts somewhat among models ( Fig. 2A) , but 217 the importation rates follow similar patterns (Fig. 2B) . Despite their differences in 218 underlying mathematical equations, there must be more nuanced differences among the 219 models that contribute their relative fit/lack-of-fit. This begs the question: "What are 220
December 16, 2019 7/15 a CC0 license.
author/funder. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available for use under
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.19.882175 doi: bioRxiv preprint these differences and are they biologically meaningful or spatially random?" 221 A cursory review of the recent literature suggests that the gravity model is the most 222 prominent in the infectious disease context, followed by the radiation model. Stouffer's 223 model and the competing destinations model have rarely been applied in this field. It is 224 therefore useful to use the gravity model fit as a baseline and explore how the 225 better-fitting models diverge from this baseline. Comparing matrices with half a million 226 entries is very difficult, so we amploy a new 'spatial likelihood contrast' (SLiC) method 227 with which to study the relative merit of the different spatial interaction models. The 228 idea is to disaggregate the overall hazard likelihoods by individual cities to study how 229 particular communities contribute to improved or diminished relative fit of each model. 230 To do so we normalize each location's contribution to the overall likelihood by the 231 number of data points each conurbation contributes to the likelihood (the accumulated 232 stretches of measles absence) and then map the model-model differences onto the 233 landscape (Fig. 3) . Inspection of Fig. 3 suggests that the poorer fit of the gravity model 234 vs. the Stouffer or competing destinations models is primarily in the urban northwest. 235 The supplement provides the full set of pairwise SLiC contrasts among all models. To 236 test whether the apparent patterns are statistically significant, we compute local 237 indicators of spatial association (LISA; [34] ) statistics for all population centers with 238 fewer than 50k inhabitants (Fig. 3A, B) . The main failure of the gravity model is in (e.g. [35, 36] ). In the case of epidemics playing out within a population center, models 254 closer to the mean field model (i.e., with uniform connectivity) often do remarkably well 255 (see [37] for a social network perspective on this). However, such models generally do 256 not capture regional spread of human infections because of the complicated patterns of 257 human movement (the third wave of the 2009 influenza pandemic across the USA, 258 perhaps, being an unusually diffusive counterexample [7] ) . To elucidate these patterns, 259 many empirical studies have been performed and a variety of abstractions have been 260 proposed. The task of choosing among these abstractions for help in help predicting 261 spread of infectious disease is complicated by the fact that (i) model formulations may 262 differ with respect to how well they actually reproduce aggregate mobility patterns, and 263 (ii) model fit may be shaped by the manner in which mobility is filtered by transmission 264 and behavior.
To advance the discussion of these issues, we considered a suite of candidate models 266 of spatial interactions and fit them to the pre-vaccination measles data from England cases. The gravity model-which has gained recent prominence in infectious disease As the gravity model has had a very rapid adoption in spatial disease epidemiology 284 during the last decade (though Noulas et al. [13] have noted the importance of 285 considering alternatives), we used this as a baseline. Despite its common usage it 286 significantly under-performs relative to the several other models. Combining SLiC and 287 LISA statistics we found that the greatest tension between the gravity model and the 288 preferred alternatives is in the Manchester-Liverpool-Leeds conurbation of northwestern 289 England, where many towns and villages commingle amongst several major cities. formulation) that spatial interaction, once discounted for higher-order interference, is 295 roughly proportional to the geometric mean of the community sizes. This curious 296 suggestion warrants future investigation.
In this study we have used actual disease incidence data, as opposed to raw 298 measurments of human movement, to further our understanding of the principles 299 governing how human hosts spread ifectious disease across a populated landscape. We 300 have found both strong evidence of higher-order interactions among population centers, 301 and that candidate models differ significantly in their ability to capture the empirical 302 patterns. We hope our findings will help stimulate a systematic, data-driven discussion 303 of the relative merit of alternative predictive models for the probable path for spatial 304 spread of infectious disease. Clearly, the prospects are good for good for further 305 refinements and improved parameterizations. Spatial interaction models predict the flux of human movements between population centers (cities, towns, villages) as a function of the distribution of the population on the landscape. In this diagram, the relative magnitude of the fluxes from a focal town to other population centers are represented by the widths of the arrows. In the widely-employed gravity models (A), interactions among cities is strictly pairwise. Thus the addition of a new town (B) has no effect on the movement to other towns. In Fotheringham's competing destinations model (C), however, competition or synergy among nearby communities can reduce or augment fluxes. Stouffer's model of intervening opportunities and the radiation model (D) posit that movement from one city to another is diminished by the presence of opportunities in communities more proximal to the source city. author/funder. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available for use under
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.19.882175 doi: bioRxiv preprint
|
