key: cord-0837001-51r3hm9m
authors: Capistran, M. A.; Capella, A.; Christen, J. A.
title: Forecasting hospital demand in metropolitan areas during the current COVID-19 pandemic and estimates of lockdown-induced 2nd waves
date: 2020-07-18
journal: nan
DOI: 10.1101/2020.07.16.20155721
sha: 72af52adce7c6f6728a193241b883416401ea8ae
doc_id: 837001
cord_uid: 51r3hm9m

We present a forecasting model aim to predict hospital occupancy in metropolitan areas during the current COVID-19 pandemic. Our SEIRD type model features asymptomatic and symptomatic infections with detailed hospital dynamics. We model explicitly branching probabilities and non-exponential residence times in each latent and infected compartments. Using both hospital admittance confirmed cases and deaths, we infer the contact rate and the initial conditions of the dynamical system, considering breakpoints to model lockdown interventions and the increase in effective population size due to lockdown relaxation. The latter features let us model lockdown-induced 2nd waves. Our Bayesian approach allows us to produce timely probabilistic forecasts of hospital demand. We have applied the model to analyze more than 70 metropolitan areas and 32 states in Mexico.

The ongoing COVID-19 pandemic has posed a major challenge to public health systems 2 in many countries with the imminent risk of saturated hospitals and patients not 3 receiving proper medical care. Although the scientific community and public health 4 authorities had insight regarding the risks and preparedness measures required at the 5 onset of a zoonotic pandemic, our knowledge of the fatality and spread rates of 6 COVID-19 remains limited [1] [2] [3] [4] . In terms of disease handling, two leading issues 7 determining the current situation are the lack of pharmaceutical treatment and our 8 inability to estimate the extent of the asymptomatic infection of COVID-19 [5] [6] [7] . 9 Under current circumstances, control measures reduce new infections by limiting the 10 number of contacts through mitigation and suppression [1] . Mitigation includes social 11 distancing, testing, tracing, and isolating infected individuals, while suppression imposes 12 temporary cancellation of non-essential activities. Mitigation and suppression pose a SARS [10] , Zika [11] , Ebola [12] , etcetera. However, surveillance data during a pandemic 21 event often suffer from serious deficiencies such as incompleteness and backlogs. 22 Another critical issue is the design of data acquisition, taking into account geographical 23 granularity [13] . Epidemic surveillance of COVID-19 is no different since there is an 24 unknown percentage of asymptomatic infections, and susceptibility is related to 25 economic vulnerability. 26 Undoubtedly, one key task during the early pandemic response efforts is using 27 epidemiological records and mathematical and statistical modeling to forecast excess 28 hospital care demand with formal quantified uncertainty. 29 In this paper, we pose a compartmental SEIRD model that considers both 30 asymptomatic and symptomatic infection, including hospital dynamics. We model the 31 residence time in each latent and infected compartments explicitly [14, 15] , and we use 32 records of daily confirmed cases and deaths to pose a statistical model that accounts for 33 data overdispersion [16, 17] . Furthermore, we use Bayesian inference to estimate the 34 initial state of the governing equations, the contact rate, and a proxy of the population 35 size to make probabilistic forecasts of the required hospital beds, including the number 36 of intensive care units. The model output has been used by Mexican public health 37 authorities to assist decision making during the COVID-19 pandemic in more than 70 38 metropolitan areas and the country's 32 states. 39 Contributions and limitations 40 • We developed a model to produce accurate midterm (several weeks) probabilistic 41 forecasting of COVID-19 hospital pressure, namely hospital beds and respiratory 42 support or mechanical ventilation demands, using confirmed records of cases at 43 hospital admittance and deaths. 44 • Our model accounts for policy changes in control measures, such as school 45 closures [18] and lockdowns, as breakpoints in the transmission rates. 46 • Assuming a given fraction of asymptomatic individuals, we infer changes in the 47 transmission rate and the effective population size before and after a given 48 lockdown-relaxation day. 49 • Inferred changes in effective population size allows us to produce a forecast of 50 lockdown-induced 2nd waves. 51 Since asymptomatic infection is not fully understood so far [19] , the fraction of 52 asymptomatic individuals is yet unknown. Therefore:

