Main

Predicting the epidemiology of the COVID-19 pandemic is a priority for guiding epidemic responses around the world. China has undergone its first epidemic wave, and, remarkably, cities across the country are now reporting few or no locally acquired cases8. Analyses have indicated that the spread of COVID-19 from Hubei to the rest of China was driven primarily by human mobility from Wuhan6,9, and that the stringent measures to restrict human movement and public gatherings within and among cities in China were associated with bringing local epidemics under control5. Key uncertainties remain as to which geographic factors drive the local transmission dynamics of COVID-19, and initial analysis suggests a limited role of climate in determining epidemic growth10.

Spatial heterogeneity in infectious disease transmission can be influenced by local differences in population or human movements, such that high local population densities might catalyze the spread of new pathogens due to higher contact rates with susceptible individuals11,12. For respiratory pathogens, the temporal clustering of cases in an epidemic (that is, the shortest period during which most cases are observed) varies with increased indoor crowding and socio-economic and climatic factors13,14,15,16,17,18. The temporal concentration of cases is minimized when incidence is spread evenly across time and increases as incidence becomes more concentrated in particular days, as has been observed for influenza13. In any given location, a higher temporal concentration of cases might require a larger surge capacity in the public health system19, especially for an emerging respiratory pathogen such as COVID-19 (ref. 20).

Results

Spatial population structure predicts the shape of epidemics of COVID-19

