First detected in December 2019 in the city of Wuhan in Hubei Province, China, the COVID-19 outbreak, caused by the newly identified coronavirus SARS-CoV-2, has spread around the globe and reached the status of a pandemic on March 11th, 2020. Due to the severity of the damages it may cause to health and its ease of transmission, a number of different strategies have been implemented by the authorities of different countries to block or reduce the spread of the virus. In some cases like China, Italy, or Spain, strict quarantine measures have been adopted [1–3]. However, strict lockdown in many cases has been impossible due to prevailing economic and social factors. In such cases, the authorities aimed for less strict mitigation policies [4], including social distancing and individual non-pharmaceutical interventions [5]. Nonetheless, an accurate description of the progression of an epidemic is of fundamental importance in helping to decide public policies to reduce its impact, especially to stay below a fixed healthcare capacity and delaying the peak of the epidemic so that the healthcare capacity can be expanded to support patients.

In Mexico, the first detected case of COVID-19 was registered on February 27th, 2020 and corresponded to an imported case from Italy. This marks the start of what was called phase 1 of the epidemic, characterized by imported cases only, with no local contagion. Mexican health authorities identified a total of three epidemiological phases, according to the degree of transmission of the disease. Phase 2 of the coronavirus pandemic was characterized by cases of the local contagion between people who have not had contact with foreigners; it was declared on March 24th, and the actions comprised primarily the suspension of certain economic activities, the suspension of lessons at schools, the restriction of mass congregations, and the recommendation of domiciliary protection for the general population. As a consequence of the evolution of confirmed cases and deaths from the disease in the country, on March 30th, a “health emergency due to force majeure” was declared, which led to the execution of additional actions for its prevention and control; the most conspicuous was the generalized voluntary quarantine of the population (the so-called Jornada Nacional de Sana Distancia). Eventually, on April 21st, phase 3, characterized by thousands of cases disseminated in all the country, was declared.

Interventional measures adopted with the intention to mitigate the spread are normally based on estimates of the progression of the outbreak. Mathematical models of infectious diseases are important tools for assessing the threat of a novel pathogen and offer the best information for mitigating an outbreak [6, 7], hence the need for epidemiological models that are able to estimate with some degree of accuracy the evolution of the outbreak to help to evaluate the impact of interventions [2, 8–11]. The paradigmatic approach traditionally used to model the dynamics of an epidemic is the well-known SIR (susceptible, infectious, and removed) compartmental model [12, 13]. In this model, the group S represents individuals who are susceptible to the disease and can become infected, the group I represents individuals who are infectious and can infect susceptible individuals, and the group R represents removed individuals who either gained life-long immunity at recovery or died; in either case, removed individuals cannot infect or be infected anymore. Although this model was successfully applied to describe the spread of an infectious agent in a well-mixed population [14], this same simplifying assumption prevents its successful application in many other cases [15], such as the recent outbreak of COVID-19 which shows an early subexponential growth [16, 17]. In a well-mixed population, a homogeneous distribution of the susceptible-infectious contacts such that any susceptible individual may be infected by any infectious individual in the whole population is assumed. However, pathogens affect populations in an uneven way [18]; there are many heterogeneities in human populations that influence virus transmission [19], for example, variability in the risk experienced by age [20], comorbidities, or other factors (e.g., behavior and nutrition); the presence of individuals that propagates the virus more efficiently (super-spreader individuals) [21]; and the limited transmission between geographically distant populations. Thus, any realistic epidemic model should take them into account to some extent.

Given that the subexponential growth seems to be a generic characteristic of the COVID-19 outbreak, independent of the suppression strategies implemented to mitigate the temporal evolution of the epidemic process, it is suggested that the existence of an underlying mechanism is responsible for this temporal behavior. The purpose of the present study is to show that the inhomogeneous transmission of the epidemics of component subpopulations may be the source of this behavior. It is also shown that the standard SIR model can be extended to include the abovementioned inhomogeneities and that the resulting model correctly captures not only the short time but also the longtime dynamics of the COVID-19 outbreaks.

A natural way to incorporate population heterogeneities or spatial structure into an infectious disease model is by means of metapopulation models [6, 22–24]. Let us assume a cross-coupled metapopulation approach [24] in which the total population is considered as if it were formed by n subpopulations or patches connected to each other with transmission lines, and no explicit mobility among subpopulations is included. In each subpopulation, an infectious agent spread is described by a standard SIR model with coupling terms. d s i d t = − s i ( t ) 1 n ∑ j = 1 n β i j i j ( t ) , (1) d i i d t = s i ( t ) 1 n ∑ j = 1 n β i j i j ( t ) − γ i i ( t ) , (2) d r i d t = γ i i ( t ) , (3)

