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

The modelling approach

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 [6, 14, 7]. 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 [15].

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 Figure 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 , 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. Our main assumption here is that only a fraction, say ∊ of the actual number of exposed cases E will experience mild to 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 [16, 17, 18]. Thus, for any practical means the compartment I_c reflects the reported confirmed infected cases.

• Download figure
• Open in new tab

Figure 1:

A schematic of the proposed compartmental SEIIRD model. The actual number of cases is unknown.

A fraction of the I_c cases that is given by the fatality ratio 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 [19]. 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 [19] 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:

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^-1) is the average per-day “effective” rate at which an exposed person becomes infectious,

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

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

• γ(d⁻¹) 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 ,

• ∊(d⁻¹) 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 [20]. With respect to the incubation period for SARS-CoV-2, a study in China [21] 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, 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 [22]. Based on the above, and on the fact that within the subset of confirmed cases the mild cases are the 81% [15], 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 [23], 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 [24]; for the case of Lombardy, the average serial number has been estimated to be 6.6 days [25]. In our model, the 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. [26]) 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 [27], 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.

Estimation of the DAY-ZERO of the outbreak, the scale of data uncertainty and the disease transmission

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 [28]) 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 [19]. 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 [5] 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₁, w₂, w₃ 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.

At this point we should note that the above optimization problem may in principle have more than one nearby optimal solutions, which may be attributed to the fact that the tuning of both DAY-ZERO and the transmission rate may in essence result in nearby values of the objective function. In order to quantify the above uncertainty in the optimization procedure, we created a grid of initial guesses within the intervals in which the optimal estimates were sought: for the DAY-ZERO (t₀) we used a step of 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, Exponential and Birnbaum-Saunders 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 [29] utilizing 6 INTEL XEON CPU X5650 cores at 2.66GHz.

Estimation of the basic reproduction number R₀ from the SEIIRD model

Here, we note that we provide an estimation of the basic reproduction number R₀ based on the estimation of the total number of (asymptomatic) infected cases in the population. Thus, it is expected that the estimated R₀ 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 Eq.(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 Eq.(2),(3),(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 [30] (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

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

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 (17) is also a sufficient condition for stability. Hence, the disease-free state is stable, if and only if, condition (17) is satisfied.

Note that in this necessary and sufficient condition (17), the fraction is the average infection time of the compartment I. Thus, the above expression reflects the basic reproduction number R0 which is qualitatively defined by . Hence, our model results in the following expression for the basic reproduction number:

Note that for ∊ = 0, the above expression simplifies to R₀ for the simple SIR model.

Forecasting

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 [31], 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 [32, 33, 34, 35]). 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 [36].

On the basis of the Google COVID-19 Community Mobility Report released on March 29 [31], 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 [36], 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 [37]. 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 [38, 39]. According to the Google COVID-19 Community Mobility Reports [31], from March 2021 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 [36], 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. For example, it has been shown that social distancing and cautiousness reduce the disease transmission rate by about 20% [34]. Thus, based on the above, it is reasonable to consider a (1-0.4)(1-0.2) (an average of 40% contribution of the reduction of the mobility and a 20% for the effect of social distancing) reduction in the effective transmission rate for the period March 8-March 19 as compared with the period February 21-March 8. For the period of March 20-21, based on the Google data [31], we considered a reduction of (1-0.65)(1-0.2) (as compared with the period February 21-March 8) reflecting the draconian measures taken then. Based on the above, we attempted a forecasting of the pandemic in Lombardy from March 8 to May 4, the first day of the relaxation of the strict isolation measures.