China and Italy provide detailed epidemiological time series for COVID-19 (refs. 2,21,22) across a wide range of geographic contexts; hence, the outbreaks in these countries provide an opportunity to evaluate the role of local factors in shaping epidemic behavior. We used daily epidemiological data from Chinese cities23,24 and Italian provinces, climate and population data and the response to local interventions as measured by human mobility data from Baidu25 and the COVID-19 Aggregated Mobility Research Dataset (https://www.google.com/covid19/mobility/) to identify drivers of transmission, with a focus on how the temporal clustering of cases differs between prefectures in China and provinces in Italy. A summary of the main findings, limitations and policy implications of our study is shown in Table 1.

Table 1 Policy summary
Full size table

We used daily incidence data of confirmed COVID-19 cases aggregated at the prefectural level (n = 293) in China (Fig. 1a) and at the province level (n = 108) in Italy. Prefectures and provinces are administrative units that typically have one urban center (Fig. 1b). We aggregated daily individual-level data collected from official government reports22. Epidemiological data in each prefecture were truncated to exclude dates before the first and after the last day of reported cases during the first epidemic. Cases reported after March 1, 2020, that were imported from outside China were excluded from the analysis. All epidemiological data from Hubei Province were excluded because of the lack of prefecture-level epidemiological data and issues with consistent reporting before January 20, 2020. The shape of epidemic curves varied among prefectures, with some showing a rapid rise and decline in reported cases and others showing more prolonged epidemics (Fig. 1a and Extended Data Fig. 1).

Fig. 1: Maps of crowding in prefectures in China.
figure 1

a, Examples of epidemic curves that are normalized to show the percentage of cases across the whole epidemic that occur at each given day. Beijing and Shanghai (red) have less peaked epidemics than Wenzhou and Zhuhai. b, Examples of prefectures in China with different levels of crowding and population size. The color scale illustrates the estimated number of inhabitants per grid cell (1 km × 1 km). c, Relationship between the Shannon index of the incidence curve and the final attack rate for prefectures in China.

Full size image

To characterize the temporal clustering of cases for each prefecture and province, we calculated the Shannon diversity index of the distribution of incident cases13. We defined the incidence distribution pij for a given city to be the proportion of COVID-19 cases during the first epidemic wave j that occurred on day i. The Shannon index of incidence for a given prefecture and year is given by \({\it{v}}_{\it{j}} = \left( { - \mathop {\sum }\limits_i p_{ij}\log p_{ij}} \right)^{ - 1}\). Because vj is a function of the disease incidence curve in each location, rather than of absolute incidence values, it is less sensitive to varying reporting rates among cities. The Shannon index is maximal when all cases occur on the same day and minimal when each day of the epidemic has the same number of incident cases (for example, ‘flat’ epidemic curves). It is highly correlated with alternative measures of epidemic peakedness, such as the proportion of cases that occur at the peak ± 1 d (Extended Data Fig. 2). The total attack rate of reported COVID-19 cases in each prefecture is strongly negatively correlated with the Shannon index in China (Fig. 1c); hence, less peaked epidemics have a larger total attack rate (Pearson’s r = −0.67, 95% confidence interval (CI), −0.73 to −0.59, P < 0.01; for Italy, R2 = 0.33, P < 0.01). We hypothesize that this variation among cities in total attack rate and the temporal clustering of cases is the result of the spatial organization of human populations.

To test this hypothesis, we used Lloyd’s index of mean crowding13,26, treating the population count of each spatial grid cell as an individual unit (Fig. 1). The term ‘mean crowding’ used here is a specific geographic metric that summarizes both population density and how density is distributed across a prefecture (that is, patchiness; Fig. 1). Higher values of Lloyd’s index suggest a spatially aggregated population structure. For example, Xi’an has high values of crowding, whereas Bozhou has a similar population density but a population that is more evenly distributed across the prefecture (Fig. 1b). We performed log-linear regression modeling to determine the association between the temporal clustering of cases with socio-economic and environmental variables, including reductions in population flows during the outbreak period (Methods).

We found that the temporal clustering of cases was significantly negatively correlated with the mean number of contacts (P < 0.01) but positively correlated with mean population density (P < 0.01) and varies widely across China and Italy (Fig. 2 and Supplementary Table 1). This observation contrasts with the expectations of simple and classical epidemiological models, which predict higher peakedness in crowded areas due to the increased availability of susceptible individuals27,28. The spatial scale at which this relationship is best explained was 10 × 10 km, but results were statistically significant at all spatial scales between 1 and 50 km2 (Extended Data Fig. 3; P < 0.01). Mean specific humidity and population mobility remained significantly negatively correlated with epidemic peakedness when included in a multivariate model with crowding (Supplementary Table 1; all P < 0.01).

Fig. 2: Crowding and the temporal clustering of transmission of COVID-19 in China.
figure 2

a, Negative association between log10 of epidemic peakedness, as measured by Shannon’s diversity index, and log population crowding, as measured by Lloyd’s mean crowding. The point sizes indicate the size of the population in each city. b, Map of epidemic peakedness in China at the prefectural level. Blue and green colors indicate lower peakedness; red and yellow colors indicate higher peakedness. Gray prefectures had either no reported cases or were not included in analyses due to potential inconsistencies in reporting of early cases (Hubei Province).

Full size image

Using weekly human mobility data, we found that within-city human mobility during the outbreak was correlated with the temporal clustering of cases—that is, prefectures that have larger reductions in mobility also have lower epidemic peakedness (Extended Data Fig. 4 and Supplementary Table 1; P < 0.01). When we combined mobility reduction in a model with crowding and humidity, we found that these variables each remained significant predictors of the temporal clustering of cases (Extended Data Table 1; P < 0.01). These results suggest that, although measures to reduce mobility can successfully lead to a flattening of the epidemic curve, population crowding is an independent contributor to the shape of epidemics in these two countries.

Our multivariate model can explain a large fraction of the variation in epidemic peakedness among Chinese cities and Italian provinces, and sensitivity analyses confirm the robustness of our results to potential noise in location-specific incidence distributions (R2 = 0.638; Extended Data Fig. 2, Supplementary Table 1 and Extended Data Fig. 5). To evaluate the out-of-sample performance of our model, we 1) performed n-fold cross validation at the prefecture level in China (Spearman’s ρ = 0.61; 95% bootstrap CI, 0.52–0.68; P < 0.01); 2) used the fitted model in China to estimate peak intensity at the corresponding administrative level 2 locations—that is, province level, in Italy (Spearman’s ρ = 0.57; 95% bootstrap CI, 0.41–0.69; P < 0.01); and 3) performed n-fold cross validation at the province level in Italy (Spearman’s ρ = 0.65; 95% bootstrap CI. 0.52–0.76; P < 0.01). These results suggest that the model is robust to both within- and between-country out-of-sample testing (Extended Data Fig. 6).