where s i ( t ) = S i ( t ) / N i , i i ( t ) = I i ( t ) / N i , and r i ( t ) = R i ( t ) / N i are the fractional representations of susceptible, infectious, and removed individuals of subpopulation i that has N i individuals and satisfies the relationship s i ( t ) + i i ( t ) + r i ( t ) = 1 . No human birth and death rates are being considered, and as stated above, neither the immigration nor emigration effects [2]. Although it is known that COVID-19 has a mean incubation period of about five days [25], in order to maintain the number of parameters and equations to a minimum, we are going to ignore it, and therefore, there is no compartment of exposed individuals. In Eqs 1–3, β i j represents the elements of a matrix that describes the transmission between and within patches, the recovery rate γ is common to all subpopulations, and, as in traditional SIR, all these quantities are time independent. The transmission rate β i j captures the rate of flow from group S i to group I j , while the recovering rate γ indicates that infectious individuals get recovered or die at a fixed average rate γ. If subpopulations are independent from each other, then β i j = β 0 δ i j , with δ i j as the Kronecker delta function. The basic reproduction number, R 0 = β 0 / γ , captures the average number of secondary infections an infected individual will cause before he or she recovers or is effectively removed from the population; in other words, it measures how many other people an infected individual will infect. However, as the number of susceptible individuals in the population declines due to a growing number of infections, the effective reproduction number over time, R i , t , is given by the product of R 0 , and the fraction of susceptible individuals in the subpopulation i is given by R i , t = R 0 s i ( t ) ; in other words, the effective reproduction number is the expected number of new infections caused by an infectious individual in a population where some individuals may no longer be susceptible. During the first few generations of disease transmission, in the absence of control interventions or reactive population behavior changes, the SIR model supports a reproduction number that is essentially invariant, for example, R i , t ≃ R 0 . Thus, in the classical SIR model with a constant transmission rate, β 0 , in a completely susceptible population, s i ( 0 ) ≃ 1 , i i ( t ) grows exponentially during the early epidemic phase [26], that is, i i ( t ) = i 0 e ( β 0 − γ ) t .

As mentioned before, it is unusual for a naturally occurring disease emergence to occur simultaneously at many locations and to propagate at a uniform rate. This means that at least during the initial phase of transmission, infectious individuals are clustered [8]. The presence of the subpopulations can be thought as clustered regions in space where a given individual spent most of his time. Clusters exchange pathogens with each other through infected or susceptible individuals traveling among them during the period of infectiousness. Thus, spatial inhomogeneities lead naturally to outbreaks that do not occur simultaneously in all subpopulations. Even if they are governed by the same dynamical equations, they are asynchronic, which means that the onset of the outbreaks in the subpopulations are not necessarily simultaneous [27–29]. In a similar way, different individuals may propagate the pathogen at different rates, and one can also consider groups of individuals with similar transmission rates as distinct subpopulations even if not geographically remote. Adding the corresponding differential equations of all subpopulations, one gets the following equations: d s d t = − 1 n ∑ i = 1 n s i ( t ) 1 n ∑ j = 1 n β i j i j ( t ) = − β 0 s e f f ( t ) i ( t ) = − β ( t ) s ( t ) i ( t ) , (4) d i d t = 1 n ∑ i = 1 n s i ( t ) 1 n ∑ j = 1 n β i j i j ( t ) − γ n ∑ i = 1 n i i ( t ) = β 0 s e f f ( t ) i ( t ) − γ i ( t ) = β ( t ) s ( t ) i ( t ) − γ i ( t ) , (5) d r d t = γ n ∑ i = 1 n i i ( t ) = γ i ( t ) , (6) where s ( t ) ≡ S ( t ) / N = ∑ s i ( t ) / n , i ( t ) ≡ I ( t ) / N = ∑ i i ( t ) / n , and r ( t ) ≡ R ( t ) / N = ∑ r i ( t ) / n are the fractional representations of susceptible, infectious, and removed individuals of the total population with N individuals, and N = n N i relates to the total population with that of each subpopulation N i which is assumed to be the same for all subpopulations. To obtain the final relations in Eqs (4–6), we have defined a time-dependent transmission rate by the following equation: β ( t ) = β 0 s e f f ( t ) s ( t ) , (7) with s e f f ( t ) ≡ ∑ i = 1 n s i ( t ) [ 1 n ∑ j = 1 n β i j i j ( t ) n β 0 i ( t ) ] . (8)

