Tracing DAY-ZERO and Forecasting the Fade out of the COVID-19 Outbreak in Lombardy, Italy: A Compartmental Modelling and Numerical Optimization Approach

Introduction

The butterfly effect in chaos theory underscores the sensitive dependence on initial conditions, highlighting the importance of even a small change in the 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 a total of 195,892 confirmed infected cases, 80,840 recovered and 7,865 deaths in 153 countries as of March 16, 2020 ([1]). The seemingly far from the epicenter, old continent became the second-most impacted region after Asia Pacific to date, mostly as a result of a dramatic divergence of the epidemic trajectory in Italy, where there have been 31,506 total confirmed infected cases, 2,941 recovered and 2,503 deaths as of March 17, 2020 ([1]).

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 [2]. 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 autonomous provinces of northern Italy placed enormous strain on the decentralized health system. Following an 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 as well as to estimate epidemiological parameters that cannot be computed directly based on clinical data, such as the transmission rate of the disease and the basic reproduction number R₀, defined as the expected number of exposed cases generated by one infected case in a population where all individuals are susceptible, many mathematical modelling studies have already appeared since the first confirmed COVID-19 case. The first models mainly focused on the estimation of the basic reproduction number R₀ using dynamic mechanistic mathematical models ([3, 4, 5, 6]), but also simple exponential growth models (see e.g. [7, 8]). 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. [6, 9]). 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 [3, 10, 11, 5], including the modelling of the influence of travel restrictions and other control measures in reducing the spread ([12].

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 that characterizes the actual number of infected cases in the total population, which is mainly due to the large percentage of asymptomatic or mild cases experiencing the disease like the common cold or the flu (see e.g. [13]), and (b) 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.

To cope with the above problems, we herein propose a novel SEIRD with two compartments, one modelling the total infected cases in the population and another modelling the confirmed cases. The proposed modelling approach is applied to Lombardy, the epicenter of the outbreak in Italy, to estimate the scale of under-reporting of the number of actual cases in the total population, the DAY-ZERO of the outbreak and for forecasting purposes. The above tasks were accomplished by the numerical solution of a mixed-integer optimization problem using the publicly available data of cumulative cases for the period February 21-March 8, the day of lockdown of all of Italy.

Methodology

The modelling approach

We address a compartmental SEIRD model that includes two categories of infected cases, namely the confirmed/reported and the unreported (unknown) cases in the total population. Based on observations and studies, our modelling hypothesis is that the confirmed cases of infected are only a (small) subset of the actual number of infected cases in the total population [5, 13, 6]. Regarding the confirmed cases of infected as of February 11, a study conducted by the Chinese CDC which was based on a total of 72,314 cases in China, about 80.9% of the cases were mild and could recover at home, 13.8% severe and 4.7% critical [14].

On the basis of the above findings, in our modelling approach, the unreported cases were considered either asymptomatic or mildly symptomatic cases that recover from the disease relatively soon and without medical care, while the confirmed cases include all the above types, but on average their recovery lasts longer than the non-confirmed, they may also be hospitalized and 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 Figure 1 for a schematic) S(t) representing the number of susceptible persons, E(t) the number of exposed, I(t) the number of unreported infected persons in the total population who are asymptomatic or experience mild symptoms and recover relatively soon without any other complications, I_c(t) the number of confirmed infected cases who may develop more severe symptoms and a part of them dies, R(t) the number of recovered persons in the total population, R_c(t) the number of confirmed recovered cases and D(t) the 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 sim422 (Inhabitants/Kmq).

• Download figure
• Open in new tab

Figure 1:

A schematic of the proposed compartmental SEIRD model. The actual number of cases are unknown.

The rate at which a susceptible (S) becomes exposed (E) to the virus is proportional to the density of infectious persons I in the total population, excluding the number of dead persons D. Our main assumption here is that upon confirmation, the infected persons I_c go into quarantine, and, thus, they don’t transmit further the disease. The proportionality constant is the “effective” disease transmission rate, say where is the average number of contacts per day and p is the probability of infection upon a contact between a susceptible and an infected.

Thus, our discrete mean field compartmental SEIRD model reads:

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): S(0) = N − 1, I(0) = 1, E(0) = 0 I_c(0) = 0, R(0) = 0, R_c(0) = 0, D(0) = 0.

The parameters of the model are:

• β(d⁻¹) is the “effective” transmission rate of the disease,

• σ(d⁻¹) is the average per-day “effective” rate at which an exposed person becomes infective,

• δ(d⁻¹) is the average per-day “effective” recovery rate within the group of unreported (asymptomatic/mild) cases in the total population,