To evaluate the potential effect of the temporal clustering of cases on the peak attack rate and total attack rate, we performed a simple linear regression (Supplementary Table 2). For locations that have a single peak, the attack rate at the peak is highest in two settings: 1) in crowded locations with high population size (prefectures that also experience high total attack rates); and 2) in locations that have lower population and lower crowding and, therefore, high temporal clustering of cases (Extended Data Fig. 7). Other prefectures that have low population and low crowding sometimes experience very short outbreaks with a small peak attack rate, suggesting local stochastic extinction possibly due to limited mixing between populations. We hypothesize that the observation that high peak attack rates can sometimes be found in low crowding areas is related to rare super-spreading events as observed in Bergamo, Italy, or Mulhouse, France.

Simulation of COVID-19 epidemics in hierarchically structured populations

We hypothesize that the mechanism underlying our central observation—that more crowded cities experience less peaked outbreaks—is that crowding enables sustained transmission among households and through a city’s population, leading incidence to be widely distributed through time. To explore this proposed mechanism, we simulated stochastic epidemic dynamics in two types of populations. Simple, well-mixed transmission models in which contact rates are high in crowded regions were not consistent with our findings, because they predict that crowded regions would have more temporally clustered outbreaks. To capture realistic contact patterns, we created hierarchically structured populations29 in which individuals had high rates of contact within their social units (which are defined broadly and could represent households, care homes, hospitals, prisons, etc); lower rates with individuals from other units but within the same neighborhoods; and relatively rare contact with other individuals in other neighborhoods within the same prefecture (Fig. 3a). These assumptions are consistent with reports that most onward transmission after lockdowns were implemented occurred in households or in other close-contact situations2,30. In this scenario, less crowded prefectures often had more peaked and shorter outbreaks that were isolated to specific neighborhoods, whereas more crowded prefectures could sustain drawn-out outbreaks of larger final size, which jumped among the more highly connected neighborhoods (Fig. 3b,c). Further, if the reproduction number of COVID-19 is over-dispersed31,32,33, then crowding could enable local outbreaks to spread more widely due to the availability of contacts34.

Fig. 3: Mechanisms generating less peaked epidemics in crowded populations.
figure 3

a, Schematic of a hierarchically structured population model consisting of households and ‘neighborhoods’ within a prefecture. Transmission is more likely among contacts connected at lower spatial levels. Crowded populations have a greater number of contacts outside the household, and interventions reduce the number of these connections in both populations (pink dotted lines). b, c, Simulated outbreak dynamics in the absence of interventions in crowded versus sparse populations. For the networks in b, blue nodes are individuals who were eventually infected by the end of the outbreak. In c, thin blue lines show individual realizations of the model, the average shown by the thick gray line. d, Simulated outbreak dynamics in the presence of strong social distancing measures in crowded versus sparse populations. The intervention was implemented at day 15 (vertical dotted line) and led to a 75% reduction in contacts, similar to observed changes in contact rates in China35,36. Mean values of median log epidemic peakedness (Shannon index) are −2.3 for low crowding and −2.8 for high crowding.

Full size image

We also simulated outbreak dynamics under extensive social distancing measures, as observed in Chinese prefectures (75% reduction in contact rates35,36). If social distancing reduces non-household contacts by the same relative amount in all locations, there will be more contacts remaining in crowded areas, because baseline contact rates are higher. Consequently, outbreaks in crowded regions could be larger and take longer to end after intervention (Fig. 3d, Fig. 1c and Extended Data Fig. 1).

