key: cord-0900478-q2d2sut8
authors: Chang, S. Y.; Pierson, E.; Koh, P. W.; Gerardin, J.; Redbird, B.; Grusky, D.; Leskovec, J.
title: Mobility network modeling explains higher SARS-CoV-2 infection rates among disadvantaged groups and informs reopening strategies
date: 2020-06-17
journal: nan
DOI: 10.1101/2020.06.15.20131979
sha: 54a105a0655b557b5aa05fa59cefcad1ad143265
doc_id: 900478
cord_uid: q2d2sut8

Fine-grained epidemiological modeling of the spread of SARS-CoV-2 -- capturing who is infected at which locations -- can aid the development of policy responses that account for heterogeneous risks of different locations as well as the disparities in infections among different demographic groups. Here, we develop a metapopulation SEIR disease model that uses dynamic mobility networks, derived from US cell phone data, to capture the hourly movements of millions of people from local neighborhoods (census block groups, or CBGs) to points of interest (POIs) such as restaurants, grocery stores, or religious establishments. We simulate the spread of SARS-CoV-2 from March 1-May 2, 2020 among a population of 105 million people in 10 of the largest US metropolitan statistical areas. We show that by integrating these mobility networks, which connect 60k CBGs to 565k POIs with a total of 5.4 billion hourly edges, even a relatively simple epidemiological model can accurately capture the case trajectory despite dramatic changes in population behavior due to the virus. Furthermore, by modeling detailed information about each POI, like visitor density and visit length, we can estimate the impacts of fine-grained reopening plans: we predict that a small minority of "superspreader" POIs account for a large majority of infections, that reopening some POI categories (like full-service restaurants) poses especially large risks, and that strategies restricting maximum occupancy at each POI are more effective than uniformly reducing mobility. Our models also predict higher infection rates among disadvantaged racial and socioeconomic groups solely from differences in mobility: disadvantaged groups have not been able to reduce mobility as sharply, and the POIs they visit (even within the same category) tend to be smaller, more crowded, and therefore more dangerous. By modeling who is infected at which locations, our model supports fine-grained analyses that can inform more effective and equitable policy responses to SARS-CoV-2.

We overlaid an SEIR disease model on the mobility network, with each CBG having its own set of SEIR compartments. New infections occur at both POIs and CBGs. The model has just three free parameters, which remain fixed over time, scaling transmission rates at POIs; transmission rates at CBGs; and the initial fraction of infected individuals. To determine the transmission rate at a given time at each POI we use the mobility network, which captures population movements as well as visit duration and the POI physical area, to estimate the density of visitors at each POI. (c) Left: To test out-of-sample prediction, we calibrated the model on data before April 15, 2020 (vertical black line). Even though its parameters remain fixed over time, the model accurately predicts the case trajectory after April 15 by using mobility data. Shaded regions denote 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations. Right: Model fit further improved when we calibrated the model on the full range of data. (d) We fit separate models to 10 of the largest US metropolitan statistical areas (MSAs), modeling a total population of 105 million people; here, we show full model fits, as in (c)-Right. While we use Washington DC as a running example throughout the paper, we include results for all other MSAs in the SI. Cumulative distribution of predicted infections over POIs Figure 2 : Assessing mobility reduction and reopening policies. (a) Counterfactual simulations (left) of the mobility reduction in March 2020-scaling its magnitude down, or shifting the timeline earlier or later-illustrate that the magnitude of mobility reduction (middle) was at least as important as its timing (right). Shaded regions denote 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations. (b) Most infections at POIs occur at a small fraction of "super-spreader" POIs: 10% of POIs account for more than 80% of the total infections that occurred at POIs in the Washington DC MSA (results for other MSAs in Extended Data Figure 3 ). (c) Left: We simulated partial reopening by clipping hourly visits if they exceeded a fraction of each POI's maximum occupancy. We plot cumulative infections at the end of one month of reopening against the fraction of visits lost by partial instead of full reopening; the annotations within the plot show the fraction of maximum occupancy used for clipping. Reopening leads to an additional 26% of the population becoming infected by the end of the month, but clipping at 20% maximum occupancy cuts down new infections by more than 80%, while only losing 40% of overall visits. Right: Compared to partially reopening by uniformly reducing visits, the clipping strategy-which disproportionately targets high-risk POIs with sustained high occupancy-always results in a smaller increase in infections for the same number of visits. The y-axis plots the relative difference between the increase in cumulative infections (from May 1 to May 31) under the clipping strategy as compared to the uniform reduction strategy. (d) We simulated reopening each POI category while keeping reduced mobility levels at all other POIs. Boxes indicate the interquartile range across parameter sets and stochastic realizations. Reopening full-service restaurants has the largest predicted impact on infections, due to the large number of restaurants as well as their high visit densities and long dwell times. (a) Across all MSAs, our model predicts that people in lower-income census block groups (CBGs) are more likely to be infected, even though they start with equal probabilities of being infected. Disparities are especially prominent in Philadelphia, which we discuss in SI Section S2. Boxes indicate the interquartile range across parameter sets and stochastic realizations. (b) Racial disparities are similar: people in non-white CBGs are typically more likely to be infected, although results are more variable. (c-e) illustrate how mobility patterns give rise to socioeconomic disparities; similar mechanisms underlie racial disparities (Extended Data Figure 6 , Table S4 ). (c) The overall disparity is driven by a few POI categories like full-service restaurants. Shaded regions denote 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations. (d) One reason for the disparities is that higher-income CBGs were able to reduce their overall mobility levels below those of lower-income CBGs. (e) Within each category, the POIs that people from lowerincome CBGs visit also tend to have higher transmission rates because they are smaller and more crowded. Thus, even if a lower-income and a higher-income person went out equally often and went to the same types of places, the lowerincome person would still have a greater risk of infection. The size of each dot indicates the total number of visits to that category. (f) We predict the effect of reopening (at different levels of clipping maximum occupancy) on different demographic groups. Reopening leads to more infections in lower-income CBGs (purple) than the overall population (blue), underscoring the need to account for disadvantaged subpopulations when assessing reopening plans.

