At the end of December of 2019 a novel coronavirus SARS-CoV-2, was detected in adults with pneumonia from the city of Wuhan, Hubei Province, China. This new virus rapidly spread from Hubei Province to all China [1] [2] [3] [4] [5] . By January 22th the virus was detected in South Korea, Taiwan, Thailand. Macao and USA. By March 12 of 2020, 128,343 cases of new virus have been detected in 116 countries and 63% of the cases has been register in Continental China. In February 26, Brazil was the first country of Latin America to detect this virus and Mexico detect his first case one day later [6] . CoVID-19 is a respiratory disease caused by the SARS-CoV-2 virus, that has been widely studied in these recents weeks and we defer its description to the appropriate sources e.g. [4, 5] . The disease has been extensively modeled as a SEIR type of epidemic. In this work we use the SEIR model framework adapting it to the particular set of problems addressed here. In each section, dependent on the model variant used, we introduce the new parameters if any and refer to the appropriate section in the Appendix for the specific model equations. The rate ω or the related time 1/ω is the mean time of application of non-pharmaceutical interventions (NPI). All the recommendations that are being advised to prevent infection affect behavior (washing hands, keeping cough/sneeze etiquette, avoiding handshakes, keeping distance from other people) and, therefore, have a learning curve and can be forgotten or relaxed or incompletely adopted in any given population; also q is the proportion of the population that is effectively put under isolation. For example, in Wuhan, about 60% of the population was put under isolation. The total population normally living in that region is of about 14 million but the population at risk was obviously larger comprising much of the Hubei region.
The incubation and infectious periods have been a most important issue in the present epidemic and estimates have varied as more information is available on the dynamics of the disease. We will use incubation period 1/γ = 6 days [3, [7] [8] [9] . The same situation occurs with the infectious period that we take 1/η = 10 days [7, 10, 11] . One of the first articles to calculate the basic reproduction number R 0 of the new coronavirus in Wuhan China, estimated (Jan 24, 2020) its value as 2.68(2.47 − 2.86) with an epidemic doubling time of 6.4 days (5.8 − 7.1) [3] . Other authors have estimated the basic reproductive number to be around 2.5, value that we use here for the Wuhan original outbreak but again there exist some variability and uncertainity in the different estimates [9] [10] [11] [12] [13] [14] . In what follows some simple models are described related to three specific situations described above and occurring during this epidemic: disease dispersion due to long range population movement, behavioral change imposed on individuals as an strategy of containment, and the effect of social distancing on the general fate of the epidemic.
We analyze the cumulative trajectory of COVID19 cases in Mexico in order to create short-term real-time forecasts of the outbreak. The data consists of daily cumulative cases from February 28, 2020 to March 26, 2020 . The data was obtained daily at 19:00 p.m. (GMT-6) from the Secretariat of Health press conferences. To create the forecasts, we use Richards model which considers that cumulative cases C(t) are given by
where r represents the growth rate of infection, K is the carrying capacity (or final epidemic size, and a is a scaling parameter. This model is an extension of the simple logistic growth model that has been recently used to predict cumulative COVID19 cases in China [15] . Parameters a, r and K must be estimated in order to properly describe the observed data and to provide accurate forecasts. We use an statistical approach through Bayesian inference. Details regarding the estimation procedure and the the results can be found on Appendix A. The growth rate r is close to 0.2 and it can range from 1.9 to 2.1 with a probability of 95%. However, the total size of the outbreak K is more difficult to estimate which is expected since Mexico is at the start of the outbreak and the current data 3 . CC-BY-NC-ND 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 March 31, 2020. shows that the contagion is accelerating. We estimate the final size of the outbreak to be 710,392 cases, ranging from 65,163 up to 1,258,296 with a probability of 95% provided no cotrol/contaiment/mitigation measures are applied. These are very long-term predictions and that is why there is too much variability in the interval. To estimate the parameters we used data up to March 24, 2020 . The last two available observations corresponding to March 25 and 26 can be seen in Figure 1 , and are closely predicted by the fitted model. Furthermore, we create 14 days predictions from March 25, 2020 to April 7, 2020 which can also be seen in Figure 1 . For example, at March 31, 2020, the model expects 1,103 cases, with a 95% probability interval of (401, 2353); and at April 7, 2020, a total of 4,202 cases are expected with a 95% probability interval of (1,440, 11,012). Table 1 shows the median posterior estimates for each parameter and 95% probability intervals. Figure 2 shows the marginal posterior distributions estimates for each parameter. Finaly, using the same data we estimated the basic reproduction number for the whole country. We use the package written in R, Estimation of R0 and Real-Time Reproduction Number from Epidemics (https://CRAN.R-project.org/package=R0) using the exponential CC-BY-NC-ND 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 March 31, 2020. . https://doi.org/10.1101/2020.03.28.20046276 doi: medRxiv preprint and maximum likelihood methods. We obtained using the exponential growth method R 0 = 2.7 (2.5, 3.2); and using the maximum likelihood R 0 = 2.3 (2.0, 2.6). Given that for the Richards curve, the growth rate for Mexico is r = 0.2 below of that of the Wuhan epidemic, see Table 1 , we use the estimate from the maximum likelihood method that renders R 0 = 2.3 with an uncertainity interval of (2.0, 2.6).
In what follows we use several mathematical models to explore the consequences of interventions on the spread of the disease, i.e., how interventions modify the worrisome picture depicted in Figure 1 .
The Mexico City airport is one of the main connections hubs of Latin America and the country, and therefore arrival of COVID19 cases has occurred. In this section, we use data on the daily number of domestic and international travelers that arrive/leave to/from the Mexico City airport [16] . Below we follow the ideas presented first in [3] . We construct an standard SEIR Kermack-McKendrick mathematical model (see Appendix B) where s, e, y a , y s and r represent susceptible, exposed, asymptomatic infectious, symptomatic infectious, and recovered (immune) individuals. N is the total population. The force on infection is (β a y a + β s y s )/N . We consider that symptomatic individuals have restricted movement since either they are too sick to travel or are immediately detected by the health services at the airport or elsewhere. This consequently decreases their overall contact rate β s compared to that of an asymptomatic individual β a . Therefore β s < β a throughout this section. Following [3] we define two rates related to the inbound number of international air passengers to Mexico City (L im ), and the inbound number of national of both airplane and terrestrial passengers to Mexico City (L sm ). The basic model is shown in Appendix B
As the main objective in this section is to study the consequences of disease importation to a large Latin American city, we assume first that a proportion of the international air passengers arrive as disease carriers. We define ρ t to be the proportion of international air passengers at time t who cannot transmit the disease, e.g., if ρ t = 0.75 then 25% of the international air passengers are infectious. Of these carriers, a proportionq arrive as exposed not yet infectious individuals and 1 −q arrive as asymptomatic infectious individuals. L ms and L sm are functions of time describing travelers either arriving to Mexico City airport of leaving Mexico City airport. These rates are time dependent since we eventually take into account the increase on the number of national travelers during holiday periods, specifically Holy Week in Mexico. For the analysis we perform a normalization in order to work with proportions and not absolute numbers of individuals (see (6) , Appendix C). One of the goals of this manuscript is to study the arrival of COVID19 to Mexico City and its plausible patterns of propagation. The timeline of our study is half a year. To estimate the daily average number of national and international inbound/outbound air and road travelers, respectively we employ information from the Secretariat of Communications and Transportation (SCT) [16] and the General Directorate for Tourism Competitiveness of Mexico City [17] . The timeline is set to include Holy Week, a short holiday which normally presents an increase of arrival/departure flights to/from Mexico City. On the other hand, the main transport means of the Mexican population is by far terrestrial by bus. Unfortunately, we only have data on the number of bus passengers outbound from Mexico City, reported by the large bus terminals Norte, Sur, Poniente and Oriente. We will assume that on average the number of bus passengers that arrive to Mexico City is greater than those that leave.
For our simulations, we employ flight data for February 2020, terrestrial data for January 2018/2019 and take February 21st as the initial date of the arrival of the first cases. This estimate is plausible since the first cases in Mexico were reported on February 28th but given that the incubation period is of 6 days they likely arrived to Mexico City on February 21st or 22nd. We also set the population size of the Greater Metropolitan Area of Mexico City at 25 210 748 inhabitants. To estimate the inbound/outbound international travelers migration rates, we take into account all flights that arrive/leave Mexico City [16] . Thus, the daily average number of international arrival and departure air travelers are L im = 22, 519 and L mi = 21, 838, respectively. See Appendix G for the procedure followed to give this estimates. To estimate the inbound/outbound national travelers migration rates we integrate air and terrestrial transportation. For national flights we arbitrarily select the thirteen more important destinations from Mexico City, i.e., those who have more passengers inbound leaving from Mexico City. On the other hand, to estimate the number of terrestrial passengers, we use data from the above mentioned four main bus terminals of Mexico City. Therefore, we calculate the averages per day of national passengers in both kinds of transport (airplane and bus). Thus L sm = 141, 664 and L ms = 141, 425 except during the Holy Week (for more details see Appendix G). For the 2020 Holy Week (April 4th to April 12th), we proceed to increase previous estimates by 25% which is a rough and conservative approximation based on data from the Secretariat of Tourism [17] . Therefore L sm = 177, 080 and L ms = 176, 781 during this period and also the information available on the main bus terminals of Mexico City (for further details see Appendix G). Finally, the incubation and the infectious period are set to γ = 1/6 and η = 1/10 per day, respectively. Also, since we do not have data on the air passenger percentage that arrives in the latent or asymptomatic stage, we propose as a plausible valueq = 0.5 as a null hypothesis from which to start our analysis.
For our first case, we assume ρ t = 0.999956 that is Mexico City airport received only one imported case per day (IPD). Figure 3 shows that the higher the proportion of asymptomatic infected arriving two things happen: the outbreak occurs earlier and also the peak of the outbreak is higher. In Table 2 , we further illustrate this case. Notice that 7 . CC-BY-NC-ND 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 March 31, 2020. . the changes in prevalence at different dates when proportions of asymptomatic infected is increased. For this model we explicitly calculated the basic reproductive number since it depends on ρ the proportion of asymptomatic infectious individuals, and explore the consequences of taking ρ = 0.30, 0.60 and 0.90, or, equivalently, R 0 = 2.0585, 2.1160 and 2.1719, respectively. These values are within the range reported by several authors and our own uncertainity interval (Section 1) [1, 10, 18] . 
In our second case, we explore the effect of we explore the effect of decreasing ρ t which is equivalent to increasing the number of IPD. Figure 4 shows that for a higher number of IPD generates an earlier outbreak but the maximum of the outbreak does not change. In Table 3 , we further illustrate this case for Mexico City. Once again we observe that prevalence changes when increasing the number of IPD. This case is set by assuming that 80% (ρ = 0.8) of the infected people are asymptomatic giving R 0 = 2.1534.
We remark that if we simulate both cases in (Figures 3-4) for a longer period of time, we can detect the existence of another epidemic outbreak as has already been noted [19] .
. CC-BY-NC-ND 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 March 31, 2020. . Figure 5 shows the expected number of exported cumulative cases (until April 1st) to several of the Mexican cities that normally receive the highest outbound air or road travel volumes. Here, as ρ = 0.8, R 0 = 2.1534. Figure 5 (a) shows the expected number of exported cumulative cases that use air travel. When considering one IPD, we see that few states receive at least one exported case from Mexico City (Cancun, Guadalajara and Monterrey). On the other hand, with three IPD, these states import 4, 3 and 3 cases from Mexico City, respectively. Figure 5 (b) shows the exported cumulative cases for terrestrial transportation. Here, if the imported cases per day (IPD) is one, then Puebla, Jalisco, Michoacan, Guerrero, and Queretaro import 4, 5, 3, 2 and 1 cases up to Aril 21st from Mexico City, respectively; while if IPD is equal to three, then Puebla, Jalisco, Michoacan, Guerrero, and Queretaro import 13, 16, 9, 7 and 4 infections from Mexico City, respectively. Figure 6 shows the expected number of exported cumulative cases to several Mexican cities. We simulate the number of cases until April 21st. We can observe the importance of the number of people who arrive infected since these can significantly increase the spread of the disease to other states. For example, Figure 6 (a) shows that for Cancun, the number of exported cases can increase by 50 cases if instead of 1 IPD, we have 3 IPD. For Monterrey the increase, under the same conditions, is approximately 40 cases. Figure 6 (b) shows a greater increase in cases for terrestrial than for air transportation. For example, we expect 9 . CC-BY-NC-ND 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 March 31, 2020. Finally, it is an obvious conclusion from our simulations that it is very important to control the number of people who arrive infected since these can significantly increase the spread of the disease to the other states. Note that the results in this section do not include any control, mitigation of containment strategy. The results shown in Tables 2-3 are some   10 . CC-BY-NC-ND 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 March 31, 2020. . Table 3 :
Mexico City prevalence (number of cases) for different days as a function of IPD (imported cases per day) under the assumption that 80% of cases end up as asymptomatic.
Once the disease arrives into a given region, city, or municipality, containment measures are at first, limited to behavioral changes called non-pharmaceutical interventions or NPIs, directed essentially to diminish transmission by reducing contact between people. In order to model this, we modify system (15) to include a variable contact rateβ k (t), and omit migratory phenomena (see Appendix D for the model equations).
Our main results in this section explore the variation of the effective contact rate under these actions. We definê
where k = a or s for asymptomatic or symptomatic individuals, respectively. Note from the formula above that we assume that behavioral change geared to diminish the contact rate from β k to q 1 β k takes a learning time θ to acquire where q 1 is the desired proportion of reduction of the initial contact rate, e.g., q 1 = 0.25 if the contact rate is reduced in 75% at the end of a time interval of length θ. Moreover, these control measures are not taken immediately after the outbreak is identified but lag a certain time T θ before being enforced or (see Figure 7 ) .
The baseline parameters for the simulations that follow are β a = 0.48, β s = 0.4, γ = 1/6, η = 1/10, ρ = 0.8,
. CC-BY-NC-ND 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 March 31, 2020. .
where, we remind the reader correspond to the contact rates of asymptomatic and symptomatic infecteds, the incubation and infectious period and the proportion of asymptomatic cases respectively. Figure 7 illustrates the behavior of the contact rates for asymptomaticsβ a (t) and for symptomaticsβ s (t). For both rates, the efficacy of education on the behavioral change is the same since we assume that education is not related to symptomatology although there is some evidence that this might not be the case. A first result is that if the time it takes to obtain the desired contact rate reduction θ increases then the reduction rate of both contact rates is slower. 9 show that if the time it takes to achieve a reduction in contact rate increases then the number of infectious cases is greater reflecting the slow rate of behavioral change. Figure 8 a) and b) show the prevalence and cumulative incidence respectively, for an scenario with a target reduction of contact rate of 25%, the initiation of the behavioral change is T θ = 15 days after 10 cases are detected and where 80% of those are asymptomatic. For comparison, a red dashed line shows the epidemic curves when no behavioral change is enforced. Here we see that although the reduction of contact rate is low (25%) (q 1 = 0.75), in the long run behavioral change helps to reduce the incidence. We also note that in addition to the reduction in cases, as the time to achieve the contact rate reduction increases, the epidemic outbreak starts earlier and is larger. Figure 9 shows the behavior of the prevalence (a) and cumulative incidence (b) when the contact rate is reduced by 50% (q 1 = 0.5). In this situation the relevance of the behavioral change to reduce the incidence is rather clear. Notice that if the time necessary 12 . CC-BY-NC-ND 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. to achieve the contact rate reduction is either θ = 30 or θ = 90 days (differing for more than two months), the outbreaks are very similar which makes us postulate the existence of a threshold value θ * such that the peak incidence is reduced if θ < θ * . Larger learning times do not have an impact on the size of the epidemic. Figures 10-12 show examples of the behavior of the cumulative incidence (at 180 days after 10 infected are detected which we use as the indicator of initiation of the outbreak in this example) as a function of certain parameter values. Figure 10 shows that if T θ is closer to the time of detection of the first cases (in this example, when 10 infected people are detected), then the cumulative incidence can be reduced more efficiently and underlines the importance of reducing the effective contact rate quickly. In our example even when the reduction is only of 50%, we see that the final cumulative incidence decreases to low levels. This observation contrasts with what is shown in Figures 8 and 9 . We will discuss this in more detail later. On the other hand, Figure 10 (b) shows that even when the time it takes to get the contact rate reduction is large (θ = 90 days), we can still recover a large reduction in the final cumulative incidence provided we have a quick response to enforce reduced contacts and a large target reduction in the effective contact rate. Figure 11 shows that when T θ = 15 days, regardless of the time θ it takes to obtain the contact rate reduction, the final cumulative incidence is very low for large reductions in contact rates (q 1 small) or high for low reductions in contact rate (q 1 large).
To finish this section, we remark that Figure 10 Figure 10 : (a) Cumulative incidence for varying θ, the learning time to achieve the target reduction and T θ , the delay to start behavioral change for a fixed reduction of contact rate of 50%, and (b) cumulative incidence for varying T θ , the delay to start behavioral change and q 1 the target proportion of reduction of the contact rate for a fixed time to obtain the target contact rate of 90 days.
reducing the contact rate for flattening the epidemic curve. On any scenario, a large reduction of the effective contact rate is necessary (Figure 12(a) ). Figures 12(b) shows an 14 . CC-BY-NC-ND 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. 
If the epidemic cannot be contained, social isolation has been implemented in many places of the world, notably in China where, given the extent and strictness of its enforcement, it has apparently been successful. We modify our model to incorporate this containment measure and let s, e, y a , y s and r as usual represent susceptible, exposed, asymptomatic and symptomatic infectious and recovered (immune) individuals in isolation which occurs when a proportion q of the general population is put on restricted movement. The population 15 . CC-BY-NC-ND 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 March 31, 2020. . https://doi.org/10.1101/2020.03.28.20046276 doi: medRxiv preprint not under isolation (i.e., a fraction 1−q of the total) is denoted by the letters c s , c e , c ya , c ys and c r for susceptible, exposed, asymptomatic and symptomatic infectious and recovered (immune) individuals. In this model isolation comprises not only physical isolation, but also the enforcement of sanitary and preventive measures like frequently washing hands, use of face masks, avoiding of public gatherings and so forth.
The actual physical isolation depends on many factors, particularly the length of the incubation period of the virus but many of the preventive measures taken by individuals and groups of individuals can last much longer. The isolation enforcement has the following modeling set up. After the first appearance of the virus a SEIR epidemic starts and develops until time T θ where quarantine and isolation are enforced and a fraction q of individuals in the population goes into isolation. We define q to be the efficacy of isolation, e.g., q = 0.9 if 90% of the population is isolated.
We then record the number of people in each of the SEIR compartments at T θ , and use these numbers as the initial conditions for the state variables of the same SEIR epidemic but under isolation starting at T θ . The subsystem that escaped or that was not put under isolation we denote by cSEIR.
We also assume that isolation measures are not perfect and individuals escape from isolation regardless of infection status at an exponential average time of 1/ω and also that individuals suspend, forget or relax sanitary measures at the rate ω. It is also assumed that the isolation time is always larger than the total time of duration of disease comprising the incubation time, 1/γ, and the infectiousness time 1/η. Therefore 1/ω > 1/γ + 1/η. We further assume that a fraction ρ of the exposed individuals develop the asymptomatic disease and that 1 − ρ become symptomatic. The contact rates before isolation is enforced are β i a and β i a for symptomatic and asymptomatic individuals. Following the ideas set in the previous section we define two time-dependent contact rates for the isolated population that can be consulted in Appendix D.
The natural mortality and recruitment rates respectively are assumed to be negligible during the time of the epidemic. With these definitions the model can is shown in Appendix E. In summary, for t < T θ the infectious disease is the standard SEIR Kermack-MacKendrick model. Once t = T θ is reached, the population is divided into the isolated and non-isolated sections with contact rates (9) and (13) respectively.
To calibrate our model, in this section we characterize the epidemic in Wuhan, China (epicenter of the first epidemic outbreak) [6] . We introduce an isolation stage and the consequent reduction in the contact rates.
For Wuhan the basic reproduction number has been estimated to be R 0 i = 2.51 [1, 3, 
. CC-BY-NC-ND 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 March 31, 2020. . https://doi.org/10.1101/2020.03.28.20046276 doi: medRxiv preprint η = 1/10, ρ = 0.80, ω = 0.01, q 1 = 0.10, q = 0.99. Figure 13 shows the contact rates for the isolated environment. For both rates, the reduction of contact is 90% (q 1 = 0.1), the day of initial isolation is January 26th. We can impose a fast behavioral change because we think that this process lasted a short time (7 days) in Wuhan. Figure 14 presents the case when the 99% of population is isolated. Observe that the peak incidence is approximately 4 935 individuals and it occurs 2 days after isolation began (January 26th, [20] ) Figure 14(b) ). Our simulations show that carrying out an isolation process promotes an earlier decrease in incidence. On the other hand, considering a total population equal to 14 million [21] the prevalence at the peak of the epidemic in Wuhan is around 19 000 cases which is a gross underestimate of the true observed value [6] . A better estimate of the at-risk population in Wuhan would, we believe, improve or estimate. Figure 15 shows a comparison of the estimated cumulative incidence (blue line) versus the data reported by the WHO (red cross) for the Hubei province [20] . Our curves overestimate the data, which may be reasonable due the known under-reporting of cases and to the fact that initially there are fewer cases making it difficult to detect them. Nonetheless we see that when the real epidemic curve begins to stabilize, the difference between both data and model decreases, which can be explained if we accept that, at this stage, there are many more cases and it is easier to identify them. Also under-reporting decreases for the same reason. The black point on the model cumulative incidence (blue curve) indicates the inflection point. The black dots over the red curve mark the range where the inflection point is located in the data.
. CC-BY-NC-ND 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 March 31, 2020. 
The mathematical model in this section is a different variant of the one we applied for Wuhan because Mexico City outbreak started with imported cases. Therefore we retake the model proposed in Section 1 and combine it with the Wuhan model of the previous section to incorporate isolation and behavioral change.
In Mexico air travel has been gradually suspended and at the time of writing (end of March) not all flights have been canceled. This implies that the flow of people from abroad is low but significative and it is arriving into a population that is not isolated but has started measures to reduce the contact rate (social distancing). The model is shown in 18 . CC-BY-NC-ND 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 March 31, 2020. Below we explore three scenarios that might occur in a large metropolitan that includes variation in the reduction of contact rates, variation in the time it takes to get the contact rate reduction θ and different dates for start de isolation T θ . As in Section 1, we consider February 21 as the initial date of infection, Figures 16-18 show the behavior of the outbreak after and before isolation.
On March 23rd authorities of Mexico City ordered the closure of cinemas, theaters, and other public places as a preventive measure. We show in Figure 16 the incidence and cummulative incidence under reductions in contact rates of 25%, 50% and 75%, and for θ = 7 days. Our simulations show that for these reductions in contact rates, a reduction of up to 91% in incidence can be obtained. However, if we take the first 35 weeks after the arrival of the first imported case, Figure 16 (b) shows that depending on the contact rate reduction, we could reach a similar situation to that of Wuhan with 69 050 total cases, (q 1 = 0.75), or a medium level situation with 10 800 total cases (q 1 = 0.50) or a more manageable situation with 4 351 total cases, (q 1 = 0.25). Figure 17 shows the scenario for different θ, the time it takes to get the desired contact rate reduction. In this case isolation also begins on March 23rd and the reduction in contact rates is 60%. We observe that by increasing the value of this parameter, the peak and cumulative incidence will be higher. For example, if the time it takes to reduce the contact rates is θ = 7 days, then the cumulative incidence 16 weeks after the arrival of the first imported case will be 6 532, while if the time is θ = 15 days, then the incidence rises to 8 909. Finally, if the time it takes to reduce contact rates is 30 days, then the cumulative incidence rises to 14 980. It is important to note that although the value of the cumulative 19 . CC-BY-NC-ND 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 March 31, 2020. Figure 16 : (a) Incidence and (b) cumulative incidence for the isolated population. The time to obtain the target contact rate reduction is 7 days. Black, red and blue lines correspond to reductions in contact rate of 25%, 50% and 75%, respectively. Horizontal axis is time in days and the vertical axis is number of individuals.
incidence increases, it does not happen with the same intensity as in the case shown in Figure 16 . 
. CC-BY-NC-ND 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 March 31, 2020. . https://doi.org/10.1101/2020.03.28.20046276 doi: medRxiv preprint
In our third case we explore variation in the time to start the isolation. Here the time to reduce the contact rates is θ = 7 days and the target reduction is 60%. The initial dates of isolation are: March 23rd (blue line), March 30th (red line) and April 6th (black line). The objective of this scenario is to asses what would happen if large cities start isolation and enforce behavioral change several days after this was done in the Mexican capital. Figure 18(a) shows the increase in incidence as the isolation start date increases. Solid lines represent the situation before isolation and the dashed lines represent the situation after the start of isolation. Figure 18(b) shows that if the onset of isolation had occurred one or two weeks after March 23rd, then by week 16, the cumulative incidence would increase by 12 058 and 45 118 new cases, respectively. This scenario underlines the importance of taking contingency measures with due anticipation, otherwise the number of infected people will notoriously increase. Figure 18 : (a) Incidence and (b) cumulative incidence for the isolated environment. Here we consider that the reduction of contact rate reduction is 60%, and the time it takes to get the contact rate reduction is equal to 7. Blue, red and black lines mean that the dates for start de isolation are 31, 38 and 45 days after the first arrival of an infected individual to Mexico City, respectively. Horizontal axis is time in days and the vertical axis is number of individuals.
The control, containment, mitigation and final resolution of the coronavirus epidemic requires careful consideration of alternative scenarios as pointed out very clearly in [19] . Models like the ones presented here can help decision-makers to have amore comprehensive understanding of the epidemic and thus a better set of criteria on which to base decisions.
. CC-BY-NC-ND 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 March 31, 2020. .
The models presented here are of the same general type that have been published and are being written elsewhere for the coronavirus epidemics. The coronavirus epidemic is an acute respiratory infection highly contagious and with a case fatality rate above of that for influenza. From the modeling point of view the Kermak-McKendrick formalism is the basic and general mathematical tool to understand the epidemic at a population level. The conclusions that we can reach from the study that we have developed in the previous sections are the following:
• Restricting population mobility has been reported as highly effective for the containment of the initial COVID-19 outbreak in China [22, 23] . We have attempted to describe the present conditions in Mexico to evaluate the consequences of delaying or advancing the enforcement of isolation and mobility restrictions putting special emphasis on the role of asymptomatic cases since these have been identified as probable important drivers of the disease [13, 24] . At the start of the epidemic in Mexico, Figure 3 shows that the higher the proportion of asymptomatic cases arriving two things happen: the outbreak occurs earlier and also the peak of the outbreak is higher. The role of asymptomatic cases has been under study in the literature particularly as their likely important role in driving the epidemic appears clearer. In our case, their impact is to generate an earlier outbreak and also a higher maximum incidence.
• Figure 4 underlines the relative insensitivity of the peak incidence on the rate of importation of cases provided this rate is relatively low (up to 5 cases per day). Nonetheless, the higher this rate the earlier the outbreak goes into exponential growth (phase 2 according to WHO criteria). We remark that running the simulations that (Figures 3-4) ilustrate for a longer period of time, we can detect the existence of a second epidemic outbreak. With respect to disease spread at a national level in Mexico, the importance of controlling the number of people who arrive infected to the larger airport hubs and bus terminals is of first priority since, not doing it, these can significantly increase the spread of the disease to the other states.
• Behavioral change is one of the main forms to reduce contact rates [19] . However, when dealing with behavioral issues we enter into the psychological realm of individuals and groups of individuals. For our simulations we assume that education or awareness and rational decisions regarding the behavioral change associated with reducing contact rates is the same for asymptomatic and symptomatic individuals although there is some evidence that this might not be the case as one can directly witness. We assume different proportions of these two types of infected individuals. This is an open area of exploration that has not been pursued here. In any case, if asymptomatic individuals or any identifiable portion of the population have more "reckless behavior" on average than other individuals infected or not, then the expected incidence will be higher in any given scenario.
. CC-BY-NC-ND 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 March 31, 2020. .
• The policies of containment and mitigation have to be applied quickly, with no hesitation. We observe from our simulations that if the time it takes to obtain the desired contact rate reduction θ increases then the reduction rate of both contact rates, for asymptomatic and symptomatic individuals, flattens the curve but more slowly. We also note that in addition to the reduction in cases, as the time to achieve the contact rate reduction increases, the epidemic outbreak starts earlier and is larger.
In the examples illustrated above, we see that if the time necessary to achieve the contact rate reduction is either θ = 30 or θ = 90 days (differing for more than two months), the outbreak sizes are similar which to us indicate the existence of a threshold value θ * such that the cumulative incidence is reduced if θ < θ * . In summary large learning times do not have an impact on the size of the epidemic. Figure 10 shows that if T θ is closer to the time of detection of the first cases (in this example, 10 infected people), then the cumulative incidence can be reduced more efficiently and underlines the importance of reducing the effective contact rate quickly. Figure 10 (b) shows that although the time it takes to get the contact rate reduction is high (θ = 90 days), we still can obtain a large reduction in the final cumulative incidence provided we have a quick response to enforce reduced contacts and a large target reduction in the effective contact rate. Figure 11 shows that when T θ = 15 days, regardless of the time θ it takes to obtain the contact rate reduction, the final cumulative incidence is very low for large reductions in contact rates (q 1 small) or high for low reductions in contact rate (q 1 large). Figure 10 (a) and 12 show the importance of reducing the contact rate for flattening the epidemic curve. On any scenario, a large reduction of the effective contact rate is necessary (Figure 12 (a)). Figures 12(b) shows an example where an insufficient reduction of the contact rate (40%), leads to a very weak reduction on the final cumulative incidence.
• Figure 14 shows the case when 99% of population is isolated. We can see that the peak incidence occurs 2 days after isolation began (January 26th) indicating that the isolation process promotes an earlier decrease in incidence. This scenario underlines the importance of taking contingency measures with due anticipation, otherwise the number of cases will notoriously increase.
• Figure 18 highlights the importance of beginning isolation early T θ = 31 (March 23rd). Our model indicates that delaying action by one or two weeks will let the outbreak grow with a significant increase in cases. Figure 16 shows, however, that isolation is not enough to lower the size of the outbreak but that a significant reduction in transmission rate is also needed. Both control measures have to be applied simultaneously to be efficient in flattening the epidemic curve. We remark that simulations in both, are for a transmission rate reduction achieved in θ = 7 days. Figure 17 shows that for longer times, the total number of cases will increase.
All of the latter leaves the message that apart from designing and enforcing control 23 . CC-BY-NC-ND 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 March 31, 2020. . measures, the full participation of the society at large in the fight against this disease is extremely important. It is necessary to reduce the contact rate to a low enough level to be able to contain it, it is necessary to isolate, trace and test individuals. These actions require from everybody to be deeply aware of what should and should not be done to stop the spreading of the disease.
Let Y j , for j = 1, 2, ...n, be the number of observed cumulative cases at time t j , with t j given in days. We assume that Y j follows a Negative Binomial distribution with mean value C(t j |a, r, K) and dispersion parameter α. Here, C(t j |a, r, K) 1 is the solution of Richards model presented in (1) . Assuming that, given the parameters, the observations Y 1 , Y 2 , . . . , Y n are conditionally independent, then
The Negative Binomial distribution allows to control the variability of the data by considering over-dispersion which is common for epidemiological data. If α = 0, then we return to the Poisson model which is often used in this context. Let θ = (a, r, K, α) be the vector of parameters to estimate. The inclusion of the parameter α, which is related to the variability of the data, not to the Richards model, is necessary since in practice this variability is unknown. Then, the likelihood function, which represent how likely is to observe the data under the Negative Binomial assumption and Richards model if we knew the parameters, is given by π(y 1 , . . . , y n |θ) = n j=1 Γ(y j + τ ) Γ(y j )Γ(τ ) τ τ + C(t j |a, r, K) τ C(t j |a, r, K) τ + C(t j |a, r, K)
Consider that parameters a, r, K and α as random variables. Assuming prior independence, the joint prior distribution for vector θ is π(θ) = π(a)π(r)π(K)π(α),
where π(a) is the probability density function (pdf) of a Uniform(0,1) distribution, π(r) is the pdf of a Gamma(shape=2, scale=1), π(K) is the pdf of a Uniform(K min , K max ), and π(α) is the pdf of a Gamma(shape=2, scale=0.1). To select the prior for parameter r, we consider that previous estimation of r are close to 0.3 [15] , and a Gamma(2,1) represent a weekly informative prior as it allows for a wide range of values of r. Also, there is no available prior information regarding the final size of the outbreak K. This is a critical parameter in the model and, in order to avoid bias, we assume a uniform prior over K min and K max . To set these last to values, we consider that the minimum number of confirmed cases is the current number of observed cases Y (t n ) times nine, i.e. K min = y n * 9. The reasoning behind this number is that it has been suggested that the observed cases are under reported by a factor of nine (REF). To set the upper bound for K, we consider a fraction of the total population K max = N * 0.01, where N is the population size of Mexico. This fraction was determined base on the observations of other countries such as Italy where the proportion of infected represents one of the worst case scenarios considering its population size. Then, the posterior distribution of the parameters of interest is π(θ|y 1 , . . . , y n ) ∝ π(y 1 , . . . , y n |θ)π(θ), and it does not have an analytical form since the likelihood function depends on the solution of the Richards model, which must be approximated numerically. We analyze the posterior distribution using an MCMC algorithm that does not require tuning called t-walk s [25] . This algorithm generates samples from the posterior distribution that can be used to estimate marginal posterior densities, mean, variance, quantiles, etc. We refer the reader to [26] for more details on MCMC methods and to [27] for an introduction to Bayesian inference with differential equations. We run the t-walk for 100,000 iterations, discard the first 25000 and use 1,000 samples to generate estimates of the parameters and short-term predictions for the next 14 days.
Equations for the model reviewed in Section 2. See main text for variables and parameter definitions.
25 . CC-BY-NC-ND 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 March 31, 2020. . https://doi.org/10.1101/2020.03.28.20046276 doi: medRxiv preprint C National spread model Equations for the model reviewed in Section 3 that includes transportation. For the analysis we standarize the model variables to work with proportions. We define the variables:
Thus, system (6) becomes
This model sets the migration functios L im = Lmi = L sm = L ms = 0 and µ = 0. Thus
with varying contact rates given as:
with k = a or s for asymptomatic and symptomatic individuals, respectively.
. CC-BY-NC-ND 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 March 31, 2020. 
After isolation is enforced, we have a model for the isolated population
and for the non-isolated population: c s = − β s (t)c ys +β a (t)c ya c s + ωx s c e = β s (t)c ys +β a (t)c ya c s − γc e + ωx e c ya = ργc e − ηc ya + ωx ya c ys = (1 − ρ)γc e − ηc ys + ωx ys c r = η(c ya + c ys ) + ωx r (12) with contact transmission for the non-isolated part of the population:
where k = a or s for asymptomatic or symptomatic individuals, respectively.q 1 is the desired proportion of reduction of the initial contact rate, e.g.,q 1 = 0.35 if the contact rate is reduced in 65% at the end of a time interval of lengthθ.
. CC-BY-NC-ND 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 March 31, 2020. 
Dynamics before isolation is given by system 7. After isolation, dynamics for the isolated population is given by
And dynamics for non-isolation is described by: 
28 . CC-BY-NC-ND 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 March 31, 2020. .
To estimate the values of L im and L mi , we calculate the average number of passengers per day arriving/departing from Mexico City from/to abroad taking the total number of passengers for each mean of transportation and divide it by the corresponding number of days. We use the data reported by the Ministry of Communications and Transportation [16] . From From Table 4 , the values of total air passengers that arrived and departure from Mexico City in February 2020 are 653050 and 633312, respectively. Therefore, Lim ≈ 22519 and Lmi ≈ 21838.
The values of L sm and L ms outside the Holy Week are estimated similarly to the previous case from data from the Ministry of Communications and Transportation (MCT) and the General Directorate for Tourism of Mexico City [16, 17] .
From Table 5 , the average number of passengers per day, outside the Holy Week air and terrestrial departing from Mexico City are 31, 173 and 110, 253, respectively. Hence L ms = 141425. Proceeding in similar manner, the value of L sm is equal to 141664. Finally, to calculate the L ms and L sm values during Holy Week, we increased the above values by 25%. We chose this value because when averaging visits to the Official Sites of the Ministry 29 . CC-BY-NC-ND 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 March 31, 2020. 
Conceptualization: MAAZ, JXVH. Formal Analysis and Investigation: MAAZ, JXVH. Performed statistical analysis: MSC; Wrote and discussed the paper: All authors.
. CC-BY-NC-ND 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 March 31, 2020. . https://doi.org/10.1101/2020.03.28.20046276 doi: medRxiv preprint