Using the fitted model from China paired with globally comprehensive covariates, we extrapolated our results to cities across the world (Fig. 4). Human mobility data from Baidu were not available for locations outside of China. Therefore, we used aggregated human mobility data from Google’s COVID Mobility Research Dataset (Methods) to capture relative differences in human mobility through time. At the global scale, cities in yellow are predicted to have concentrated and peaked epidemics, whereas cities in blue are predicted to have more prolonged outbreaks (Fig. 4b; a full list is provided in the Supplementary Information). In general, the epidemics in coastal cities were less peaked and were larger and more prolonged, which could be attributable to high levels of population crowding in coastal cities. These predictions rely on fitted relationships of the first epidemic curves from Chinese and Italian cities and, therefore, should be interpreted very cautiously when generalizing to other settings.

Fig. 4: Predicted epidemic peakedness across the world.
figure 4

a, Maps of cities and their population densities at a 1 × 1-km scale. Madrid, Spain, and Colombo, Sri Lanka, have low predicted peakedness, whereas Novosibirsk, Russia, and Ulaanbaatar, Mongolia, have high predicted peakedness. b, Map of predicted epidemic peakedness for 310 cities across the world for which both human population data and mobility data were available for the study period.

Full size image

Discussion

Our findings confirm previous work on the peakedness of epidemics transmission for influenza in cities13. Our work provides empirical support for the role of spatial organization in determining infectious disease dynamics29,37 and, specifically, spatial variability in transmission parameters38. Furthermore, with lower total incidence in small cities compared to larger cities, the risk of resurgence could be elevated owing to lower population immunity after the first wave of the epidemic. Higher seroprevalence for COVID-19 in urban areas39 provides initial data to support these findings; however, there remains an urgent need to expand serological data collection and provide a full picture of attack rates across cities worldwide40. Even though our model does not account for over-dispersion in COVID-19 transmission, there is a theoretical link between the reproduction number in heterogeneous environments and Lloyd’s crowding index of aggregation41, such that the reproduction number increases with higher aggregation34. We report that, in dense cities, reductions in mobility tend to be larger, which potentially elevates the effectiveness of non-pharmaceutical interventions in dense cities42. However, assessing the effect of within-city connectivity and its spatial heterogeneity on disease dynamics will be critical to further our understanding of how COVID-19 spreads in urban areas. We found that there is an association between climatic factors and the peakedness of epidemics, but particular caution will need to be applied in interpreting these relationships outside the two studied countries (Italy and China). More work is needed to provide causal evidence for the effect of climatic factors on transmission dynamics of COVID-19 during the pandemic and post-pandemic phases10.

Currently, non-pharmaceutical interventions are the primary control strategy for COVID-19. As a result, public health measures are often focused on ‘flattening the curve’ to lower the risk of essential services running out of capacity. We show that spatial context, especially crowding, are important factors for assessing the shape of epidemic curves. Therefore, it will be critical to view non-pharmaceutical interventions through the perspective of crowding—that is, how does an intervention reduce the circle of contacts of an average individual—in cities across the world.

Methods

Epidemiological data

