On the relation of the COVID-19 reproduction number to the explosive timescales: the case of Italy
OOn the relation of the COVID-19 reproduction number to the explosivetimescales: the case of Italy
Dimitris G. Patsatzis ∗ School of Chemical Engineering, National Technical University of Athens, 15780 Athens, Greece
Abstract
A great issue of discussion of an infectious disease is its basic reproduction number R , which provides anestimation of the contagiousness of the disease. When R >
1, the disease spread will potentially lead toan outbreak, such that of the ongoing COVID-19 pandemics. During the evolution of an outbreak, variousnon-pharmaceutical interventions are employed, the impact of which is frequently assessed by the reductionthat they introduce to the effective reproduction number R t ; reduction below 1 is an indication of eventualdying out of the disease spread. Motivated by the fact that R essentially expresses the stability of thedisease-free equilibrium, in this work, R t was examined in the view of timescale analysis. It was shown thatduring the evolution of the COVID-19 outbreak in Italy, when various interventions were in place, R t hada clear relation with the explosive timescale characterizing the dynamics of the outbreak. In particular, it isshown that the existence of an explosive timescale during the progression of the epidemics implies R t > R t <
1. In addition, as this timescale converges/diverges with the immediatelyslowest one, R t approaches to/withdraws from its threshold value 1. These results suggest that timescaleanalysis can be utilized for the assessment of the impact of various interventions, since it reflects the insightprovided by the effective reproduction number, without being hindered by the selection of the populationmodel, nor the parameter estimation process followed for model calibration. Keywords:
COVID-19, reproduction number, timescale analysis, population dynamics
1. Introduction
As of March 11, 2020, the novel coronavirus disease (COVID-19) was declared a pandemic by WorldHealth Organization (WHO) [1]. By January 15, 2021, the COVID-19 pandemics has been spread to morethan 219 countries and territories, reporting more than 93 million infected cases and 2 million deaths [2].A great issue of discussion of the COVID-19 pandemics is its basic reproduction number, estimations forwhich were provided in numerous early studies; see Refs within [3]. The basic reproduction number, R , is theaverage number of secondary infections produced by an infectious individual in a population where everyoneis considered susceptible [4, 5]. Being dependent on human behavior and the biological characteristics ofthe pathogen, R provides an estimation of the contagiousness of the infectious disease [4] and serves as athreshold parameter; when R > R < R is constant in time, it cannot monitor the effect of the undertaken measures;instead, the time-varying effective reproduction number R t is utilized, that estimates the secondary infectionsproduced by an infectious individual during the course of an outbreak, thus, in a population where noteveryone is considered susceptible. As a result, during the evolution of the epidemics, the undertaken controlmeasures affect R t , since they influence (i) the duration of contagiousness, (ii) the likelihood of infection percontact and (iii) the contact rate of the infection [4, 6]. The impact of various interventions (case isolation, ∗ Corresponding author
Email address: [email protected] (Dimitris G. Patsatzis)
Preprint submitted to arXiv.org January 18, 2021 a r X i v : . [ q - b i o . P E ] J a n ontact tracing, travel restrictions, etc.) in R t has been assessed in a number of studies to provide guidelinesin decision-making policies [7–12].The use of R as a threshold parameter is related to the stability of the disease-free equilibrium (DFE)of the epidemiological model under consideration [4, 5, 13], which is locally assessed by the existence ofpositive eigenvalues. During the evolution of the system, the local dynamics is characterized by timescalesof dissipative/explosive nature - associated with positive/negative eigenvalues - the action of which tends todrive the system towards to/away from equilibrium [14, 15]. Timescale analysis has been frequently employedto address the dynamical properties of systems arising from reactive flows [16, 17], systems biology [18, 19],pharmacokinetics [20], etc, but, to my knowledge, it hasn’t been widely applied to population dynamics.Motivated by the fact that R mathematically expresses the stability of the DFE; i.e., the existence ofpositive eigenvalues at day zero, here the relation of R t to the explosive timescales (positive eigenvalues) dur-ing the course of COVID-19 outbreak in Italy was investigated. It is shown that the existence of an explosivetimescale implies R t >
1, while its absence implies R t <
1. In addition, as this timescale converges/divergeswith the immediately slowest one, R t was shown to approaches to/withdraws from its threshold value 1. Fi-nally, by performing the analysis in 4 different population dynamics models, it is demonstrated that timescaleanalysis is a robust methodology to monitor the progression of the epidemics, since it directly reflects thevariations in R t , without being hindered by the complexity of the selected model.
2. Materials and Methods
Compartmental modeling is widely used for the analysis of various infectious diseases [21–24], among whichCOVID-19 pandemics [8, 11, 25, 26]. Four population dynamics compartmental models in the framework ofthe SIR model [27] were analyzed here. The effective reproduction number R t was calculated on the basisof these models and conclusions were drawn on its relation to the timescales characterizing the dynamics ofeach model.The four compartmental models are presented in Section 2.1, followed by the parameter estimation processconsidered for their calibration against the data of Italy in Section 2.2. The methodology to calculate theeffective reproduction number R t and the timescales τ i on the basis of each model is presented in Sections 2.3and 2.4, respectively. The SIR model formulates the transmission of an infectious disease among three population groups,namely the susceptible , infected and recovered individuals [27]. In this framework, four population dynamicsmodels were considered, the SIRD, SEIRD, SEInsRD and SIDARTHE models, the governing equations ofwhich can be written in the ODE form: ddt y = g ( y ) (1)where y is the N-dim. column state vector , which includes the fraction of each population group over thetotal population and g ( y ) is the N-dim. column vector field , which incorporates the transition rates fromone population group to another.The simplest compartmental model to capture COVID-19 pandemics is the SIRD model, which essentiallyis the SIR model with the addition of a compartment accounting for the dead individuals. Denoting S , I , R and D the fraction of susceptible, infected, recovered and dead individuals respectively, over the totalpopulation N , the SIRD model is written in the form of Eq. (1) as: ddt SIRD = − βSIβSI − ( γ + µ ) IγIµI (2)where β is the transmission ratio, γ the recovery ratio, which also expresses the inverse of the infection periodof the disease, and µ the fatality ratio.A more realistic assumption for COVID-19 infection is the existence of an incubation (latency) period,during which an individual is infected but yet not infectious [28, 29]. Such an assumption can be incorporated2n the SIRD model with the addition of a compartment accounting for exposed individuals. Denoting E theirfraction over the total population, the resulting SEIRD model is written in the form of Eq. (1) as: ddt SEIRD = − βSIβSI − σEσE − ( γ + µ ) IγIµI (3)where σ is the transition ratio from exposed to infected individuals, expressing the inverse of the incubationperiod of the disease.In addition, it has been shown that the COVID-19 infected individuals have symptoms of different severity,varying from mild to severe [30]. Since the severely infected individuals are in need of immediate health care,a more biologically realistic assumption for COVID-19 infection is the distinction between normally infected and severely infected individuals. Such an assumption can be incorporated in the SEIRD model, by dividingthe infected compartment in two sub-compartments. Denoting IN and IS the fraction of normally andseverely infected individuals over the total population, the resulting SEInsRD model in the form of Eq. (1)reads: ddt SEINISRD = − β N S.IN − β S S.IS − µ T P Sβ N S.IN + β S S.IS − σE − µ T P E (1 − ss ) σE − γIN − µ N INssσE − γIS − µ S ISγ ( IN + IS ) − µ T P Rµ N IN + µ S IS (4)where the subscripts N and S indicate the normally and severely infected transmission β and fatality µ ratios, ss denotes the fraction of severely over normally infected individuals and µ T P is the physiologicaldeath ratio.Finally, a more detailed compartmental model was considered, accounting for susceptible ( S ), asymp-tomatic detected and undetected infected ( I and D ), symptomatic detected and undetected infected ( A and R ), severely symptomatic ( T ), healed ( H ) and extinct ( E ) individuals, namely the SIDARTHE model[25]. Here, the SIDARTHE model was considered for validation purposes and thus, only a brief descriptionof the model is provided in Appendix A; details can be found in [25]. Note that the SIDARTHE model isalso written in the form of Eq. (1); see Eq. (A.13) and Eqs. (1-8) in Methods section in [25].
In this study only the SEIRD and SEInsRD models in Eqs. (3, 4) respectively, were calibrated, since (i)the relation of R t with the timescales can be reached analytically on the basis of the SIRD model in Eq. (2)and (ii) the parameter values of the SIDARTHE model are provided in [25].The SEIRD and SEInsRD models in Eqs. (3, 4) were calibrated to the daily reported data of infected,recovered and dead individuals in Italy, as reported by John Hopkins database [31]. The parameter estimationprocess was performed in a weekly basis accounting for the data from February 26 (week 0) to September 30(week 30). February 26 was selected as starting day in order to minimize early data distortion, since morethan 400 infected individuals were reported at that date. Initially, given the reported fraction of infected,recovered and dead population groups at week 0 - and the susceptible one, through conservation of thetotal population - the fraction of the exposed, normally infected and severely infected population groups wasestimated. In the following, a parameter estimation process was performed in a weekly basis, given these 3reported data sets, through a genetic algorithm provided by the open-source COPASI software [32, 33]. Theinitial conditions at day 0 of each week were the predicted values at day 7 of the previous week, in orderto preserve continuity in the solution. The resulting parameter sets are depicted for SEIRD and SEInsRDmodels in Fig. B.1 of Appendix B.A schematic representation of the SEIRD and SEInsRD models is provided in the left panels of Figs. 1and 2, respectively. The profiles of infected, recovered and dead individuals, resulted from the aforementionedparameter estimation process, are in very good agreement with the reported data, as shown in the middlepanels of Figs. 1 and 2 for the SEIRD and SEInsRD models, respectively. Note that in the case of SEInsRD3 E DR I weeks × × InfectedRecoveredDead weeks SEIRD
Figure 1: The SEIRD model (left), its fitting against the reported data for Italy in circles (middle) and the profiles of all thepopulation groups (right). The parameter estimation process was performed in a weekly basis from the 26th of February (week0) to the 30th of September (week 30), accounting for the reported data of infected, recovered and dead individuals.
SEDR ISIN weeks × × InfectedRecoveredDead weeks SEINISRDI
Figure 2: The SEInsRD model (left), its fitting against the reported data for Italy in circles (middle) and the profiles of all thepopulation groups (right). The parameter estimation process was performed in a weekly basis from the 26th of February (week0) to the 30th of September (week 30), accounting for the reported data of infected, recovered and dead individuals. model, the sum of the normally and severely infected individuals I = IN + IS was fitted against the reporteddata set of the infected individuals. The very good agreement of the model parameters to the reported datacan be also demonstrated by the R values of both fittings shown in Table B.1 of Appendix B, combinedwith the respective p-values, which in all cases are p (cid:28) .
05. Finally, the profiles of all the population groupsare displayed in the right panels of Figs. 1 and 2, in which great agreement of the population profiles betweenthe two models is reported.Finally, the SIDARTHE model parameters were directly adopted by [25], following a slightly differentapproach than the one followed here for SEIRD and SEInsRD models. First, due to availability of data,here the SEIRD and SEInsRD models were calibrated for Italy from February 26 to September 30, whilethe SIDARTHE model was calibrated for Italy from February 20 to April 5. Second, here the SEIRD andSEInsRD models were calibrated in a constant, 7-days long, time frame, while the SIDARTHE model wascalibrated in varying time frames, depending on the interventions undertaken in Italy [25]. Last but not least,the reported data of infected, recovered and dead individuals were considered in our analysis, while in [25]only the ones of infected and recovered (not fitted to the healed compartment of the model) individuals.
The basic reproduction number, R , is a constant biological parameter that provides an estimation of thecontagiousness of the infectious disease. It also serves as a threshold parameter; when R >
1, one infectedindividual can trigger an outbreak, while when R <
1, the infection will not spread in the population [4, 5].When various non-pharmaceutical interventions (NPI) are in place, the effective reproduction number R t is utilized, instead of R , to monitor the reproduction number during the evolution of the outbreak. R t R t for COVID-19 pandemics in Italy using the Next Generation Matrix (NGM) approach[13, 34, 35], which yields in the following expressions for the SIRD, SEIRD, SEInsRD and SIDARTHEmodels: R SIRDt = βγ + µ = R SEIRDt R SEInsRDt = σσ + µ T P (cid:18) (1 − ss ) β N γ + µ N + ssβ S γ + µ S (cid:19) (5) R SIDART HEt = αr + β(cid:15)r r + γζr r + δθζr r r + δ(cid:15)ηr r r where α, β, γ, δ, (cid:15), ζ, η, θ are model parameters of the SIDARTHE model and r = (cid:15) + λ + ζ , r = η + ρ , r = κ + µ + θ annd r = ν + ξ . Note that the expression of R t for SIDARTHE model estimated here via theNGM approach is the same with the one derived in [25].A brief discussion on NGM approach is provided in Appendix A, along with details on the calculationof R t on the basis of the four population dynamics models. Given a system of ODEs in the matrix form of Eq. (1), the timescales are calculated as the inversemodulus of the eigenvalues of the N × N Jacobian matrix J ( y ) = ∇ y ( g ( y )) [14, 15]. The timescales are ofdissipative/explosive nature, i.e., the components of the system that generate them tend to drive the systemtowards to/away from equilibrium, when the respective eigenvalue has negative/positive real part.When a complex mathematical model in the form of Eq. (1) is encountered, it is usually impossible tocalculate analytic expressions for its eigenvalues and thus its timescales. This is the case of the SEIRD,SEInsRD and SIDARTHE models, for which the timescales were calculated numerically. However, in the caseof the SIRD model, the non-zero eigenvalues can be calculated analytically as: λ , = 12 (cid:16) X ± (cid:112) X − Y (cid:17) X = − γ − βI − µ + βS Y = βI ( γ + µ ) (6)Therefore, the related timescales are of explosive nature (either real or complex λ , ) if and only if: X > ⇒ β ( S − I ) > γ + µ ⇒ β ( S − I ) γ + µ > R t .
3. Results
The impact of the undertaken NPIs in COVID-19 pandemics is assessed by the effect that they introducein the reproduction number [7–12]. Here, we show that the insights provided by the utilization of the effectivereproduction number R t during the progression of the COVID-19 pandemics can be deduced by timescaleanalysis. In particular, it is shown that:i) the existence of an explosive timescale during the progression of COVID-19 epidemics implies R t > R t <
1, andii) the tendency of this timescale to converge/diverge with the immediately slowest one, implies that R t tends to approach to/withdraw from its threshold value 1.These results are reached on the basis of the four population dynamics models discussed in Section 2.1, forthe case of Italy. 5 .1. The explosive timescales in relation to the reproduction number The first indication on the relation of the explosive timescales to the reproduction number is providedby the analysis of the SIRD model in Eq. (2), that is the simplest model to describe the progression ofCOVID-19 epidemics. In contrast to more complicated models, the timescales of the SIRD model can becalculated analytically. According to Section 2.4, the evolution of the SIRD model is characterized by theaction of two timescales τ , = 1 / | λ , | ; the expressions of λ , were derived in Eq. (6). Both τ , are ofexplosive/dissipative nature when the condition in Eq. (7) holds/is violated. Given the expression of R t forSIRD model in Eq. (5) and that S − I < S (0) = 1, Eq. (7) yields: Re ( λ , ) > ⇔ X > ⇔ β ( S − I ) γ + µ > ⇒ βS (0) γ + µ > ⇒ R t > R t >
1, while their absence implies R t <
1. Note that this outcome, holds true not only for COVID-19 pandemics, but also for any infectiousdisease, since it was derived by analytical means on the basis of the SIRD model.Next, the relation of the explosive timescales with R t was examined using reported data for COVID-19pandemics in Italy. The SEIRD model in Eq. (3) was adopted and fitted against the reported data sets ofinfected, recovered and dead individuals in Italy from February 26 to September 30. In order to account forthe NPIs undertaken, the SEIRD model was calibrated in a weekly basis following the parameter estimationprocess described in detail in Section 2.2. The resulting solution is in great agreement with the reported data,as shown in Fig. 1.The timescales and R t , estimated on the basis of the SEIRD model in Eq. (5), are displayed in Fig. 3from week 0 (starting in Feb. 26) to week 30 (ending in Sep. 30). As shown in the left panel of Fig. 3, theevolution of the SEIRD model is characterized by three timescales τ , , , the fastest of which, τ , is alwaysdissipative in nature, while τ , are either dissipative or explosive. In particular, during weeks 0-6 and 20-26, τ , are explosive, as indicated by the shaded background in Fig. 3. The values of R t are depicted in the rightpanel of Fig. 3, in which the red dashed horizontal line indicates the threshold value R t = 1. As indicated bythe shaded background, the time periods when the explosive nature of timescales τ , is reported coincideswith the ones that R t > τ , are of dissipative nature, R t < τ , , and vice-versa, is immediate, since model calibration is performed in a weekly basis.Comparison of the explosive timescales and R t in Fig. 3 reveals the following trend: as the gap between τ and τ decreases/increases, R t approaches to/withdraws from its threshold value unity. This is particularlyclear during the first wave of COVID-19 pandemics in Italy (weeks 0-12). During weeks 0-6, where τ and τ are explosive, their gap tends to decrease, so that R t decreases, approaching close to unity values. Atweek 7, τ and τ become dissipative and R t attains values below 1. From this point on and up to week 12,the gap of τ and τ increases, so that R t continues to decrease, this time withdrawing from its threshold1. This behaviour is additionally supported by the fact that during weeks 4-8, 17, 19, 21 and 30, in which R t weeks i DissipativeExplosive weeks R t Figure 3: The timescales (left) and the effective reproduction number R t (right) estimated on the basis of the solution of theSEIRD model shown in Fig. 1. τ and τ is small; to the point where τ = τ in week 19, in which R t = 0 . In order to demonstrate the robustness of the relation of the explosive timescales to the reproductionnumber, a more complicated population dynamics model was considered, the SEInsRD model. The SEInsRDmodel in Eq. (4) was adopted and calibrated to the same reported data sets with SEIRD model in Section 3.1,corresponding to infected, recovered and dead individuals in Italy from February 26 to September 30. Similarlyto SEIRD model, the SEInsRD model calibration was performed in a weekly basis following the processdescribed in Section 2.2 and the resulting solution is in great agreement with the reported data, as shown inFig. 2. weeks i DissipativeExplosive weeks R t Figure 4: The timescales (left) and the reproduction number R t (right) calculated on the basis of the solution of the SEInsRDmodel shown in Fig. 2. The timescales and R t , estimated on the basis of the SEInsRD model in Eq. (5), are displayed in Fig. 4from week 0 (starting in Feb. 26) to week 30 (ending in Sep. 30). As shown in the left panel of Fig. 4, theevolution of SEInsRD model is characterized by 5 timescales: three of which are always of dissipative natureand the remaining ones are either explosive or dissipative; denoting τ exp,f the fast explosive timescale and τ exp,s the slow one. In particular, during weeks 0-6 and 20-26, τ exp,f and τ exp,s are explosive, as indicatedby the shaded background in Fig. 4. The right panel of Fig. 4 displays the values of R t in comparison tothe threshold value R t = 1 indicated by the red dashed horizontal line. Similarly to the SEIRD model, it isshown by the shaded background that R t > τ exp,f and τ exp,s are explosive (weeks 0-6 and 20-26),while R t < τ exp,f and τ exp,s is again reflected in R t approaching to/withdrawingfrom its threshold value 1. In particular, it is shown that the closer the values of R t to 1, (weeks 4-8, 17,19, 21 and 30), the smaller the gap between τ exp,f and τ exp,s ; to the point where R t = 0 .
97 in week 30, inwhich τ exp,f ≈ τ exp,s . In summary, the qualitative results on the relation of the explosive timescales to R t are maintained on the basis of the SEInsRD model. In order to validate the qualitative results, reached on the basis of the SEIRD and SEInsRD models,regarding to the relation of the explosive timescales to R t , a more complicated SIDARTHE model wasconsidered [25], as briefly discussed in Section 2.1.The profiles of the population groups accounted for in the SIDARTHE model were reproduced, adoptingthe model parameters in [25]. On the basis of the SIDARTHE solution, the timescales were calculated and R t was estimated according to the expression in Eq. (5). The resulting values are displayed in Fig. 5 startingfrom day -6 (Feb 20) and ending in day 40 (Apr 5); day 0 was chosen to be Feb 26 for comparison withFigs. 3 and 4. As shown in the left panel of Fig. 5, the evolution of SIDARTHE model is characterized by sixtimescales, four among which are always dissipative in nature, while the remaining two are either dissipativeor explosive; denoted as τ exp,f and τ exp,s . In particular, τ exp,f and τ exp,s are explosive from day -6 to day 22,7
14 28 days i DissipativeExplosive days R t Figure 5: The timescales and reproduction number R t calculated on the basis of SIDARTHE model, that was calibrated forItaly data in [25]. as indicated by the shaded background in Fig. 5. The values of R t are depicted in the right panel of Fig. 5,in which the red dashed horizontal line indicates the threshold value R t = 1. As indicated by the shadedbackground, the explosive nature of timescales τ exp,f and τ exp,s implies R t > R t < τ exp,f and τ exp,s becomes smaller, R t approaches its threshold value unity, to the point when R t = 0 .
986 during days 23-33, in which τ exp,f = τ exp,s .It should be noted here that, the evolution of timescales τ exp,f and τ exp,s and the reproduction number R t , calculated on the basis of SIDARTHE model, is different in comparison to those calculated on the basisof SEIRD and SEInsRD models. In particular, on the basis of the SIDARTHE model, the timescales areexplosive in nature and R t > R t , deduced on the basis of SEIRD andSEInsRD models, is validated by the analysis with the SIDARTHE model.
4. Conclusions
The progression of an infectious disease spread like COVID-19 pandemics is frequently examined by pop-ulation dynamics models [8, 11, 21–26]. Their evolution as dynamical systems is characterized by timescalesthat are either of dissipative or explosive nature; i.e., their action tends to drive the system either towards toor away from equilibirum [14, 15]. The basic reproduction number R as a threshold parameter provides suchan intuition, in the view that when R < R > R ≈ − R t , [7–12]; ideally making R t <
1, which indicatesthat the disease spread will eventually die out. In this work, the relation of the effective reproduction num-ber R t with the timescales characterizing the evolution of the epidemic spread was examined in the case ofCOVID-19 pandemics in Italy from February 26 to September 30.In particular, it was demonstrated analytically on the basis of the SIRD model and numerically on thebasis of the SEIRD model in Section 3.1, that when two of the timescales characterizing the evolution ofthe epidemic spread are of explosive nature, the effective reproduction number is above its threshold value;i.e., R t >
1. On the contrary, when all the timescales are of dissipative nature it is implied that R t <
1. Inaddition, the following trending behaviour was revealed: as the gap between the two explosive timescales in-creases/decreases, R t approaches to/withdraws from its threshold value 1, as shown in Fig. 3. These outcomessuggest that the insights provided by the utilization of R t as a threshold parameter can be also obtained bytimescale analysis. 8his work additionally suggests that timescale analysis is a robust methodology to assess the progressionof the epidemic spread, since it is not hindered by the complexity of the selected model, nor the calibrationprocess followed to fit the model against the reported data. Following the same model calibration procedure tothe SEInsRD model, resulted in timescales that are almost equal to the ones of the SEIRD model; see Figs. 3and 4. Such a result indicates that the relation of the explosive timescales to R t is not affected by modelselection, as discussed in Section 3.2. In addition, this relation is not affected by the parameter estimationprocess either, as demonstrated through the analysis of the SIDARTHE model, the calibration of which in[25] had significant differences with the one followed here for SEIRD and SEInsRD models; see Section 2.2.In conclusion, timescale analysis is a rigorous mathematical methodology to assess the progression of anepidemic spread, since it can effectively provide the insight obtained by the reproduction number. Timescaleanalysis is not hindered by model selection in contrast to the reproduction number that is highly dependableon the structure of the selected model [4]. In addition, the expression of the reproduction number becomesmore complex as the detail of the model increases, as shown in Eq. (5); compare for example R t fo SEIRDand SIDARTHE models. In contrast, timescale analysis can be performed in an algorithmic fashion, utilizingthe diagnostic tools of Computational Singular Perturbation [14, 36] that have been effectively employed toaddress the dynamical properties of systems arising from a wide variety of fields [16–20]. More importantly,the use of timescale analysis for the assessment of various NPIs is promising, since it can determine via itsalgorithmic tools the factors that play the most significant role on the control of ongoing COVID-19 outbreak.
5. Acknowledgements
This publication is based upon work supported by the Khalifa University of Science and Technology,under Award No. CPRA-2020-Goussis.
References [1] World Health Organization. WHO Director-General’s opening remarks at the mediabriefing on COVID-19, March 11, 2020, , accessed: 2020-03-11.[2] Worldometer. COVID-19 Coronavirus Pandemic, ,accessed: 2021-1-15.[3] Y. Liu, A. A. Gayle, A. Wilder-Smith, J. Rockl¨ov, The reproductive number of covid-19 is highercompared to sars coronavirus, Journal of travel medicine.[4] P. L. Delamater, E. J. Street, T. F. Leslie, Y. T. Yang, K. H. Jacobsen, Complexity of the basicreproduction number (r0), Emerging infectious diseases 25 (1) (2019) 1.[5] O. Diekmann, J. A. P. Heesterbeek, J. A. Metz, On the definition and the computation of the basic repro-duction ratio r 0 in models for infectious diseases in heterogeneous populations, Journal of mathematicalbiology 28 (4) (1990) 365–382.[6] G. Viceconte, N. Petrosillo, Covid-19 r0: Magic number or conundrum?, Infectious disease reports 12 (1).[7] J. Hellewell, S. Abbott, A. Gimma, N. I. Bosse, C. I. Jarvis, T. W. Russell, J. D. Munday, A. J.Kucharski, W. J. Edmunds, F. Sun, et al., Feasibility of controlling covid-19 outbreaks by isolation ofcases and contacts, The Lancet Global Health.[8] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, F. Sun, M. Jit,J. D. Munday, et al., Early dynamics of transmission and control of covid-19: a mathematical modellingstudy, The lancet infectious diseases.[9] A. Pan, L. Liu, C. Wang, H. Guo, X. Hao, Q. Wang, J. Huang, N. He, H. Yu, X. Lin, et al., Associationof public health interventions with the epidemiology of the covid-19 outbreak in wuhan, china, Jama323 (19) (2020) 1915–1923. 910] B. J. Cowling, S. T. Ali, T. W. Ng, T. K. Tsang, J. C. Li, M. W. Fong, Q. Liao, M. Y. Kwan, S. L. Lee,S. S. Chiu, et al., Impact assessment of non-pharmaceutical interventions against coronavirus disease2019 and influenza in hong kong: an observational study, The Lancet Public Health.[11] B. Tang, X. Wang, Q. Li, N. L. Bragazzi, S. Tang, Y. Xiao, J. Wu, Estimation of the transmission riskof the 2019-ncov and its implication for public health interventions, Journal of clinical medicine 9 (2)(2020) 462.[12] B. Tang, N. L. Bragazzi, Q. Li, S. Tang, Y. Xiao, J. Wu, An updated estimation of the risk of transmissionof the novel coronavirus (2019-ncov), Infectious disease modelling 5 (2020) 248–255.[13] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria forcompartmental models of disease transmission, Mathematical biosciences 180 (1-2) (2002) 29–48.[14] S. Lam, D. Goussis, Understanding complex chemical kinetics with computational singular perturbation,in: Symposium (International) on Combustion, Vol. 22, Elsevier, 1989, pp. 931–941.[15] U. Maas, S. B. Pope, Simplifying chemical kinetics: intrinsic low-dimensional manifolds in compositionspace, Combustion and flame 88 (3-4) (1992) 239–264.[16] D. M. Manias, E. A. Tingas, C. E. Frouzakis, K. Boulouchos, D. A. Goussis, The mechanism by whichch2o and h2o2 additives affect the autoignition of ch4/air mixtures, Combustion and Flame 164 (2016)111–125.[17] E. A. Tingas, D. C. Kyritsis, D. A. Goussis, Autoignition dynamics of dme/air and etoh/air homogeneousmixtures, Combustion and Flame 162 (9) (2015) 3263–3276.[18] P. D. Kourdis, R. Steuer, D. A. Goussis, Physical understanding of complex multiscale biochemicalmodels via algorithmic simplification: Glycolysis in saccharomyces cerevisiae, Physica D: NonlinearPhenomena 239 (18) (2010) 1798–1817.[19] D. G. Patsatzis, D. A. Goussis, A new michaelis-menten equation valid everywhere multi-scale dynamicsprevails, Mathematical biosciences 315 (2019) 108220.[20] D. G. Patsatzis, D. T. Maris, D. A. Goussis, Asymptotic analysis of a target-mediated drug dispositionmodel: algorithmic and traditional approaches, Bulletin of mathematical biology 78 (6) (2016) 1121–1161.[21] S. Gao, L. Chen, Z. Teng, Pulse vaccination of an seir epidemic model with time delay, NonlinearAnalysis: Real World Applications 9 (2) (2008) 599–607.[22] L. Canini, F. Carrat, Population modeling of influenza a/h1n1 virus kinetics and symptom dynamics,Journal of virology 85 (6) (2011) 2764–2770.[23] L. Esteva, C. Vargas, Coexistence of different serotypes of dengue virus, Journal of mathematical biology46 (1) (2003) 31–47.[24] A. Stegeman, A. R. Elbers, J. Smak, M. C. de Jong, Quantification of the transmission of classical swinefever virus between herds during the 1997–1998 epidemic in the netherlands, Preventive veterinarymedicine 42 (3-4) (1999) 219–234.[25] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Matteo, M. Colaneri, Modellingthe covid-19 epidemic and implementation of population-wide interventions in italy, Nature Medicine(2020) 1–6.[26] L. Russo, C. Anastassopoulou, A. Tsakris, G. N. Bifulco, E. F. Campana, G. Toraldo, C. Siettos, Tracingday-zero and forecasting the fade out of the covid-19 outbreak in lombardy, italy: A compartmentalmodelling and numerical optimization approach., medRxiv.1027] W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedingsof the royal society of london. Series A, Containing papers of a mathematical and physical character115 (772) (1927) 700–721.[28] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich,J. Lessler, The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmedcases: estimation and application, Annals of internal medicine 172 (9) (2020) 577–582.[29] Q. Li, X. Guan, P. Wu, X. Wang, L. Zhou, Y. Tong, R. Ren, K. S. Leung, E. H. Lau, J. Y. Wong, et al.,Early transmission dynamics in wuhan, china, of novel coronavirus–infected pneumonia, New EnglandJournal of Medicine.[30] V. Surveillances, The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases(covid-19)—china, 2020, China CDC Weekly 2 (8) (2020) 113–122.[31] COVID-19 Data Repositoy by the Center for Systems Science and Engineering (CSSE) at Johns HopkinsUniversity, https://github.com/CSSEGISandData/COVID-19 , accessed: 2020-06-26.[32] T. P. Runarsson, X. Yao, Stochastic ranking for constrained evolutionary optimization, IEEE Transac-tions on evolutionary computation 4 (3) (2000) 284–294.[33] S. Hoops, S. Sahle, R. Gauges, C. Lee, J. Pahle, N. Simus, M. Singhal, L. Xu, P. Mendes, U. Kummer,Copasi—a complex pathway simulator, Bioinformatics 22 (24) (2006) 3067–3074.[34] J. M. Heffernan, R. J. Smith, L. M. Wahl, Perspectives on the basic reproductive ratio, Journal of theRoyal Society Interface 2 (4) (2005) 281–293.[35] P. van den Driessche, Reproduction numbers of infectious disease models, Infectious Disease Modelling2 (3) (2017) 288–303.[36] S. Lam, D. Goussis, The csp method for simplifying kinetics, International journal of chemical kinetics26 (4) (1994) 461–486.
Appendix A. Derivation of the effective reproduction number
The Next Generation Matrix (NGM) approach is utilized for the calculation of the basic reproductionnumber R [13, 34, 35]. Given a system of ODEs in the form of Eq. (1), let y j be the j = 1 , . . . , m infectedpopulation groups among all the y i populations groups of the i = 1 , . . . , n compartments in y . In turn, let F i ( y ) be the rate of appearance of new infections in the i -th compartment and V i ( y ) = V − i ( y ) − V + i ( y ) thetransition rates out of ( V − ) and into ( V + ) the i -th compartment. By definition, it is implied that: dy i dt = F i ( y ) − V i ( y ) = F i ( y ) + V + i ( y ) − V − i ( y ) (A.1)Let the matrices F and V be: F = (cid:20) ∂F i ( y ∗ ) ∂y j (cid:21) and V = (cid:20) ∂V i ( y ∗ ) ∂y j (cid:21) (A.2)where y ∗ is the disease-free equilibrium and i, j = 1 , . . . , m . According to the NGM approach, the basicreproduction number R is the spectral radius (largest eigenvalue) of the matrix F · V − ; i.e., R = ρ ( F · V − )[13, 34, 35]. However, since the model parameters vary in time (different parameter values in each week), theNGM approach utilization results in the calculation of the effective reproduction number R t . In the following,the analytical expressions of R t for SIRD, SEIRD, SEInsRD and SIDARTHE models in Eq. (5) are derived.The SIRD mathematical model in Eq. (2) can be written in the form of Eq. (A.1) as: ddt SIRD = βSI + γIµI − βSI ( γ + µ ) I = F i ( y ) + V + i ( y ) − V − i ( y ) (A.3)11he disease-free equilibrium is y ∗ = ( S (0) , , , F = βS (0) and V = γ + µ (A.4)Given that S (0) = 1 as fraction of the total population, the effective reproduction number for the SIRDmodel is: R t = ρ ( F · V − ) = βγ + µ (A.5)The SEIRD mathematical model in Eq. (3) can be written in the form of Eq. (A.1) as: ddt SEIRD = βSI + σEγIµI − βSIσEγI + µI = F i ( y ) + V + i ( y ) − V − i ( y ) (A.6)The disease-free equilibrium is y ∗ = ( S (0) , , , , F = (cid:20) βS (0)0 0 (cid:21) and V = (cid:20) σ − σ γ + µ (cid:21) (A.7)Given that S (0) = 1, the effective reproduction number for the SEIRD model is: R t = ρ ( F · V − ) = βγ + µ (A.8)Note that the R t of SEIRD model is the same to that of SIRD model in Eq. (A.5).The SEInsRD mathematical model in Eq. (4) can be written in the form of Eq. (A.1) as: ddt SEINISRD = β N S.IN + β S S.IS + − ss ) σEssσEγ ( IN + IS ) µ N IN + µ S IS − β N S.IN + β S S.IS + µ T P
SσE + µ T P
EγIN + µ N INγIS + µ S ISµ
T P R = F i ( y )+ V + i ( y ) − V − i ( y )(A.9)The disease-free equilibrium is y ∗ = ( S (0) , , , , , F = β N S (0) β S S (0)0 0 00 0 0 and V = σ + µ T P − (1 − ss ) σ γ + µ N − ssσ γ + µ S (A.10)Given that S (0) = 1, the effective reproduction number for the SEInsRD model is: R t = ρ ( F · V − ) = σσ + µ (cid:18) (1 − ss ) β N γ + µ N + ssβ S γ + µ S (cid:19) (A.11)Note that when considering the µ T P (cid:28) σ limit, R t of SEInsRD model in Eq. (A.11) is simplified to: R t µ TP (cid:28) σ = (cid:18) (1 − ss ) β N γ + µ N + ssβ S γ + µ S (cid:19) (A.12)which is similar to that of SIRD and SEIRD models in Eqs. (A.5, A.8) when setting ss = 0; i.e., whenneglecting the severely infected individuals from the model.12inally, the SIDARTHE mathematical model in [25] can be written in the form of Eq. (A.1) as: ddt SIDARTHE = − S ( αI − βD − γA − δR ) S ( αI + βD + γA + δR ) − ( (cid:15) + ζ + λ ) I(cid:15)I − ( η + ρ ) DζI − ( θ + µ + κ ) AηD + θA − ( ν + ξ ) RµA + νR − ( σ + τ ) TλI + ρD + κA + ξR + σTτ T = S ( αI + βD + γA + δR )000000 ++ (cid:15)IζIηD + θAµA + νRλI + ρD + κA + ξR + σTτ T − S ( αI − βD − γA − δR )( (cid:15) + ζ + λ ) I ( η + ρ ) D ( θ + µ + κ ) A ( ν + ξ ) R ( σ + τ ) T = F i ( y ) + V + i ( y ) − V − i ( y ) (A.13)where the parameter notation is explained in detail in [25]. The disease-free equilibrium is y ∗ = ( S (0) , , , , , , , F = αS (0) βS (0) γS (0) δS (0) 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0 and V = (cid:15) + λ + ζ − (cid:15) η + ρ − ζ κ + µ + θ − η − θ ν + ξ
00 0 µ ν σ + τ (A.14)Given that S (0) = 1, the effective reproduction number for the SIDARTHE model is: R t = ρ ( F · V − ) = αr + β(cid:15)r r + γζr r + δθζr r r + δ(cid:15)ηr r r (A.15)where r = (cid:15) + λ + ζ , r = η + ρ , r = κ + µ + θ and r = ν + ξ . Note that the expression in Eq. (A.15)derived here in the context of NGM approach is the same with the one in Eq. (18) derived in [25] using adifferent approach. Appendix B. The SEIRD and SEInsRD model parameters
The parameter estimation process described in Section 2.2, that was followed to fit the reported datasets of infected, recovered and dead individuals of Italy from February 26 to September 30 in a weekly basis,resulted in the model parameters shown in Fig. B.1. The left panel shows the distribution of the SEIRDmodel parameters β , σ , γ and µ and the right panel shows the ones of SEInsRD model β N , β S , σ , γ , µ N , µ S and ss . The values of parameter µ T P of the SEInsRD model are not shown, since they are smaller than10 − .Figure B.1 indicates that the parameters expressing the transition from a population group to anotherattain similar values in both models: transmission rate ( β and β N , β S ), incubation period (1 /σ ), recoveryrate ( γ ) and fatality rate ( µ and µ N , µ S ) constants.As shown in Fig. B.1, the following trends in the parameter values are indicated: • the transmission rate constant β attains high/low values in the periods where explosive timescale arepresent/absent. The values of β tend to decrease during the transition from an explosive to a dissipativeregion and vice-versa. • the rate constant σ (inverse of incubation period) tend to increases during the explosive regions.13 weeks -4 -2 βσγ µ weeks -4 -2 NS σγ µ N µ S ss Figure B.1: The parameters estimated for the SEIRD (left) and SEInsRD (right) models. The shaded regions indicate the weeksfor which R t > • the recovery rates γ are almost constant • the fatality rates µ tend to decrease, despite the explosive/dissipative region transition. They tend toincrease only in the last few weeks. • the normally to severely infected ratio ss is almost constant.population group SEIRD SEInsRDinfected, I . . R . . D . . Table B.1: R values of the solution acquired on the basis of the SEIRD and SEInsRD models with the parameter distributionshown in Fig. B.1, with reference to the reported data for infected, recovered and dead individuals in Italy.values of the solution acquired on the basis of the SEIRD and SEInsRD models with the parameter distributionshown in Fig. B.1, with reference to the reported data for infected, recovered and dead individuals in Italy.