• The effective population size is only a proxy, and its absolute value is not 54 meaningful, but only its relative value before and after a relaxation day.

• Without serological studies in the open population -ideally after an outbreak-it 56 is impossible to forecast the population fraction that will be in contact with the 57 virus by the end of the current outbreak.

• At this point, our model does not address next pandemic outbreaks beyond 59 lockdown-induced 2nd waves.

• Finally, although the model does not account explicitly for biases due to 61 behavioral changes [20, 21] , population clustering and super spreading events [22] , 62 we argue that our approach to lockdowns and relaxation events is a proxy model 63 of these more general events. Other COVID-19 models have been used to explore exit strategies [29, 30] , the role of 83 recovered individuals as human shields [31] , digital contact tracing [32] , break points in 84 the contact rate to account for changes in suppression and mitigation policies [18] and 85 lockdown-induced 2nd COVID waves [33] under the assumption that population is 86 temporally geographically isolated.

"Models should not be presented as scientific truth" [34] . Indeed, models are tools 89 intended to serve a specific purpose, evaluate or forecast particular aspects of 90 phenomena and ideally should be developed following the processes of predictive 91 science [35] . Namely, identify the quantities of interest (QoI), verify the computational 92 and mathematical approximation error, including their implication in the estimation of 93 QoI, and calibrate the parameters to adjust the model in light of data to bring it closer 94 to experimental observation. When considering uncertainty, Bayesian inference may be 95 used to calibrate some key model features given data. Finally, a validation process must 96 be used to build confidence in the accuracy of the QoI predictions. Our model is built 97 out of three interrelated components; a law for dynamics, a law for uncertainty, and the 98 choice of parameters.

Dynamical model 100 As a proxy of hospital pressure, the quantities of interest in our model are the evolving 101 demand of ICU/respiratory-support beds and no-ICU hospital beds. We developed a infectious and a proportion f of them become sufficiently severe symptomatic cases (I S ) to approach a hospital, while the remaining cases become mild-symptomatic to asymptomatic (I A ). The asymptomatic/mild-symptomatic cases remain infectious a mean time of 1/γ 1 days and eventually recover. For the symptomatic cases (I S ) we assume that after an average time of 1/σ 2 days a proportion g of infected individuals will need hospitalization (H 1 ), while the rest (I C ) will receive ambulatory care, recovering after an average convalescent time of 1/γ 2 days in quarantine. The hospitalized patients (H 1 ) remain an average time of 1/σ 3 days until a fraction h will need assisting respiratory measures or ICU care such as mechanical ventilation (U 1 ). The remaining fraction 1 − h of hospitalized patients (H 2 ) will recover after 1/γ 3 days in average. Respiratory-assisted/ICU patients (U 1 ) remain in that state an average of 1/σ 4 days until a critical day is reached when a proportion i of them will die (D) and the remaining proportion 1 − i will recover (H 3 ) after an average period of 1/γ 4 days. Similar models have been proposed by [2, 31, 32, 36] . For the infection force (λ) we assume that only mild-symptomatic/asymptomatic (I A ) and symptomatic (I S ) individuals spread the infection, that is

where β A and β S are the contact rates of asymptomatic/mild-symptomatic and 109 symptomatic individuals, respectively.

Parameters and observational model

The model has two kinds of parameters that have to be calibrated or inferred; the ones 112 related to COVID-19 disease and hospitalization dynamics (such as residence times and 113 proportions of individuals that split at each bifurcation of the model) and those 114 associated with the public response to mitigation measures such as the contact rates β's 115 and the effective population size N ef f during the outbreak. Some of these parameters 116 can be estimated from hospital records or found in recent literature or inferred from 117 reported cases and deaths, but some remain mostly unknown. In the latter category, we 118 have N ef f and the fraction 1 − f of asymptomatic/mild-symptomatic infections.