No officially reported line list was available for cases in China. We used a standardized protocol43 to extract individual-level data from December 1, 2019, to March 30, 2020. Sources were mainly official reports from provincial, municipal or national health governments. Data included basic demographics (age and sex), travel histories and key dates (dates of onset of symptoms, hospitalization and confirmation). Data were entered by a team of data curators on a rolling basis, and technical validation and geo-positioning protocols were applied continuously to ensure validity. A detailed description of the methodology is available22. Lastly, total numbers were matched with officially reported data from China and other government reports. Daily case counts from Italian provinces (n = 107) were extracted from the Presidenza del Consiglio dei Ministri Dipartimento della Protezione Civile (https://github.com/pcm-dpc/COVID-19).

Estimating epidemic peakedness

Epidemic peakedness was estimated for each prefecture by calculating the inverse Shannon entropy of the distribution of COVID-19 cases. Inverse Shannon entropy was used to fit time series of other respiratory infections (influenza)13. The inverse Shannon entropy of incidence for a given prefecture in 2020 is then given by \(v_j = \left( { - \mathop {\sum}\nolimits_i {p_{ij}\log p_{ij}} } \right)^{ - 1}\). Because vj is a function of incidence distribution in each location rather than raw incidence, it is invariant under differences in overall reporting rates between cities or attack rates. We then assessed how peakedness \({\it{v}} \propto \mathop {\sum}\nolimits_j {v_j}\) varied across geographic areas in China. As an alternative measure of temporal clustering of cases, we estimated the proportion of cases at the peak ± 1 d (Extended Data Fig. 2).

Proxies for COVID-19 interventions using within-city human mobility data from China

Estimates of within-city reductions of human mobility between the period before and after the lockdown was implemented on January 23, 2020, were extracted from Lai et al.36. Daily measures of human mobility were extracted from the Baidu Qianxi web platform to estimate the proportion of daily movement within prefectures in China. Relative mobility volume was available from January 2, 2020, to January 25, 2020. For each city, change in relative mobility was defined by mi=mil(lockdown)/mib(baseline), where mi is defined as mobility in prefecture i. Baidu’s mapping service is estimated to have a 30% market share in China, and more data can be found5,6.

Data on drivers of transmission of COVID-19

Prefecture-specific population counts and densities were derived from the 2020 Gridded Population of The World, a modeled continuous surface of population estimated from national census data and the United Nations World Population Prospectus44. Population counts are defined at a 30-arc-second resolution (approximately 1 km × 1 km at the equator) and extracted within administrative 2 level cartographic boundaries defined by the National Bureau of Statistics of China. Lloyd’s mean crowding, \(\frac{{\left[ {\mathop {\sum }\nolimits_i \left( {q_i - 1} \right)q_i} \right]}}{{\mathop {\sum }\nolimits_i q_i}}\), was estimated for each prefecture, where qi represents the population count of each non-zero pixel within a prefecture’s boundary and the resulting value estimates an individual’s mean number of expected neighbors13,45. When fitting the models, we consider the numerator \(\left[ {\mathop {\sum}\nolimits_i {\left( {q_i - 1} \right)q_i} } \right]\), which we refer to as ‘contacts’, and the denominator \(\mathop {\sum}\nolimits_i {q_i}\) (that is, population size) as separate predictors. We note that a negative slope for ‘contacts’ and a positive slope for ‘population’ support a negative coefficient for Lloyd’s mean crowding.

Daily temperature (°F), relative humidity (%) and atmospheric pressure (Pa) at the centroid of each prefecture was provided by The Dark Sky Company via the Dark Sky API and aggregated across a variety of data sources. Specific humidity (kg/kg) was then calculated using the R package humidity16. Meteorological variables for each prefecture were then averaged across the entirety of the study period.

Statistical analysis

We normalized the values of epidemic peakedness between 0 and 1 and, for all non-zero values, fit a generalized linear model of the form

$${\mathrm{log}}\left( {Y_j} \right)\sim \beta _0 + \beta _1{\mathrm{log}}\left( {C_j} \right) + \beta _2{\mathrm{log}}\left({h_j}\right) + \beta _3{\mathrm{log}}\left( {P_j} \right) + \beta _4{\mathrm{log}}\left( {f_j} \right) + \beta _5{\mathrm{log}}\left( {t_j} \right)$$

where, for each prefecture j, Y is the scaled inverse Shannon entropy measure of epidemic peakedness derived from the COVID-19 time series; C is the mean number of contacts26,46; h is the mean specific humidity over the reporting period in kg/kg; P is the estimated population density; f is the relative change in population flows within each prefecture; and t is daily mean temperature.

Projecting epidemic peakedness in cities around the world

We selected 310 urban centers from the European Commission Global Human Settlement Urban Centre Database and their included cartographic boundaries47. To ensure global coverage, up to the five most populous cities in each country were selected from the 1,000 most populous urban centers recorded in the database. Population count, crowding and meteorological variables were then estimated following identical procedures used to calculate these variables in the Chinese prefectures. Weather measurements were averaged over the 2-month period starting on February 1, 2020.

The parameters from the model of epidemic peakedness predicted by humidity, crowding and population size (Supplementary Table 1, model 6) were used to estimate relative peakedness in the 310 urban centers. A full list of predicted epidemic peakedness values can be found in Supplementary Table 3.

Global human mobility data

We used the Google COVID-19 Aggregated Mobility Research Dataset, which contains anonymized relative mobility flows aggregated over users who have turned on the Location History setting, which is off by default. This is similar to the data used to show how busy certain types of places are in Google Maps, helping identify when a local business tends to be the most crowded. The mobility flux is aggregated per week, between pairs of approximately 5-km2 cells worldwide, and for the purpose of this study aggregated for 310 cities worldwide. We calculated both mobility within each city’s shapefile and mobility coming into each city. For each city, change in relative mobility was defined by mi = mil(April)/mib(December), where mi is defined as mobility in city i.

To produce this data set, machine learning was applied to log data to automatically segment it into semantic trips48. To provide strong privacy guarantees, all trips were anonymized and aggregated using a differentially private mechanism49 to aggregate flows over time (https://policies.google.com/technologies/anonymization). This research is done on the resulting heavily aggregated and differentially private data. No individual user data were ever manually inspected; only heavily aggregated flows of large populations were handled.

All anonymized trips were processed in aggregate to extract their origin and destination location and time. For example, if users traveled from location a to location b within time interval t, the corresponding cell (a, b, t) in the tensor would be n ± err, where err is Laplacian noise. The automated Laplace mechanism adds random noise drawn from a zero-mean Laplace distribution and yields (𝜖, δ)-differential privacy guarantee of 𝜖 = 0.66 and δ = 2.1 × 10−29. The parameter 𝜖 controls the noise intensity in terms of its variance, whereas δ represents the deviation from pure 𝜖-privacy. The closer they are to zero, the stronger the privacy guarantees. Each user contributes, at most, one increment to each partition. If they go from a region a to another region b multiple times in the same week, they contribute only once to the aggregation count.

These results should be interpreted in light of several important limitations. First, the Google mobility data are limited to smartphone users who have opted in to Google’s Location History feature, which is off by default. These data might not be representative of the population as whole, and, furthermore, their representativeness might vary by location. Importantly, these limited data are viewed only through the lens of differential privacy algorithms, specifically designed to protect user anonymity and obscure fine detail. Moreover, comparisons across, rather than within, locations are descriptive only because these regions can differ in substantial ways.

Simulating epidemic dynamics

We simulated a simple stochastic SIR model of infection spread on weighted networks created to represent hierarchically structured populations. Individuals were first assigned to households using the distribution of household sizes in China (data from the United Nations Population Division; mean, 3.4 individuals). Households were then assigned to ‘neighborhoods’ of ~100 individuals, and all neighborhood members were connected with a lower weight. A randomly chosen 10% of individuals were given ‘external’ connections to individuals outside the neighborhood. The total population size was n = 1,000. Simulations were run for 300 d, and averages were taken over 20 iterations. The SIR model used a per-contact transmission rate of 𝛽 = 0.15 per day and recovery rate 𝛾 = 0.1 per day. For the simulations without interventions, the weights were wHH = 1, wNH = 0.01 and wEX = 0.001 for the crowded prefecture and wEX = 0.0001 for the less crowded prefecture. For the simulations with interventions, the household and neighborhood weights were the same, but we used wEX = 0.01 for the crowded prefecture and wEX = 0.001 for the ‘sparse’ prefecture. The intervention reduced the weight of all connections outside the household by 75%.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.