Spatial dynamics of the 1918 influenza pandemic in England, Wales and the United States

There is still limited understanding of key determinants of spatial spread of influenza. The 1918 pandemic provides an opportunity to elucidate spatial determinants of spread on a large scale. To better characterize the spread of the 1918 major wave, we fitted a range of city-to-city transmission models to mortality data collected for 246 population centres in England and Wales and 47 cities in the US. Using a gravity model for city-to-city contacts, we explored the effect of population size and distance on the spread of disease and tested assumptions regarding density dependence in connectivity between cities. We employed Bayesian Markov Chain Monte Carlo methods to estimate parameters of the model for population, infectivity, distance and density dependence. We inferred the most likely transmission trees for both countries. For England and Wales, a model that estimated the degree of density dependence in connectivity between cities was preferable by deviance information criterion comparison. Early in the major wave, long distance infective interactions predominated, with local infection events more likely as the epidemic became widespread. For the US, with fewer more widely dispersed cities, statistical power was lacking to estimate population size dependence or the degree of density dependence, with the preferred model depending on distance only. We find that parameters estimated from the England and Wales dataset can be applied to the US data with no likelihood penalty.


INTRODUCTION
Spatially explicit models are critical to understanding the spread of infectious diseases through populations and to better inform policy aimed at controlling that spread. Indeed, recent outbreaks of communicable diseases in human populations have triggered a series of studies addressing the spread of directly transmissible infections at a country level [1 -5]. Identifying a possible backbone of high probability transmission paths through populations may underpin the development of effective interventions to curtail spread on the population network [6]. For example, in human diseases, spatial models and microsimulations can quantify the possible role of border control, quarantine or transport reductions in curtailing local and international spread [2 -4,6 -11]. Spatial microsimulation models like these are critical to making effective policy decisions. Spatial models also allow limited control resources to be used where they might be most effective. For instance, during the initial stages of an outbreak of a new virus, disease incidence tends to occur in spatial clusters with occasional long-range infection events [10]. As case numbers increase in the start location, the frequency of long-range infection events increases. This spatial pattern was seen during the early stages of the 2009 influenza pandemic in Mexico, leading to local foci seeded by long-range interactions to other countries [12,13].
However, those epidemic models rely on a set of structural assumptions that need to be validated from data. A basic assumption of many spatially explicit transmission models is that flows between urban centres are a function of the distance between them and their attributes, most notably residential or worker population sizes [14,15], resulting in a so-called gravity model. However, for human diseases, little work has been done to validate the underlying assumption that human travel patterns are predictive of the spatial spread of diseases.
Early models of spatial coupling in ecology assumed that connectivity between populations was inversely related to the distance between them [16]. For people (and most animals) distance-based coupling is too simplistic an assumption. Movement between large population centres is disproportionately more frequent than between smaller ones [17]. Xia et al. [18] found that a distance-only model for spread of measles in the UK was a poor fit to weekly measles data from England and Wales from 1944 to 1967. The gravity model used in that study measures connectivity between population centres as a function of distance and a function of the population sizes of the origin and destination cities. However, measles is a childhood disease, and so the spatial dependencies of the host are different than for infections that affect both adults and children, such as pandemic influenza. Gravity models and other spatial interaction models allow understanding of the movement of populations from one location to another in the absence of movement data. Models, once validated, can predict modifications in connectivity when populations grow or shrink, when workflows vary owing to economic changes or if restrictions are imposed on one city and not others. This is in contrast to movement surveys, which are context-specific and provide a snapshot of the movement habits of a population.
The strength of connection between cities may be density independent, that is, the sum of connectivity of a city to all its neighbours does not depend on the number of neighbours that city has. In contrast, density-dependent connectivity links two cities at a strength solely determined by the sizes of those cities and their distance apart, so that the total connectivity of any one city scales with the number of close neighbours. Density-independent transmission gives a total force of infection, which is independent of the remoteness of the population, whereas density-dependent transmission will cause populations with many neighbours to experience a higher force of infection than those cities that have few neighbours. Thus, a densitydependent model will predict that isolated populations are less likely to become infected than populations with many neighbours, or few very large neighbours. The concepts of density dependence/independence have not only been used to model interactions between cities; they have also been used extensively in individual based models of disease spread. Most past studies tend to assume either density dependence (e.g. for animal epidemics; [19,20]), or density independence (e.g. for most human diseases) [3,4]. Here, we explore the extent to which city-to-city (rather than individualto-individual) contacts are density independent by constructing a model that can capture intermediate levels of density dependence.
In this paper, we analyse mortality datasets from England and Wales and the United States from 1918 to 1919 to examine the pattern of spatio-temporal spread and the extent to which gravity models can reproduce observed trends. The 1918 pandemic constitutes a rare example of a well-documented epidemic in a largely susceptible human population, where the high mortality gives a clear incidence signal, and is therefore a rare opportunity to validate models of epidemiologically relevant geographic coupling. We examine the effect of city-specific characteristics (e.g. location, distance from other cities, population size and the number of influenza-related deaths) on the pattern of spread seen. We also investigate the impact of the distribution of cities in each country. The analysis provides further insight on the spatial variation in the spread of the 1918 influenza pandemic at a country level, much of which remained unexplained in past studies [18,[21][22][23].

Data
The 1918 pandemic H1N1 virus appears to have entered the general population of the UK and US in the spring of 1918 causing a reportedly mild disease [24]. This early wave was associated with increased mortality, but was probably only noticed because influenza is rare in summer. This epidemic waned later in the summer, but infection reappeared in the autumn with much increased mortality. It is unclear whether the viruses causing the spring and summer waves were closely related, but there is increasing evidence that the spring wave gave immunological protection against the autumn wave at a population level [25,26]. By September 1918, the pandemic was a prominent global phenomenon. The autumn wave was virtually universal, albeit with some variation between countries in precise timing. In both the UK and the US, a third wave of influenza occurred in early 1919 although with greater heterogeneity in mortality rates between cities [24,27,28]. US cities had more variation in the severity of the major wave than the UK, probably in part because some enacted more stringent non-pharmaceutical interventions to mitigate the epidemic [27,28]. The third wave was less pronounced in US cities than in England and Wales, again perhaps partly because of the effect of interventions.

England and Wales
The England and Wales dataset shown in figure 1a was published in the Supplement to the 81st Registrar General Report [29]. It provides weekly death counts and annualized mortality rates per 1000 from 83 county boroughs, 84 municipal boroughs, 71 urban districts and three unclassified urban centres in England and Wales, for a 46 week interval, 29 June 1918 -10 May 1919. The first 10 weeks are designated as wave 1, the next 19 as wave 2 and the last 17 as wave 3. Our analysis focuses on the second wave because it occurred in all cities (unlike wave 3) and because recording of mortality had begun in all cities before its arrival (unlike wave 1). In addition, reporting of influenza and influenza-related mortality changed between the first and second waves. Point locations of all urban centres were determined from the current or historical records, with Euclidean distance used to quantify inter-centre separation. Further information is given in the electronic supplementary material.

United States of America
We compiled a US city dataset (figure 1b) from five publications reporting the Weekly Health Index as collated by the Bureau of the Census [24,[30][31][32][33]. It covers the period 14 September 1918 to 15 March 1919 and contains weekly pneumonia and influenza death counts for 47 cities in the US. The US data therefore covers a period of two waves, with not all cities experiencing the later wave. There is very good agreement between different sources where they overlap. We used the Euclidean distance between cities (accounting for curvature of the Earth) to measure separation. Further information is given in the electronic supplementary material.

Reconstructing the week of infection
The analysis requires an estimate of when each city became infected to allow potential sources of that infection to be identified. For each city infected in week t, the candidate infectors are those infected in any week before t. We define the infection week of city i, t i , to be the first that meets a set of conditions on mortality in weeks t i þ 1, t i þ 2 and t i þ 3. We use mortality values ahead of t i to include the time from infection to death. A week could be designated the infection week if either (or both) of two sets of criteria were met for mortality in the following weeks. The first set of criteria required the mortality rate in week t i þ 1 to be above a certain threshold, to have increased in t i þ 2 and to be above a higher threshold in t i þ 3. These criteria are intended to ensure that the epidemic in that city is patently increasing. The second set of criteria was designed to capture cities where there was a rapid onset of increased influenza-related mortality. They therefore used a higher threshold on mortality in week t i þ 1, but less strict conditions on rate of growth in the following two weeks. The week of infection determined was found to be relatively robust to the precise choice of thresholds used. For further details on the algorithm and a spatio-temporal display of the result, see the electronic supplementary material.

Spatial models
In formulating our inter-city transmission model, we take the city as our unit of study. Each of N cities, i, has an infection time t i , an invariant population size P i and a time-varying mortality rate, r i,t at time t. Infected city i is separated from susceptible city j by distance d ij . Each week, each city can be in one of the three disease states: Susceptible, Latent or Infectious. We assume that all cities are susceptible at the start of a wave, they are latently infected for one week on infection, and that a city becomes infectious the week after it becomes infected. We assume that all transmission is endogenous to England and Wales or US after external seeding to the first infected city in each territory. If a city becomes infected in week t i , the candidate infectors are only those cities that are infectious in week t i . We assume that the transmission parameters are constant through time.
The model formulation aims to capture the effect of distance and population size on the connectivity of cities. Three modes of spatial transmission are considered: density-independent connectivity, densitydependent connectivity and an intermediate form where the degree of density dependence is estimated. The model also examines the different assumptions regarding a city's infectivity over time.
The force of infection, l is the hazard of infection from one city to another. From infected city i on susceptible city j at time t, it is: where source city population and distance are normalized together. n and m are estimated parameters on source and destination population sizes, respectively. d ij represents the distance between cities with power parameter g to be estimated. w is an estimated parameter relating the infectivity of a city to its mortality rate. When w ¼ 1, the infectiousness of a city at time t i is proportional to the death rate in that city at time t i þ 1. We use one week as a lag from infection to death [24,[34][35][36]. A value of w ¼ 0 gives a flat infectiousness profile, independent of the death rate in the source city. Intermediate values of w give variation in infectiousness, which scales sub-linearly with weekly mortality. Estimating w allows us to assess whether mortality rate is a good proxy for infectiousness in an infected city. b is a time-invariant estimated infectivity term.
Parameter 1 describes the strength of connection of a susceptible city to all possible infectors. 1 ¼ 0 gives the density-dependent model and 1 ¼ 1 gives the density independent model. By allowing 1 to vary, we allow the model to estimate the degree of density dependence in connectedness between the cities.
The total force of infection on city j at time t is given by: where Since the force of infection is a hazard, the probability that a susceptible city j is infected in a week t j is given by: t¼0 Àl j;t ! ð1 À expðÀl j;t j ÞÞ:

Estimation of parameters
We used a Bayesian framework for statistical inference.
The log-likelihood is given by: ln P t j À Á : We explore the joint posterior distribution of parameters by Markov Chain Monte Carlo (MCMC) sampling [37,38]. We sampled parameters on a logscale using a random walk update scheme. b and g were jointly updated, while 1, m and n were updated singly. Five MCMC chains were started from a variety of start points within a credible range to assess convergence. Convergence was achieved within 100 000 iterations for all models from all starting parameter values. For each model, the chain was run for 500 000 iterations including a burn in of 100 000. Parameter estimates and equal-tailed 95 per cent credible intervals were obtained from the posterior distribution of 80 000 values thinned from the last 400 000 samples of the MCMC chains.

Model variants
To investigate which components are most important for describing the spread of influenza, we consider a set of simplified variants of the model presented above. In those variants, each parameter can be either fixed at 0, at 1 or be estimated by MCMC. For a full comparison of each component of the model, see the electronic supplementary material.
The Deviance Information Criterion (DIC) is used to compare models [39]. This is calculated using the median parameter values owing to non-normality in the likelihood [39,40]. Lower values are preferable and a difference of around 5 units is considered important [41].

Epidemic trees
The model is used to generate epidemic trees [42,43]. We sample 1000 parameter sets from the joint posterior distribution and calculate the probability of infection for each potential infector city. The most likely tree for each parameter set is generated by calculating which 'infector' city has the highest probability of infecting each 'infectee' city. The distance to this infector, the probability of the infector-infectee pair and the number of infectees each infector creates are calculated for each parameter set. Mean values are weighted by the frequency of infector -infectee pairs from 1000 trees.

Validation
We examine the ability of the models to recreate the observed epidemic by simulation. We use 1000 parameter sets sampled from the joint posterior distribution and for each set, we simulate an epidemic using the first infected city as a source of infection. Once a city is infected, the observed mortality curve is used to model the infectiousness of that city through time.
We also calculate the probability distribution of the week of infection for each city conditional upon the observed epidemic up to that time point. We use 1000 parameter sets sampled from the joint posterior distribution and for each set, we calculate the probability of infection each week for each city given the epidemic observed up to that time point. We use Welch's two-tailed t-test to differentiate outlying groups.

Model components
Comparison of spatial and non-spatial population-independent models shows that inclusion of distance substantially improves model fit for both England and Wales and the US (D DIC ¼ 64.7 and D DIC ¼ 33.3, respectively). DIC and parameter estimates for the distance-only model are given in column 1 of tables 1 and 2.
Previous formulations of the gravity kernel in the literature have considered either density-dependent (1 ¼ 0) or density-independent transmission (1 ¼ 1). Figure 2 compares the fit (expressed by the posterior deviance) of these two formulations and with that from the model where the degree of density dependence, 1, is estimated. This comparison is made for models assuming no linear or a fitted power-dependence of spatial coupling on both source and destination city population size. In figure 2a, w is estimated, whereas in b it is fixed at 1. See model components in the electronic supplementary material for further comparisons.
For England and Wales (figure 2a), in each population context the variant that estimates the degree of density dependence (the lightest curve of each colour) gives a slightly better fit than models with no density dependence, with pure density dependence fitting substantially less well. The comparison also shows that the models which estimate the effect of origin and destination city population sizes on the connectivity of cities are much better than either the populationindependent or linear population size-dependent models.
The same set of comparisons is made for the US in figure 2b. The situation is more complex with the posterior distribution of many model variants lying in the same area. Comparisons by DIC value cannot distinguish these models. Unlike in England and Wales, there is no density-dependence variant which has lower deviance for all three of the population sizedependence variants examined. Inclusion of nonlinear population size-dependence does not penalize the fit of the US model, and so cannot be definitively excluded as being consistent with the data. The models presented in columns 5 and 6 in table 2 have different population relationships, but the same DIC score. The credible intervals on the population parameters of the densitydependent population with infectivity model (column 5) are very wide suggesting that little information is added by the inclusion of these parameters.
In England and Wales, the lowest DIC model is one where the degree of density dependence is estimated and the effect of population is also estimated. This is in contrast to the US where population-independent models either with density-dependent or estimated density dependence spatial interaction terms are indistinguishable.

Impact of infectivity profile
We tested models with three types of infectiousness profile through time: constant infectivity, a linear relationship between infectivity and mortality in the week ahead, and an estimated power-law relationship between mortality and infectivity. Mixing was poor when estimating w with the US data so we only compare the first two models in that setting.
In England and Wales, the linear infectivity model has a DIC value of more than 25 above either the constant or estimated infectivity model. Parameter estimates for the constant-infectivity and estimated-infectivity model variants are shown in columns 4 and 5 of table 1 using the density-dependent population-dependent framework from the previous comparison in England and Wales. These two models are indistinguishable by DIC (D DIC ¼ 0.1). Estimates for all other parameters are very comparable between these two models.
The estimated relationship includes two inputs from the infected city: the mortality rate and the population of the city. It can be more difficult to estimate parameters regarding infectivity, so we tested a model which takes only one piece of information from the infected city. The final column in table 1 shows a model which takes the mortality rate from the infected city into account but does not include the population size of that city. There is an improvement in the DIC score for this model of 4.4 over the constant infectivity model. Table 2 shows parameter estimates for models in the US. The difference in DIC score between a constant infectivity model and one with a linear relationship between mortality and infectivity is negligible in either a distance-only model framework (columns 1 and 6) or a population-dependent framework (columns 2 and 5). Adding infectivity information does not improve the fit of the model.
In the England and Wales dataset, the lowest DIC model is the single infected city parameter model in column 6 of table 1. The model is dependent on the destination population size, has estimated dependence of infectivity on mortality and an estimated intermediate degree of density dependence. The distance power g was estimated as 1.18 (0.96, 1.39). A lower value was found for models in the US, where for the most parsimonious low DIC model (density dependent, population independent), g was estimated as 0.79 (0.54, 1.00). Figure 3 shows the distance kernels for the two datasets. The credible intervals for g overlap for the two datasets.
The power parameter on the destination city population, m, was estimated at 0.40 (0.25, 0.54) in England and Wales. The credible intervals exclude 1, demonstrating that as population size increases, the susceptibility of the city increased more slowly.

Comparison between datasets
We used the posterior median parameter estimates fitted to the England and Wales dataset to calculate a likelihood value in the US dataset. By likelihood ratio test, this value was not different from the most parsimonious low DIC US model (255.23, 257.72, p . 0.97). We therefore cannot reject the assumption that spread had the same characteristics in the US and England and Wales, though clearly the smaller size of the US dataset reduces inferential power. the epidemic during which each city was infected. Inferred city-to-city infection events more frequent than 70 per cent (in 1000 trees) are shown in black, events of lower frequency are shown in grey. Interactions in weeks 0 -3 are longer range than those in weeks 4 -7 ( p , 0.01), which are in turn longer range than those in weeks 8 -10 ( p , 0.01) (figure 4d ). The probability that the most likely infector was responsible for each infection falls as the epidemic progresses because there are many more potential infectors available later (figure 4e). Cities infected early give rise to more infections than those infected late in the wave, as expected, but the range is large, with some early cities giving rise to no new infections (figure 4f ). Figure 5 shows results for two models using the US data. We compare the most likely infection trees for the distance-only constant infectivity model with parameters inferred from the US data (figure 5a) with a model where parameters used to generate the trees are taken from the England and Wales single-infected city parameter model (figure 5c). In the distance-only constant infectivity model, the nearest infected city is always the most likely infector. In contrast, with the England and Wales parameters, some links between cities are high frequency, while other cities have several potential infectors of intermediate frequency ( figure 5b). As in England and Wales, infection events inferred early in the epidemic have a higher support than those later in the epidemic. There are some exceptions owing to the distribution of cities in the US dataset-Oakland and San Francisco are distant from all other cities but very close to each other. In the distance-only constant infectivity model, some cities may give rise to a large number of new infections (e.g. Pittsburgh gives rise to nearly a quarter of infections) (figure 5d ). The effect of a city acting as a hub of infection is reduced in the more complex model, as the risk of infection from one city to another is the combined effect of several factors including distance.

Validation
For England and Wales, there is a relatively good agreement between observed and simulated epidemic curves (figure 6b). The observed epidemic curve rises more steeply than the simulation curves in the early stages of the epidemic, and peaks one week earlier than the simulation mean. This suggests that the model may underestimate the external infection pressure early in the wave. We calculated the probability that a city was infected in each week given the observed behaviour of all other cities up to that time. In England and Wales, 245 of 246 cities lie within the 95 per cent interval of their expected distribution. Figure 6a shows the cities which the observed infection week lies outside the stricter inter quartile interval. For further information see the electronic supplementary material. There are no population size ( p ¼ 0.36) or density trends ( p ¼ 0.11) in these cities, which are typically infected later in the epidemic ( p , 0.01 for difference in infection week). In the US, all cities lie within the 95 per cent probability interval and all but three lie within the inter quartile interval. Those three outlier cities are smaller than other cities ( p ¼ 0.01) but equally distributed in space ( p ¼ 0.88) and time ( p ¼ 0.06).
We have tested the effect on parameter estimates in England and Wales of relaxing the single-introduction assumption inherent in the model. We re-estimated the parameters conditioning on infections that occurred from week 3 of the epidemic onward. There is a small increase in the kernel power parameter estimate, which causes the kernel to decay more rapidly with distance (electronic supplementary material, figure S9). This suggests that the very long-range interactions, which are forced to occur early in the epidemic impact the shape of the kernel. However, the credible intervals largely overlap which indicates this assumption does not affect the fit of the model to a large degree.
In the US the simulated curves for the distance-only constant infectivity model are shown in figure 5e and for the England and Wales parameters in figure 5f. In both cases, the mean simulated and observed curve are very comparable, with the distance-only constant infectivity model giving peak incidence in the same week as observed. Figure 5g shows the observed week of infection against the simulated week of infection for all 1000 simulated epidemics. There is good correlation between the observed and simulated weeks of infection for both parametrizations.
We have tested the effect of thinning the England and Wales dataset so that it more closely resembles the US dataset to determine if the differences in formulation between the best models for each dataset are owing to the smaller number of cities in the US dataset. We removed all cities with fewer than 90 000 inhabitants in England and Wales leaving 46 cities distributed quite evenly in England and Wales as shown in the electronic supplementary material, figure S10. There were identifiability problems in estimating the density-dependence parameter 1 using the thinned dataset. The best model by DIC comparison gave a distance-only interaction (no dependence on population size) with infectivity scaling linearly with mortality in a density-independent framework. As we found with the US data, it is difficult to disentangle the effects of population and infectivity parameters because these feature in different combinations in comparable DIC models.

DISCUSSION
We have presented a statistical analysis of the spatiotemporal spread of the 1918 influenza pandemic between cities in England and Wales and the US. The results demonstrate that for England and Wales, a model with intermediate levels of density dependence in the connectivity between cities gives the best fit to the observed pattern. For the US dataset, where there are few, large and widely spaced population centres, estimating the degree of density dependence does not improve the fit. In both contexts, city population size affects inter-city coupling sub-linearly. Parameter estimates and model formulation inferred from the data of England and Wales explain the US dataset well. Gravity model parameter estimates generated in this study are comparable with values found in studies describing the spread of seasonal influenza [5,44].
Our analysis demonstrates the degree of spatial locality in the large-scale geographical spread of influenza in both England and Wales and the US in 1918. However, it is difficult to directly compare the kernel power estimate from this study with those from other studies owing to differences in the functional forms used. For instance, Viboud et al. [5] estimate two power parameters above and below a given distance threshold when modelling the spread of seasonal influenza in the US. Gravity models used to describe the spread of measles in the UK by Xia et al. [18] assumed a kernel power of 1, rather than fitting this parameter. The distance power estimates we found for England and Wales and the US are quite different from each other. It is not surprising that there is a disparity in the distance kernel in England and Wales and the US, as the spatial scale in the US is much larger than in England and Wales. In comparing the US and UK, it should be noted that the mean distance between cities is of course much larger in the US (see electronic supplementary material, table S1). In theory, this gives better resolution for estimating the kernel shape, as the range of inter-city separations is an order of magnitude larger than for the UK. However, this is counterbalanced by the smaller size of the US dataset, which reduces inferential power. The low kernel power parameter estimates we have found in both England and Wales and the US suggests that long-distance interactions were important in spreading influenza between distant cities in both countries. At the start of the major autumn wave in 1918, the armistice was more than two months away and it is likely that travel relating to the war effort, including troop movements, might have enhanced the frequency of long-distance movements. Density-dependent gravity models are frequently used to explain the connectivity of urban centres for human diseases [5,18]. There is good evidence from the estimates of 1 in this analysis for England and Wales that the density-dependent model underestimates the total force of infection on remote cities. The 95 per cent credible interval for 1 includes 1, which indicates that the density-independent model formulation cannot be definitively excluded as an explanation for the data. There was limited statistical power to estimate the degree of density dependence in the US context, but use of the England and Wales best-fit model to describe the US data gave a very similar DIC to the best-fit US model. Hence, it is not clear if the difference in the estimated density dependence found between the US and the England and Wales is because of the large differences in the degree of population coverage between the US and England and Wales datasets. Our results using a  subset of the England and Wales dataset suggest that the degree of density dependence is a difficult parameter to estimate when coverage is low. The low power on destination city size found in fitting the gravity model to the England and Wales dataset shows that connectivity of a city increases sub-linearly with population increase. When modelling spread of influenza in the US, Viboud et al. [5] found very comparable low values of the population exponents with the infectious city lower than the susceptible. Differing results come from the analysis of measles data in Great Britain with a power coefficient on infectious populations estimated at approximately 1.5 [18]. We found the best-fit model in England and Wales does not include the population size of the infector (origin) city, a result which needs further examination in future work. Differences in our population parameter estimates and those from studies on contemporary populations are likely to differ owing to changes in human mobility patterns since 1918. Our estimate for England and Wales that the infectivity of a city is sub-linearly related to mortality, suggests that the rate of death in a city is not as important to infectivity as the presence or absence of disease. Other studies have used constant infectivity terms for the analysis of human seasonal influenza [21]. However, our estimates do support some level of mortality dependence, suggesting that cities with a very high influenza burden, usually later in the epidemic wave, are more infectious than newly infected cities. In the US dataset, the best-fit model gave constant infectiousness, but again this may be due to a lack of power to estimate such parameters from the US dataset. It may also be caused by non-uniform infection pressure from cities not in the dataset, which could mask an infectivity relationship for the cities that are given.
Future possible extensions of this work include relaxing the assumption that all cities were equally susceptible at the start of the autumn wave of the pandemic. The variation in the onset of infection in cities may, in part, be due to the differing susceptibility of each city owing to differing attack rates experienced in the spring -summer wave, or population-level immunity from the 1890 pandemic or seasonal strains. However, the low case fatality of the first wave and age-specificity of infection between waves need to be understood before spatial heterogeneity in susceptibility can be discerned. There are varying reports on the magnitude and mechanism of the effect of infection during the first wave on attack rates in subsequent waves [24][25][26][45][46][47]. Further analysis of the datasets considered here may provide an opportunity to disentangle these effects.