• δ_c(d⁻¹) is the average per-day “effective” recovery rate within the subset of confirmed infected cases

• γ(d⁻¹) is the average per-day “effective” mortality rate within the subset of confirmed infected cases,

• ϵ(d⁻¹) is the per-day rate of the all cases of infected in the total population that get confirmed. This proportionality rate quantifies the uncertainty in the actual number of unreported cases in the total population.

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 [15]. With respect to the incubation period for SARS-CoV-2, a study in China [16] 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, 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 .

Regarding the recovery period, the WHO-China Joint Mission in a study that is based on 55,924 laboratory-confirmed cases 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 [17]. Based on the above and on the fact that within the subset of confirmed cases the mild cases are the 81% [14], 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 [18], i.e. we have set δ = 1/7.

Note, that per-day mortality rate in the model does not coincide with the case fatality (mortality) ratio (CFR) which is usually reported and computed as the ratio between the reported deaths and the cumulative number of infected. Furthermore, the transmission rate cannot be obtained by clinical studies, but only by mathematical models.

Finally, regarding DAY-ZERO in Lombardy, what has been 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.

Results

As discussed in the Methodology, we used five different values of E (0.01, 0.05, 0.1, 0.15, 0.2) to assess the actual number of cases in the total population. Thus, for our computations, we run 50 times the numerical optimization procedure and for further analysis we kept the value of E that gave the smaller fitting error over all runs. In particular, by using the genetic algorithm as described in the Methodology, the best fitting to the reported data was obtained with E=0.05 (best fitting residual norm of the objective function: 301,231). For this scale of under-reporting, we present the results for DAY-ZERO, the “effective” per-day transmission and mortality rates as computed with the reported data from February 21 to March 8. The solution of the mixed-integer optimization problem for the DAY-ZERO and the “effective” per-day transmission and mortality rates of the compartmental SEIRD model resulted in the following values (with w₁=1, w₂=4, w₃=16 to balance for the different scales of the number of infected vs. the number of recovered and dead).

For the computed “effective” per-day transmission and mortality rates, we also report the corresponding 90% confidence intervals instead of the more standard 95% CI because of the small size of the data. Under this scenario, the DAY-ZERO for the outbreak in Lombardy was found to be the 21st of January. The “effective” per-day transmission rate was found to be β = 0.779 (90% CI: 0.777-0.781) and the “effective” per-day mortality rate for the confirmed cases was found to be γ = 0.0173 (90% CI: 0.0154-0.0192). Based on the derived value of the “effective” per-day disease transmission rate, the basic reproduction number was found to be R₀= 4.04 (90% CI: 4.03-4.05).

Thus, based on the computed expression for R₀, the “effective” per-day transmission rate should have been below ∼0.19, implying that the average contacts per person should have dropped by at least ∼75% for the outbreak to fade out at the initial stage.

Using these estimated values for the epidemiological parameters, we ran the simulator from DAY-ZERO (21st of January) to March 8. On March 8, simulations resulted in the following numbers for the cumulative cases: ∑ ∆E(t) = 66, 016, ∑ ∆I(t) = 37, 691, ∑ ∆I_c(t) = 4, 232, R(t) = 12, 123, R_c(t) = 623, D(t) = 220. The reported cumulative numbers for the day of lockdown were ∑ ∆I_c(t) = 4, 189, R_c = 550, D(t) = 267. Figures (2),(3),(4) depict the simulation results based on the optimal estimates, starting from the 21st of January to the 8th of March. As shown, the predictions of the model are quite close to the reported number of confirmed cases for that period.

• Download figure
• Open in new tab

Figure 2:

Cumulative number of (confirmed) infected cases resulting from simulations from DAY-ZERO (January 21) until the 16th of March. The DAY-ZERO, β = 0.779, γ = 0.0173 were computed by solving the mixed-integer optimization problem for the period DAY-ZERO to March 8. The validation of the model was performed using the reported data of confirmed cases from March 9 to March 16 (shaded area) by taking a 90% reduction in the β due to the lockdown of March 8. Dots correspond to the reported data of confirmed cases.

• Download figure
• Open in new tab

Figure 3:

Cumulative number of (confirmed) recovered cases resulting from simulations from DAY-ZERO (January 21) until the 16th of March. The DAY-ZERO, β = 0.779, γ = 0.0173 were computed by solving the mixed-integer optimization problem for the period DAY-ZERO to March 8. The validation of the model was performed using the reported data of confirmed cases from March 9 to March 16 (shaded area) by taking a 90% reduction in the β due to the lckdown of March 8. Dots correspond to the reported data of confirmed cases.

