key: cord-347317-qcghtkk0
authors: Russo, Lucia; Anastassopoulou, Cleo; Tsakris, Athanasios; Bifulco, Gennaro Nicola; Campana, Emilio Fortunato; Toraldo, Gerardo; Siettos, Constantinos
title: Tracing day-zero and forecasting the COVID-19 outbreak in Lombardy, Italy: A compartmental modelling and numerical optimization approach
date: 2020-10-30
journal: PLoS One
DOI: 10.1371/journal.pone.0240649
sha: 
doc_id: 347317
cord_uid: qcghtkk0

INTRODUCTION: Italy became the second epicenter of the novel coronavirus disease 2019 (COVID-19) pandemic after China, surpassing by far China’s death toll. The disease swept through Lombardy, which remained in lockdown for about two months, starting from the 8th of March. As of that day, the isolation measures taken in Lombardy were extended to the entire country. Here, assuming that effectively there was one case “zero” that introduced the virus to the region, we provide estimates for: (a) the day-zero of the outbreak in Lombardy, Italy; (b) the actual number of asymptomatic infected cases in the total population until March 8; (c) the basic (R(0))and the effective reproduction number (R(e)) based on the estimation of the actual number of infected cases. To demonstrate the efficiency of the model and approach, we also provide a tentative forecast two months ahead of time, i.e. until May 4, the date on which relaxation of the measures commenced, on the basis of the COVID-19 Community Mobility Reports released by Google on March 29. METHODS: To deal with the uncertainty in the number of the actual asymptomatic infected cases in the total population Volpert et al. (2020), we address a modified compartmental Susceptible/ Exposed/ Infectious Asymptomatic/ Infected Symptomatic/ Recovered/ Dead (SEIIRD) model with two compartments of infectious persons: one modelling the cases in the population that are asymptomatic or experience very mild symptoms and another modelling the infected cases with mild to severe symptoms. The parameters of the model corresponding to the recovery period, the time from the onset of symptoms to death and the time from exposure to the time that an individual starts to be infectious, have been set as reported from clinical studies on COVID-19. For the estimation of the day-zero of the outbreak in Lombardy, as well as of the “effective” per-day transmission rate for which no clinical data are available, we have used the proposed SEIIRD simulator to fit the numbers of new daily cases from February 21 to the 8th of March. This was accomplished by solving a mixed-integer optimization problem. Based on the computed parameters, we also provide an estimation of the basic reproduction number R(0) and the evolution of the effective reproduction number R(e). To examine the efficiency of the model and approach, we ran the simulator to “forecast” the epidemic two months ahead of time, i.e. from March 8 to May 4. For this purpose, we considered the reduction in mobility in Lombardy as released on March 29 by Google COVID-19 Community Mobility Reports, and the effects of social distancing and of the very strict measures taken by the government on March 20 and March 21, 2020. RESULTS: Based on the proposed methodological procedure, we estimated that the expected day-zero was January 14 (min-max rage: January 5 to January 23, interquartile range: January 11 to January 18). The actual cumulative number of asymptomatic infected cases in the total population in Lombardy on March 8 was of the order of 15 times the confirmed cumulative number of infected cases, while the expected value of the basic reproduction number R(0) was found to be 4.53 (min-max range: 4.40- 4.65). On May 4, the date on which relaxation of the measures commenced the effective reproduction number was found to be 0.987 (interquartiles: 0.857, 1.133). The model approximated adequately two months ahead of time the evolution of reported cases of infected until May 4, the day on which the phase I of the relaxation of measures was implemented over all of Italy. Furthermore the model predicted that until May 4, around 20% of the population in Lombardy has recovered (interquartile range: ∼10% to ∼30%).

Italy became the second epicenter of the novel coronavirus disease 2019 (COVID-19) pandemic after China, surpassing by far China's death toll. The disease swept through Lombardy, which remained in lockdown for about two months, starting from the 8th of March. As of that day, the isolation measures taken in Lombardy were extended to the entire country. Here, assuming that effectively there was one case "zero" that introduced the virus to the region, we provide estimates for: (a) the day-zero of the outbreak in Lombardy, Italy; (b) the actual number of asymptomatic infected cases in the total population until March 8; (c) the basic (R 0 )and the effective reproduction number (R e ) based on the estimation of the actual number of infected cases. To demonstrate the efficiency of the model and approach, we also provide a tentative forecast two months ahead of time, i.e. until May 4, the date on which relaxation of the measures commenced, on the basis of the COVID-19 Community Mobility Reports released by Google on March 29.

To deal with the uncertainty in the number of the actual asymptomatic infected cases in the total population Volpert et al. (2020), we address a modified compartmental Susceptible/ Exposed/ Infectious Asymptomatic/ Infected Symptomatic/ Recovered/ Dead (SEIIRD) model with two compartments of infectious persons: one modelling the cases in the population that are asymptomatic or experience very mild symptoms and another modelling the infected cases with mild to severe symptoms. The parameters of the model corresponding to the recovery period, the time from the onset of symptoms to death and the time from a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 exposure to the time that an individual starts to be infectious, have been set as reported from clinical studies on COVID-19. For the estimation of the day-zero of the outbreak in Lombardy, as well as of the "effective" per-day transmission rate for which no clinical data are available, we have used the proposed SEIIRD simulator to fit the numbers of new daily cases from February 21 to the 8th of March. This was accomplished by solving a mixed-integer optimization problem. Based on the computed parameters, we also provide an estimation of the basic reproduction number R 0 and the evolution of the effective reproduction number R e . To examine the efficiency of the model and approach, we ran the simulator to "forecast" the epidemic two months ahead of time, i.e. from March 8 to May 4. For this purpose, we considered the reduction in mobility in Lombardy as released on March 29 by Google COVID-19 Community Mobility Reports, and the effects of social distancing and of the very strict measures taken by the government on March 20 and March 21, 2020.