Reported values of the proportion of asymptomatic/ mild-symptomatic infections cases 120 1 − f range from 10% to 75%, and even 95% in children population [6, 7, 37] interventions, a breakpoint is established at which β = β 1 before and β = β 2 after the 145 intervention day. This creates a non-linear time-dependent β(t) [18, 38] . In the same 146 fashion, further intervention days may be included by adding more change points and β 147 parameters. These additional parameters are then included in the inference.

Assuming that effective population size N ef f and transmission rates are fixed,

SEIRD type models converge to the attractor E = 0, I = 0, i.e. the system models an 150 epidemic that dies out after one single acme. Even for sensible non-constant 151 transmission rates, these kinds of models are only able to produce single epidemic 152 outbreaks. Therefore, in order to be able to estimate secondary outbreak waves after 153 lockdown relaxation measures, one necessarily needs to estimate a different N ef f before 154 and after relaxation days. An increase of N ef f is aimed as a proxy to the rise in the 155 average connectivity between population clusters after relaxation days. Our approach 156 here is as follows: we include the new parameter(s) ω i ∈ (0, 1) and set

where N is the total population of the metropolitan area or region under study. Next, 158 we postulate a value of f fixed and estimate ω 1 and ω 2 before and after the relaxation 159 day, respectively. To add flexibility to our method and model possible changes in 160 population behavior, we also force a shift from β 1 to β 2 before and after a relaxation 161 day. With this procedure, given f we are able to estimate N ef f in terms of the ω's.

Nevertheless, due to the confounding effect of the product

are still unable to estimate the real N ef f 's until the real value of f is known. Relaxation 164 days may include both lockdown relaxation day or mayor changes in the public response 165 to mitigation measures (see SM for a methodology on how to choose relaxation days). 166

Setting lockdown and relaxation days 167 As explained before, we model interventions and relaxation days as discontinuities in β 168 and (β, ω), respectively. We consider a lockdown day on 22 March 2020, where a 169 country-wide lockdown started. In Mexico City, we include a second intervention day to 170 model a further local intervention in early April. The methodology to set the relaxation 171 days is as follows: we computed R t following [39] and look for days where local 172 minimums occur. We choose the relaxation days judiciously, keeping them at least three 173 weeks apart to have enough data to perform the inference. A rise in R t can be a 174 consequence of an increase in contact rates, effective population size, or both. In either 175 case, our model captures these changes, that serve as a model for the connectivity 176 changes in the outbreak cluster structure.

Observational model and data 178 To make our inferences, we use both confirm cases and deceased counts. In some 179 regions, sub reporting of COVID-19 related deaths may become relevant, especially in 180 places hit by a severe outbreak [40] . We consider daily deaths counts d i and its theoretical expectation that is estimated in terms of the dynamical model as µ D (t i ) = D(t i ) − D(t i−1 ) for the metropolitan area or region being analyzed. Analogously, we consider daily cases c i and its corresponding µ c (t i ) given by the daily flux entering the H 1 compartment, which may be calculated as in [17] , namely

where I S m (t) is the last state variable in the I S Erlang series. We calculate the above 194 integral using a simple trapezoidal rule with 10 points (1/10 day).

Bayesian Inference