In response to the SARS-CoV-2 crisis, numerous stay-at-home orders were enacted across the 2 United States in order to reduce contact between individuals and slow the spread of the virus. 1 3 As of May 2020, these orders are being relaxed, businesses are beginning to reopen, and mobility 4 is increasing, causing concern among public officials about the potential resurgence of cases. 2 5 Epidemiological models that can capture the effects of changes in mobility on virus spread are 6 a powerful tool for evaluating the effectiveness and equity of various strategies for reopening or 7 responding to a resurgence. In particular, findings of SARS-CoV-2 "super-spreader" events 3-7 8 motivate models that can reflect the heterogeneous risks of visiting different locations, while well- 9 reported racial and socioeconomic disparities in infection rates [8] [9] [10] [11] [12] [13] [14] require models that can explain 10 the disproportionate impact of the virus on disadvantaged demographic groups. 11 To address these needs, we construct a mobility network using US cell phone data from 12 March 1-May 2, 2020 that captures the hourly movements of millions of people from census 13 block groups (CBGs), which are geographical units that typically contain 600-3,000 people, to 14 points of interest (POIs) such as restaurants, grocery stores, or religious establishments. On top of 15 this dynamic bipartite network, we overlay a metapopulation SEIR disease model that tracks the 16 infection trajectories of each CBG over time as well as the POIs at which these infections are likely 17 to have occurred. The key idea is that combining even a relatively simple epidemiological model 18 with our fine-grained, dynamic mobility network allows us to not only accurately model the case 19 trajectory, but also identify the most risky POIs; the most at-risk populations; and the impacts of 20 different reopening policies. This builds upon prior work that models disease spread using mobility 21 data, which has used aggregate 15-21 , historical [22] [23] [24] , or synthetic [25] [26] [27] mobility data; separately, other 22 work has directly analyzed mobility data and the effects of mobility reductions in the context of 23 SARS-CoV-2, but without an underlying epidemiological model of disease spread. 28 -33 24 We use our model to simulate the spread of SARS-CoV-2 within 10 of the largest metropoli- 25 tan statistical areas (MSAs) in the US, starting from a low, homogeneous prevalence of SARS- 26 CoV-2 across CBGs. For each MSA, we examine the infection risks at individual POIs, the ef- 27 fects of past stay-at-home policies, and the effects of reopening strategies that target specific types 28 of POIs. We also analyze disparities in infection rates across racial and socioeconomic groups, 29 identify mobility-related mechanisms driving these disparities, and assess the disparate impacts of 30 reopening policies on disadvantaged groups. 31 als at each POI to determine its transmission rate. The model has only three free parameters, which information about the case trajectory or social distancing measures. In contrast, a baseline SEIR 60 6 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint model that does not use the mobility network has considerably worse out-of-sample fit (Extended 61 Data Figures 1b and 2b ). All subsequent results were generated using the models calibrated on the 62 entire range of case counts from March 8-May 9, 2020. 63 Evaluating mobility reduction and reopening policies 64 We can estimate the impact of a wide range of mobility reduction and reopening policies by apply- 65 ing our model to a modified mobility network that reflects the expected effects of a hypothetical 66 policy. We start by studying the effect of the magnitude and timing of mobility reduction poli- 67 cies from March 2020. We then assess several fine-grained reopening plans, such as placing a 68 maximum occupancy cap or only reopening certain categories of POIs, by leveraging the detailed 69 information that the mobility network contains on each POI, like its average visit length and visitor 70 density at each hour. 71 The magnitude of mobility reduction is as important as its timing. US population mobility 72 dropped sharply in March 2020 in response to SARS-CoV-2; for example, overall mobility in the 73 Washington DC MSA fell by 58.5% between the first week of March and the first week of April 74 2020. We constructed counterfactual mobility networks by scaling the magnitude of mobility compared to a less than 2× increase had people begun reducing their mobility one full week later. 83 We observe similar trends across other MSAs (Tables S1 and S2) .

A minority of POIs account for a majority of infections. Since overall mobility reduction re-85 duces infections, we next investigated if how we reduce mobility-i.e., to which POIs-matters. 86 infections at POIs (Figure 2b ; Extended Data Figure 3 shows similar results across MSAs). Note 90 that infections at POIs represent a majority, but not all, of the total infections, since we also model 91 infections within CBGs; across MSAs, the median proportion of total infections that occur at POIs 92 is 73%. These "superspreader" POIs are smaller and more densely occupied, and their occupants 93 stay longer, suggesting that it is especially important to reduce mobility at these high-risk POIs.

In the DC MSA, the median number of hourly visitors per square foot was 4.6× higher for the 95 riskiest 10% of POIs than for the remaining POIs; the median dwell time was 2.3× higher.

Reducing mobility by clipping maximum occupancy. We simulated the effects of two reopen-97 ing strategies, implemented beginning on May 1, on the increase in infections by the end of May.

First, we evaluated a "clipping" reopening strategy, in which hourly visits to each POI return to 99 those in the first week of March (prior to widespread adoption of stay-at-home measures), but are 100 capped if they exceed a fraction of the POI's maximum occupancy, 36 which we estimated as the 101 maximum hourly number of visitors ever recorded at that POI. A full return to early March mobil-102 ity levels without clipping produces a spike in predicted infections: in the DC MSA, we project that 103 an additional 26% of the population will be infected within a month ( Figure 2c ). However, clipping 104 substantially reduces risk without sharply reducing overall mobility: clipping at 20% maximum oc-105 cupancy in the DC MSA cuts down new infections by more than 80% but only loses 40% of overall 106 visits, and we observe similar trends across other MSAs (Extended Data Figure 4 ). This highlights 107 the non-linearity of infections as a function of visits: one can achieve a disproportionately large 108 reduction in infections with a small reduction in visits. 109 We also compared the clipping strategy to a baseline that uniformly reduces visits to each 110 POI from their levels in early March. Clipping always results in fewer infections for the same total 111 number of visits: e.g., clipping at 20% maximum occupancy reduces new infections by more than 112 25% compared to the uniform baseline for the same total number of visits in the Washington DC 113 MSA ( Figure 2c ). This is because clipping takes advantage of the heterogeneous risks across POIs, 114 disproportionately reducing visits at high-risk POIs with sustained high occupancy, but allowing 115 lower-risk POIs to return fully to prior mobility levels.

Relative risk of reopening different categories of POIs. We assessed the relative risk of re-117 opening different categories of POIs by reopening each category in turn on May 1 (and returning 118 its mobility patterns to early March levels) while keeping mobility patterns at all other POIs at 119 8 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint their reduced, stay-at-home levels ( Figure 2d ). We find a large variation in reopening risks: on 120 average across the 10 evaluated MSAs (Extended Data Figure 5 ), full-service restaurants, cafes, 121 gyms, limited-service restaurants, and religious establishments produce the largest increases in in-122 fections when reopened. Reopening full-service restaurants is particularly risky: in the Washington 123 DC MSA, we predict an additional 296k infections by the end of May, more than double the next 124 riskiest POI category. These risks are the total risks summed over all POIs in the category, but the 125 relative risks after normalizing by the number of POIs are broadly similar, with restaurants, gyms, 126 cafes, and religious establishments predicted to be the most dangerous on average per individual 127 POI. These categories are more dangerous because their POIs tend to have higher visit densities 128 and/or visitors stay there longer ( Figures S4-S13 ).

Infection disparities between socioeconomic and racial groups 130 We characterize the differential spread of SARS-CoV-2 along demographic lines by using US Cen-131 sus data to annotate each CBG with its racial composition and median income, then tracking how 132 infection disparities arise across groups. We use this approach to study the mobility mechanisms 133 behind disparities and to quantify how different reopening strategies impact disadvantaged groups.