Based on the proposed methodological procedure, we estimated that the expected dayzero was January 14 (min-max rage: January 5 to January 23, interquartile range: January 11 to January 18). The actual cumulative number of asymptomatic infected cases in the total population in Lombardy on March 8 was of the order of 15 times the confirmed cumulative number of infected cases, while the expected value of the basic reproduction number R 0 was found to be 4.53 (min-max range: 4.40-4.65). On May 4, the date on which relaxation of the measures commenced the effective reproduction number was found to be 0.987 (interquartiles: 0.857, 1.133). The model approximated adequately two months ahead of time the evolution of reported cases of infected until May 4, the day on which the phase I of the relaxation of measures was implemented over all of Italy. Furthermore the model predicted that until May 4, around 20% of the population in Lombardy has recovered (interquartile range: *10% to *30%).

The butterfly effect in chaos theory underscores the sensitive dependence on initial conditions, highlighting the importance of even a small change in the initial state of a nonlinear system. The emergence of a novel coronavirus, SARS-CoV-2, that caused a viral pneumonia outbreak in Wuhan, Hubei province, China in early December 2019 has evolved into the COVID-19 acute respiratory disease pandemic due to its alarming levels of spread and severity, with more than 3.5 million cases and 250,000 deaths globally, as of May 7, 2020 ( [1, 2] ). The seemingly far from the epicenter, old continent became the second-most impacted region after Asia Pacific, mostly as a result of a dramatic divergence of the epidemic trajectory in Italy first, where there have been 214,457 total confirmed infected cases and 29,684 deaths, and then in Spain where there have been 220,325 total confirmed infected cases and 25,857 deaths, as of May 7, 2020 ( [1, 2] ).

The second largest outbreak outside of mainland China officially started on January 31, 2020, after two Chinese visitors staying at a central hotel in Rome tested positive for SARS-CoV-2; the couple remained in isolation and was declared recovered on February 26 [3] . A 38-year-old man repatriated back to Italy from Wuhan who was admitted to the hospital in Codogno, Lombardy on February 21 was the first secondary infection case ("patient 1"). "Patient 0" was never identified by tracing the first Italian citizen's movements and contacts. In less than a week, the explosive increase in the number of cases in several bordering regions and in the two autonomous provinces of Trento and Bolzano (the northerner in Italy) placed enormous strain on the decentralized health system. Following a dramatic spike in deaths from COVID-19, Italy transformed into a "red zone", and the movement restrictions were expanded to the entire country on the 8th of March. All public gatherings were cancelled and school and university closures were extended through at least the next month.