In order to carry out a likelihood-based analysis, we assume that epidemic data has 197 more variation than implied by a standard Poisson process, as is the case in other 198 ecological studies. Following [16] , we postulate that the number of both the confirmed 199 cases and deaths follows a negative binomial distribution N B. Denoting the mean and 200 variance as µ and σ 2 , and requiring that σ 2 = ωµ + θµ 2 > µ we enforce overdispersion 201 for suitable chosen parameters ω and θ. For data y i , namely c i and d i , we reparametrize 202 the negative binomial distribution and let y i ∼ N B(pµ(t i ), ω, θ), with fixed values for 203 the overdispersion parameters ω, θ and an additional reporting probability p (see SM for 204 further details). 205 We assume conditional independence in the data, and therefore from the NB model, 206 we obtain a likelihood. Our parameters are the contact rate parameter β's, the ω's and 207 July 17, 2020 6/27

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. . https://doi.org/10.1101/2020.07.16.20155721 doi: medRxiv preprint crucially we also infer the initial conditions E(0), I A (0), I S (0). Letting 208 S(0) = N − (E(0) + I A (0) + I S (0)) and setting the rest of the parameters to zero, we 209 have all initial conditions defined and the model may be solved numerically to obtain 210 µ D and µ c to evaluate our likelihood. 211 Finally, regarding the elicitation of the parameters prior distribution, we use Gamma 212 distributions with scale 1 and shape parameter 10 to model the initial conditions 213 E(0), I A (0), I S (0) of the community transmission. Following the argumentation of Cori 214 et al. [39] , we assume that local transmission starts when there are 10 confirmed cases. 215 The rationale for using Gamma distribution priors is that we can specify the 216 distribution by prescribing its first two moments, and the resulting distribution verifies 217 a maximum entropy condition. Namely, we obtain the less informative distribution that 218 has the prescribed mean and (log) variance [41] . The prior for the first transmission 219 rate β 0 , is a long tail, log Normal with σ 2 = 1 and scale parameter 1; that is 220 log(β 0 ) ∼ N (0, 1). For the subsequent β's, we use autoregressive priors to impose some 221 coherence from one change point to the next with log(β i ) ∼ N (log(β i−1 ), 1). The prior 222 on ω i is a Beta(1 + 1/6, 1 + 1/3) restricted to ω i > ω i−1 . This beta distribution is a 223 fairly flat near uniform density in [0, 1], touches zero in 0 and 1, and is slightly skewed 224 to lower values. We model here the unlikely values ω i = 0, 1 and that under current 225 social distancing measures we expect smaller rather than larger N ef f . Otherwise, the 226 prior for the ω i 's is rather diffuse and non-informative.

To sample from the posterior we resort to MCMC using the "t-walk" generic 228 sampler [42] . The MCMC runs semi-automatic, with consistent performances in most 229 data sets. For any state variable V , the MCMC allows us to sample from the posterior 230 predictive distribution for V (t i ). By plotting some of its quantiles sequentially, we may 231 produce predictions with quantified probabilistic uncertainty, as seen in Figure 3 . Notice that the forecast begins after June 20th, and an uncertainty cone opens to the 249 right the next 3 to 4 weeks. However, the attractor of the dynamical system closes the 250 cone for longer times, and the predictive power of our forecast decreases.

The estimate of total hospital bed occupancy corresponds to the daily integral of H 1 252 and H 2 in the model and the ICU occupancy corresponds to the U 1 daily integral. We 253 consistently overestimated the total number of hospital beds and ICU units and shifted 254 to earlier times. We calibrate residence times from reports on daily demand of hospital 255 beds and intensive care unit records from Instituto Mexicano del Seguro Social or (which was not certified by peer review) 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 July 18, 2020. All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. . https://doi.org/10.1101/2020.07.16.20155721 doi: medRxiv preprint 

We present a compartmental SEIRD model to make probabilistic forecasts of hospital 277 pressure during COVID-19 outbreaks in metropolitan areas. In this work, we also 278 consider lockdowns and lockdown-relaxations as two different kinds of interventions. We 279 model the former as a change in transmission rates while the latter also allows for 280 changes in the effective population size. These changes in effective population size and 281