If all the subpopulations were synchronic and had the same transmission rate, then β i j = β 0 , and we can recover the standard SIR result β ( t ) = β 0 and s e f f ( t ) = s ( t ) . However, in general, β ( t ) < β 0 for inhomogeneous systems with a distribution of onset times and transmission rates. This means that s e f f ( t ) < s ( t ) . In other words, the fact that the population is not homogeneous (in contrast to the well-mixed assumption of standard SIR) implies that d s / d t is not proportional to the product of the total number of susceptible individuals in the population times the total number of infectious individuals as in the traditional SIR model. On the contrary, these inhomogeneities (the population is not well-mixed) make d s / d t proportional (with the same proportionality constant β 0 ) to the product of the number of infectious individuals multiplied by a smaller number of susceptible individuals, s e f f ( t ) . Such effective number would represent a subset of the susceptible individuals who are able to mix with the infectious population. The exact expression for s i ( t ) appearing in Eq. (8) in principle could be obtained from the solution of Eqs (1, 2); however, the quantities β i j are unknown. Here, we propose to approximate s i ( t ) in Eq. (8) by the solution of a standard SIR model, s S I R ( t + t i ) , evaluated at time t + t i that depends on the subpopulation i. The onset times t = − t i of the local epidemics can be different (they are asynchronic) from one subpopulation to the other, and transmission rates β i j can also be different. The initial conditions s S I R ( 0 ) = 1 − 1 / N i , i S I R ( 0 ) = 1 / N i , and r S I R ( 0 ) = 0 are assumed. Then, s i ( t ) ≃ s S I R ( t + t i ) can be directly substituted into Eq. (8) to obtain s e f f ( t ) ≃ ∑ i = 1 n s S I R ( t + t i ) [ 1 n ∑ j = 1 n β i j i j ( t ) n β 0 i ( t ) ] . (9)

Summarizing, inhomogeneous transmission has been incorporated in an SIR framework through a time-dependent transmission rate, β ( t ) . As we will show below, a time-dependent transmission rate may be used to explain the early subexponential growth of the spread, even in the absence of susceptible depletion or interventional measures. Let us stress that other SIR models also consider time-dependent transmission rates, but in those cases, it is introduced externally and mainly to model the reactive behavior of the population in response to containment measures [14, 18, 30, 31]. In contrast, in the present model, β ( t ) is obtained as part of the solution of the dynamic equations, and its time-dependence appears even if the population contact dynamics is uniform and the depletion of susceptible individuals is negligible.

The calculation of s e f f ( t ) in Eq. 9 would depend on the distribution of subpopulation onset times, − t i , and transmission rates, β i j , and is influenced by the amount of contacts among the individuals of different populations. The effective fractional susceptibility s e f f ( t ) can be thought as a weighted average over subpopulations, with the quantity in square brackets in Eq. 9 playing the role of a weight factor. Assuming that we can transform the average over subpopulations into a time average and making a continuous approximation, we propose to transform Eq. 9 into the following equation: s e f f ( t ) ≃ ∫ − ∞ ∞ s S I R ( t + t ′ ) ρ ( t , t ′ ) d t ′ ∫ − ∞ ∞ ρ ( t , t ′ ) d t ′ , (10) with ρ ( t , t ′ ) as a weight function that plays the role of the quantity in square brackets in Eq. 9. Lacking additional information about the distribution of onset times, population sizes, transmission rates, or mixture rules, here we make the parsimonious assumption that ρ ( t − t ′ ) = Θ ( a t − t ′ ) Θ ( t ′ − t 0 + t ) , with Θ ( x ) as the Heaviside step function and a as a constant. In other words, the effective susceptible population, intervening in the transmission rate given by the average in Eq. 10, has been approximated by a set of independent SIR subsystems with a uniform distribution of onset times, between an initial time t 0 − t and a final time a t : s e f f ( t ) ≃ 1 a t − ( t 0 − t ) ∫ t 0 − t a t s S I R ( t + t ′ ) d t ′ . (11)