Mobility patterns contribute to disparities in infection rates. Despite only having access to 135 mobility data and no other demographic information, our models correctly predicted higher risks of 136 infection among disadvantaged racial and socioeconomic groups. [8] [9] [10] [11] [12] [13] [14] Across all MSAs, individuals 137 from CBGs in the bottom decile for income were substantially likelier to have been infected by the 138 end of the simulation, even though all individuals began with equal likelihoods of infection in our 139 simulation ( Figure 3a ). This overall disparity was driven primarily by a few POI categories (e.g., 140 full-service restaurants), which infected far larger proportions of lower-income CBGs than higher-141 income CBGs (Figure 3c ; similar trends hold across all MSAs in Figure S1 ). We similarly found 142 that CBGs with fewer white residents had higher relative risks of infection, although results were 143 more variable (Figure 3b ). Our models also recapitulated known associations between population 144 density and infection risk 37 (median Spearman correlation between CBG density and cumulative 145 incidence proportion, 0.42 across MSAs), despite not being given any information on population 146 density. In SI Section S2, we confirm that the magnitude of the disparities our model predicts 147 are generally consistent with real-world disparities and explore the large predicted disparities in race. In the analysis below, we focus on the mechanisms producing higher relative risks of infection 150 among lower-income CBGs, and we show in Extended Data Figure 6 and Table S4 that similar   151 results hold for racial disparities as well.

Lower-income CBGs saw smaller reductions in mobility. Across all MSAs, we found that 153 lower-income CBGs were not able to reduce their mobility as sharply in the first few weeks of 154 March 2020, and had higher mobility than higher-income CBGs for most of March through May 155 (Figure 3d , Extended Data Figure 6 ). For example, over the month of April, lower-income CBGs in 156 the Washington DC MSA had 17% more visits per capita than higher-income CBGs. Differences 157 in mobility patterns within categories partially explained the within-category infection dispari-158 ties: e.g., lower-income CBGs made substantially more visits per capita to full-service restaurants 159 than did higher-income CBGs, and consequently experienced more infections at that category (Ex-160 tended Data Figure 7 ).

POIs visited by lower-income CBGs tend to be more dangerous. Differences in the number of 162 visits per capita between lower-and higher-income CBGs do not fully explain the infection dispar-163 ities: for example, in the DC MSA, grocery stores were visited more often by higher-income CBGs 164 but still caused more predicted infections among lower-income CBGs. We found that even within a 165 POI category, the transmission rate at POIs frequented by people from lower-income CBGs tended 166 to be higher than the corresponding rate for higher-income CBGs (Figure 3e; Table S3 ), because 167 these POIs tended to be smaller and more crowded. It follows that, even if a lower-income and 168 higher-income person had the same mobility patterns and went to the same types of places, the 169 lower-income person would still have a greater risk of infection.

As a case study, we examine grocery stores in further detail. Across all MSAs but Dallas, 171 visitors from lower-income CBGs encountered more dangerous grocery stores than those from 172 higher-income CBGs (median transmission rate ratio of 2.11, Table S3 ). Why was one visit to the 173 grocery store twice as dangerous for a lower-income individual? Taking medians across MSAs, 174 we found that the average grocery store visited by lower-income individuals had 45% more hourly 175 visitors per square foot, and their visitors stayed 27% longer on average. These findings highlight 176 how fine-grained differences in mobility patterns-how often people go out, which categories of 177 places they go to, which POIs they choose within those categories-can ultimately contribute to 178 dramatic disparities in infection outcomes. 179 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

Reopening plans must account for disparate impact. Because disadvantaged groups suffer a 180 larger burden of infection, it is critical to not just consider the overall impact of reopening plans 181 but also their disparate impact on disadvantaged groups specifically. For example, our model 182 predicted that full reopening in the Washington DC MSA would result in an additional 35% of 183 the population of CBGs in the bottom income decile being infected within a month, compared 184 to 26% of the overall population (Figure 3f ; results for all MSAs in Extended Data Figure 4 ).

Similarly, Extended Data Figure 8 illustrates that reopening individual POI categories tends to 186 have a larger impact on the bottom income decile. More conservative reopening plans produce 187 smaller absolute disparities in infections-e.g., we predict that clipping visits at 20% occupancy 

We model the spread of SARS-CoV-2 using a dynamic mobility network that encodes the hourly 192 movements of millions of people between 60k neighborhoods (census block groups, or CBGs) 193 and 565k points of interest (POIs). Because our data contains detailed information on each POI, 194 like visit length and visitor density, we can estimate the impacts of fine-grained reopening plans-

predicting that a small minority of "superspreader" POIs account for a large majority of infections, 196 that reopening some POI categories (like full-service restaurants) poses especially large risks, and 197 that strategies that restrict the maximum occupancy at each POI are more effective than uniformly 198 reducing mobility. Because we model infections in each CBG, we can infer the approximate de-199 mographics of the infected population, and thereby assess the disparate socioeconomic and racial 200 impacts of SARS-CoV-2. Our model correctly predicts that disadvantaged groups are more likely 201 to become infected, and also illuminates two mechanisms that drive these disparities: (1) dis-202 advantaged groups have not been able to reduce their mobility as dramatically (consistent with 203 previously-reported data, and likely in part because lower-income individuals are more likely to 204 have to leave their homes to work 10 ) and (2) when they do go out, they visit POIs which, even 205 within the same category, are smaller, more crowded, and therefore more dangerous.

The cell phone mobility dataset we use has limitations: it does not cover all populations (e.g., 207 prisoners, young children), does not contain all POIs (e.g., nursing homes), and cannot capture 208 sub-CBG heterogeneity in demographics. These limitations notwithstanding, cell phone mobil-209 ity data in general and SafeGraph data in particular have been instrumental and widely used in 210 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. MSAs, differentials in susceptibility (due to pre-existing conditions or access to care), various 213 transmission-reducing behaviors (e.g., hand-washing, mask-wearing), as well as POI-specific risk 214 factors (e.g., ventilation). Although our model recovers case trajectories and known infection dis-215 parities even without incorporating these processes, we caution that this predictive accuracy does 216 not mean that our predictions should be interpreted in a narrow causal sense, and that it is impor-217 tant to recognize that certain types of POIs or subpopulations may disproportionately select for 218 certain types of omitted processes. However, the predictive accuracy of our model suggests that 219 it broadly captures the relationship between mobility and transmission, and we thus expect our 220 broad conclusions-e.g., that lower-income CBGs have higher infection rates in part because they 221 have not been able to reduce mobility by as much, and because they tend to visit smaller, denser 222 POIs-to hold robustly.

Our results can guide policymakers seeking to assess competing approaches to reopening 224 and tamping down post-reopening resurgence. Despite growing concern about racial and socioe- As society reopens and we face the possibility of a resurgence in cases, it is critical to build 239 models which allow for fine-grained assessments of the effects of reopening policies. We hope 240 that our approach, by capturing heterogeneity across POIs, demographic groups, and cities, helps 241 address this need. 242 12 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. 15 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint

The methods section is structured as follows. We describe the datasets we use in Methods M1 244 and the mobility network that we derive from these datasets in Methods M2. In Methods M3, we 245 discuss the SEIR model we overlay on the mobility network, and in Methods M4, we describe how 246 we calibrate this model and quantify uncertainty in its predictions. In Methods M5, we provide 247 details on the experimental procedures used for our analysis of physical distancing, reopening, and 248 demographic disparities. Finally, in Methods M6, we elaborate on how we estimate the mobility 249 network from the raw mobility data. Table 1 ). We chose these MSAs by taking a random subset of the SafeGraph Patterns data and 269 picking the 10 MSAs with the most POIs in the data. Our methods in this paper can be straightfor- CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. percentile of median dwell times in that period. Summary statistics of the post-processed data are 286 in Extended Data Table 1 . Overall, we analyze over 59,000 CBGs from the 10 MSAs, and over 287 250M visits from these CBGs to over 565,000 POIs. 288 SafeGraph data has been used to study consumer preferences 46 and political polarization. 47

More recently, it has been used as one of the primary sources of mobility data in the US for tracking 290 the effects of the SARS-CoV-2 pandemic. 28, 30, [48] [49] [50] In SI Section S1, we show that aggregate trends CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. From US Census data, each CBG c i is labeled with its population N c i , income distribution, 314 and racial and age demographics. From SafeGraph data, each POI p j is similarly labeled with its 315 category (e.g., restaurant, grocery store, or religious organization), its physical size in square feet 316 a p j , and the median dwell time d p j of visitors to p j .

The central technical challenge in constructing this network is estimating the network weights 318

ij } from SafeGraph data, since this visit matrix is not directly available from the data.

Because the estimation procedure is involved, we defer describing it in detail until Methods M6; 320 in Methods M3-M5, we will assume that we already have the network weights. We use a SEIR model with susceptible (S), exposed (E), infectious (I), and removed (R) 328 compartments. Susceptible individuals have never been infected, but can acquire the virus through 329 contact with infectious individuals, which may happen at POIs or in their home CBG. They then 330 enter the exposed state, during which they have been infected but are not infectious yet. Individuals 331 18 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint transition from exposed to infectious at a rate inversely proportional to the mean latency period.

Finally, they transition into the removed state at a rate inversely proportional to the mean infectious 333 period. The removed state represents individuals who cannot infect others, because they have 334 recovered, self-isolated, or died.

Each CBG c i maintains its own SEIR instantiation, with S (t)

At each hour t, we sample the transitions between states as follows:

where λ

ij , the ij-th entry of the visit matrix from 339 the mobility network (Methods M2), is the number of visitors from CBG c i to POI p j at time t;

c i is the base rate of infection that is independent of visiting POIs; δ E is the mean latency period; 341 and δ I is the mean infectious period. 342 We then update each state to reflect these transitions. Let ∆S (t)

M3. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint New exposures from visiting POIs. We assume that any susceptible visitor to POI p j at time t 349 has the same independent probability λ (t) p j of being infected and transitioning from the susceptible 350 (S) to the exposed (E) state. Since there are w 

ij . 355 We model the infection rate at POI p j at time t, λ (t)

p j , as the product of its 

is the total number of visitors to p j at time t, 358 We model the transmission rate at POI p j at time t as

where a p j is the physical area of p j , and ψ is a transmission constant (shared across all POIs) that 360 we fit to data. The inverse scaling of transmission rate with area a p j is a standard simplifying 

For sufficiently large values of ψ and a sufficiently large proportion of infected individuals, the 368 expression above can sometimes exceed 1. To address this, we simply clip the infection rate to 1.

However, this occurs very rarely for the parameter settings and simulation duration that we use. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. .

c k /N c k in that CBG, although we note that the scaling factor ψ can 373 account for differences in the ratio of infectious individuals who visit POIs. This gives 374 at POIs that are absent from the SafeGraph data. We assume that at each hour, every susceptible

c i probability of becoming infected and transitioning to the exposed

is proportional to the infection density at CBG c i , and β base is a constant that we fit to data. 

new infections from visiting POIs

384 We model exposed individuals as becoming infectious at a rate inversely proportional to the mean 385 latency period δ E . At each time step t, we assume that each exposed individual has a constant, 386 time-independent probability of becoming infectious, with

21 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

We estimate both δ E and δ I from prior literature; see Methods M4.

In our experiments, t = 0 is the first hour of March 1, 2020. We approximate the infectious I and removed R compartments at t = 0 as initially empty, with all infected individuals in the exposed 393 E compartment. We further assume the same expected initial prevalence p 0 in every CBG c i . At 394 t = 0, every individual in the MSA has the same independent probability p 0 of being exposed E 395 instead of susceptible S. We thus initialize the model state by setting

397 Most of our model parameters can either be estimated from SafeGraph and US Census data, or 398 taken from prior work (see Extended Data Table 2 for a summary). This leaves 3 model parameters 399 that do not have direct analogues in the literature, and that we therefore need to calibrate with data: 400 1. The transmission constant in POIs, ψ (Equation (9)) 401 2. The base transmission rate, β base (Equation (11)) 402 3. The initial proportion of exposed individuals at time t = 0, p 0 (Equation (16)).

In this section, we describe how we fit these parameters to published numbers of confirmed cases, 404 as reported by The New York Times. We fit models for each MSA separately. In Methods M4.4,

we show that the resulting models can accurately predict the number of confirmed cases in out-of-406 sample data that was not used for model fitting. 407 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. Transmission rate factors ψ and β base . We select parameter ranges for the transmission rate fac-409 tors ψ and β base by checking if the model outputs match plausible ranges of the basic reproduction 410 number R 0 pre-lockdown, since R 0 has been the study of substantial prior work on SARS-CoV-411 2. 55 Under our model, we can decompose R 0 = R base + R POI , where R POI describes transmission 412 due to POIs and R base describes the remaining transmission (as in Equation (12)). We first establish 413 plausible ranges for R base and R POI before translating these into plausible ranges for β base and ψ. 414 We assume that R base ranges from approximately 0.1-1. R base models transmission that is 415 not correlated with POI activity, which includes within-household transmission. We chose the 416 lower limit of 0.1 because beyond that point, base transmission would only contribute minimally 417 to overall R, whereas previous work suggests that within-household transmission is a substantial 418 contributor to overall transmission. 56, 57 However, household transmission alone is not estimated to 419 be sufficient to tip overall R 0 above 1; for example, a single infected individual has been estimated 420 to cause an average of 0.32 (0.22, 0.42) secondary within-household infections. 56 We therefore 421 chose an upper limit of 1, corresponding to the assumption that R 0 < 1 when there is no POI 422 activity whatsoever (i.e., R POI = 0).

The plausible range for R POI is then determined by combining R POI = R 0 − R base with an 424 overall range, estimated from prior work, for pre-lockdown R 0 of 2-3. 55 Thus, R POI pre-lockdown 425 plausibly ranges from roughly 1-3.

To determine the values of R base and R POI that a given pair of β base and ψ imply, we seeded a 427 fraction of index cases and then ran the model on looped mobility data from the first week of March 428 to capture pre-lockdown conditions. We initialized the model by setting p 0 , the initial proportion 429 of exposed individuals at time t = 0, to p 0 = 10 −4 , and then sampling in accordance with Equation 430 (16) . Let N 0 be the number of initial exposed individuals sampled. We computed the number of 431 individuals that these N 0 index cases went on to infect through base transmission, N base , and POI 432 transmission, N POI , which gives