July 17, 2020 9/27 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. . https://doi.org/10.1101/2020.07.16.20155721 doi: medRxiv preprint Daily R t 's calculated as in [39] for Cancun metropolitan area. Two relaxation days, marked with red vertical lines, were judiciously included at local minimums, allowing for a minimum gap of 3 weeks. transmission rates are used as a proxy of behavioral changes, changes in the 282 connectivity between population clusters and super spreading events. 283 We in-line monitor exogenous changes in SEIRD parameters to inform the model of 284 interventions and relaxation days, in the same spirit of our inferences in the 285 transmission rate and initial conditions. Our analysis showed that we account for the 286 outbreak evolution in many metropolitan areas by setting one lockdown on 22 May 287 2020, the beginning of country-wide lockdown, and two more relaxation days. One 288 around 10 May, mother's day, a popular celebration in Mexican culture where families 289 gather together. The other relaxation day is set around four to eight days after 1 June, 290 the announced date for the country-wide lockdown relaxation measures. In some cases, 291 we also impose somewhat different relaxation days to account for local changes, such as 292 the opening of tourist activity. Of note, all relaxation days where set exclusive by our 293 R t based methodology. The interpretation, as above, came afterward. In a prior version of this ama model -before reaching the first acme-we did not infer 295 the effective population size, namely ω. Instead, by modeling several countries that 296 already had passed the acme, we postulated ω = 1 and f = 0.05. This yields some 297 controversy that even reached the news media [43] . With the current model we can 298 show that the product ω × f , that is the quantity that affects the estimates, takes 299 July 17, 2020 10/27

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. . https://doi.org/10.1101/2020.07.16.20155721 doi: medRxiv preprint maximum a posteriori values between 0.04 and 0.12 in most of our studied cases, see 300 Figure 6 . This fact explains the reason why our forecast weeks ahead on the acme have 301 been rather accurate. It is also important to notice that this is evidence of a 302 confounding effect between ω, , and the transmission rates β on the early stages of the 303 outbreak. Many forecasting models that fail to recognize this fact will consequently fail 304 on their acme estimates. 305 We note that after the acme SEIRD type models converge to the attractor 306 E = 0, I = 0, and long term estimates tend to overestimate the decline of an outbreak. 307 The underlying assumption that social contact and other conditions remain constant is 308 not reasonable for most societies over long periods (e.g., more than eight weeks).

Conservative short to mid-term forecasts most be preferred and change points added 310 when necessary.

The age-independent model has proven to be adequate to produce accurate forecasts 312 for the hospitalization dynamics during current outbreaks. Therefore the added 313 complexity of the age-structure model may not justify its use at this point. However, we 314 continue to work in our age structure model.

Since our QoI are related to the hospital pressure, we choose all parameters 316 conservatively. However, Erlang densities can not correctly approximate residence times 317 of some cases real distributions and more general distributions should be considered.

Moreover, as health professionals learn to treat the disease, hospital residence times also 319 change. Both these effects should also be considered to obtain more accurate estimates. 320

Our observation model is designed to integrate data after the nonlinear term in the 321 flow diagram of the dynamic model (see Figure 1 ) and the rest of the dynamics is 322 proportional to the hospital occupancy curves, therefore the model forecasts can be used 323 as a proxy of the full outbreak. We follow this idea to predict the date when an 324 outbreak will reach its first peak as well as lockdown induced second waves.

The confounding effect between the population size, namely ω, and the fraction of All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. July 17, 2020 13/27

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. (which was not certified by peer review) 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 July 18, 2020. (which was not certified by peer review) 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 July 18, 2020. . https://doi.org/10.1101/2020.07.16.20155721 doi: medRxiv preprint S1.4 Acapulco All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. . https://doi.org/10.1101/2020.07.16.20155721 doi: medRxiv preprint S1.5 Culiacan 

We developed a dynamic transmission compartmental model to simulate the spread of the novel coronavirus SARS-CoV-2. A definition of the state variables is given in Table S1 . Additionally, an "Erlang series" is included for most of these state variables to account for non-exponential residence times. The model may be described conceptually with the graph in Figure 1 . Without showing the Erlang series for the July 17, 2020 17/27

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. . Hospitalized patients recovering after ICU or respiratory support R Recovered D Deceased sub-compartments the system of equations in the model is as follows:

In Table S2 , we give a brief description of all the parameters in the model. For the infection force (λ) we assume that individuals that spread the infection correspond to the mild-symptomatic/asymptomatic (I A ) and symptomatic individuals (I S ) before they get contact with the health care system or doctor, e.g.

We compute the basic reproductive number R 0 of the epidemic by the next generation matrix method [44] and obtain