Approximation (Eq. 11), together with Eqs 4–6, pretends to describe the effects that the inhomogeneous transmission between subpopulations have on the epidemic propagation. The lower limit of the integral appearing in Eq. 11 is chosen so that in the average, there will always be contributions from subpopulations that evolve more slowly than the nominal subpopulation, defined as the one with t ′ = 0 , that is, s S I R ( t ) . On the other hand, the upper limit is chosen to allow subpopulations that evolve more rapidly. More precisely, the interval ( t 0 , − t ) corresponds to subpopulations that have not started their local outbreaks at time t, as can be the case of segments of the population or regions where the pathogen is not present at time t. The interval ( − t , 0 ) corresponds to subpopulations with delayed outbreaks, and finally, the interval ( 0 , a t ) corresponds to outbreaks apparently in advance to the nominal one. Since the upper limit a t is chosen to consider subpopulations where the pathogen propagates more rapidly than the average, this would imply that for those subpopulations, their local transmission rates β i j would be larger than the nominal value β 0 and their evolution equations should be written in terms of those transmission rates. Correspondingly, the number of susceptible individuals at time t has evolved more rapidly than that in subpopulations with smaller values of β i j . Here, however, we have described the evolution of all the subpopulations using the same nominal value β 0 for the transmission rates. This means that the evolution of those populations with β i j larger than β 0 has been approximated by populations evolving with β 0 but evaluated at latter times ( a > 0 ). As a consequence of the distribution of transmission rates, the time dispersion of the local outbreaks increases with time. This is reflected in the time dependence of the limits of the integral in Eq. 11. The value of the parameters t 0 and a is obtained from the fitting to the data points of the epidemic curve for a given population. The parameter a is a measure of the strength of the asynchronicity in the sense that the larger the distribution of transmission rates, the larger the parameter a.

To examine the validity of the proposed approximation, we apply it to the COVID-19 epidemics in different countries starting with the case of Mexico. Figure 1 (top panel) shows a semilogarithmic plot of s ( t ) and i ( t ) , as obtained from Eqs (4, 5), and compares it with the observed data (by confirmation date) [32] for both the total number of infected (cumulative cases), 1 − s ( t ) , and infectious individuals (active cases), i ( t ) , for the case of Mexico. In practice, active cases correspond to individuals whose symptoms started within the previous 14 days from the day a given data point was released. Let us stress that identifying the measured active cases with the number of infectious individuals is only an approximation. There is some uncertainty on the number of days an infected individual remains infectious, and also, there is some variability in the degree of infectiousness at different days. Nonetheless, here we are going to consider that the active cases correspond to infectious individuals. An initial basic reproduction number R 0 = 2.2 and recovery rate γ = 1 / 6 were used as parameters for phase 1, that is, before containment measures were adopted. After t = 20 days, once containment measures were adopted and their effects started to be noticeable, the initial reproduction number was changed to R 0 = 1.7 . The data have been plotted from March 9th, which corresponds to t = 0, where the data points show the beginning of a regular behavior (see Materials and Methods). The bottom panel shows the empirical new daily cases and the model fit. Remarkably, the model is able to correctly reproduce the empirical infectious cases (yellow line) and the new daily cases data, considering that the fitting was performed only for the total number of infected (blue line), and no independent fittings for each curve were required. We did not exclude the possibility that the other set of parameters (using a different value for γ) corresponds better to real values even if the yellow curve would not be as accurately fit. The inset in the first panel shows that traditional SIR using the same set of parameters fails to reproduce the epidemic curves. The exponential growth and characteristics of traditional SIR makes the comparison with real data, which shows algebraic growth, to eventually differ largely. Figure 2 shows the early growth of the total number of cases. After a first short exponential growth, an algebraic dependence with scaling law ≃ t 3.4 follows. Let us stress that, in contrast to other approaches [16, 17], this scaling law appears exclusively due to the time-dependence of the transmission rate β ( t ) , arising from the assumption of inhomogeneous transmission even without reactive population behavioral changes and before the susceptible depletion sets in. Figure 3 shows with lines the behavior of β ( t ) as obtained from the model, Eq. 7 with the approximation given by Eq. 11. The first section of this curve (blue line) corresponds to the beginning of the outbreak, when propagation was free, without containment measures. The second section of the curve (yellow line) corresponds to days after containment measures were implemented. In both cases, the transmission rate decreases with time. The data points are obtained solving for β ( t ) from Eq. 4, employing real data for Mexico [32]. We used the fact that the measured new daily cases equal − d s / d t , that the cumulative cases equal 1 − s ( t ) , and that the infectious cases are the measured active cases. The gray vertical line signals the end of the generalized voluntary quarantine (“Jornada Nacional de Sana Distancia”), and mitigation measures started to be released locally, depending on the strength of the epidemic in each state of the country. Accordingly, rebounds are verified after this date. Other approaches [14, 31, 34] attributed the time-dependence of β ( t ) to the reactive behavior of the population in response to the non-pharmaceutical mitigation measures and are introduced externally to the model. Even if no centralized information is provided about the presence of a disease, the impact of information diffusion, through first-hand observation and word of mouth, on the epidemic propagation can also be incorporated in a time-dependent basic reproductive number [35–37]. In contrast, in the present model, the temporal variation of β ( t ) is attributed to the inhomogeneous transmission arising from a distribution of transmission rates and of onset times of the local outbreaks, even if the behavior of the population remains unaltered and is obtained from the model dynamical equations. The right axis shows the corresponding values of the time-dependent basic reproduction number defined as [30], R 0 ( t ) ≡ β ( t ) / γ .