We averaged these quantities over 20 stochastic realizations per MSA. Figure S2 shows that, as 434 23 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint expected, R base is linear in β base and R POI is linear in ψ. R base lies in the plausible range when 435 β base ranges from approximately 0.001-0.012, and R POI lies in the plausible range (for at least one 436 MSA) when ψ ranges from approximately 1,000-10,000, so these are the parameter ranges we 437 consider when fitting the model. As described in Methods M4.2, we verified that case count data 438 for all MSAs can be fit using parameter settings for β base and ψ within these ranges. Using the parameter ranges above, we grid searched over ψ, β base , and p 0 to find the models that 446 best fit the number of confirmed cases reported by The New York Times (NYT). 35 For each of the 10 447 MSAs studied, we tested 1,260 different combinations of ψ, β base , and p 0 in the parameter ranges 448 specified above, with parameters linearly spaced for ψ and β base and logarithmically spread for p 0 .

In Methods M3, we directly model the number of infections but not the number of confirmed 450 cases. To estimate the number of confirmed cases, we assume that an r c = 0.1 proportion of in-451 fections will be confirmed, and moreover that they will confirmed exactly δ c = 168 hours (7 days) 452 after becoming infectious. We assume that these parameters are time-invariant. As a sensitivity 453 analysis, we alternatively stochastically sampled the number of confirmed cases and the confirma-454 tion delay from distributions with mean r c and δ c , but found that this did not change predictions 455 noticeably. We estimated these parameters, r c and δ c , from prior work (Extended Data Table 2 ).

From these assumptions, we can calculate the predicted number of newly confirmed cases 457 across all CBGs in the MSA on day d,

where for convenience we define N . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint of cases that our model predicts by computing the root-mean-squared-error (RMSE) 59 over the 462 D = T /24 days of our simulations,

For each combination of model parameters and for each MSA, we quantify model fit with the NYT 464 data by running 20 stochastic realizations and averaging their RMSE.

Our simulation spans March 1 to May 2, 2020, and we use mobility data from that period.