Since our QoI are related to the hospital pressure we choose all parameters 369 conservatively. For each metropolitan area, we assume that N corresponds to its full 370 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) 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 July 18, 2020. To make the intrinsic generation-interval of the renewal equation in each compartment 374 more realistic we divide each compartment of the model into m equal sub-compartments 375 to generate an Erlang-distributed waiting time [14] . The Erlang distributions of each 376 compartment is calibrated by two parameters: the rate λ E and the shape m, a positive 377 integer that corresponds to the number of sub-compartments on the model. In terms of 378 these parameters the mean of the Erlang distribution is m/λ E , this mean correspond to 379 the average times in the dynamic model. 380 We use recent publications and information generously shared by the Instituto In our methodology, we aim to infer the force of the infection λ. This parameter is 387 defined in terms of contact rate of asymptomatic/mild-symptomatic individuals β A and 388 contact rate of symptomatic individuals β S . Due to the functional dependence of λ in 389 these parameters, there is a lack of identifiability between β A and β S that can not be 390 resolved without further assumptions. We assume that the relative strength between the 391 transmission rate of asymptomatic/mild-symptomatic and symptomatic is modeled as a 392 fixed ratio κ. We model the value of κ directly as the ratio of the viral load of 393 symptomatic and asymptomatic/mild-symptomatic patients [48, 49] Poisson distribution may be needed to correctly, and safely, model these types of data. 401 Following [16] (see main paper) the NB distribution is re parametrized in terms of its To sample from the posterior we resort to MCMC using the "t-walk" generic 434 sampler [42] . The MCMC runs semi-automatic, with a fairly consistent burn-in of 1,000 435 iterations (sampling initial values from the prior). We perform subsampling using the Adding group ages is straightforward in these types of models. The number a of age 450 groups is decided, and our model is repeated a times. Different residence times may be 451 included [50] but we preferred to concentrate on the different transition probabilities 452 g, h that vary nearly two orders of magnitude using age groups 453 [0, 25], (25, 50] , (50, 65] , (65, 100] [1] . The age structure is used to divide the initial 454 infectious and susceptible population proportional to each age group size. We infer the 455 same number of parameters, using a single β, with an optional weighting contact 456 matrix [51] to model different contact rates among age groups in specific regions. Using 457 a sufficiently flexible software design, progressing to an age-structure model is not 458 complex; nonetheless, the MCMC may run substantially slower. We have experimented 459 with our age-structured model using census data from Mexico and both uniform and 460 non-uniform contact matrices. However, we do not report any of these results here, 461 given that our non-age structure model has sufficiently enough predictive power, as 462 already discussed. To explain the confounding effect of N ef f × f we have two observations. First, if we let f =f /α for some α ∈ (0,f ) then differential equations for the variables I c , H 1 , H 2 , H 3 , U 1 , U 2 and D remain invariant and the equations for I s becomes dI S dt =f α σ 1 E − σ 2 I S .

By lettingẼ = E/α the equation for I s is also invariant with the substitution of E bỹ E. Now, the equation forẼ is given by

By lettingÑ ef f = αN ef f the latter equations becomes also invariant under the 465 substitution of E byẼ. Therefore for the lower branch in the model (see Figure ? ?) the 466 system of equations is invariant under the change of f and N ef f byf andÑ ef f 467 providedÑ ef f ×f = N ef f × f holds. We need to adapt the equations for S, I A , and R 468 to get a consistent system of equations. (which was not certified by peer review) 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 July 18, 2020. . Fig S6. Outbreak analysis for Mexico city metropolitan area, data until 15 May. Total number of recovered R; with N ef f = N , total population, and f = 0.05, R(∞) reaches approximately 17.5 · 10 6 while with N ef f = N/3, f = 0.05 · 3, R(∞) only reaches roughly 5.5 · 10 6 . However, the fit for cases and deaths and the predictive curves for hospital demand are identical (results not shown).

andf is considered. There is a range of validity for α where the inference of β does not 475 change, but we do not explore this property further. 476 We also present numerical simulations to confirm this confounding effect (see 477 Figure S6 ), but until the asymptomatic infection is fully described, it is not possible to 478 resolve this issue.