Since the first outbreak of COVID-19 in Mexico remains in progress (at the time of writing this manuscript), it is interesting to evaluate the performance of the model in countries where the first outbreak has nearly ended. We have chosen the eight European countries with the largest number of COVID-19 cases to further validate the model [38]. Figure 4 shows the epidemic curves for Belgium, Italy, France, Germany, Netherlands, Spain, Sweden, and United Kingdom. The model-fits capture the epidemic progression surprisingly well in all cases in spite of the different mitigation strategies applied by each country [39, 40]. Note that the fit starts to fail around t > 80 particularly in the case of Sweden where a second wave is already apparent. This is due to the fact that we have not taken into consideration the progressive release of the containment measures once the first outbreak has mostly finished. If a larger value for R 0 is used around t = 80 , then an increase of the epidemic curve can be obtained again, giving rise to a second wave. In Figure 4, we have used the present model to fit only the first wave in all cases.

Table 1 shows the parameters used to fit the model in Figure 4. In most of these countries, containment measures were taken at about 20–40 days after time t = 0 days, and thus, the epidemic curve can be divided into two sections, one before and one after the containment measures were adopted. The negative values of t 0 found in the first section, before mitigation measures were adopted, reflect the fact that a fraction of the subpopulations whose local outbreaks had not started contributes to the average in Eq. 11. On the other hand, the positive values of t 0 in the second section reflect the fact that when confinement measures were adopted, there were already a considerable number of infected individuals. Since the local epidemics in this second section evolve with an assumed attenuated value R 0 , then, in order to produce the same number of initial cases (which in reality were obtained with the original larger value R 0 ), larger times ( t 0 > 0 ) are needed in the SIR subsystems appearing in Eq. 11. In contrast, the case of Mexico shows a negative value of t 0 after the containment measures were adopted, reflecting the fact that the measures were taken very soon in the epidemic progression, when very few cases were present. This is also consistent with the fact that the evolution of the outbreak is taking considerably longer in Mexico than in most of the other countries analyzed. In most cases, the value γ = 1 / 6 was used as a parameter and the rest of the parameters were obtained by a least-squares fit. We cannot disregard that other choices for the parameter γ could also produce good fits but with different values for the rest of the fit parameters. Additionally, the fact that we are ignoring the incubation time of the disease may give rise to modified fit parameters; thus, they should be considered only as rough estimates.

The effect of the strength of asynchronicity is shown in Figure 5, which shows the evolution of the total number of infected, considering different values for the parameter a, for a given R 0 . We observe that the longtime total number of infected individuals will be highly dependent on a. The blue line corresponds to the fitted value of a for the case of Mexico, while the yellow and green lines are the predictions for a smaller (smaller a) or a larger (larger a) asynchronicity.

Summarizing, we have shown that the spread of the local outbreak onset times and transmission rates can be incorporated in an SIR formulation through the use of a time-dependent basic reproduction number. This quantity is obtained as a solution of the dynamical equations and not introduced externally. Thus, its time-dependence does not arise as a result of changes in the social behavior in response to containment measures or of the depletion of susceptible individuals. Instead, it is the result of an inhomogeneous transmission expressed as a time-dependent average of asynchronic SIR subpopulations (Eq. 11). We have shown that a simple assumption for the distribution of onset times and transmission rates can be at the origin of the algebraic growth observed at the early stages of the COVID-19 outbreak. This contradicts the common assumption that the early growth phase should be exponential in the absence of susceptible depletion or interventional measures. In the present model, containment measures contribute by decreasing the initial reproduction number or equivalently, the initial transmission rate. Other epidemic outbreaks also show early subexponential growth, and a number of potential mechanisms have been proposed to explain it. Among them are the spatial heterogeneity and clustering of contacts arising from the fact that the number of noninfected individuals in the immediate neighborhood of infecting agents is strongly constrained [41, 42]. Reactive population behavior has also been proposed to explain the changes that can gradually mitigate the transmission rate [16, 26, 31, 35–37, 43]. Related to these mechanisms, a range of mechanistic models that can reproduce the subexponential growth dynamics before susceptible depletion sets in have been proposed. These include models with gradually declining contact rate over time [14] and spatially structured models such as household-community networks [26], among others. However, for real epidemics, the underlying mechanisms governing the subexponential growth have been difficult to disentangle, and the matter remains debated [26, 41, 44].