However, because we assume that cases will be confirmed δ c = 7 days after individuals become 467 infectious (Extended Data parameters ψ, β base , and p 0 . The latter is equivalent to assuming that the posterior probability over 485 the true parameters is uniformly spread among all parameter sets within the 20% threshold. 486 25 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint M4.4 Model validation on out-of-sample data 487 We validate our models by showing that they predict the number of confirmed cases and deaths on 488 out-of-sample data when we have access to corresponding mobility data. We then confirm that the 489 mobility data used as input in the model improves the fit to case and death data by comparing to a 490 model that does not use mobility data.

Out-of-sample prediction of the number of cases (Extended Data Figure 1 ). For each MSA, 492 we split the available NYT dataset into a training set (spanning March 8, 2020 to April 14, 2020) 493 and a test set (spanning April 15, 2020 to May 9, 2020). We fit the model parameters ψ, β base , and 494 p 0 , as described in Methods M4.2, but only using the training set. We then evaluate the predictive 495 accuracy of the resulting model on the test set. When running our models on the test set, we 496 still use mobility data from the test period. Thus, this is an evaluation of whether the models can 497 accurately predict the number of cases, given mobility data, in a time period that was not used for 498 model calibration. Extended Data Figure 1a shows that the models fit the out-of-sample case data 499 fairly well, demonstrating that they can extrapolate beyond the training set to future time periods.

Note that we only use this train/test split to evaluate out-of-sample model accuracy. All 501 other results are generated using parameter sets that best fit the entire dataset, as described in 502 Methods M4.2.

Out-of-sample prediction of the number of deaths (Extended Data Figure 2 ). In addition to 504 the number of confirmed cases, the NYT data also contains the daily reported number of deaths 505 due to COVID-19 by county. We use this death data as an additional source of validation. To 506 estimate the number of deaths N deaths , we use a similar process as for the number of cases N cases , 507 except that we replace r c with r d = 0.66%, the infection fatality rate for COVID-19, and δ c with 508 δ d = 432 hours (18 days), the number of days between becoming infectious and passing away 509 (Extended Data Table 2 ). This gives

Because we assume that deaths occur δ d = 18 days after individuals become infectious, we com-511 pare with NYT death data starting on March 19, 2020. 512 26 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint Extended Data Figure 2a demonstrates that the calibrated models also fit death counts sur-513 prisingly well, even though their parameters are selected to minimize RMSE in predicting cases, 514 not deaths. In some MSAs, the model fits the death data less well; this is unsurprising, because 515 our case and death count predictions assume constant case detection rates and fatality rates across 516 MSAs.

Comparison to baseline that does not use mobility data. To determine whether mobility data 518 aids in modeling case and death counts, we compare to a baseline SLIR model that does not 519 use mobility data and simply assumes that all individuals within an MSA mix uniformly. In this 520 baseline, an individual's risk of being infected and transitioning to the exposed state at time t is 521

where I (t) is the total number of infectious individuals at time t, and N is the total population size 522 of the MSA. As above, we performed a grid search over β base and p 0 , and calibrated the models would have happened if we changed the magnitude or timing of mobility reduction, we modify the 534 real mobility networks from March 1-May 2, 2020, and then run our models on the hypothetical 535 data. In Figure 2a , we report the cumulative incidence proportion at the end of the simulation (May 536 2, 2020), i.e., the total fraction of people in the exposed, infectious, and removed states at that time.

To simulate a smaller magnitude of mobility reduction, we interpolate between the mobility 538 27 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint network from the first week of simulation (March 1-7, 2020), which we use to represent typical 539 mobility levels (prior to mobility reduction measures), and the actual observed mobility network 540 for each week. Let W (t) represent the observed visit matrix at the t-th hour of simulation, and let 541 f (t) = t mod 168 map t to its corresponding hour in the first week of simulation, since there are 542 168 hours in a week. To represent the scenario where people had committed to α ∈ [0, 1] times 543 the actual observed reduction in mobility, we construct a visit matrixW

If α is 1, thenW

, and we use the actual observed mobility network for the simulation.

On the other hand, if α = 0, thenW 

If d is positive, this corresponds to starting mobility reduction d days later; if we imagine time on p j (Equation (9)). Then, we count the total expected number of infections per POI 565 by summing over hours. In Figure 2b , we sort the POIs by their expected number of infections and 566 report the proportion of all infections caused by the top x% of POIs.

Reducing mobility by clipping maximum occupancy (Figure 2c , Extended Data Figure 4) . 568 We implemented two partial reopening strategies: one that uniformly reduced visits at POIs to a 569 fraction of full activity, and the other that "clipped" each POI's hourly visits to a fraction of the 

where, as above, f (t) = t mod 168. This function leaves t unchanged if there is observed mo-579 bility data at time t, and otherwise maps t to the corresponding hour in the first week of our 580 simulation.

To simulate a reopening strategy that uniformly reduced visits to an γ-fraction of their origi-582 nal level, where γ ∈ [0, 1], we constructed the visit matrix

where R represents the first hour of reopening (May 1, 2020, 12AM) . In other words, we use the 584 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. β whose i, j-th entry is

This corresponds to the following procedure: for each POI p j and time t, we first check if t < R 593 (reopening has not started) or if V This reopening analysis is similar to the above analysis on clipping vs. uniform reopening.

As above, we set the reopening time R to May 1, 2020, 12AM. To simulate reopening a POI 610 category, we take the set of POIs in that category, V, and set their activity levels after reopening to 611 30 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 

As in the above reopening analysis, f (t) maps t to the corresponding hour in the first week of Average transmission rate of a POI category (Figure 3e) . We compute the average transmis- 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. Notation. We use a hat to denote quantities that we read directly from SafeGraph data, and r 639 instead of t to denote time periods longer than an hour. 

Estimating the POI marginals V (t) . We estimate the POI marginals V (t) ∈ R n , whose j-th

p j represents our estimate of the number of visitors at POI p j (from any CBG) at time t.

The number of visitors recorded at POI p j at hour t in the SafeGraph data,V 

33 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . This correction factor is approximately 7, using population data from the most recent 1-year ACS 684 (2018).

Estimating the CBG marginals U (t) . Next, we estimate the CBG marginals U (t) ∈ R m . Here, 686 the i-th element U (t) c i represents our estimate of the number of visitors leaving CBG c i (to visit 687 any POI) at time t. We will also use N c i ; recall that N c i is the total population of c i , which is 688 independent of t. 689 We first use the POI marginals V (t) to calculate the total number of people who are out 690 visiting any POI from any CBG at time t, 691 N (t)

Since the total number of people leaving any CBG to visit a POI must equal the total number of 692 people at all the POIs, we have that N (t)

Next, we estimate the number of people from each CBG c i who are not at home at time t as 694ĥ (t) c i N c i . In general, the total number of people who are not at home in their CBGs, iĥ To correct for this discrepancy, we assume that the relative proportions of POI visitors com-700 ing from each CBG follows the relative proportions of people who are not at home in each CBG. 701 We thus estimate U 

where N c i is the total population of CBG i, as derived from US Census data. This construction 704 ensures that the POI and CBG marginals match, i.e., N (t)

Iterative proportional fitting procedure (IPFP). IPFP is a classic statistical method 34 for ad-706 justing joint distributions to match pre-specified marginal distributions, and it is also known in the 707 34 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint literature as biproportional fitting, the RAS algorithm, or raking. 62 In the social sciences, it has 708 been widely used to infer the characteristics of local subpopulations (e.g., within each CBG) from 709 aggregate data. 63-65 710 We estimate the visit matrix W (t) by running IPFP on the aggregate visit matrixW , the 711 CBG marginals U (t) , and the POI marginals V (t) constructed above. Our goal is to construct a 712 non-negative matrix W (t) ∈ R m×n whose rows sum up to the CBG marginals U (t) ,

and whose columns sum up to the POI marginals V (t)

but whose distribution is otherwise "as similar as possible", in the sense of Kullback-Leibler di-715 vergence, to the distribution over visits induced by the aggregate visit matrixW .

Algorithm 1: Iterative proportional fitting procedure to estimate visit matrix W (t) Input: Aggregate visitsW ∈ R m×n CBG marginals U (t) ∈ R m ; POI marginals V (t) ∈ R n Number of iterations τ max Initialize W (t,0) =W for τ = 1, . . . , τ max do if τ is odd then for i = 1, . . . , m do

// Compute scaling factor for row i W (t,τ ) i,:

// Compute scaling factor for col j W 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint IPFP is an iterative algorithm that alternates between scaling each row to match the row 717 (CBG) marginals U (t) and scaling each column to match the column (POI) marginals V (t) . We 718 provide pseudocode in Algorithm 1. For each value of t used in our simulation, we run IPFP 719 separately for τ max = 100 iterations. Note that IPFP is invariant to scaling the absolute magnitude 720 of the entries inW , since the total number of visits it returns is fixed by the sum of the marginals; 721 instead, its output depends only on the distribution over visits inW .

The notion of similarity invoked above has a maximum likelihood interpretation: if IPFP 723 converges, then it returns a visit matrix W (t) whose induced distribution minimizes the Kullback-724 Leibler divergence to the distribution induced byW . 66 We further discuss the convergence of IPFP 725 in our setting in SI Section S3.

Here, we describe how we estimate the dwell time d p j , which we use to estimate the hourly mean 728 occupancy at each POI p j . For simplicity, we treat d p j as independent of t, i.e., we average across 729 all times t in our data for each POI p j .

Quantities from SafeGraph data. To compute the average timeδ p j spent at each POI p j , we 731 average the values in the median dwell field in the Patterns datasets from 2020.δ p j is measured 732 to minute-level resolution and expressed in units of hours, e.g.,δ If a visit straddles multiple hour boundaries at a POI, SafeGraph treats it as multiple visits in 735 each of those hours for the purposes of computing the visit counts that we use in Methods M6.1.

As a hypothetical example, consider a POI p 1 which has 1 new visitor come at the start of every 737 hour and stay for exactly 1 hour, and another POI p 2 which also has 1 new visitor come at the start 738 of every hour, except that visitors to p 2 stay for exactly 2 hours. SafeGraph data will reflect this 739 difference in visit times, withδ p 1 = 1 andδ p 2 = 2. However, SafeGraph will also record that p 2 740 has twice as many visitors at every hour than p 1 -in the notation of Methods M6.1, we would have 741 that V CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint Estimating the dwell time fraction d p j . Our goal is to estimate a correction factor for each 747 POI p j that corrects for the average visit duration of visitors to p j . As the example above shows, 748 setting the correction factor to be directly proportional to the average visit timeδ p j would overcount 749 visitors who stay across multiple hour boundaries. Instead, we define the dwell time fraction 750 d p j ∈ [0, 1] as the average fraction of an hour that a visitor to POI p j at any hour will spend 751 there. In other words, conditioned on a visitor being at p j at some time within an hour t, d p j is the 752 expected fraction of the hour t that the visitor physically spends at p j .

To estimate d p j , we make two assumptions: first, that every visitor to p j stays for exactly 754δ p j hours, and second, that the arrival times of visitors are uniformly distributed over all possible 755 arrival times. Concretely, if a visitor to p j stays forδ p j hours and is recorded as present at hour t, 756 then we assume that they are equally likely to have arrived at any time from [t −δ p j , t + 1). From 757 these assumptions, we can calculate

We truncate the departure time at t + 1 because any time spent after t + 1 does not count towards 759 the hour t. Similarly, we truncate the arrival time at t. This expression for d p j simplifies into

See SI Section S4 for the derivation. 37 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint Extended Data Figure 3 : A small fraction of POIs account for a large fraction of the predicted infections at POIs. We ran our models on the observed mobility data from March 1-May 2, 2020 and recorded the number of infections that occurred at each POI. Shaded regions denote 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations. See Methods M5 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. Extended Data Figure 4 : The predicted increase in infections with the "clipping" partial reopening strategy. We simulate reopening starting on May 1, 2020 and run the simulation until the end of the month. Each dot represents a degree of clipping, e.g., clipping at 50% maximum occupancy, at 20% maximum occupancy, etc. The y-coordinate of each dot represents the predicted number of new infections incurred after reopening (per 100k population) and its x-coordinate represents the fraction of visits lost from partial reopening compared to full reopening. In 6 MSAs, the lower-income CBGs incur substantially more infections from reopening. Only in New York City (NYC) is this trend reversed; this is because such a high fraction of lower-income CBGs in NYC had been infected before reopening (62%) that after reopening, there was only a minority of the lower-income population that is still susceptible. In comparison, none of the other MSAs saw such a high incidence proportion among lower-income CBGs before reopening; for example, the second highest was 31% for Philadelphia, and the rest ranged from 2-13%. Shaded regions denote 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations. See Methods M5 for reopening details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 48 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . # visitors from CBG c i to POI p j at time t Variable (SafeGraph) a p j area of POI p j in square feet Variable (SafeGraph) p 0 initial proportion of latent population 

Variable (Estimated) S (0) c i initial susceptible population in CBG c i (1 − p 0 )N c i E (0) c i initial latent population in CBG c i p 0 N c i I (0) c i initial infected population in CBG c i 0 R (0) c i initial removed population in CBG c i 0 Extended Data

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. To assess the reliability of the SafeGraph datasets, we measured the correlation between mobility 765 trends according to SafeGraph versus Google. 51 Google provides a high-level picture of mobility 766 changes around the world for several categories of places, such as grocery stores or restaurants. 767 We analyzed four of the six categories defined by Google: Retail & recreation (e.g., restaurants, 768 shopping centers, movie theaters), Grocery & pharmacy (e.g., grocery stores, farmers markets, 769 pharmacies), Parks (e.g., local parks, national parks, public beaches), and Residential (i.e. places 770 of residence). We omitted Transit stations because they are not well-covered by SafeGraph POIs, 771 and Workplaces because we do not model whether people are at work. To account for the first three 772 categories, we used POI visits in the SafeGraph Patterns datasets, identifying POIs in each category 773 based on their 6-digit NAICS codes (Table S5) . For the Residential category, we used SafeGraph 774 Social Distancing Metrics, which provides daily counts of the number of people in each CBG who 775 stayed at home for the entire day.

For each US region and category, Google tracks how the number of visits to the category 777 has changed over the last few months, compared to typical levels of activity before SARS-CoV-2.

To set this baseline, they compute the median number of visits to the category for each day of the 779 week, over a 5-week span from January 3-February 6, 2020. For a given day of interest, they then 780 compute the relative change in number of visits seen on this day compared to the baseline for the 781 corresponding day of week. We replicated this procedure on SafeGraph data, and compared the 782 results to Google's trends for Washington DC and 14 states that appear in the MSAs that we model.

For each region and category, we measured the Pearson correlation between the relative change in Table S6 . The Pearson correlations are high for all categories aside from Parks. Since CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint only covers 2.5%. Overall, these results indicate that SafeGraph mobility trends recapitulate those 794 in Google data, providing a validation of the reliability of SafeGraph data.

S2 Plausibility of predicted racial/socioeconomic disparities 796 To assess the plausibility of the predicted disparities in infection rates in Figure 3 , we compared 797 the model's predicted racial disparities to observed racial disparities in mortality rates. (Data on 798 socioeconomic disparities in mortality was not systematically available on a national level.) The 799 racial disparities in Figure 3 are generally of the same magnitude as reported racial disparities in 800 mortality rates-for example, the overall reported black mortality rate is 2.4× higher than the white 801 mortality rate, 71 which is approximately the same as the median racial disparity across MSAs of 802 2.4× that our model predicts (Figure 3b ). However, we note that this is an imperfect comparison 803 because many factors besides mobility contribute to racial disparities in death rates.

In addition, we observed that our model predicted unusually large socioeconomic and racial 805 disparities in infection rates in the Philadelphia MSA. To understand why the model predicted 806 such large disparities, we inspected the mobility factors discussed in the main text; namely, how 807 much each group was able to reduce their mobility, and whether disadvantaged groups encountered 808 higher transmission rates at POIs.

First, we find in Philadelphia that higher-income CBGs were able to reduce their mobility 810 substantially more than lower-income CBGs (Extended Data Figure 6 left). The CBGs with the 811 greatest percentage of white residents were also able to reduce their mobility more than the CBGs 812 with the lowest percentage of white residents (Extended Data Figure 6 right). These gaps are 813 noticeable, but not obviously larger than those in other MSAs. The key to Philadelphia's outlier 814 status seems to lie in the comparison of transmission rates. Within the same category of POI-815 e.g., full-service restaurants-individuals from lower-income CBGs tend to visit POIs with higher 816 transmission rates than individuals from high-income CBGs (Table S3 ). This is particularly true for Table S4 .) The 822 transmission rates encountered by individuals from lower-income CBGs in Philadelphia are often 823 dramatically higher than those encountered by higher-income CBGs; for example, up to 11.8× 824 52 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint higher for grocery stores. Digging deeper, this is because the average grocery store visited by 825 lower-income CBGs has 4× the number of hourly visitors per square foot, and visitors tend to stay 826 twice as long. Furthermore, Philadelphia's large discrepancy in density between lower-income and 827 higher-income POIs in SafeGraph data is consistent with Census data, which shows that the corre-828 lation between population density and median household income is larger in Philadelphia than in 829 any of the other MSAs that we examine (Spearman correlation 0.55 in Philadelphia, as compared 830 to a median of 0.31 across MSAs).

Since there are many other factors of disparity that we do not model, we do not place too 832 much weight on our model's prediction that Philadelphia's disparities will be larger than those 

Error in row marginals

Error in column marginals ,

which sums up the errors in the row (CBG) and column (POI) marginals of the visit matrix W (t,τ ) 842 from the τ -th iteration of IPFP. Each iteration of IPFP monotonically reduces this L 1 -error E (t,τ ) , 843 i.e., E (t,τ ) ≥ E (t,τ +1) for all τ ≥ 0. 72 In other words, the row and column sums of W (t,τ ) (which 844 is initialized as W (t,0) =W ) progressively get closer to (or technically, no further from) the target 845 marginals as the iteration number τ increases. Moreover, IPFP maintains the cross-product ratios 846 of the aggregate matrixW , i.e.,

for all matrix entries indexed by i, j, k, , for all t, and for all iterations τ . 848 53 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint IPFP converges to a unique solution, in the sense that W (t) = lim τ →∞ W (t,τ ) , if there exists 849 a matrix W (t) that fits the row and column marginals while maintaining the sparsity pattern (i.e., 850 location of zeroes) ofW . 72 If IPFP converges, then the L 1 -error also converges to 0 as τ → ∞, 72 851 and W (t) is the maximum likelihood solution in the following sense. For a visit matrix W = {w ij }, 852 let P W represent a multinomial distribution over the mn entries of W with probability proportional 853 to w ij , and define U (t) ⊆ R m×n + and V (t) ⊆ R m×n as the set of non-negative matrices whose row 854 and column marginals match U (t) and V (t) respectively. Then, if IPFP converges, 855

where KL (p q) is the Kullback-Leibler divergence KL (p q) = E p log p(x) q(x) . In other words, IPFP 856 returns a visit matrix W (t) whose induced distribution P W (t) is the I-projection of the aggregate 857 visit distribution PW on the set of distributions with compatible row and column marginals. 66 In 

ij .

Note that in this scenario, IPFP still monotonically decreases the L 1 -error. 72

In our implementation (Algorithm 1), we take τ max = 100, so IPFP ends by fitting the column 869 (POI) marginals. This ensures that our visit matrix W (t) is fully compatible with the POI marginals

while still minimizing the L 1 -error E (t,τ ) with respect to the CBG marginals U (t) . Empirically, we 872 find that τ max = 100 iterations of IPFP are sufficient to converge to this oscillatory regime.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

This expression for d p j simplifies into

To see this, first consider the case whereδ p j ≤ 1. To keep notation simple, let t = 0 mark the start 877 of the hour being considered. Then,

• Visitors who arrive at −δ p j will spend 0 time during the hour being considered. This time 879 increases linearly from 0 toδ p j as the arrival time increases from −δ p j to 0.

• Visitors who arrive during [0, 1−δ p j ] will spend the fullδ p j time in the hour being considered.

• This time decreases linearly fromδ p j to 0 as the arrival time increases from 1 −δ p j to 1.

In total, we thus have

Similarly, consider the case whereδ p j > 1. Then,

• Visitors who arrive at −δ p j will spend 0 time during the hour being considered. This time 885 increases linearly from 0 to 1 as the arrival time increases from −δ p j to 1 −δ p j .

• Visitors who arrive during [1 −δ p j , 0] will spend the full hour.

• This time decreases linearly from 1 to 0 as the arrival time increases from 0 to 1.

In total, we likewise have

55 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. Table S2 : What if the magnitude of mobility reduction changed? Each column represents a counterfactual scenario where the magnitude of mobility reduction is only a some percentage of the observed mobility reduction, i.e., 0% corresponds to no mobility reduction, and 100% corresponds to the real, observed level of mobility reduction. We report the expected ratio of the number of infections predicted under the counterfactual to the number of infections predicted using observed mobility data; a ratio lower than 1 means that fewer infections occurred under the counterfactual. The numbers in parentheses indicate the 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations. See Methods M5 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. 

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Table S7 : Model parameters used for each MSA. # sets counts the number of parameter sets that are within 20% of the RMSE of the best-fit parameter set, as described in Section M4. For each of β base , ψ, and p 0 , we show the best-fit parameter set and, in parentheses, the corresponding minimum and maximum within the 20% threshold.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint

Cumulative infections (per 100k) at category Figure S1 : For each POI category, we plot the predicted cumulative number of infections (per 100k population) that occurred at that category for CBGs in the bottom-(purple) and top-(gold) income deciles. Shaded regions denote 2.5th and 97.5th percentiles across sampled parameters and stochastic realizations.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Completely at home (SafeGraph Data) NEW YORK Figure S3 : Google versus SafeGraph mobility trends for New York state. The y-axis represents mobility levels compared to baseline activity in January and February 2020. For the categories from left to right, the Pearson correlations between the datasets are 0.96, 0.76, 0.57, and 0.91. See SI Section S1 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. Figure S5 : POI attributes in Chicago. See Figure S4 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 17, 2020. Figure S6 : POI attributes in Dallas. See Figure S4 for details. Figure S7 : POI attributes in Houston. See Figure S4 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. Figure S8 : POI attributes in Los Angeles. See Figure S4 for details. Figure S9 : POI attributes in Miami. See Figure S4 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. Figure S10 : POI attributes in New York. See Figure S4 for details. Figure S11 : POI attributes in Philadelphia. See Figure S4 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. Figure S12 : POI attributes in San Francisco. See Figure S4 for details. Figure S13 : POI attributes in Washington DC. See Figure S4 for details.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted June 17, 2020. . https://doi.org/10.1101/2020.06.15.20131979 doi: medRxiv preprint

Stay-at-home orders across the country

As states reopen, governors balance existing risks with new ones

Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China

Clustering and superspreading potential of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infections in Hong Kong

Full genome viral sequences inform patterns of SARS-CoV-2 spread into and within Israel. medRxiv (2020)

Coronavirus Disease Outbreak in Call Center

Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, China: a retrospective cohort study

These Graphs Show How COVID-19 Is Ravaging New York City's Low-Income Neighborhoods

Hospitalization Rates and Characteristics of Patients Hospitalized with Laboratory-Confirmed Coronavirus Disease 2019 -COVID-NET, 14 States

Class and COVID: How the less affluent face double risks. The Brookings Institution (2020)

Ethnicity and COVID-19: an urgent public health research priority

COVID-19 exacerbating inequalities in the US

COVID-19 and African Americans

Racial Health Disparities and Covid-19-Caution and Context

Projections for first-wave COVID-19 deaths across the US using social-distancing measures derived from mobile phones. medRxiv (2020)

A cell phone data driven time use analysis of the COVID-19 epidemic. medRxiv (2020)

State-level tracking of COVID-19 in the United States (2020)

Population flow drives spatio-temporal distribution of COVID-19 in China

Differential Effects of Intervention Timing on COVID-19 Spread in the United States. medRxiv (2020)

Effect of non-pharmaceutical interventions to contain COVID-19 in China

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2)