S7 Data sources 480 We take metropolitan areas delimitation and population from [52] and [53], respectively. Official records of COVID-19 confirmed cases and deaths are reported in [54] .

Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand. London: Imperial College COVID-19 Response Team

Estimates of the severity of COVID-19 disease. medRxiv

Zhonghua liu xing bing xue za zhi= Zhonghua liuxingbingxue zazhi, Number = 2, Pages = 145, Title = The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases (COVID-19) in China

Follow-up of asymptomatic patients with SARS-CoV-2 infection

Asymptomatic Transmission, the Achilles' Heel of Current Strategies to Control Covid-19

Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship

Covid-19: identifying and isolating asymptomatic people helped eliminate virus in Italian village

The 2009 pandemic in Mexico: Experience and lessons regarding national preparedness policies for seasonal and epidemic influenza

A lesson learned from the MERS outbreak in South Korea in 2015

The lessons of SARS

Delays in global disease outbreak responses: lessons from H1N1, Ebola, and Zika

The next epidemic-lessons from Ebola

Modeling and predicting human infectious diseases

Equivalence of the Erlang-distributed SEIR epidemic model and the renewal equation

Appropriate Models for the Management of Infectious Diseases

Using the negative binomial distribution to model overdispersion in ecological count data

Model selection for seasonal influenza forecasting. Infectious Disease Modelling

Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

Reconciling early-outbreak estimates of the basic reproductive number and its uncertainty: framework and applications to the novel coronavirus (SARS-CoV-2) outbreak. medRxiv

Systematic biases in disease forecasting-The role of behavior change

Moving Beyond a Peak Mentality: Plateaus, Shoulders, Oscillations and Other'Anomalous' Behavior-Driven Shapes in COVID-19 Outbreaks. medRxiv

Secondary attack rate and superspreading events for SARS-CoV-2. The Lancet

Projection of COVID-19 Cases and Deaths in the US as Individual States Re-open

Covid Act Now

Short-term forecasts of COVID-19 deaths in multiple countries

Forecasting the impact of the first wave of the COVID-19 pandemic on hospital demand and deaths for the USA and European Economic Area countries. medRxiv

Projecting hospital utilization during the COVID-19 outbreaks in the United States

Estimating the burden of SARS-CoV-2 in France

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

Expected impact of reopening schools after lockdown on COVID-19 epidemic in Ile-de-France

Modeling shield immunity to reduce COVID-19 epidemic spread

Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing

Modelling lockdown-induced secondary COVID waves in France. medRxiv

Predictive Mathematical Models of the COVID-19 Pandemic: Underlying Principles and Value of Projections

Computer predictions with quantified uncertainty, part I. SIAM News

Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The Lancet Infectious Diseases

Estimating cases for COVID-19 in South Africa Update

Epidemic models with nonlinear infection forces

A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics

Preliminary Estimate of Excess Mortality During the COVID-19 Outbreak

Derivation of some frequency distributions using the principle of maximum entropy (POME)

A general purpose sampling algorithm for continuous distributions (the t -walk)

As Official Toll Ignores Reality

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Updated rapid risk assessment from ECDC on the novel coronavirus disease 2019 (COVID-19) pandemic: increased transmission in the EU/EEA and the UK

Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside Hubei province, China: a descriptive and modelling study. The Lancet Infectious Diseases

Modellierung von Beispielszenarien der SARS-CoV-2-Ausbreitung und Schwere in Deutschland

SARS-CoV-2 viral load in upper respiratory specimens of infected patients

Temporal dynamics in viral shedding and transmissibility of COVID-19

Projecting hospital utilization during the COVID-19 outbreaks in the United States

Projecting social contact matrices in 152 countries using contact surveys and demographic data

Delimitación de las zonas metropolitanas de México

Datos Abiertos -Dirección General de Epidemiología

The authors wish to thank ElenaÁlvarez-Buylla, Paola Villarreal (CONACYT) and 339 Hugo López-Gatell (Secretaría de Salud) for their support and comments for improving 340 this forecast model and also Instituto Mexicano del Seguro Social (IMSS) for sharing