• Download figure
• Open in new tab

Figure 4:

Cumulative number of (confirmed) deaths resulting from simulations from DAY-ZERO (January 21) until the 16th of March. The DAY-ZERO, β = 0.779, γ = 0.0173 were computed by solving the mixed-integer optimization problem for the period DAY-ZERO to March 8. The validation of the model was performed using the reported data of confirmed cases from March 9 to March 16 (shaded area) by taking a 90% reduction in the β due to the lockdown of March 8. Dots correspond to the reported data of confirmed cases.

Thus, according to the above results, on the 8th of March, the actual cumulative number of infected cases in the total population (taking into account the exposed cases to the virus) was of the order of 15 times more the confirmed cumulative number of infected cases.

To validate the model with respect to the reported data of confirmed cases from March 9 to March 16, we have set as β = 0.078, (i.e. taking a 90% reduction in the value of β found with the optimization procedure) and as initial conditions the values resulting from the simulation on March 8. Based on the above, the model resulted in the following numbers for the cumulative cases for March 16: ∑ ∆E(t) = 81, 619, ∑ ∆I(t) = 75, 149, ∑ ∆I_c(t) = 14, 317, R(t) = 40, 946, R_c(t) = 3, 158, D(t) = 1, 139 (the reported confirmed cumulative numbers on March 16 were ∑ ∆I_c(t) = 14, 649, R_c = 2, 368, D(t) = 1, 420). Thus, the model predicted fairly well the period from March 9 to March 16.

As discussed in the Methodology, we also attempted to forecast the evolution of the outbreak based on our analysis. To do so, we have considered a 90% reduction in the effective transmission rates starting on March 8, the day of lockdown not only of Lombardy, but of all Italy. The result of our forecast is depicted in Figure 5. As predicted by simulations, if the strict isolation measures continue to hold, the outbreak in Lombardy will fade out by the end of May, 2020.

• Download figure
• Open in new tab

Figure 5:

Estimated number of infected cases in the total population at each day resulting from simulations from March 8 to May 31 (the day of the lockdown of all Italy) using β = 0.779 − 90%0.779, γ = 0.0173 as computed by solving the mixed-integer optimization problem for the period DAY-ZERO (21st January) to March 8.

Discussion

The crucial questions about an outbreak is how, when (DAY-ZERO), why it started and when it will end. Answers to these important questions would add critical knowledge in 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. No doubt, 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 in 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.

This study aimed exactly at shedding more light into this problem, taking advantage of state-of-the-art tools of mathematical modelling and numerical analysis/optimization tools. To achieve this goal, we addressed a new compartmental SEIRD with two infectious compartments in order to bridge the gap between the number of reported cases and the actual number of cases in the total population. By following the proposed methodological framework, we found that the DAY-ZERO in Lombardy was the 21th of January, a date that precedes by one month the fate of the first confirmed case in Lombardy. Furthermore, our analysis revealed that the actual cumulative number of infected cases in the total population is around 15 times more the cumulative number of confirmed infected cases. Importantly, based on our simulations, we predict that the fade-out of the outbreak in Lombardy will be by the end of May, if the strict isolation measures continue to hold.

To this end, we would like to make a final comment with respect to the basic reproduction number R₀, the significance and meaning of which are very often misinterpreted and misused, thereby leading to erroneous conclusions. Here, we found an R₀ ∼ 4, which is similar to the values reported by many studies in China. For example, Zhao et al. estimated R₀ 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 [7]. Similar estimates, were obtained for R₀ by Imai et al. 2.6 (95% CI: 1.5-3.5) [4], Li et al. [21], Wu et al. 2.68 (95% CI: 2.47–2.86), as well as by Anastassopoulou et al. recently 3.1 (90% CI: 2.5-3.7) [6].

However, we would like to stress that R₀ 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, the demographics as well as, importantly, by the social behavior of the population (see for example the discussion in [22]). Thus, a value for R₀ that is found in a part of world (and even in a region of the same country) cannot be generalized as a global biological constant for other parts of 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.

We hope that the results of our analysis help to mitigate some of the severe consequences of the currently uncontrolled pandemic.

Funding

We did not receive any specific funding for this study.

Data Availability

The data used in this paper are given in the Supporting information.

Supporting information

All the relevant data used in this paper are publicy available and accessible at https://lab.gedidigital.it/gedi-visual/2020/coronavirus-i-contagi-in-italia/. In Table S1 are given the reported cumulative numbers from Febrary 21 to March 16. The data from February 21 to March 8 have been used for the calibration of the model parameters and the data from March 9 to March 16 have been used for the validation of the model.

S1 Table..