Lockdown Measures and their Impact on Single- and Two-age-structured Epidemic Model for the COVID-19 Outbreak in Mexico

Lockdown Measures and their Impact on Single- and Two-age-structured Epidemic Model for the COVID-19 Outbreak in Mexico J. Cuevas-Maraver Grupo de F´ ısica No Lineal, Departamento de F´ ısica Aplicada I, Universidad de Sevilla. Escuela Polit´ ecnica Superior, C/ Virgen de ´ Africa, 7, 41011-Sevilla, Spain Instituto de Matem´ aticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain P. G. Kevrekidis, Q. Y. Chen, and G. A. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-4515, USA V´ ıctor Villalobos-Daniel National Center of Disease Prevention and Control Programs - CENAPRECE, Avenida Benjam´ ın Franklin, 132, 11800-Ciudad de M´ exico, CDMX Z. Rapti Department of Mathematics and Carl R.Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign Y. Drossinos European Commission, Joint Research Centre, I-21027 Ispra (VA), Italy The role of lockdown measures in mitigating COVID-19 in Mexico is investigated using a comprehensive nonlinear ODE model. The model includes both asymptomatic and presymptomatic populations with the latter leading to sickness (with recovery, hospitalization and death possibilities). We consider the situation involving the imposed application of partial social distancing measures in the time series of interest and find optimal parametric fits to the time series of deaths (only), as well as to that of deaths and cumulative infections. We discuss the merits and disadvantages of each approach, we interpret the parameters of the model and assess the realistic nature of the parameters resulting from the optimization procedure. Importantly, we explore a model involving two sub-populations (younger and older than a specific age), to more accurately reflect the observed impact as concerns symptoms and behavior in different age groups. For definitiveness and to separate people that are (typically) in the active workforce, our partition of population is with respect to members younger vs. older than the age of 65. The basic reproductive number of the model is computed for both the single- and the two-population variant. Finally, we consider what would be the impact on the number of deaths and cumulative infections upon imposition of partial lockdown (involving only the older population) and full lockdown (involving the entire population). 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 2 I. INTRODUCTION COVID-19, the disease caused by the novel coronavirus SARS-CoV-2 has, as of August 2020, affected 216 countries and changed the daily lives of billions of people [1]. It has at the same time been the focus of numerous studies, both clinical and mathematical in nature. The study of compartmental models that address the spreading of such epidemics has a time-honored history since the seminal contribution of [2], which by now has been summarized in various reviews [3] and books [4–6]. In recent years, variations to such models focusing on the particularities of coronaviruses have been incorporated, such as the role of asymptomatic carriers of the virus, both as regards earlier CoV examples, such as MERS (see for a related example the work of [7]), and lately in the case of COVID-19 (see for a related example the work of [8]). Following on some of these more recent developments, the present study focuses on a comprehensive compartmental epidemiological model that takes into account some of the intricacies of COVID-19, while considering its applicability to an urgent and important case example, namely the country of Mexico. More specifically, the model is an extension of the standard SEIR (Susceptible, Exposed, Infectious, Recovered) model that includes a presymptomatic stage, during which a person experiences no symptoms, but is nevertheless infectious [9]. The proposed mathematical setup also accounts for both the asymptomatic infectious and symptomatic infectious individuals. Asymptomatic infectious cases have been found in numerous studies and reports argue that they may be significantly under-reported [8, 10], which may complicate mitigation efforts such as contact tracing and self-isolation. Those with severe disease symptoms may require lengthy hospitalization, which has strained the health system of many countries [11]. In light of that, the model also includes a compartment describing the hospitalizations. Another distinctive feature of the disease is the nonhomogeneity with which it manifests in different age groups, especially as it pertains to symptom severity and mortality risk [12, 13]. Other factors, such as preexisting conditions and intergenerational contacts may also play a role [14]. While population age-structure may be often averaged out and deemphasized in numerous modeling attempts [15], in our model we choose to consider both a single age-group and a two age-group version of the model. The rationale behind this choice is the multifold inhomogeneity in the population of various countries (including our example of interest). Firstly, as mentioned above, the severity in younger people (especially children [16]) is smaller than that in adults. Secondly, older and more vulnerable people may shed more viral particles, thus being more infectious [17]. Thirdly, contacts per day [18] and the contact network itself of older people is different from those of younger people. The partition especially between professionally active (i.e., non-retired) individuals and retirees is important in connection to the above two points, both as regards the differential in average number of contacts of these two groups, and as regards the potential vulnerabilities thereof. As a case study, we focus on the COVID-19 outbreak in Mexico. While studies for Mexico based on mathematical models exist, they differ from ours in several significant ways. Some ignore social-distancing and other mitigation measures [19], others focus on the estimation of R0and infections, using a Bayesian hierarchical model [20], and yet others have since become outdated [21]. Mexico faces a unique challenge, due to the prevalence of COVID-19 risk factors, such as obesity, diabetes and hypertension, among its population [22]. This is being reflected in the reported data and our predictive results, which show almost as many fatalities in the <65 years old group as in the >65 years old group. We feel that these particular features of the Mexican population in conjunction with the large number of infections and especially of deaths in the country warrant an examination through the prism of different age-structured models (e.g., single-age vs. two-age models; in future studies, possibly further partitioning may be of interest) and an assessment of the potential impact of lockdown measures in the immediate future. Following the formulation of the single-population model and the results obtained through it in section II, we continue with the two age-group model in section III. In each case, we obtain the optimal parameters of the model in matching the available data regarding deaths, which are considered to be the single most reliable piece of available information. We do discuss the advantages and disadvantages of potentially matching the number of cumulative infections (and the number of deaths). Once the optimal fitting parameters are obtained we assess the impact in both deaths (but also cumulative infections) of immediate lockdown measures in either the case of the entire population or in that of just the older age group. The prediction of the model is that thousands of lives may be saved in just the following month alone, should such measures be imposed effective immediately. In section IV we summarize our conclusions and present ideas for future investigation. II. SINGLE-POPULATION MODEL A. Equations In the model presented herein, we modify somewhat the setup of the earlier work of a subset of the present authors [23], by incorporating the effect of presymptomatic individuals. More concretely, we start with a susceptible (S) population that can become exposed (E) to the SARS-CoV-2 virus upon interaction with three categories of already infected individuals: (a) the presymptomatic (P), individuals who are infected, infectious, and eventually will develop symptoms; (b) the asymptomatic (A), individuals who are infected, infectious, and will not develop (clinical) symptoms; and (c) the symptomatically infected/sick (I) population members carrying the virus (infectious). Upon such interaction, the susceptible become exposed to the virus. 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 3 Notice that in the susceptible population, there is the potential for “removal” (i.e., death) by other causes, via the last term of the corresponding ODE proportional to µ, however this term is practically irrelevant for our purposes. For this reason, we set µ= 0. Once a member of the population becomes exposed (E), a latent period (τl= 1/σ1) of the virus follows (expected to be in the vicinity of 3 days [24]), during which the exposed population is infected but not infectious. After this period, we assume that the host can naturally be partitioned to either asymptomatic (A) or presymptomatic (P). The fraction of the former is φ, while of the latter 1−φ. While both A and P play a role (along with the infected I) in further transmitting the virus, and indeed A have been argued to play a crucial role [9, 10], it is only P that will present symptoms after an additional time scale, the preclinical period τp= 1/σ2. The incubation period, i.e., the period from infection to the development of (clinical) symptoms, for this partition τinc =τl+τp= 1/σ1+ 1/σ2is of the order of 5 days [24], and represents the time till the onset of symptoms. From there on, the asymptomatics A continue as if nothing happened, given that they have minimal or no symptoms. Their path is only towards recovery with a characteristic rate MAR, or with a time scale representing the asymptomatic infectious period τA inf = 1/MAR that is typically expected to be on the order of 7days. This is why in most countries quarantine is expected to last around 14 days (7 days of infectiousness till recovery and another 7 till recovery of any other person that may be infected among the ones in close contact with the person of interest). In order to distinguish those recovered from asymptomatics (which cannot be directly monitored, unless extensive testing is performed in a community) from those coming from a path of symptoms/sickness (which can be –at least partially– monitored), we denote those recovered from A as AR. The path of the presymptomatics P is more complicated. Indeed, these may still recover (R) without the need for hospitalization and without severe manifestation of symptoms (going through class I). However, in a fraction γof the cases hospitalization is needed. Recovery or hospitalization in the model is associated with a time scale 1/M, the symptomatic infectious period τI inf = 1/M. Subsequent steps involve a fraction ωof the hospitalized that die, over a time scale 1/ψ and a fraction 1−ω that recover over a time scale 1/χ. The above offers an, in principle, complete description of the modeled quantities within our system. We should note that βrefer to the coefficients of interaction between A (or P) and S, as well as I and S, leading to new infections; these are, respectively, βSA and βSI . We note here, that we somewhat abuse notation as far as β’s are concerned. Namely, it is relevant consider the transmission rate β=˜ β N, where Nis the total population. What we report in the tables that follow is actually ˜ β. In terms of the relevant time scales, the infectious period associated with symptomatically infected is τI inf = 1/M, with asymptomatic individuals τA inf = 1/MAR, whereas the infectious period associated with presymptomatic individuals is τP inf =τp+τI inf = 1/σ2+ 1/M. The transcription of the above steps in equations leads to the following ODEs: dS dt=−βSA(t)S(A+P)−βSI (t)SI −µS dE dt=βSA(t)S(A+P) + βSI (t)SI −σ1E dP dt= (1 −φ)σ1E−σ2P dA dt=φσ1E−MARA dAR dt=MARA dI dt=σ2P−MI dH dt=γMI −(1 −ω)χH −ωψH dR dt= (1 −γ)MI + (1 −ω)χH dD dt=ωψH A schematic diagram of the model is shown in figure 1. In the practical aspects of what follows, we will consider the data for Mexico, with a total population of 127,575,528 people (in 2019) and more than 430,000 confirmed cases and 47,000 deaths by the beginning of August 2020. Data were taken from “Direcci´ on General de Epidemiolog´ ıa (DGE)” of “Gobierno de M´ exico” [28]. These data have the particularity of including each clinical case, from which we retrieve three basic pieces of information: the date when symptoms start, the death date (if applicable) and the age. The data are updated on a daily basis. As the way of measuring always implies an underestimation of the number of cases and deaths in the days close the report’s (especially because of the delay in death communications to the DGE), we have performed fits up to dates about 20 days from the report date (i.e. the report date is July 29 and fits are performed until July 10). 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 4 SEP A AR I H R D 𝜔𝜔𝜓𝜓 MAR 𝝋𝝋𝝈𝝈1 (1-𝝋𝝋)𝝈𝝈1𝝈𝝈2(1-𝝲𝝲)M 𝝲𝝲M (1-𝜔𝜔)𝟀𝟀 𝞫𝞫SII 𝞫𝞫SAA 𝞫𝞫SPP EAAR Latent period 1/𝝈𝝈1 Infectious period 1/MAR EPIR Infectious period 1/𝝈𝝈2+1/M Incubation period 1/𝝈𝝈1+ 1/𝝈𝝈2 Preclinical period 1/𝝈𝝈2 FIG. 1. Schematic diagram of the single-population model (left panel) and a diagram of the two disease progression pathways: the asymptomatic and the symptomatic one (right panel, based on a relevant variation of [8] adapted to the specific compartments and time scales of the present model). Our principal diagnostic quantities in order to obtain the optimal parameters for the model will be the time series for deaths D(t), but we will also monitor the cumulative infections. The latter in the realm of the present model amounts to C(t) = I(t) + H(t) + R(t) + D(t), i.e., the sum of the individuals going through the part of the network involving the symptomatically infected. We take t= 0 as March 22, and fit data to July 10. This reflects our effort to be (in terms of the total numbers of both diagnostics) well within the “well mixed” regime from the beginning of the infection where the ODEs of interest are expected to be relevant. We are particularly interested in the effect of non-pharmaceutical intervention strategies on the spreading and development of the disease. Such strategies render the transmission rates time dependent. On April 21 (tq= 30) a light form of social distancing was enforced in Mexico. In addition, we will consider a number of scenarios according to which a more severe lockdown may be imposed on August 10 (tL= 141). We are particularly interested in how non-pharmaceutical intervention strategies modify the spreading and development of the disease. Such strategies render the transmission rates time dependent. On April 21 (tq= 30) a light form of social distancing was enforced in Mexico. In addition, we will consider a number of scenarios according to which a more severe lockdown may have been imposed on August 10 (tL= 141). The effect of interventions on the overall transmission rate βmay be estimated by considering biological and physical properties of expelled respiratory droplets, which are the carriers of the pathogens and specifically of SARS-CoV-2. The transmission rate is usually written as β=cp, with cthe number of contacts per day a susceptible individual has and pthe transmission probability. For pathogen transmission via infectious respiratory droplets (of diameter d), be they airborne (airborne transmission) or settled (contact transmission), the transmission rate has been expressed as, cf. Refs. [25] and [26] (for airborne transmission) β=βd κd αd ,(1) where κdis the respiratory droplet emission rate (viral shedding) by e.g. breathing, speaking, coughing, sneezing, αdis the droplet effective removal rate, by e.g. gravitational settling, ambient airflow, pathogen inactivation, and βdthe transmission rate per deposited respiratory droplet. The subscript “d” refers to a specific droplet size. As we will use it Eq. (1) to estimate how the overall transmission rate changes we will neglect its complex dependence on d. Lastly, the droplet transmission rate βddepends on the number of effective contacts a susceptible has with other individuals, the number of pathogens contained in an infectious droplet (its pathogen load) , and the virus transmission probability per inhaled/deposited droplet. Social distancing, or other lockdown measures that restrict human mobility, decreases the average number of daily contacts, and possibly the average duration of contact, thereby decreasing βd(and consequently β). In what follows we denote this fractional decrease as ηSA and ηSI . Another common intervention measure is the use of surgical face masks. Milton et al. (2013) [27] argued that their use produced a 3.4-fold reduction in viral aerosol shedding. Face masks also render physical contact of infected hands with susceptible areas on an individual’s face (mouth, eyes, nose) more difficult. They, also, modify the expelled air flow, with consequential effects on droplet transport and dispersion in the environment and their airborne lifetime (thus their removal rate). We suggest that the combined effect of wearing surgical face masks is primarily conditioned by the decrease in viral shedding. We, thus, estimate that their use (if the whole population used them continuously and correctly fitted) would decrease the overall transmission rate to ∼0.2β. The effect of lockdown measures is modelled herein by the parameter ζ, Eq. (2): we set the after lockdown transmission rate to be ζβ. This parameter incorporates the effect of all lockdown measures, including the requirement that face mask be worn. As the previous estimate leads to a considerable decrease in the transmission rate, we decided to be conservative and we varied ζfrom 1.0 to 0.5. 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 5 These features should be kept in mind, as we aim not only to capture the current trend of the pandemic, but also to suggest mitigation strategies that may reduce the cumulative infections, as well as the fatalities as a result of COVID-19 in the time series of interest, namely in the case of Mexico. In what follows, we suggest possible scenarios of social distancing, and we accordingly evaluate their impact as regards the potential reduction induced in the number of deaths and cumulative infections for the time series of interest that is obtained from the DGE data. To explore the impact of further lockdown effects to mitigate the spread of the infection, we will consider the following time dependence for the β’s: βSI (t) =βSI ηSI + (1 −ηSI )1−tanh[2(t−tq)] 2 βSA(t) =βSA ηSA + (1 −ηSA)1−tanh[2(t−tq)] 2+ηSA(ζ−1)1 + tanh[2(t−tL)] 2(2) with 0< ζ ≤1, as argued, and tq,tLpreviously specified. For ζ= 1, Eqs. (2) reduce to the equations modelling the decrease in the number of personal contacts when a light form of social distancing is imposed (no lockdown, no requirement to wear face masks). Figure 2 shows the time-dependence of the β’s. The first jump in the value of βSI and βSA (i.e., the transmission rates to infected and presymptomatic-asymptomatic virus carriers) occurs at tq. The next jump occurs at tLand is assumed to take place only for the presymptomatic or asymptomatic carriers of the virus through the imposition of lockdown restrictions at t=tL. As our data will show below, the growth of fatalities as well as infections in the country is considerable and thus our model suggests the relevance of the application of significant lockdown restrictions so that the progress of the infection be substantially curbed. In the parameters of the single-population model, we assume the following constraints: •ηSI ≤1 •ηSA ≤1 •5≤1/σ1+ 1/σ2≤6 The first two are rather trivial (and without loss of generality) implying that we go from a value of the transmission rate β, to a lower value β×ηin both the interactions of S with A (or P) and in those of S with I. The third constraint is based on observations associated with SARS-CoV-2 [24] positioning the incubation period associated with this virus between roughly 5and 6days from the start of exposure. FIG. 2. Single-population model. Transmission-rate time dependence (βSI and βSA). For the asymptomatic carriers of the virus the transmission rate, which initially decreases due to social distancing (reflected in ηSA) is further assumed to decrease by a factor ζ(shown for different cases), reflecting the effect of obligatory wearing of (effective) face masks or other lockdown-based decrease of contacts. For the infected population instead, we assume that βSI decreases only once due to social distancing, (ηSI ), but it does not further decrease, given the self-isolation of such individuals due to the presence of symptoms. 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 6 B. Results We now present the results determined by an optimization procedure that identifies under the above minimal constraints the optimal parameters of the model in comparison to the time series of data obtained from the Government of Mexico [28]. Table I shows the optimal parameters found from minimizing the norm N=X iα1|log(Cnum(ti)) −log(Cobs(ti))|+α2|log(Dnum(ti)) −log(Dobs(ti))|(3) with α1= 0,α2= 1 as displayed in Fig. 3. We will compare these results shortly with the case of α1=α2= 0.5shown in Fig. 4. The rationalization of these two choices is as follows. In the former case, we take the view that the only “ground truth” data is that of the deaths; in fact, even those can be under-estimated (deaths with “suspicion” of COVID-19, but no definitive test) or –perhaps less likely– over-estimated (deaths attributed to COVID-19 without explicit testing), but here we will assume that this is the most well-defined piece of data, as is generally expected to be the case. On the other hand, infections are broadly expected to be under-reported. This is because many of the cases with symptoms do not get to be serious enough to lead to hospitalization or to be reported. In that light, it is expected that assuming the deaths as ground truth, we should expect to find a significant over-estimation of the number of infections (we return to this point below). If, on the other hand, we ”force” the model to match the current reporting of infections, then we will end up with a better approximation ”on average” to both curves but with a potential adverse by-product in the resulting number of deaths that we will discuss below. This will be due, in addition, to issues of measurement errors in the setting of cumulative incidences, as discussed, e.g., in [29]. Our observation in Fig. 3 is that the model can provide an excellent fit to the total number of deaths, but in that case, there is a considerable over-prediction of the number of infections, presumably because numerous of the incurred infections are not reported in the official data. The results can be compared with Fig. 4 where both deaths and cumulative infections are attempted to be fitted, i.e., α1=α2in the minimization procedure above. Here, we see that while the model is capable of doing a very adequate job in following C(t), it also does a reasonable job of capturing D(t).Nevertheless, there is a caveat to the latter. A closer observation of the semilogarithmic scale of the graph will lead the astute reader to observe that as the model is trying to juggle the optimization of both time-series, it slightly over-predicts deaths early on, slightly under-predicts them at the middle of the time series and eventually slightly overpredicts again at longer times, likely predicting a much more catastrophic scenario (with multiple hundreds of thousands of deaths at the end of the examined evolution) than is warranted by the data trends. For this reason, we will stick to the consideration of the former case of Fig. 3 hereafter. Before we discuss the implications of mitigation strategies, let us briefly comment on the optimal parameter values identified by the model, as illustrated in Table I. The latent period 1/σ1is indeed found to be in the vicinity of 3 days (2.8765), while the total incubation period is comfortably within the prescribed interval of 5-6 days (1/σ1+ 1/σ2≈5.3836). The time scale of recovery for asymptomatics is a little under 7 days as expected (1/MAR = 6.0983), while the time scale of going from symptomatic infected to hospitalization is close to 4 days which is also fairly reasonable (1/M = 3.6221). The model predicts a large fraction of asymptomatics (φ= 0.8134) which is in line with discussions such as that of [10]. However, such results should be taken with a grain of salt. This is because of issues associated with the notion of identifiability [30]. In particular, a systematic analysis of the model in connection to identifiability (the formal mathematical details of which are outside the scope of the present work) suggest that φitself will not end up being an identifiable parameter, but rather the product of φwith β’s will be on such. As a result, we do not significantly focus on the larger value obtained for φ, but we do note it. Interestingly, the model predicts that roughly 30% of those symptomatically infected need hospitalization, while the rest recover. Of those hospitalized, nearly 30% result in fatalities, leading to a death percentage of about 10% among those that present symptoms. Indeed, this is a rather significant fatality percentage accounting for the large number of deaths in the population. The average time scale of recovery upon hospitalization is about 10.7days (1/χ), while that of death (1/ψ) is found to optimally be near 11.5 days. However, here too we should highlight that ωand ψare not independently identifiable (only their product is) and neither is (1 −ω)and χ(again only their product is). Hence, the relevant optimization parameter values should be considered as coming together with the corresponding caveats (and the relevant products as being the genuine model derived quantities), despite their realistic individual values. It is relevant to note that at the initial time of the model, it is not known how many members of the population are exposed, asymptomatic or presymptomatic. We thus optimize those parameters too, obtaining the last 3 entries of the Table. Our aim is now to explore the impact of mitigating measures reducing the obtained optimal case by a factor of ζ < 1(ζ= 1 is the current optimization result without additional measures). Table II summarizes the results of Fig. 5, in which the predictions for ζ= 0.9,0.7,0.5are displayed. It can be seen that a decrease of βby a factor between 0.9and 0.5could have a catalytic result as concerns the predictions of the model both for the deaths and as regards the cumulative infections. In particular, a decrease of ζby a factor of even as little as 0.9will reduce the number of deaths by more than 800, only within the time frame between August 10 and September 10, while a decrease by a factor of 1/2is predicted within this single population model to save close to 4000 lives in this interval alone. The corresponding effect to the number of infections is perhaps even more striking (also accounting for the fact that many of these infections could lead to fatalities at a later stage). In particular, a decrease of ζ to 0.9leads to about 24000 less infections, while to a factor of 1/2would lead to about 100000 less infections again just within 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 7 TABLE I. Optimal parameters for the single population model. For a discussion thereof, see the text. Parameter Symbol Optimal value Parameter Symbol Optimal value Transmission rate [per day] βSI 1.4789 Infectivity period (A) [days] τA inf 1/MAR 6.0983 Transmission rate [per day] βSA 0.1982 Conversion fraction (I to H, R) γ0.2991 Social distancing effect ηSI 0.5806 Conversion fraction (H to R, D) ω0.3095 Social distancing effect ηSA 0.4585 Recovery period(H to R) [days] 1/χ 10.6927 Latent period [days] τl1/σ12.8765 H to D period [days] 1/ψ 11.5331 Preclinical period [days] τp1/σ22.5071 Initial population fraction (E) E(0)/I(0) 2.0851 A/P partitioning φ0.8134 Initial population fraction (A) A(0)/I(0) 2.0532 Infectivity period (I) [days] τI inf 1/M 3.6221 Initial population fraction (P) P(0)/I(0) 0.5462 FIG. 3. Single-population model. Number of cases C(t)(left) and of deaths D(t)(right) found by minimizing norm (3) with α1= 0 and α2= 1. To produce the fit we have used the data until July 10. The prediction of the optimized model is shown by the solid line, while the official time series [28] is given by dots. FIG. 4. As in the previous figure, but for α1=α2. The number of cases C(t)is captured significantly better, but the lower accuracy in capturing D(t)and its implications are further discussed in the text. this 31-day period. One can clearly see the significant potential impact of further lockdown measures, a feature that may be worthwhile to factor into further public health considerations. 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 8 TABLE II. Single-population model. Predictions for Cand Dat September 10 if lockdown measures had been applied on August 10. ζ= 1 ζ= 0.9ζ= 0.7ζ= 0.5 C(total) 1090578 1066436 1024199 988931 C(from Aug. 10) 292666 268524 226289 191022 D(total) 83593 82744 81203 79848 D(from Aug. 10) 24362 23514 21972 20617 FIG. 5. Single-population model. Evolution of the number of cases (left) and of deaths (right) for different ζif lockdown had been applied on August 10 (t= 141). Notice the significant curbing of the pandemic as a result of such intervention measures, especially so at the level of long term effects in the case of ζ= 0.7and ζ= 0.5. As our final comment about the single population variant of the model, we note that an important quantity in epidemiological models of this general type is the notion of the basic reproductive number R0; see, e.g., [3]. The number effectively represents the expected new infections (so-called secondary infections) from a single infection in a population where all subjects are susceptible. Using the so-called next-generation approach [31], we can determine the basic reproductive number on the basis of the parameters of the model and the susceptible population [see the Appendix for details], according to the expression: R0=βSA(1 −φ) σ2 +βSAφ MAR +βSI (1 −φ) MS∗,(4) with S∗the initial susceptible population, namely S∗= 1 as parameters βare normalized by N, as discussed above. For the parameters in Tab. I, one obtains a value of R0= 2.0755. The effective reproductive number at the beginning of social distancing measures (t=tq) is Re= 1.5745 which is unable to mitigate the pandemic effects. In order that Re<1, it is needed that ζis smaller than 0.8506. That is why among the case examples that we considered, those with ζ= 0.7and 0.5present a significant decrease in the number of deaths not only imminently (i.e., within the interval of August 10 to September 10) but also over the longer scale prediction of Fig. 3. III. TWO-POPULATION MODEL A. Equations We now turn to the two-population variant of the model. Recall that due to the different structural characteristics of the two populations as considered herein, namely below and above 65 years, we expect that this model will be more adequate in capturing both the deaths and the cumulative infections of the full population. This is because on the one hand, the younger population in our considerations is more active (belonging typically in the workforce), hence carries a different number of contacts. On the other hand, the older population has its own vulnerabilities to the impact of the virus SARS-CoV-2 and the associated disease, 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 9 namely COVID-19. On the other hand, as explained in the Introduction, the prevalence of various risk factors within the Mexican population [22] render this partition even more relevant to consider towards capturing the detailed data trends. Firstly, we discuss the mathematical structure of the model. Here, we basically assume that each of the populations has its own set of parameters. Thus, the superscript ywill denote the parameters associated with the younger population, while the superscript owill be connected to the older population. Naturally, the number of variables (S, E, A, P, I, H, R, D)now doubles with each part having a younger and an older component. It is worthwhile that for some quantities that are associated with the virus, such as σ1and σ2, characterizing the latent and incubation time thereof, we assume these to be independent of age. Lastly, we explore a mildly anisotropic variant of the model where in terms of the interactions βo≡βoy =βyo =βoo and βy≡βyy. That is to say, we assume that the older population has a different interaction within itself and with the younger individuals, than the younger members of the population between themselves [18]. This is a reasonable assumption under the present conditions where the more sensitive older members of the population are advised to reduce their interactions. While, in principle, we could have used a fully anisotropic variant of the model, we find it relevant to attempt to reduce the overall number of model parameters, hence the above choice. dSy dt=−βyy SI (t)SyIy−βyy SA(t)Sy(Ay+Py)−βyo SI (t)SyIo−βyo SA(t)Sy(Ao+Po) dEy dt=−σ1Ey+βyy SI (t)SyIy+βyy SA(t)Sy(Ay+Py) + βyo SI (t)SyIo+βyo SA(t)Sy(Ao+Po) dPy dt= (1 −φy)σ1Ey−σ2Py dAy dt=φyσ1Ey−My ARAy dAy R dt=My ARAy dIy dt=σ2Py−MyIy dHy dt=γyMyIy−(1 −ωy)χyHy−ωyψyHy dRy dt= (1 −γy)MyIy+ (1 −ωy)χyHy dDy dt=ωyψyHy dSo dt=−βoo SI (t)SoIo−βoo SA(t)So(Ao+Po)−βoy SI (t)SoIy−βoy SA(t)So(Ay+Py) dEo dt=−σ1Eo+βoo SI (t)SoIo+βoo SA(t)So(Ao+Po) + βoy SI (t)SoIy+βoy SA(t)So(Ay+Py) dPo dt= (1 −φo)σ1Eo−σ2Po dAo dt=φoσ1Eo−Mo ARAo dAo R dt=Mo ARAo dIo dt=σ2Po−MoIo dHo dt=γoMoIo−(1 −ωo)χoHo−ωoψoHo dRo dt= (1 −γo)MoIo+ (1 −ωo)χoHo dDo dt=ωoψoHo (5) In a natural extension of what we assumed for the fitting of the single population model, we assume that as a result of the 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 16 chance to be hospitalized. In other settings where there is heavy testing involved, asymptomatics may be counted in the reported infections: while this is not an explicit assumption of the model, we did not count any fraction of asymptomatics in the cumulative infections when presenting the relevant comparison. Additionally, also, some of the hospitalized individuals may transmit the virus (e.g. to medical personnel) despite the much more substantial health and safety protocols applicable within hospitals. In any event, we consider these features to be the exception rather than the rule and hence have excluded them from our more mainstream considerations. Within the realm of the model, we have exposed the meaning of the relevant parameters (e.g., transmission rates, latent and incubation times, fractions of asymptomatics vs. presymptomatics, of hospitalized vs. directly recovered, and of recovered vs. dying individuals at the hospital; also the time scales of the latter partitions were considered). I.e., we have attempted to assign an epidemiological meaning to our different parameters and to examine the associated results of the optimal fit of these parameters to the data from a specific time series to illustrate the “reasonable” nature of the findings. Notice that in addition to explaining the fitting process (to either deaths or deaths and cumulative infections), we have taken the approach of using a minimal number of assumptions to avoid constraining the system to the degree possible. We have also illustrated how the model can be partitioned to different age groups, based on the data that may be available for the country or region of interest. Here, we have opted to consider the simplest partition to 2-age models, yet while tedious, it is structurally straightforward (and of some interest in its own right) to generalize considerations to many age group models. As our prototypical illustration of choice, we have used the case of data from Mexico which have been available through [28]. This is a case where a significant number of cases has arisen and the consideration of potential further lockdown measures is an important topic of ongoing debate. Indeed, our findings suggest that the present measures appear not to be sufficient to mitigate the catastrophic consequences of the pandemic, since the current value of the basic reproductive number of the epidemic is R0>1and hence the situation appears to need further mitigation measures and strategies to avoid significant loss of life. In that vein, we have discussed in the realm of the model what consequences different measures would have at the level of deaths and of cumulative infections. It was found, e.g., that a reduction of transmission rates by a factor of around 1/2via social distancing or related measures would lead to a nontrivial curbing of the rampant growth of both C(t)and D(t). The relevant reductions could of the order of 100000 in terms of infections and of more than 2000 in terms of deaths in the interval of the next 30 days alone. While our results are only suggestive (and relevant with the context of the model), we hope that they may be of some value towards policy considerations in the near future. Naturally, there are numerous directions that are worthwhile to consider towards extensions of the present work. A natural feature of many of the models (e.g., associated with the US [15], but also elsewhere) is the incorporation of uncertainty. We are currently in the process of buidling into the formulation an uncertainty quantification framework on the basis of polynomial chaos considerations [34] and the use of suitable distributions for quantities such as the transmission rates. Constructing a robust such framework would be of considerable value to models such as the one proposed herein. In addition, as indicated in numerous cases, there are data broken down by age groups (e.g., not only for Mexico, but also for other countries such as Portugal [35], etc.). Clearly, a generalization of the model that considers the data by decade would offer a more complete and systematic picture of the impact of COVID-19 to different sub-populations and hence their potential risk. This would offer, in turn, a clearer picture of which age groups to attempt to protect and would be worthwhile (even if somewhat cumbersome). Lastly, as different countries are delving into a re-opening exercise, the formulation of a meta-population model with different hubs and a quantification of the transport coefficients between these [36] would be central towards going beyond the well-mixed assumption and factoring in a spatial structure and transportation features within the model. Such directions are currently under active investigation and will be reported in future work. Appendix A: R0calculations 1. One-population model We use the next generation matrix approach to find R0[31]. We define the relevant vectors: F=                 βSAS(A+P) + βSI SI 0 0 0 0 0 0 0 0                 ,V=                 σ1E −(1 −φ)σ1E+σ2P −φσ1E+MARA −σ2P+MI −γMI + (1 −ω)χH +ωψH βSAS(A+P) + βSI SI +µS −MARA −(1 −γ)MI −(1 −ω)χH −ωψH                 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 17 We then focus on the 5infectious/infected compartments (E,P,A,I,H) and ignore the rest (S,AR,R,D). We find the Jacobians of F,Vwith respect to E, P, A, I, H in the order in which they appear. This will yield two 5×5matrices: F=       0βSAS∗βSAS∗βSI S∗0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0        (A1) V=       σ10 0 0 0 −(1 −φ)σ1σ20 0 0 −φσ10MAR 0 0 0−σ20M0 0 0 0 −γM (1 −ω)χ+ωψ        (A2) The basic reproductive number is the spectral radius of FV −1which in our case is R0= (1 −φ)βSAS∗ σ2 +φβSAS∗ MAR + (1 −φ)βSI S∗ M.(A3) The first term is due to the presymptomatic hosts P, the second due to the asymptomatic hosts A, and the last one due to the symptomatic infectious hosts I. In each term, the numerator yields the rate of new infections βSAS∗,βSI S∗and this is then multiplied with the average duration of the stay in that infectious stage 1 σ2,1 MAR,1 M. Each term is multiplied with the fraction of hosts in that stage/state φfor asymptomatics and 1−φfor symptomatics. 2. Two-population model Using again the next generation matrix approach results in the following two matrices: F=                   0βyy SASyβyy SASyβyy SI Sy0 0 βyo SASyβyo SASyβyo SI Sy0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0βoy SASoβoy SASoβoy SI So0 0 βoo SASoβoo SASoβoo SI So0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                   (A4) V=                   σ10 0 0 0 0 0 0 0 0 −(1 −φy)σ1σ20 0 0 0 0 0 0 0 −φyσ10My AR 0 0 0 0 0 0 0 0−σ20My0 0 0 0 0 0 0 0 0 −γyMy(1 −ωy)χy+ωyψy0 0 0 0 0 0 0 0 0 0 σ10 0 0 0 0 0 0 0 0 −(1 −φo)σ1σ20 0 0 0 0 0 0 0 −φoσ10Mo AR 0 0 0 0 0 0 0 0 −σ20Mo0 0 0 0 0 0 0 0 0 −γoMo(1 −ωo)χo+ωoψo                   . (A5) Since Vis a block matrix it holds V= A0 0D!⇒V−1= A−10 0D−1!, 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 18 from which is readily follows that the spectral radius of FV −1is given by R0=Ryy 0+Roo 0 2+p(Ryy 0−Roo 0)2+ 4Ryo 0Roy 0 2,where Ryy 0= (1 −φy)βyy SASy σ2 +φyβyy SASy My AR + (1 −φy)βyy SI Sy My Roo 0= (1 −φo)βoo SASo σ2 +φoβoo SASo Mo AR + (1 −φo)βoo SI So Mo Ryo 0= (1 −φo)βyo SASy σ2 +φoβyo SASy Mo AR + (1 −φo)βyo SI Sy Mo Roy 0= (1 −φy)βoy SASo σ2 +φyβoy SASo My AR + (1 −φy)βoy SI So My. Ryy 0and Roo 0are the basic reproductive numbers in the young and old age groups, respectively, if they were completely isolated of each other. Ryo 0is the basic reproductive number if susceptible young hosts come into contact with only old infectious hosts, and Roy 0is the basic reproductive number if old hosts come into contact with only young infectious hosts. Acknowledgments This material is based upon work supported by the US National Science Foundation under Grants No. DMS-1815764 (ZR), PHY-1602994, and DMS-1809074 (PGK). PGK also acknowledges support from the Leverhulme Trust via a Visiting Fellowship and thanks the Mathematical Institute of the University of Oxford for its hospitality during this work. Disclaimer The views expressed in this manuscript are purely those of the authors and may not, under any circumstances, be regarded as an official position of the European Commission. [1] https://www.who.int/emergencies/diseases/novel-coronavirus-2019 [2] W.O. Kermack, A.G. McKendrick, Proc. Roy. Soc. A 115, 700 (1927). [3] H. Hethcote, SIAM Rev. 42, 599 (2000). [4] N.T.J. Bailey, The mathematical theory of infectious diseases and its applications, Griffin (London, 1975). [5] R.M. May, R.M Anderson, Infectious diseases of humans: dynamics and control, Oxford University Press (Oxford, 1991). [6] F. Brauer, C. Castillo-Ch´ avez, Mathematical Models in Population Biology and Epidemiology, Springer-Verlag (New York, 2001). [7] Z.-Q. Xia, J. Zhang, Y.-K. Xue, G.-Q. Sun, Z. Jin, PLoS One 10, e0144778 (2015). [8] X. Hao, et al., Nature, https://doi.org/10.1038/s41586-020-2554-8 (2020). [9] M. M. Arons et al., Presymptomatic SARS-CoV-2 infections and transmission in a skilled nursing facility, NEJM, DOI: 10.1056/NEJ- Moa2008457. [10] M. Peirlinck, K. Linka, F.S. Costabal, J. Bhattacharya, E. Bendavid, J.P.A. Ioannidis, E. Kuhl https://doi.org/10.1101/2020.05.23.20111419 [11] E. J. Emanuel et al., N Engl J Med 382, 2049-2055 (2020). [12] J. B. Dowd et al., PNAS 117, 9696-9698 (2020). [13] Z. Zheng et al., J. of Infection 81, e16-e25 (2020). [14] M.R. Nepomuceno et al., PNAS 117, 13881-13883 (2020). [15] For instance in the models reported in: https://www.cdc.gov/coronavirus/2019-ncov/covid-data/forecasting-us.html for modeling the US, only one (Northeastern) appears to be focusing on age-structured modeling [see, e.g., https://covid19.gleamproject.org/ and references therein]. [16] J. E. Ludvigsson, Acta Paediatrica 109 1088–1095 (2020). [17] K. Xu et al., Clinical Infectious Diseases 71, 799-806 (2020). [18] J. Mossong et al., PLoS Med 5, 374 (2008). [19] O. Torrealba-Rodriguez, R. A. Conde-Gutti´ errez, A. L. Hern´ andez-Javier, Chaos, Solitons and Fractals 138, 109946 (2020). [20] M. Anzarut et al., https://arxiv.org/abs/2007.09117 [21] U. Avila-Ponce de Le´ on, A. G. C. P´ erez, E. Avila-Vales, https://arxiv.org/abs/2004.08288 [22] J. P. Gutierrez and S. M. Bertozzi, https://www.medrxiv.org/content/10.1101/2020.05.27.20115204v2 [23] P.G. Kevrekidis, J. Cuevas-Maraver, Y. Drossinos, Z. Rapti, G.A. Kevrekidis, arXiv:2005.04527. [24] See, e.g., https://elifesciences.org/articles/57309 [25] N.I. Stilianakis and Y. Drossinos, J. R. Soc. Interface 7, 1355 (2010). [26] Y. Drossinos and N.I. Stilianakis, Aerosol Sci. Technol. 54, 639 (2020). [27] D.K. Milton, M. Patricia Fabian, B.J. Cowling, M.L. Grantham, J.J. McDevitt, PLoS Pathog. 9, e1003205 (2013). [28] https://www.gob.mx/salud/documentos/datos-abiertos [29] A.A. King, M. Domenech de Cell´ es, F.M.G. Magpantay, P. Rohani, Proc. Roy. Soc. B 282, 20150347 (2015). [30] M.C. Eisenberg, S.L. Robertson, J.H. Tien, J. Theor. Biol. 324, 84 (2013). 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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint 19 [31] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz, Journal of Mathematical Biology 28, 365–382 (1990). [32] A.S. Fokas, J. Cuevas-Maraver, P.G. Kevrekidis, https://doi.org/10.1101/2020.05.08.20095380. [33] J. M. Hyman and J. Li, SIAM J Appl. Math. 57 1082-1094 (1997). [34] D. Xiu and G. Em. Karniadakis, SIAM J. Sci. Comput. 24, 619 (2002). [35] https://github.com/dssg-pt/covid19pt-data [36] R. Li et al., Science 368, 489 (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 preprintthis version posted August 14, 2020. .https://doi.org/10.1101/2020.08.11.20172833doi: medRxiv preprint