In an attempt to assess the dynamics of the outbreak for forecasting purposes, it is important to estimate epidemiological parameters that cannot be computed directly based on clinical data, such as the transmission rate (or as otherwise called "effective contact rate" of the disease and the basic reproduction number, R 0 . The transmission rate is defined as the product of the probability of transmitting the virus given a contact between a susceptible and an infected individual and the average rate of contacts between susceptibles and infected. R 0 is defined as the expected number of exposed cases generated by one infected case in a population where all individuals are susceptible [4] .

Since the first confirmed COVID-19 case many mathematical modelling studies have already appeared. The first models mainly focused on the estimation of the basic reproduction number R 0 using dynamic mechanistic mathematical models ( [5] [6] [7] [8] ), but also simple exponential growth models (see e.g. [9, 10] ). Compartmental epidemiological models like SIR, SIRD, SEIR and SEIRD have been proposed to estimate other important epidemiological parameters, such as the transmission rate and for forecasting purposes (see e.g. [8, 11] ). Other studies have used metapopulation models, which include data of human mobility between cities and/or regions to forecast the evolution of the outbreak in other regions/countries far from the original epicenter in China [5, 7, 12, 13] , including the modelling of the influence of travel restrictions and other control measures in reducing the spread [14] .

Among the perplexing problems that mathematical models face when they are used to estimate epidemiological parameters and to forecast the evolution of the outbreak, two stand out: (a) the uncertainty regarding the day-zero of the outbreak, the knowledge of which is crucial to assess the stage and dynamics of the epidemic, especially during the first growth period, and (b) the uncertainty that characterizes the actual number of the asymptomatic infected cases in the total population (see e.g. [15, 16] ).

At this point we should note that what is done until now with dynamical epidemiological models is the investigation of several scenarios including different "days-zero" or just fixing the day-zero and run different levels of asymptomatic cases e.t.c. To cope with the above problems, we herein address a methodological framework that provides estimates for the day-zero of the outbreak and the number of the asymptomatic cases in the total population in a systematic way. Towards this goal, and for our demonstrations, we address a conceptually simple SEIRD model with a total of five compartments, with one of them modelling the asymptomatic infected cases in the population and another one modelling the part of the infected cases that will experience mild to severe symptoms, a significant share of which will be hospitalized, admitted to intensive care units (ICUs) or die from the disease. The proposed approach is applied to Lombardy, the epicenter of the outbreak in Italy. Furthermore, we provide a twomonths ahead of time forecast from March 8 (the day of lockdown of all Italy) to May 4 (the first day of the relaxation of the strict isolation measures). The above tasks were accomplished by the numerical solution of a mixed-integer optimization problem using the publicly available data of daily new cases for the period February 21-March 8, and the COVID-19 Community Mobility Reports released by Google on March 29.

We address a compartmental SEIIRD model that includes two categories of infected cases, namely the asymptomatic (unknown) cases in the total population and the cases that develop mild to more severe symptoms, a significant share of which are hospitalized, admitted to ICUs and a part of them dies. In agreement with other studies and observations, our modelling hypothesis is that the confirmed cases of infected are only a (small) subset of the actual number of asymptomatic infected cases in the total population [7, 8, 16] . Regarding the confirmed cases of infected as of February 21, a study conducted by the Chinese CDC which was based on a total of 72,314 cases in China, found that about 80.9% of the cases were mild and could recover at home, 13.8% severe and 4.7% critical [17] .

On the basis of the above findings, in our modelling approach, it is assumed that the asymptomatic or very mildly symptomatic cases recover from the disease relatively soon and without medical care, while for the other category of infected, on average their recovery lasts longer than the non-confirmed, they may also be hospitalized, admitted to ICUs or die from the disease.

Based on the above, let us consider a well-mixed population of size N. The state of the system at time t, is described by (see also Fig 1 for a schematic) S(t) representing the number of susceptible persons, E(t) the number of exposed, I(t) the number of asymptomatic infected persons in the total population who experience very mild or no symptoms and recover relatively soon without any other complications, I c (t) the number of infected cases who may develop mild to more severe symptoms and a significant part of them is hospitalized, admitted to ICUs or dies, R(t) the number of asymptomatic cases in the total population that recover, R c (t) the number of the recovered cases that come from the compartment of I c and D(t) the reported number of deaths. For our analysis, and for such a short period, we assume that the total number of the population remains constant. Based on demographic data, the total population of Lombardy is N = 10m; its surface area is 23,863.09 kmq and the population density is *422 (Inhabitants/Kmq).

The rate at which a susceptible (S) becomes exposed (E) to the virus is proportional to the density of infectious persons I. The proportionality constant is the "effective" disease transmission rate, say b ¼ � cp, where � c is the average number of contacts per day and p is the probability of infection upon a contact between a susceptible and an infected. Our main assumption here is that only a fraction, say � of the actual number of exposed cases E will experience mild to 

Tracing DAY-ZERO and forecasting the COVID-19 outbreak in Lombardy more severe symptoms denoted as I c (t) and a significant part of them will be hospitalized, admitted to ICUs or die. Thus, we assume that the infected persons that belong to the compartment I c go into quarantine at home or they are hospitalized, and, thus, it is assumed that for any practical means they don't transmit further the disease. Here, it should be noted that a wide testing policy may also result in the identification of asymptomatic cases belonging to the compartment I that would then be assigned to compartment I c . However, as a generally reported rule in Italy, tests were conducted only for those who presented for treatment with symptoms like fever and coughing. Thus, people who did not seek medical attention were tested very scarcely [18] [19] [20] . Thus, for any practical means the compartment I c reflects the reported confirmed infected cases.

A fraction of the I c cases that is given by the fatality ratio f = D(t)/(I c (t) + R c (t) + D(t)) dies with a mortality rate δ the inverse of which is the average time from the onset of symptoms to death, while the remaining part ((1 − f)) of the I c compartment recovers with a rate γ c , the inverse of which corresponds to the average time from the onset of symptoms to full recovery.

We note that while more compartments could conceptually be included, we aimed at keeping a low level of complexity in order to avoid the introduction of more parameters, and thus a model that would suffer from the "curse of dimensionality".

At this point we should note that on March 8, the date of the general lockdown, the number of confirmed infected cases was 3,372, the number of cases in ICUs was 399 and the number of hospitalized persons was 2,616 [21] . That is, until March 8, the number of confirmed cases was approximately equal to the number of hospitalized cases and the cases that were admitted to ICUs. Therefore, until March 8, any difference between the asymptomatic cases as represented by our model by the compartment I and the compartment I c would approximately reflect a level of under-reporting of the actual asymptomatic cases in the total population. We should also note that in the available data of reported cases [21] there is no distinction between cases that recovered at home and those that recovered at and were dismissed from hospitals. Thus, in the absence of such information, if one were to consider as a separate category the cases that are hospitalized, an extra parameter would have to be introduced (the fraction of recovered cases dismissed from hospitals). On one hand, such a piece of information is not available, and, on the other, such an attempt would add an extra degree of freedom that would need calibration or to be fixed at a certain value; however, due to the small size of the data and the "curse of dimensionality", this would also introduce unnecessary computational burden and further modelling uncertainty.

Thus, our discrete mean field compartmental SEIIRD model reads:

Sðt À 1ÞIðt À 1Þ À sEðt À 1Þ ð2Þ

The above system is defined in discrete time points t = 1, 2, . . ., with the corresponding initial condition at the very start of the outbreak (day-zero):

The the term fδI c (t − 1) in Eq 4 represents the fraction f of the I c cases that dies with a mortality rate γ and the term (1 − f)γ c I c (t − 1) represents the complementary part ((1 − f)) of the I c cases that recovers with a rate γ c .

The parameters of the model are:

is the average per-day "effective" rate at which an exposed person becomes infectious,

is the average per-day "effective" recovery rate within the group of asymptomatic cases in the total population,

is the average per-day "effective" recovery rate within the subset of the I c infected cases that finally recover,

• δ(d −1 ) is the average per-day "effective" mortality rate within the subset of I c infected cases that finally die,

• f is the probability that a I c case will die. Here, this, is given by the "emergent" case fatality ratio, computed as f

is the fraction of the actual (all) cases of exposed in the total population that enter to the compartment I c .

Here, we should note the following: as new cases of recovered and dead at each time t appear with a time delay (which is generally unknown but an estimate can be obtained by clinical studies) with respect to the corresponding infected cases, the above per-day rates are not the actual ones; thus, they are denoted as "effective/apparent" rates.

The values of the epidemiological parameters σ, γ, γ c , δ that were fixed in the proposed model were chosen based on clinical studies.

In particular, in many studies that use SEIRD models, the parameter σ is set equal to the inverse of the mean incubation period (time from exposure to the development of symptoms) of a virus. However, the incubation period does not generally coincide with the time from exposure to the time that someone starts to be infectious. Regarding COVID-19, it has been suggested that an exposed person can be infectious well before the development of symptoms [22] . With respect to the incubation period for SARS-CoV-2, a study in China [23] suggests that it may range from 2-14 days, with a median of 5.2 days. Another study in China, using data from 1,099 patients with laboratory-confirmed 2019-nCoV ARD from 552 hospitals in 31 provinces/provincial municipalities suggested that the median incubation period is 4 days (interquartile range: 2 to 7). In our model, as explained above, 1/σ represents the period from exposure to the onset of the contagious period. Thus, based on the above clinical studies, for our simulations, we have set 1/σ = 3.

Regarding the recovery period, in a study that is based on 55,924 laboratory-confirmed cases, the WHO-China Joint Mission has reported a median time of 2 weeks from onset to clinical recovery for mild cases, and 3-6 weeks for severe or critical cases [24] . Based on the above, and on the fact that within the subset of confirmed cases the mild cases are the 81% [17] , we have set the recovery period for the confirmed cases' compartment to be δ c = 1/21 in order to balance the recovery period with the corresponding characterization of the cases (mild, severe/critical). The average recovery period of the unreported/non-confirmed part of the infected population, which in our assumptions experiences the disease like the flu or a common cold, is set equal to one week [25], i.e. we have set δ = 1/7. This choice is based also on reports on the serial interval of COVID-19. The serial interval of COVID-19 is defined as the time duration between a primary case-patient (infector) having symptoms and a secondary case-patient having again symptoms. For example, it has been reported that the serial interval for COVID is estimated at 4.4-7.5 days [26] ; for the case of Lombardy, the average serial number has been estimated to be 6.6 days [27] . In our model, the 1/σ = 3 period refers to the period from exposure to the onset of the contagiousness. In this period, obviously there are no symptoms. Thus, the serial interval in our model is 7 days (this is the average number of days in which an infectious becomes recovered and no longer transmits the disease). Importantly, there are studies (see e.g. Nushiura et al. [28] ) suggesting that a substantial proportion of secondary transmission may occur prior to illness onset. Thus, the 7 days period that we have taken as the average period that an infectious person can transmit the disease before he/she recovers, reflects exactly this period; it refers to the serial interval for the cases that are asymptomatic and for cases with mild symptoms.

Finally, the median time from the onset of symptoms until death for Italy has been reported to be eight days [29] , thus in our model we have set γ = 1/8.

We have set f = 11% for the optimization. For the forecasting (i.e. for the period March 9 to May 4), the value of f was not fixed but it was computed dynamically each day t through the model simulations as

The transmission rate β, as it cannot be obtained in general by clinical studies, but only by mathematical models, was estimated through the optimization process.

Regarding day-zero in Lombardy, that is also unknown and estimated by the optimization process, what has been officially reported is just the date on which the first infected person was confirmed to be positive for SARS-CoV-2. That day was February 21, 2020, which is the starting date of public data release of confirmed cases.

The day-zero of the outbreak, the per-day "effective" transmission rate β, and the ratio � were computed by the numerical solution of a mixed-integer optimization problem with the aid of genetic algorithms to fit the reported data of daily new cases (see the discussion in [30] ) from February 21 to March 8, the day of the lockdown of Lombardy.

As already mentioned, on March 8, the number of confirmed infected cases in the population was 3,372, the number of cases in ICUs was 399 and the number of hospitalized persons was 2,616 [21] . That is, until March 8 the number of mild and severe cases that were hospitalized and admitted to ICUs was approximately equal to the number of confirmed infected cases. Thus, for the period of calibration, it is reasonable to assume that for any practical means the number of confirmed cases was approximately the same as for those that experienced mild to more severe symptoms and were admitted for medical care. Thus, until March 8, the parameter � reflects also the level of under-reporting of the asymptomatic cases in the total population.

Here, for our computations, we have used the genetic algorithm "ga" provided by the Global Optimization Toolbox of Matlab [6] to minimize the following objective function:

where,

ΔX SEIRD (t), (X = I, R, D) are the new cases resulting from the SEIRD simulator at time t.

The weights w 1 , w 2 , w 3 correspond to scalars serving in the general case as weights to the relevant functions for balancing the different scales between the number of infected, recovered cases and deaths. The convergence tolerance was set to 1.E −06 , the population size (distributed between the lower and upper bounds) was selected as 40 resulting from min(max (10 � nvars,40),100) for mixed-integer problems (here nvars = 3), ceil(0.05 � PopulationSize) individuals are guaranteed to survive to the next generation, the migration fraction was set to 0.2, while the number of generations that take place between migrations of individuals between subpopulations was set to 20 [6] .

At this point we should note that the above optimization problem may in principle have multiple nearby optimal solutions (MNOS). Finding and assessing the information contained from MNOS (known also as niching) is a particular challenging problem [31] . Here, we created a grid of initial guesses within the intervals in which the optimal estimates were sought: for the day-zero (t 0 ) we used a step of 2 days within the interval December 27, 2019 until the 5th of February, 2020 i.e. ±20 days around the 16th of January, for β we used a step of 0.05 within the interval (0.3, 0.9) and for � we used a step of 0.02 within the interval (0.01, 0.29). The numerical optimization procedure was repeated 48 times for each combination of initial guesses. For our computations, we kept the best fitting outcome for each combination of initial guesses.

Next, in order to reveal structured patterns of distributions vs. uniformly random distributions, we fitted the resulting probability distributions of the optimal values using several functions, including the Normal, Log-normal, Weibull, Beta, Gamma, Burr [32] , Exponential and Birnbaum-Saunders [33] distributions and kept the one resulting in the maximum Log-likelihood (see in the Supporting Information for more details). For the computed parameters of the corresponding best distributions, we also provide the corresponding 95% confidence intervals. Such fitting can demonstrate that the obtained values are not uniformly distributed.

Finally, we have run simulations based on all obtained values to assess the efficiency of the model and obtained results and the forecasting uncertainty until May 4.

For our computations, we used the parallel computing toolbox of Matlab 2020a [34] utilizing 6 INTEL XEON CPU X5650 cores at 2.66GHz.

Estimation of the basic reproduction number. Here, we note that we provide an estimation of the basic reproduction number R 0 based on the estimation of the total number of (asymptomatic) infected cases in the population. Thus, it is expected that the estimated R 0 will be larger than the ones reported using just the confirmed number of cases; the latter may underestimate the actual R0.

Initially, when the spread of the epidemic starts, all the population is considered to be susceptible, i.e. S � N. On the basis of this assumption, we computed the basic reproduction number based on the estimates of the epidemiological parameters computed using the data from the 21st of February to the 8th of March with the aid of the SEIIRD model given by Eqs (1)-(7) as follows.

Note that there are three infected compartments, namely E, I, I c and two of them (E,I) determine the outbreak. Thus, considering the corresponding equations given by Eqs (2), (3) and (4) , and that at the very first days of the epidemic S � N and D � 0, the Jacobian of the system as evaluated at the disease-free state reads:

The eigenvalues (that is the roots of the characteristic polynomial of the Jacobian matrix) dictate if the disease-free equilibrium is stable or not, that is if an emerging infectious disease can spread in the population. In particular, the disease-free state is stable, meaning that an infectious disease will not result in an outbreak, if and only if all the norms of the eigenvalues of the Jacobian J of the discrete time system are bounded by one. Jury's stability criterion [35] (the analogue of Routh-Hurwitz criterion for discrete-time systems) can be used to determine the stability of the linearized discrete time system by analysis of the coefficients of its characteristic polynomial. The characteristic polynomial of the Jacobian matrix reads:

where a 2 ¼ 1

The necessary conditions for stability read:

The sufficient conditions for stability are given by the following two inequalities:

The first inequality (13) results in the necessary condition bð1 À �Þ g < 1:

It can be shown that the second necessary condition (14) and the sufficient condition (15) are always satisfied for the range of values of the epidemiological parameters considered here.

Thus, the necessary condition (16) is also a sufficient condition for stability. Hence, the disease-free state is stable, if and only if, condition (16) is satisfied.

Note that in this necessary and sufficient condition (16) , the fraction (1 − �)/δ is the average infection time of the compartment I. Thus, the above expression reflects the basic reproduction number R 0 which is qualitatively defined by R 0 = β � 1/infection time. Hence, our model results in the following expression for the basic reproduction number:

Note that for � = 0, the above expression simplifies to R 0 for the simple SIR model. Estimation of effective reproduction number. For the estimation of the effective reproduction number R e , representing the average number of secondary infections from an infectious individual when in the population exist already non-susceptible individuals. For the calculation of R e , we now use the next generation matrix approach [36] .

The next generation matrix G with elements = g ij is formed by the average number of secondary infections of type i from an infected individual of type j. Formally, it is constructed by:

where F is the vector containing the transmission rates from the model, and V is the vector containing the transition rates between the infected compartments. In our model:

The effective reproduction number R e is the spectral radius, i.e. the dominant eigenvalue of G. Thus, R e ðtÞ ¼ bð1 À �Þ g Sðt À 1Þ N À Dðt À 1Þ À R c ðt À 1Þ À I c ðt À 1Þ

: ð20Þ

As discussed, we used the proposed approach to forecast the evolution of the pandemic in Lombardy from March 8 to May 4, i.e. from the first day of lockdown to the first day of the relaxation of the social isolation. Our estimation regarding the as of March 8 reduction of the "effective" transmission rate was based on the combined effects of prevention efforts and behavioral changes. In particular, our estimation was based on (a) the COVID-19 Community Mobility Reports released by Google on March 29 [37] , and, (b) an assessment of the synergistic effects of such control measures as the implementation of preventive containment in workplaces, stringent "social distancing", and the ban on social gatherings, as well as the public awareness campaign prompting people to adopt cautious behaviors to reduce the risk of disease transmission (see also [38] [39] [40] [41] ). The effect of the distribution of contacts at home, work, when travelling, and during leisure activities can be also assessed. For example, based on an analysis for the social contacts and mixing patterns relevant to the spread of infectious diseases that was conducted in various countries, it has been found that for Italy, around 20-23% of all physical social contacts during a day are attributed to workplaces, around 17-18% to schools, 2-3% to transportation, 20% to leisure activities, 15-18% to home, 15% to other activities (contacts made at locations other than home, work, school, travel, or leisure) and a 7-8% to contacts made at multiple other locations during the day, not just at a single location [42] .

On the basis of the Google COVID-19 Community Mobility Report released on March 29 [37] , the average reduction in the mobility in Lombardy during the period February 16-March 8, compared to the period before February 16, was *15% in retail & recreation activities, *20% in transit stations, *12% in workplaces, while it was increased in parks by * 11% and was almost the same in groceries and pharmacies. In the period March 9 to March 20, the mobility was reduced by an average of *73% in retail & recreation activities, by *75% in transit stations, by *55% in workplaces, by *58% in parks and by * 32% in groceries and pharmacies. Thus, taking into account the coarse effect of different activities in the physical contact [42] , the average reduction was of the order of *40% in mobility when compared to the period February 21-March 8. In fact, on March 17, based on the release of mobile phone data, the vice-president of Lombardy, announced that the average mobility in the region (for distances more than 500 meters) had been reduced by a *50% with respect to the period before February 20 [43] . On March 20, the government announced the implementation of even stricter measures that included the closure of all public and private offices, closing all parks, walking only around the residency and not even in pairs, and the prohibition of mobility to second houses [44, 45] . According to the Google COVID-19 Community Mobility Reports [37] , from March 20-21 until April 30, activities were reduced by an average of *87% in retail & recreation activities, by *84% in transit stations, by *70% in workplaces, by *80% in Parks and by * 50% in groceries and pharmacies. Thus, taking into account the coarse effect of different activities in the physical contact [42] , the average reduction was of the order of *65% in the mobility when compared to the period of February 21-March 8.

A further reduction may be attributed to behavioral changes [46] . For example, it has been shown that social distancing and cautiousness reduce the disease transmission rate by about 20% [40] . Thus, based on the above, it is reasonable to consider a (1-0. 

As discussed, for our computations we ran 48 times the numerical optimization procedure for each combination of initial guesses based on the daily reported new cases from February 21 to For all the near-optimal points obtained using the genetic algorithm optimization, the residuals were of the order of * 4,750,000. Regarding the values of the optimal parameters, we fitted their cumulative probability distributions using several functions, including the Normal, Log-normal, Weibull, Beta, Gamma, Burr, Exponential and Birnbaum-Saunders functions, and kept the one yielding the maximum Log-likelihood (see in the S1 File). Table 1 summarizes the mean values of the optimal parameters and their interquartiles, for the day-zero, the transmission rate (β) and the fraction (�) of the actual cases of exposed in the total population that enter to the compartment I c , and also the information on the values of the parameters of the best fitting distributions.

Note that the optimal values of day-zero were between January 5-January 23 (interquartile range: January 11 to January 18) (see S1 Fig in S1 File) , the optimal values of β were between 0.636 and 0.86 (interquartile range: 0.665 to 0.719) (see S2 Fig in S1 File) , and the optimal values of � were between 0.01 and 0.249 (interquartile range: 0.023 to 0.1) (see S3 Fig in S1 File) . The best fit to the distribution of optimal values of the day-zero was obtained using a Normal CDF with mean 1.67 (i.e. sim 2 days before the 16th of January) (95% CI:1.53, 1.80) and variance 4.43 (95% CI: 4.33, 4.52); thus taking the round value at 2 days, the expected day-zero corresponds to January 14 (interquartile range: January 11 to January 18).

The best fit to the distribution of the optimal values of β, was given by fitting a Burr CDF with α = 0.654 (95% CI: 0.653, 0.655), c = 276.061 (95% CI: 248.911, 306.172), k = 0.0537 (95% CI: 0.048, 0.060) having a mean value of 0.7 (interquartile range: 0.665 to 0.719). Finally, the best fit to the distribution of the optimal values of � was given by fitting a Birnbaum-Saunders CDF with parameters μ = 0.0492 (95% CI: 0.048, 0.0505) (scale parameter) and α = 0.934 (95% CI: 0.914, 0.954) (shape parameter), resulting in an expected value 0.07 (interquartile range: 0.023 to 0.1) (see in S1 File).

Thus, based on the derived values of the "effective" per-day disease transmission rate, the basic reproduction number R 0 is 4.53 (min-max range: 4.40-4.65).

Finally, we ran the simulator for all values of the optimal triplets from the corresponding distributions as found by the solution of the optimization problem using the data from February 21 until March 8, then from March 9 to March 20 for validation purposes and from March 29 to May 4 for forecasting purposes.

Figs (2)-(4) depict the simulation results based on the optimal estimates, until the 19th of March. To validate the model with respect to the reported data of confirmed cases from March 9 to March 19, we have considered an (1-0.4)(1-0.2) reduction in the "effective" transmission rate and as initial conditions the values resulting from the simulation on March 8 as described Table 1 . Optimal parameter values and interquartiles for day-zero, transmission rate (β) and fraction (�) of the cases of exposed individuals in the total population that enter to the compartment I c . The resulting value of R 0 along with the minimum and maximum value is also given. 

Tracing DAY-ZERO and forecasting the COVID-19 outbreak in Lombardy in the methodology. Thus, the model approximated fairly well the dynamics of the pandemic in the period from February 21 to March 19. As discussed in the Methodology, we also attempted based on our methodology and modelling approach to forecast the evolution of the outbreak until May 4, the first day of the relaxation of the measures. To do so, as described in the methodology, we have considered a (1-0.65) (1-0.2) reduction in the "effective" transmission rate starting on March 20 (compared to the period February 21-March 8), the day of announcement of even stricter measures in the region of Lombardy (see in Methodology). The result of our forecast is depicted in Fig 5. As shown, the model predicts fairly the evolution of the epidemic two months ahead of March 8 (the model parameters and the day-zero were estimated using the reported data from February 21 to March 8). Note that, regarding the confirmed cases, the mean values of I c (t) from over all simulations almost coincide with the reported values of the confirmed cases (see top panel of Fig 5) .

The reported recovered cases are in the lower part of the model predictions (R c (t)) (see medium panel of Fig 5) . Also, the model predicts fairly the total number of deaths until May 4 (see bottom panel of Fig 5. However, at May 4 there is a significant difference between the mean value of deaths as obtained from simulations and the actual number of deaths. This is due to the following reasons. First, the difference that is observed with respect to the reported number of deaths and model forecasts in the period from March 18 to April 15 can be attributed to facts that the model did not take into account, such as the saturation of ICUs in that period, which could potentially lead to a larger number of deaths. Indeed, on March 8, 400 new deaths were reported. However, we note that the model was calibrated with the data reported until March 8. Until March 8, the case fatality ratio was of the order of *9-10% and the number of cases admitted to ICUs was relatively small. After March 20 when there were many cases admitted to ICUs, the death toll increased significantly, it almost doubled, thus reaching the *17-18% of the confirmed cases. Second, regarding the difference that is observed between the model predictions and actual number of deaths at May 4, this is due to the fact that we did not calibrate the model parameters to take into account the changes in the death and recovery rates. It is expected that when the evolution of the epidemic is abrupt, as it was in Lombardy in the period March 12 until the early of April, due to the saturation of the 

Tracing DAY-ZERO and forecasting the COVID-19 outbreak in Lombardy health system and ICUs the mortality rates would be higher and the recovery rates longer. Indeed, during this period, it has been reported, that the median time from the onset of symptoms until death in Lombardy was eight days [29] which a very short period. For example, in other studies it has been reported that the time between symptom onset and death ranged from about 2 weeks to 8 weeks which from two to seven times longer than the one reported in Lombardy [24, 47] . So, this difference explains the difference between simulations and reported number of deaths at May 4. Here, for our demonstrations we have used the information about the epidemiological parameters, that was available at the early phase of the epidemic. Thus, we must underline the very good agreement of the model predictions with the actual number of confirmed infected cases for the entire period March 8-May 4. This number, in contrast to the number of recovered and number of deaths is not shaped/biased by the capacity of ICUs and this is the number on which policy makers should focus in order to decide ahead of time of the necessary resources (e.g. necessary number of beds and ICUs) needed to keep the fatality ratio as low as possible.

The model predicts that until May 4, an average of 20% of the population in Lombardy has already recovered (interquartile range: *10% to *30%) (see Fig 6) .

Finally, in Fig 7, we report the evolution of the effective reproduction number R e until May 4.

On May 4 the estimated mean value of R e was 0.987 (interquartile range: 0.857 to 1.133), thus marking the onset of a critical date for the post-lockdown period. 

Tracing DAY-ZERO and forecasting the COVID-19 outbreak in Lombardy

The crucial questions about an outbreak is how, when (day-zero), why it started, and if and when it will end. Answers to these important questions would add critical knowledge to our arsenal to combat the pandemic. The tracing of day-zero, in particular, is of outmost importance. It is well known, that minor perturbations in the initial conditions of a complex system, such as the ones of an outbreak, may result in major changes in the observed dynamics. Undoubtedly, a high level of uncertainty for day-zero, as well as the uncertainty in the actual numbers of exposed people in the total population, raise several barriers to our ability to correctly assess the state and dynamics of the outbreak and to forecast its evolution and its end. Such pieces of information would lower the barriers and help public health authorities respond fast and efficiently to the emergency. For example, an over-or under-estimation of day-zero would result in an under-or over-estimation of the the transmission rate β and therefore of the basic reproduction number R 0 and consequently the effective reproduction number R e . Furthermore, the correct estimation of day-zero is important for the assessment of the number of asymptomatic and actual recovered cases in the total population. This is turn will bias the assessment and ultimately the design of efficient control policies in real time.

This study aimed exactly at shedding more light into this problem. To achieve this goal, we addressed a conceptually simple compartmental SEIIRD model with two infectious compartments in order to bridge the gap between the number of asymptomatic cases in the total population and the cases that will experience mild to more severe symptoms.

What is done until now with mathematical epidemiological models is the investigation of several scenarios, by changing e.g. the (or assuming a fixed) initial day (day-zero), the level of asymptomatic cases etc. Our work, is the first that introduces a methodological framework, to estimate the day-zero as well as the level of asymptomatic cases in the total population in a 

Tracing DAY-ZERO and forecasting the COVID-19 outbreak in Lombardy systematic way. Following the proposed methodological framework, we found that the dayzero in Lombardy was around the middle of January, a period that precedes by one month the fate of the first confirmed case in the hardest-hit northern Italian region of Lombardy. Interestingly enough, when we submitted a preprint of the work at MEDRXIV on March 20 2020 ( [48] , another study that was also submitted at the same day at MEDRXIV, that was based on genomic and phylogenetic data analysis, reports the same time period, between the second half of January and early February, 2020, as the time when the novel coronavirus SARS-CoV-2 entered northern Italy [49] . Our analysis further revealed that the actual number of asymptomatic infected cases in the total population in the period until March 8 was around 15 times the number of confirmed infected cases, which until March 8 was also approximately equal to the number of cases that were hospitalized and admitted to ICUs.

Our model and methodological approach assume that there was one effective "zero" infected case that introduced the virus to the region; one could certainly argue that there were more than one cases that introduced the virus to the region on the same day; such scenarios can be investigated in a straightforward manner based on our proposed methodological approach. Furthermore, the proposed approach could be used for the quantification of the uncertainty of the evolving dynamics, taking into account the reported, from clinical studies, distributions of the epidemiological parameters rather than their expected values. A critical point that is connected with the above is that with such a small number of infectious at the initial stage of the simulations, a stochastic model or a hybrid stochastic model could be more realistic in which uncertainty could be modelled in the form of realistic perturbations (see for example [50] ). Furthermore, we did not consider the effect of an ongoing sampling strategy in the total population for the estimation of the level of under-reporting (represented in our model with the variable �). As mentioned in the methodology, within the first period of the outbreak the number of confirmed cases was approximately equal to the number of hospitalized cases, i.e. there was no sampling strategy. Furthermore, as also reported for the later period, tests were conducted only for those who seeked medical care and had symptoms like fever and coughing. Thus, people who did not seek for medical attention were tested very scarcely [18] [19] [20] . We will consider such type of modelling and analysis in a future work.

Regarding the forecasting in Lombardy from March 8 until May 4 (the first day of relaxation of the measures), we have taken into account the very latest facts on the drop of human mobility, as released by Google [37] until April 30 for the region of Lombardy; these were shaped by the very strict measures announced on March 20-21 that included the closure of all parks, public and private offices and the prohibition of any pedestrian activity, even individually [45] . Our modelling approach approximated fairly well the reported number of infected cases in Lombardy two months ahead of time. The mean value of the evolution of the compartment that in our model reflects the confirmed cases, almost coincides with the reported cases for the entire period from March 8 to May 4. The differences that are observed between the reported number of deaths and simulations as discussed in the results should be attributed to the very short times from the onset of symptoms to death that have been reported in Lombardy at the early phases of the pandemic in the region (8 days instead of 2 to 8 weeks that have been reported in other studies for China that is linked with the saturation of the ICUs [47, 51] . Furthermore another important factor that is missing from the data used is that the number of deaths is largely affected by the criteria of death notification. Indeed as it has been reported, the global coronavirus death toll could be 60% higher than the confirmed [52] . However, this fact did not affect the model predictions with respect to the confirmed cases, but the rate with which the confirmed cases die; the later is also dependent on the rate of the outbreak and the relevant capacity of the ICUs.

To this end, we would like to make a final comment with respect to the basic reproduction number R 0 , the significance and meaning of which are very often misinterpreted and misused, thereby leading to erroneous conclusions. Here, we found a R 0 * 4.5, which is higher compared to the values reported by many studies in China, and also in Italy, and in Lombardy in particular. For example, Zhao et al. estimated R 0 to range between 2.24 (95% CI: 1.96, 2.55) and 3.58 (95% CI: 2.89, 4.39) in the early phase of the outbreak [9] . Similar estimates were obtained for R 0 by Imai et al., 2.6 (95% CI: 1.5, 3.5) [ [8] .

Regarding Italy, D'Arienzo and Coniglio [54] used a SIR model to fit the reported data in nine Italian cities and found that R 0 ranged from 2.43 to 3.10. In another study, the authors provided an estimate of the basic reproduction number by analyzing the first 5,830 laboratoryconfirmed cases. By doing so they estimated the basic reproduction number at 3.1 [27] .

First, we would like to stress that R 0 is not a biological constant for a disease as it is affected not only by the pathogen, but also by many other factors, such as environmental conditions, demographics, as well as, importantly, by the social behavior of the population (see for example the discussion in [4] ). Thus, a value for R 0 that is found in one part of the world (e.g. in China), and even in a region of the same country, e.g. in Tuscany, Italy, cannot be generalized as a global biological constant for other parts of the world, or even for other regions of the same country. Obviously, the environmental factors and social behavior of the population in Lombardy are different from the ones, for example, prevailing in Hubei.

Second, most of the studies that provide estimates of the basic reproduction number are based solely on the reported cases, thus the actual number of infected cases in the total population, that may be asymptomatic but transmit the disease, is not considered; this fact may lead to an underestimation of the basic reproduction number. Moreover, in our approach as compared to clinical studies, the computation of R 0 comes out as the necessary and sufficient condition of the stability as derived from the proposed model whose parameters are computed based on the available reported data, thus with a delay between the first actual case (see also the discussion in [55] .

Our relatively simple conceptually model and approach do not aspire to accurately describe the complexity of the emergent dynamics, which in any case is an overwhelmingly difficult, if not impossible, task in the long run, even with the use of detailed agent-based models. We tried to keep the structure of the model as simple as possible in order to be able to model (in a coarse way) the uncertainty in both the "day-zero" and the number of asymptomatic actual cases in the total population using as few parameters as possible. The results of our analysis have indeed proved that the modelling approach succeeded in providing fair predictions of the evolution of the epidemic two months ahead of time. Such an early assessment would help authorities to evaluate the required measures to control the epidemic, such as the scale of diagnostic tests that have to be performed and the number of ICU beds required. While more complicated models can, in principle, be constructed to take into account more detailed information, such as the number of hospitalized patients and patients in ICUs, for any practical means such an approach would suffer from the "curse of dimensionality" as it would introduce many more parameters that would need calibration based on a relatively small size of data especially at the beginning of an outbreak. An attempt to compute the values of some of these additional parameters, which can only be roughly estimated by clinical studies at the early stages of an emerging novel infectious disease, would introduce additional uncertainty, thereby further complicating matters rather than solving the problem.

To this end, we hope that our conceptually simple, but pragmatic, modelling approach and methodological framework help to provide improved insights into the currently uncontrolled pandemic and to contribute to the mitigation of some of its severe consequences.

Author Contributions

Coronavirus disease 2019 (COVID-19). Situation report

Coronavirus COVID-19 Global Cases by

Covid-19: preparedness, decentralisation, and the hunt for patient zero

Complexity of the basic reproduction number (r0)

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in wuhan, china: a modelling study

Transmissibility of 2019-ncov

Estimating the scale of covid-19 epidemic in the united states: Simulations based on air traffic directly from wuhan, china

Data-based analysis, modelling and forecasting of the covid-19 outbreak

Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in china, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak

COVID-19 and italy: what next?

Breaking down of the healthcare system: Mathematical modelling for controlling the novel coronavirus (2019-nCoV) outbreak in wuhan

Estimating the risk on outbreak spreading of 2019-ncov in china using transportation data

Using predicted imports of 2019-ncov cases to determine locations that may not be identifying all imported cases

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

Coronavirus-scientific insights and societal aspects

Current US coronavirus cases are "just the tip of the iceberg,-former USAID director says

The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases (covid-19) in china

In one Italian town, we showed mass testing could eradicate the coronavirus

Perch& x00E9; non stiamo facendo pi& x00F9; tamponi?

Covid-19-Situazione in Italia

Prevention, How COVID-19 Spreads

Early transmission dynamics in wuhan, china, of novel coronavirus infected pneumonia

Report of the WHO-China Joint Mission on Coronavirus Disease

How will country-based mitigation measures influence the course of the COVID-19 epidemic?

The early phase of the covid-19 outbreak in lombardy, italy (2020)

Serial interval of novel coronavirus (COVID-19) infections

Superiore di Sanit Characteristics of COVID-19 patients dying in Italy Report based on available data on

Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to ebola

Niching methods and multimodal optimization performance

Cumulative frequency functions

A new family of life distributions

Inners and stability of dynamic systems

The construction of next-generation matrices for compartmental epidemic models

Advancing the right to health: the vital role of law, World Health Organization

Case investigations of infectious diseases occurring in workplaces, united states

Quantifying social distancing arising from pandemic influenza

Nonpharmaceutical measures for pandemic influenza in nonhealthcare settings-social distancing measures, Emerging Infectious Diseases

Social contacts and mixing patterns relevant to the spread of infectious diseases

Il 40% ancora si sposta

A Milano militari per le strade. Posti di blocco a Roma

Coronavirus, ordinanza Lombardia con nuove limitazioni

Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases

Real estimates of mortality following COVID-19 infection

Tracing day-zero and forecasting the covid-19 outbreak in lombardy, italy: A compartmental modelling and numerical optimization approach

Genomic characterisation and phylogenetic analysis of sars-cov-2 in italy

Bounded Noises in Physics, Biology, and Engineering

WHO Statement Regarding Cluster of Pneumonia Cases in Wuhan

Global coronavirus death toll could be 60% higher than reported | Free to read

Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus Infected Pneumonia

Assessment of the SARS-CoV-2 basic reproduction number, r0, based on the early phase of COVID-19 outbreak in italy

A new framework and software to estimate timevarying reproduction numbers during epidemics