The model predicts a final number of cumulated cases that is substantially smaller than that predicted for homogenous well-mixed populations in agreement with models that predict smaller disease-induced herd immunity when population heterogeneity is taken into account [19].

The present model was validated using the case of the COVID-19 outbreak in different countries, with an emphasis in the case of Mexico. As of the origin of the inhomogeneities in the transmission rates and onset times, one can propose different underlying mechanisms, including the transmission between geographical dispersal subpopulations as individuals travel among them [6, 8, 15]. Another possible source of inhomogeneities is the existence of different social cohorts, with transmission rates between them that are lower than those between individuals of the same cohort [19, 45]. In general, the overall connectivity between subpopulations will determine to a large extent the rate of propagation of the viral agent and the final number of cases. Subpopulations more efficiently connected between them will have smaller dispersion of onset times and transmission rates.

The model is consistent with mitigation strategies, consisting in the design of containment mechanisms oriented to increase inhomogeneities, for example, imposing travel restrictions for long-distance routes, partially isolating subpopulations from each other [46, 47], or by increasing the time an individual takes to move from one subpopulation to the other. Furthermore, travel restrictions could be targeted to highly connected individuals. The reason is that the effective number of susceptible individuals, s e f f ( t ) , that are the responsible for the rate of change of s ( t ) , decreases with the implementation of such mitigation mechanisms. For countries where a strict quarantine is impracticable, this could be a more realistic alternative. We believe that the insight obtained from the present model may be useful for planning non-pharmaceutical responses to better mitigate or block the overall spread of an epidemic outbreak.

Identifying the first infectious case is a difficult task, so considering t = 0 days as the day the first case was reported, does not seem to be appropriate. For this reason, before trying to fit the original data points, we first shifted them in time so that the time t = 0 days corresponds to the point where a discernible regular behavior starts to appear. This is shown in Supplementary Figure S1 where the data for Mexico that were disregarded for fitting purposes correspond to dates ranging from February 27th, the day of the first detected case, until March 9th. Those data show an irregular behavior probably due to the fewer number of cases and large relative fluctuations.

From a practical point of view, it is numerically easier to approximate s S I R ( t ) in Eq. 11 by the function s S I R ( t ) ≃ 1 − 1 ( c + q e − t / τ ) 1 / ν , (12) where the parameters c, q, ν, and τ are obtained from the fitting to the exact form for s S I R ( t ) once β 0 and γ are known. The initial conditions s S I R ( 0 ) = 1 − 1 / N i , i S I R ( 0 ) = 1 / N i , and r S I R ( 0 ) = 0 were assumed. Supplementary Figure S2 shows an example for s S I R ( t ) as obtained from the standard SIR model (solid line) and its approximation Eq. 12 (dashed line). By incorporating Eq. 12 into Eq. 11 and the last result into Eqs 4–7, the evolution of s ( t ) , i ( t ) , and r ( t ) for the whole population is finally obtained.

The fitting of the model to the data, before containment measures were taken, was done choosing the value of a spline function that softens the data taken arbitrarily at time t = 4 days as an initial condition, as shown in Supplementary Figure S3. In this way, possible misfitting due to fluctuations of the data at the very first data points is minimized. The fitting of the second section of the model, that is, once containment measures were taken at day 20 in the case of Mexico, considered as initial values of the model, the values obtained from the first section evaluated at day 20.

For the case of Mexico, the parameters used were R 0 = 2.2 and γ = 1 / 6 , values generally accepted for the initial regime without containment measures, and R 0 = 1.7 for the second regime after containment measures were adopted. The first fitting process considered the first data points up to t = 20 days with obtained fitting parameters t 0 ≃ − 26.8 days and a ≃ 6.5 . Analogously, the fitting for the second regime consisted of data from t = 21 days to the last data point available. For this second regime, the obtained fitting values were t 0 ≃ − 896 days and a ≃ 11 . The results were largely insensitive to variations in the exact day of transition between the two regimes.