Initial Simulation of SARS-CoV2 Spread and Intervention Effects in the Continental US. medRxiv (2020)

Modeling the impact of social distancing, testing, contact tracing and household quarantine on second-wave scenarios of the COVID-19 epidemic. medRxiv (2020)

COVID-19: How to Relax Social Distancing If You Must. medRxiv (2020)

Social network-based distancing strategies to flatten the COVID-19 curve in a post-lockdown world

Adaptive cyclic exit strategies from lockdown to suppress COVID-19 and allow economic activity. medRxiv (2020)

Mapping county-level mobility pattern changes in the united states in response to covid-19

Available at networkscienceinstitute.org/publications/assessing-changes-in-commuting-and-individualmobility-in-major-metropolitan-areas-in-the

Rationing social contact during the COVID-19 pandemic: Transmission risk and social benefits of US locations

Human Mobility in Response to COVID-19 in

SafeGraph. Measuring and Correcting Sampling Bias in Safegraph Patterns for More Accurate Demographic Analysis

Places Manual (2020)

Discrete multivariate analysis

Synthesis-a synthetic spatial information system for urban and regional analysis: methods and examples

The reliability of using the iterative proportional fitting procedure. The Professional Geographer

Combining sample and census data in small area estimates: Iterative proportional fitting with standard software

I-divergence geometry of probability distributions and minimization problems. The Annals of Probability

Estimating unobserved SARS-CoV-2 infections in the United States. medRxiv (2020)

Average detection rate of SARS-CoV-2 infections has improved since our last estimates but is still as low as nine percent on March

The unseen and pervasive threat of COVID-19 throughout the US. medRxiv (2020)

Estimates of the severity of coronavirus disease 2019: a model-based analysis

The color of coronavirus: COVID-19 deaths by race and ethnicity in the

Biproportional scaling of matrices and the iterative proportional fitting procedure