Limits of Compartmental Models and New Opportunities for Machine Learning: A Case Study to Forecast the Second Wave of COVID-19 Hospitalizations in Lombardy, Italy

: Compartmental models have long been used in epidemiological studies for predicting disease spread. However, a major issue when using compartmental mathematical models concerns the time-invariant formulation of hyper-parameters that prevent the model from following the evolution over time of the epidemiological phenomenon under investigation. In order to cope with this problem, the present work suggests an alternative hybrid approach based on Machine Learning that avoids recalculation of hyper-parameters and only uses an initial set. This study shows that the proposed hybrid approach makes it possible to correct the expected loss of accuracy observed in the compartmental model when the considered time horizon increases. As a case study, a basic compartmental model has been designed and tested to forecast COVID-19 hospitalizations during the ﬁrst and the second pandemic waves in Lombardy, Italy. The model is based on an extended formulation of the contact function that allows modelling of the trend of personal contacts throughout the reference period. Moreover, the scenario analysis proposed in this work can help policy-makers select the most appropriate containment measures to reduce hospitalizations and relieve pressure on the health system, but also to limit any negative impact on the economic and social systems.


Introduction
In early December 2019, the first COronaVIrus Disease-2019 (COVID- 19) case was identified in Wuhan, China [1], which started a pandemic that spread across the globe in just a few months [2]. COVID-19 is a disease caused by a new form of coronavirus named Severe Acute Respiratory Syndrome CoronaVirus type 2 (SARS-CoV-2) [1,3]. The most common symptoms are fever, fatigue, and dry cough, which can evolve, in severe cases, into pneumonia, severe acute respiratory syndrome, and death [4]. The disease has been spreading fast and more than 65 million people have been infected as of fifth December 2020, with more than 1.5 million deaths around the world [2]. In Italy, the first COVID-19 case was announced on 20 February 2020 in Codogno, Lombardy [5]. In a few days, the contagion started spreading throughout the country and, for this reason, the Italian government announced a national lockdown on 9 March 2020 to contain the virus [6]. The lockdown period ended on 4 May 2020, but the circulation of the virus has never really stopped in Italy, with few cases in Summer, but a rapid increase of the pandemic spread since September reaching a total of about 1.7 million cumulative cases on 6 December [7]. According to scientific literature, compartmental mathematical models have been widely used to emulate the behavior of infectious diseases and to map specific phenomenon-based characteristics within the model. Compartmental models were first introduced by Kermack and McKendrick [8] in 1991 with the Susceptible, Infected, and Recovered (SIR) formulation [9]. Several other studies using compartmental models have also been published, especially aimed at modelling COVID-19 spread in different countries by focusing on multiple variables. In particular, Chen and colleagues [10] developed a multilevel transmission network (called BHRP: bats, host, reservoir, people) based on the Susceptible, Exposed, Infected and Recovered (SEIR) paradigm in order to emulate COVID-19 transmission in China. Starting from the initial outbreak among bats, the method by Runge-Kutta [11] was adopted to estimate the model equations hyper-parameters based on Chinese official data. Moreover, an age-structured SEIR model was proposed by Prem and colleagues [12] to evaluate the effect of social distancing measures on the outcomes of the COVID-19 epidemic in Wuhan, China. Kucharski and colleagues [13] conducted a mathematical study based on a stochastic SEIR model, with the goal of modelling the early dynamics of COVID-19 pandemic transmission in Wuhan, also considering the impact of travelling on the virus spread outside Wuhan. A modified SEIR model was used by Yang and colleagues [14] to derive the COVID-19 epidemic curve in China. A Machine Learning approach, trained on the 2003 SARS data, was also used to predict the epidemic spread. An Extended SEIR model was proposed by Tang and colleagues [15,16] to predict the trend of COVID-19 hospitalizations in China, which included some timevariant hyper-parameters to implement a dynamic setting of the model. The same model was used for the Italian case [17] to forecast COVID-19 related variables in two Italian regions. Bollon and colleagues [18] used compartmental mathematical models to assess the role of containment measures in curbing the demand for hospitalization. Rǎdulescu and colleagues [19] adapted a traditional SEIR model based on four age groups to analyze the effect of containment measures strategies in USA. Moreover, Giordano and colleagues [20] developed a new mathematical model based on eight compartments to forecast COVID-19 trend in Italy. Ahmad and colleagues [21] proposed a fractional SEIR model to investigate the spread of the COVID-19 in Pakistan during the first pandemic wave using, also, the next generation matrix algorithm to compute the basic reproduction number R 0 . Kuzdeuov and colleagues [22] developed a COVID-19 simulator for Kazakhistan through a stochastic SEIR model. The simulator exploits a variety of data such as epidemiological, population demographics and mobility. Moreover, Gaglione and colleagues [23] used a Bayesian sequential and adaptive dynamic estimation method for computing the infection and recovery parameters, and to track and predict the epidemiological curve in Lombardy, Italy, and in USA. A SEIR model has also been exploited by Bherwani and colleagues [24] to emulate the spread of the COVID-19 epidemic in India. Battineni et al. [25] also developed a SEIR model explaining the infection growth across Italy and presenting epidemic rates after and before country lockdown. A study related to the COVID-19 second wave in Italy has been proposed by Chintalapudi et al. [26] evaluating the basic reproduction number.
Moreover, Zisad et al. [27] proposed an integrated Neural Network and SEIR Model to predict COVID-19 spread in Bangladesh. Furthermore, Alanazi et al. [28] developed a framework based on SIR model and Machine Learning for measuring and preventing the continued spread of COVID-19 in the Kingdom of Saudi Arabia (KSA).
Growth models have been also used to understand COVID-19 epidemic outbreaks by Tovissodé et al. [29], as well as in [30] to fit the number of cumulative confirmed cases for revealing the patterns and the rate of COVID-19 spread in different countries. A second wave of the COVID-19 pandemic has started in Europe since the end of Summer, which achieved a peak in November [31]. However, all the mathematical models based on the first wave of the pandemic were not able to provide an accurate forecast of the second COVID-19 pandemic curve, mainly due to substantial parametric changes from the first to the second wave. In general, a major issue when using compartmental mathematical models relates to the gradual loss of accuracy coming from the time-invariant formulation of hyper-parameters, which actually prevent the model from following the evolution over time of the epidemiological phenomenon under investigation. Such hyper-parameters are generally calculated on the basis of appropriate estimation algorithms (e.g., MCMC, Maximum Likelihood) that do not provide time-variant formulations though. For this reason, each hyper-parameter should be defined with a time-variant formulation that would allow the model to adapt to the different boundary conditions observed throughout the evolution of the investigated phenomenon. The problem is exacerbated in the case of complex compartmental models like the model in this work, based on nine compartments and 29 hyper-parameters. In order to cope with this problem, the present work suggests an alternative approach, based on Machine Learning, that avoids recalculation of hyperparameters and only uses an initial set. As the time frame increases, the expected loss of accuracy in the model is corrected through an approach based on machine learning. Indeed, in the present study, a Hybrid SEIR Model (from now on called HSM) framework, has been designed by coupling a Basic SEIR Model (from now on called BSM) and Machine Learning (ML), which led to a novel hybrid framework that takes advantage of both approaches. The developed HSM has been used here to produce accurate predictions about COVID-19 hospitalizations for the second epidemic wave, to simulate different scenarios based on different containment measures and to assess the impact of such measures on the epidemiological curve, using the HSM framework.

Population and Data Sources
To avoid model uncertainty due to inter-country and inter-regional variability, we restricted the application of the model to the Lombardy region, North Italy, which is one of the most affected regions in Italy and Europe, with more than 400,000 cases and 20,000 deaths in 10 million inhabitants on 6 December, with a cumulative incidence of about 4% [32]. The number of hospitalizations was chosen as the target variable because it reflects the number of severe cases and does not depend on the number of swabs, which has varied substantially over time in Italy. Keeping track of hospitalizations is strategic for policy-makers as it allows predicting the potential saturation of the hospital system, and when this might happen. Moreover, data concerning COVID-19 hospitalizations, deaths, isolated infected patients, infected patients, number of swabs and recovered patients in the period from 9 March to 6 January 2021 were collected from the Italian Civil Protection Department (ICPD) repository [33]. The average recovery period of hospitalized patients has been estimated to be equal to 14 days by the Regional Health Agency (ARS) of Tuscany, Italy [34]. This period corresponds to the worst case identified by ARS among the five trajectories related to hospitalization duration according to different severity scenarios. For the BSM modelling phase, other parameters regarding the COVID-19 epidemic spread have been collected from Tang and colleagues [15,16] as well as from Reno and colleagues [17].

The Basic SEIR Model
The Basic SEIR Model (BSM) is a closed-population compartmental model. Based on the SEIR model first introduced by Tang and colleagues [15,16] and subsequently adapted by Reno and colleagues [17] to the Italian case, BSM has been proposed by redefining some hyper-parameters related to the underlying differential equations. As shown in Figure 1, BSM takes into account nine compartments. some hyper-parameters related to the underlying differential equations. As shown in Figure 1, BSM takes into account nine compartments. Figure 1. The BSM architecture to forecast the hospitalizations trend in Lombardy (Italy). According to BSM, the population is divided into nine compartments: susceptible (S), exposed (E), infected (I), recovered (R), quarantined susceptible (Sq), asymptomatic (A), hospitalized (H), quarantined exposed (Eq), isolated infected (L). Interactions among compartments are depicted as directed segments, also reporting the related transition rates.
Each compartment is connected to one or more compartments through directed segments, and their interaction is regulated by the rate associated with each segment. Formally, BSM is defined by the following system of ordinary differential equations that allows modelling important epidemiological characteristics related to the COVID-19 spread in Lombardy: The solution of the system above is obtained by integrating differential equations over a time interval starting from a set of initial conditions, reported in Table A1 of Appendix A, and related to compartments. The time evolution of each compartment depends on a series of hyper-parameters whose values have been estimated based on the available data and information about the first pandemic wave. The full set of hyper-parameters is reported in Table A2 of Appendix A. Hyper-parameters play an important role in modelling region-specific outbreak dynamics, since the simulated trends strongly depend on the parameter's choice. Most hyper-parameters values are those proposed by Tang and colleagues [16] as well as by Reno and colleagues [17] because some parallelism seems to exist between the outbreak dynamics observed in Lombardy (Italy) and in Wuhan (China), even though with different population characteristics. Further hyper-parameters have been added whereas others have been redefined in order to improve the BSM forecasting performance in Lombardy. As the virus is mainly transmitted through contact, the modeling of an appropriate "contact function" was necessary to allow for the impact of containment measures on the gradual trend in the number of contacts. Information about the contact rate was not available at the time of writing. Hence, an intuitive approach has been applied using mobility data provided by Apple [35], as better explained in Appendix A. Therefore, the contact rate ( ) has been defined as: Figure 1. The BSM architecture to forecast the hospitalizations trend in Lombardy (Italy). According to BSM, the population is divided into nine compartments: susceptible (S), exposed (E), infected (I), recovered (R), quarantined susceptible (S q ), asymptomatic (A), hospitalized (H), quarantined exposed (E q ), isolated infected (L). Interactions among compartments are depicted as directed segments, also reporting the related transition rates.
Each compartment is connected to one or more compartments through directed segments, and their interaction is regulated by the rate associated with each segment. Formally, BSM is defined by the following system of ordinary differential equations that allows modelling important epidemiological characteristics related to the COVID-19 spread in Lombardy: The solution of the system above is obtained by integrating differential equations over a time interval starting from a set of initial conditions, reported in Table A1 of Appendix A, and related to compartments. The time evolution of each compartment depends on a series of hyper-parameters whose values have been estimated based on the available data and information about the first pandemic wave. The full set of hyper-parameters is reported in Table A2 of Appendix A. Hyper-parameters play an important role in modelling regionspecific outbreak dynamics, since the simulated trends strongly depend on the parameter's choice. Most hyper-parameters values are those proposed by Tang and colleagues [16] as well as by Reno and colleagues [17] because some parallelism seems to exist between the outbreak dynamics observed in Lombardy (Italy) and in Wuhan (China), even though with different population characteristics. Further hyper-parameters have been added whereas others have been redefined in order to improve the BSM forecasting performance in Lombardy. As the virus is mainly transmitted through contact, the modeling of an appropriate "contact function" was necessary to allow for the impact of containment measures on the gradual trend in the number of contacts. Information about the contact rate was not available at the time of writing. Hence, an intuitive approach has been applied using mobility data provided by Apple [35], as better explained in Appendix A. Therefore, the contact rate c(t) has been defined as: where: and with The corresponding hyper-parameters values are defined in Table A2 in Appendix A. As shown in Equation (2), the contact function c(t) has been designed as a linear composition of two terms c w 1 (t) and c w 2 (t), which model the trend of personal contacts during the first and the second wave respectively. The contact function c(t) estimates the trend of contacts between people, while also mapping the adoption of various containment measures aimed at curbing the epidemic spread (i.e., lockdown). Specifically, c w 1 (t) is expressed as a function of two time-varying terms: the first term c 1 (t), shown in Equation (4), has been derived by Tang and colleagues [16] to model the exponentially decreasing rate of contacts upon the adoption of containment measures for the first pandemic wave (9 March-18 May 2020) in the present analysis. The second term c 2 (t), shown in Equation (5), has been used to model the increased rate of contacts after gradually lifting containment measures at the end of the first wave lockdown. Moreover, as shown in Equation (6), c w 2 (t) is expressed as a function of the term c 3 (t), defined in Equation (7) and relates to the decreasing trend of contacts that was observed after the new containment measures put in place to control the pandemic spread during the second wave (started on 16 October 2020). In Equations (3) and (6) suitable sigmoid functions have been used for simulating the gradual decrease or increase in contacts when entering or exiting lockdown, respectively. Furthermore, the c e parameter in the c 3 (t) term of the c w 2 (t) contact function allows modelling the asymptotic presumed average rate of personal contacts after the second wave lockdown imposed in Lombardy on 6 November 2020. Hence, c e allows estimating different hospitalization scenarios based on the pandemic evolution in Lombardy. The solution of the BSM mathematical model strongly depends on the initial conditions (given on 9 March 2020) and on the compartment hyper-parameters, that have been set to model hospitalizations during the first pandemic wave. Moreover, it has been possible to achieve a better setting of BSM through the fine-tuning of hyper-parameters, by means of a trial and error procedure that compares BSM predictions with real hospitalization data. Specifically, the BSM best setting has been achieved by setting the exponential decreasing rate of diagnose rate r 2 and the recovery rate of asymptomatic infected individuals γ A to 0.1 and 0.17529, respectively. However, as the simulation horizon increases over time, BSM accuracy progressively decreases since BSM settings (hyper-parameters and initial conditions) should be constantly updated according to the new time window under consideration. BSM is therefore not properly suited for estimating hospitalizations during the second COVID-19 pandemic wave. For this reason, the hybrid framework (HSM) has been developed.

The Hybrid SEIR Model
The use of BSM to accurately forecast hospitalizations during the second pandemic wave, with the first containment measures that were introduced on 13 October 2020 [36] would have meant redefining the whole set of initial conditions and hyper-parameters. Therefore, a hybrid framework has been designed to develop the HSM. The aforementioned framework derives from using the two different predictive approaches obtained through both the BSM mathematical model and the Machine Learning (ML) module. The latter takes care of fixing mistakes in the long-term forecasts produced by the BSM and it relies on polynomial models that have been properly trained on hospitalization data in Lombardy from 15 October to 5 November 2020 (date when the second wave containment measures became effective) [37]. Once trained, the ML module is able to make long-term predictions (starting from 6 November 2020). Polynomial models of different orders have been considered for assessing their accuracy against real hospitalization data. The best fitting has been achieved by means of second order polynomial models.

The Machine Learning Module
The Hybrid SEIR Model (HSM) framework, introduced in the present study, exploits the features of compartmental SEIR models and Machine Learning (ML) in order to provide accurate forecasts of hospitalizations in Lombardy during the second wave lockdown. Hence, HSM relies on a compartmental model (BSM), that is unable to provide long-term accurate estimations, and on a ML module that aims to learn and correct the estimation error produced by BSM. Figure 2 shows the HSM training and inference phases. Equation (8) was used in the supervised training phase to define polynomial hypotheses h d θ (d) , x of different degree d, that have been trained with the hospitalizations observed in Lombardy from 15 October to 5 November 2020 as target and with BSM estimations for the same period (backward lagged by 1 day) as input.
where θ 0 , . . . , θ d are unknown regression coefficients which have been estimated using the Normal Equation algorithm. Once trained, each polynomial hypothesis allows the mapping between the BSM inaccurate estimations and the corresponding observed hospitalizations, thus providing improved forecasts. Moreover, the HSM training phase provides as an outcome the trained polynomial hypotheses to be used during the inference phase. In the inference phase, the already trained polynomial models are exploited to correct BSM estimations and to forecast the long-term trend of future hospitalizations, which has not been observed yet. Indeed, inaccurate BSM estimates, starting from 6 November, have been backward lagged by 1 day in this phase, and used as input of each polynomial hypothesis previously trained to produce hospitalizations forecasts starting from 6 November, when the second wave lockdown became effective [37]. The different polynomial hypotheses have been assessed in the model selection block, shown in Figure 2, to select the optimal model which produces the most accurate hospitalizations forecasts. The Root Mean Squared Error (RMSE), defined by (9), has been selected as the evaluation metric for the error produced by each polynomial hypothesis in the prediction of hospitalizations from 6 November to 8 December 2020.
whereŷ i is the i-th prediction, y i is the i-th observed sample and n is the number of samples. The Root Mean Squared Error (RMSE), defined by (9), has been selected as the evaluation metric for the error produced by each polynomial hypothesis in the prediction of hospitalizations from 6 November to 8 December 2020.
where is the i-th prediction, is the i-th observed sample and is the number of samples.

What-If Analysis: Scenarios during the First and the Second Lockdown
The what-if analysis allows the assessment of different scenarios related to the evolution of the pandemic, by taking advantage of time variability in the ( ) contact function. Two scenario analyses have been carried out with regard to the first and the second wave of the pandemic, respectively. The first analysis was aimed at estimating the trend of hospitalizations in Lombardy under various scenarios with a different lockdown duration during the first wave of the pandemic. In particular, this was achieved by running BSM multiple times with different values of the variable, which represents the number of days of simulated lockdown.
The second analysis focused instead on assessing the effect of potential containment measures to keep the rate of personal contacts unaltered during the second wave lockdown. The HSM framework allowed minimizing the long-term error made by BSM in estimating second wave lockdown hospitalizations. Different values, that define the minimum average rate of contacts during the second wave lockdown (started on 6 November 2020), have been assessed through several sensitivity analyses. For each sensitivity analysis all the hyper-parameters have been left unchanged except for in order to assess its impact on future hospitalization trend scenarios, comparing HSM forecasts with

What-If Analysis: Scenarios during the First and the Second Lockdown
The what-if analysis allows the assessment of different scenarios related to the evolution of the pandemic, by taking advantage of time variability in the c(t) contact function. Two scenario analyses have been carried out with regard to the first and the second wave of the pandemic, respectively. The first analysis was aimed at estimating the trend of hospitalizations in Lombardy under various scenarios with a different lockdown duration during the first wave of the pandemic. In particular, this was achieved by running BSM multiple times with different values of the L 1 variable, which represents the number of days of simulated lockdown.
The second analysis focused instead on assessing the effect of potential containment measures to keep the rate of personal contacts unaltered during the second wave lockdown. The HSM framework allowed minimizing the long-term error made by BSM in estimating second wave lockdown hospitalizations. Different c e values, that define the minimum average rate of contacts during the second wave lockdown (started on 6 November 2020), have been assessed through several sensitivity analyses. For each sensitivity analysis all the hyper-parameters have been left unchanged except for c e in order to assess its impact on future hospitalization trend scenarios, comparing HSM forecasts with the actual available hospitalizations data. Several c e values have been considered for modelling potential scenarios during the second epidemic wave, when containment measures were imposed on the population. In particular, the scenario with the highest value of contacts equal to 6 has been considered first, as strongly suggested by the decree of the Italian Prime Minister on 13 October 2020 [36]. Additionally, this value was also imposed by the UK government on 14 September 2020 [38] to slow down the pandemic spread by limiting the number of personal contacts. Moreover, further c e parameter values have been considered for the evaluation of more restrictive scenarios (average number of contacts equal to 3) and of other less stringent scenarios (average number of contacts equal to 10 and 12). For each scenario, the average percentage error ( % ) between HSM hospitalization predictions and actual hospitalizations data has also been computed.

The Daily Reproduction Number Estimation
The daily reproduction number R d (t) has been estimated by using the next-generation matrix method [39] (see Appendix B): Since R d (t) depends on the contact function c(t), by modifying the variable L 1 in Equation (3), different daily reproduction number trends could be assessed in Lombardy during the first wave of the pandemic, with respect to the critical threshold R d (t) = 1, (i.e., when the disease is endemic). When computed at the beginning of the pandemic (t = 0), R d (0) represents the basic reproduction number (known as R 0 in literature), usually employed for evaluating the potential infectivity of a disease. Herd immunity (i.e., the minimum threshold allowing a population to become self-immune to a virus) has also been evaluated based on 9 March 2020 data in Lombardy, using the following equation proposed by Gorkin [40]: where i represents the portion of the population that needs to be immunized in order to reach the herd immunity threshold. Figure 3 reports the results of the simulations performed using BSM. For both panels, the x-axis reports the number of days from 9 March to 9 July 2020, whereas the y-axis reports the absolute number of hospitalizations. Real hospitalized people are depicted by a red line, whereas BSM predictions are shown in green through a dashed line. BSM hyperparameters have been first initialized to the values reported in Appendix A (Table A2), and the corresponding predictions are shown in Figure 3a. BSM hyper-parameters have then been further tuned for the best setting, leading to a better fitting with respect to the actual hospitalizations trend observed (Figure 3b).  Moreover, the BSM best setting has been used as the starting point to perform a what-if analysis on the first pandemic wave, by taking advantage of the time-variant behavior of the contact function c(t). The what-if analysis has been performed for assessing the effect of different lockdown durations on the hospitalizations trend ( Figure 4a) and on the corresponding daily reproduction number R d (t), leading to several scenarios (Figure 4b).

Results
The daily reproduction number was further evaluated on 9 March 2020 (t = 0) upon the adoption of containment measures for the first COVID-19 wave. The results reported in Figure 4b, show a value of R d (0) = R 0 = 2.625. This value has been used for estimating the immunity threshold for Lombardy using Equation (11):     Figure 5 shows the results obtained by means of the second order HSM hybrid model to predict hospitalization evolution during the second pandemic wave. As described in Section 2, the c e parameter in the c 3 (t) term of the contact function allows modelling the asymptotic presumed rate of personal contacts after the second wave lockdown imposed in Lombardy on 6 November 2020. Four different average rates of contact equal to 3, 6, 10, and 12 have been considered to assess the impact of different containment measures on the hospitalizations trend. Panel (a) reports the estimates of contacts evolution from 9 March 2020 to 31 December 2020, under four different scenarios based on the c e parameter. Panel (b) shows the results of the BSM prediction during the same observation period under the four different scenarios based on the c e parameter. An increasing loss of accuracy can be observed in the BSM with wider time windows, starting from 17 June 2020, producing an overestimate of hospitalizations in the period from the end of the first wave lockdown to the introduction of new containment measures during the second wave (15 October 2020). Panel (c) shows the results using HSM and four different rates of personal contacts during the lockdown period imposed in Lombardy on 6 November 2020: the lower the contact rate, and the earlier and lower the peak of the hospitalization curve. The comparison between BSM and HSM applications shows an improvement of more than 94% of HSM with respect to BSM ( Table 1). The average percentage error ( % ) was also calculated for each contact scenario and, as it turned out, the c e = 10 scenario recorded the lowest error (2.39%) among all scenarios considered in the analysis (Figure 5c). To emphasize such result, Figure 6 reports the hospitalizations trend as predicted by HSM compared with the actual number of hospitalizations in Lombardy, from 6 November to 8 December 2020. As shown, HSM predictions exhibited a very close match with actual hospitalizations data, falling into the confidence interval of prediction (CI: 95%) throughout the observed period.  To emphasize such result, Figure 6 reports the hospitalizations trend as predicted by HSM compared with the actual number of hospitalizations in Lombardy, from 6 November to 8 December 2020. As shown, HSM predictions exhibited a very close match with actual hospitalizations data, falling into the confidence interval of prediction (CI: 95%) throughout the observed period.

Discussion
The present work suggests a hybrid approach based on Machine Learning that avoids recalculation of compartmental models' hyper-parameters. This study shows that the proposed hybrid approach makes it possible to correct the expected loss of accuracy observed in the compartmental model when the considered time horizon increases.
Various mathematical models have been applied to obtain the most accurate predictions of COVID-19 hospitalizations in the first and second pandemic waves, according to social contact restrictions for lockdown, using the Lombardy region as a case study. As regards the first pandemic wave, from March to May 2020, a better setting of the classical BSM has been found through the fine-tuning of hyper-parameters, by means of a trial and error procedure that compares BSM predictions with real hospitalization data. This procedure allowed an improvement in hospitalizations predictions by about 26% (Figure 3b), with respect to forecasts achieved by the initial BSM version (Figure 3a). However, for the second wave, from October to December 2020, an HSM framework based on the classical BSM and a ML module, which relies on polynomial models that have been properly trained on real data, provided more accurate long-term predictions of COVID-19 hospitalizations than the traditional BSM. The application of both modelling approaches confirmed that containment measures, such as the lockdown and consequently contact reduction through social distancing, can heavily contribute to control the pandemic spread, depending on the magnitude of reduction of individual contacts and the length of the period: the lower the contact number and the longer the lockdown period, the lower the peak of hospitalizations and the farther the point the peak is reached (Figure 4a). The analysis of the different scenarios obtained through BSM may well support the selection of the most

Discussion
The present work suggests a hybrid approach based on Machine Learning that avoids recalculation of compartmental models' hyper-parameters. This study shows that the proposed hybrid approach makes it possible to correct the expected loss of accuracy observed in the compartmental model when the considered time horizon increases.
Various mathematical models have been applied to obtain the most accurate predictions of COVID-19 hospitalizations in the first and second pandemic waves, according to social contact restrictions for lockdown, using the Lombardy region as a case study. As regards the first pandemic wave, from March to May 2020, a better setting of the classical BSM has been found through the fine-tuning of hyper-parameters, by means of a trial and error procedure that compares BSM predictions with real hospitalization data. This procedure allowed an improvement in hospitalizations predictions by about 26% (Figure 3b), with respect to forecasts achieved by the initial BSM version (Figure 3a). However, for the second wave, from October to December 2020, an HSM framework based on the classical BSM and a ML module, which relies on polynomial models that have been properly trained on real data, provided more accurate long-term predictions of COVID-19 hospitalizations than the traditional BSM. The application of both modelling approaches confirmed that containment measures, such as the lockdown and consequently contact reduction through social distancing, can heavily contribute to control the pandemic spread, depending on the magnitude of reduction of individual contacts and the length of the period: the lower the contact number and the longer the lockdown period, the lower the peak of hospitalizations and the farther the point the peak is reached (Figure 4a). The analysis of the different scenarios obtained through BSM may well support the selection of the most appropriate lockdown duration to reduce the number of hospitalized patients and relieve pressure or stress from the health system, while also limiting any negative impact on the economic and social systems. Results actually show that a 60-, 90-or 120-day lockdown would produce a very similar trend in terms of hospitalizations. The results achieved by the BSM best setting confirmed also the impact of lockdown duration on the daily reproduction number when it is above the critical threshold, meaning that the disease is spreading. Indeed, as reported in Figure 4b, a shorter lockdown causes R d (t) to persist above the critical threshold for more days, as compared with a longer lockdown. Furthermore, for each considered scenario, the results reported in Figure 4 show a dependence between the day in which R d (t) achieves the endemic critical threshold (R d (t) = 1) and the day when the hospitalization peak is reached. The time interval between these two days might be related to the recovery period, estimated as 14 days according to Italian figures [34] and to the incubation period. Finally, the herd-immunity threshold of 61.9% is an optimistic result, since R 0 was computed on 9 March 2020 when some containment measures had already been introduced in Lombardy to reduce the COVID-19 spread. The adopted SEIR approach proved to properly model the trend of hospitalizations in Lombardy during the first pandemic wave. Unfortunately, due to the strong dependence of the BSM mathematical model on the initial conditions and on the hyper-parameters values that date back to 9 March 2020, an increasing loss of accuracy can be observed in the BSM with wider time window, as shown in Figure 5b, particularly, from the end of the first wave lockdown to the introduction of new containment measures during the second wave (15 October 2020). For this reason, the HSM framework has been applied. The appropriate variation of the c e parameter in the contact function that defines the asymptotic average value of contacts throughout the entire second wave, made it possible to outline four scenarios based on a different rate of personal contacts during the lockdown period imposed in Lombardy on 6 November 2020 [37]. Modelling the gradual decrease in the rate of personal contacts throughout the entire lockdown required using sigmoid functions in the contact function. Therefore, the results of applying the HSM framework depended on the contact function values during the transition period, that started on 15 October 2020, with the first Italian Prime Minister's Decree [36]. The HSM framework proved very close to the hospitalization trend observed in Lombardy during the second wave of the pandemic (see Figure 5c), as it allowed to reduce the estimation error of more than 94% with respect to BSM in the hospitalization trend from 15 October to 5 November 2020, for all the considered scenarios. For the c e = 10 scenario, the hospitalizations trend predicted by HSM is very close to the real hospitalizations trend in Lombardy, from 6 November to 8 December 2020, with an average percentage error of 2.39%, as shown in Figure 6. Indeed, all the proposed scenario analyses depend on the choice of the contact function. The modelling procedures that have been proposed to describe the trend of personal contacts have proved appropriate for the estimation of different plausible scenarios about the trend of hospitalizations in Lombardy, both in the first and the second phase. In addition, the different scenarios related to the trend of contacts during the second wave suggest that the adoption of stricter containment measures with respect to the allowed rate of personal contacts might promote a more moderate trend of hospitalizations, with a lower peak and earlier in time. Under the most plausible scenario that sets the rate of allowed personal contacts to ten throughout the entire lockdown period in the second wave, predictions show a peak at around 9300 hospitalized patients, followed by a rapidly decreasing trend that should end on 1 January 2021. The main strength of the study lies in the development of a novel hybrid framework based on ML to reduce the errors that are usually introduced by mathematical compartmental models in the long term. The analysis was also stretched out to include the period from 8 December 2020 to 6 January 2021 by forecasting the contacts trend during the considered period. The results of the extended analysis are reported in Figure 7. As clearly expected, HSM was not able to produce accurate forecasts in the period of interest ( % = 30.84%), even though forecasts fit in the 95% confidence interval. The reason for this depends on a number of factors related to the period of interest, including the quick alternation of different containment measures. As a matter of fact, the Prime Minister's Decree of 3 November 2020 [37] imposed new restrictions on the country based on three risk scenarios related to increased critical levels (yellow, orange and red) in the various regions. The Decree provides for the assessment of restrictions based on the index value observed in each Italian region. The definition of the critical level per region will take place on a weekly basis and it will last 15 days. Each region shall then comply with the assigned critical level measures for at least two weeks. For this reason, the Italian situation is extremely dynamic and difficult to model. In addition, the vaccination campaign against COVID-19 was launched in Europe on 27 December 2020. Proper modelling should also account for the impact that vaccines are supposed to have on the epidemic curve. However, the developed framework is extremely general and although the focus of the present work was on the hospitalized people compartment, the framework could also have been efficiently used for other compartments of the SEIR mathematical model, thus enabling the accurate prediction of other COVID-19 variables as well. The modelling proposed for the contact function allows the assessment of different, more or less stringent containment measures based on a variable average number of contacts. Such measures help outline different hospitalization trend scenarios that are useful to monitor their impact on the national health system. For these reasons, the framework represents a valid and accurate decision support tool for policy-makers. Yet, as already said, the use of compartmental models strongly depends on the initial conditions and on a set of hyper-parameters that govern the evolution of the various compartments over time. The use of a model based on different time scenarios (i.e., phase two and any further phase of the pandemic) requires the definition of time variant hyper-parameters, for the model to adapt to changes in working conditions also in the long term. Clearly, in the case of complex compartmental models such as As clearly expected, HSM was not able to produce accurate forecasts in the period of interest ( % = 30.84%), even though forecasts fit in the 95% confidence interval. The reason for this depends on a number of factors related to the period of interest, including the quick alternation of different containment measures. As a matter of fact, the Prime Minister's Decree of 3 November 2020 [37] imposed new restrictions on the country based on three risk scenarios related to increased critical levels (yellow, orange and red) in the various regions. The Decree provides for the assessment of restrictions based on the R t index value observed in each Italian region. The definition of the critical level per region will take place on a weekly basis and it will last 15 days. Each region shall then comply with the assigned critical level measures for at least two weeks. For this reason, the Italian situation is extremely dynamic and difficult to model. In addition, the vaccination campaign against COVID-19 was launched in Europe on 27 December 2020. Proper modelling should also account for the impact that vaccines are supposed to have on the epidemic curve. However, the developed framework is extremely general and although the focus of the present work was on the hospitalized people compartment, the framework could also have been efficiently used for other compartments of the SEIR mathematical model, thus enabling the accurate prediction of other COVID-19 variables as well. The modelling proposed for the contact function allows the assessment of different, more or less stringent containment measures based on a variable average number of contacts. Such measures help outline different hospitalization trend scenarios that are useful to monitor their impact on the national health system. For these reasons, the framework represents a valid and accurate decision support tool for policy-makers. Yet, as already said, the use of compartmental models strongly depends on the initial conditions and on a set of hyper-parameters that govern the evolution of the various compartments over time. The use of a model based on different time scenarios (i.e., phase two and any further phase of the pandemic) requires the definition of time variant hyper-parameters, for the model to adapt to changes in working conditions also in the long term. Clearly, in the case of complex compartmental models such as the BSM model with nine compartments used in the present study, redefining all hyper-parameters would mean having extremely finegrained information and epidemiological point data available both in space and in time, along with using the appropriate estimation algorithms, as for instance those based on Maximum Likelihood. Consequently, the impossibility to redefine all hyper-parameters for very long-time frames causes a loss of accuracy in long-term predictions. To overcome this limitation, the authors exploited the Machine Learning. Moreover, as the contact function has a strong impact on the model predictions, the best definition of such function would require knowing the number of personal contacts at a regional level for the entire time frame under observation. The lack of such information led the authors to the assumption of a contact trend that is governed by the mobility trend provided by Apple. The choice of Lombardy as the case study for the proposed modelling approach can be a possible limit of the study. However, Lombardy has a larger population than various European countries, was highly affected by COVID-19 in both first and second wave and has large industrial and commercial activities determining a wide network of social contacts. However, the Italian regions have different public health systems and adopted nonidentical approaches for contrasting the COVID-19 pandemic, making a prediction model not suitable for any region, as somewhat different parameters should be included in the model. In conclusion, various mathematical forecasting models for the COVID-19 pandemic have been tried so far, but the choice of the best model depends on various assumptions and parameters which are valid in a specific context. Among them, the average number of daily contacts is of the utmost relevance to assess the effectiveness of confinement measures, due to the high transmission of SARS-CoV-2 through human contacts. However, the duration and strictness of containment measures have a substantial impact on the socio-economic fabric and the health system, therefore the estimate of different scenarios of the pandemic spread according to daily contacts and other parameters can be a valuable tool for choosing the most appropriate compromise solution to adopt. Data Availability Statement: Publicly available datasets were used in this study. As reported throughout the manuscript, these datasets can be found at the following links: https://github.com/ pcm-dpc/COVID-19 (accessed on 21 July 2021); https://covid19.apple.com/mobility (accessed on 21 July 2021); http://dati.istat.it/?lang=en (accessed on 21 July 2021).

Acknowledgments:
The authors would like to acknowledge Antonio Aloisio for his editing and proofreading work on this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Compartmental Model Specification.  Quarantined rate of exposed individuals 1.887 × 10 −7 [17] σ Transition rate of individuals exposed to the infected class 1/14 [43] λ Rate at which the quarantined uninfected contacts were released into the wider community 1/14 [16] ρ Probability of symptoms among infected individuals 0.86834 [16] δ I0 Initial transition rate of symptomatic infected individuals to the quarantined infected class 0.3266 [16]   BSM refers to the parameters that have been updated through the "trial and error" approach illustrated in the manuscript.
Formally, the Basic SEIR Model (BSM) used in our analysis is defined by the system of ordinary differential equations described by (1) in the manuscript and also used by Reno and colleagues [17]. The solution to the differential equation system is obtained by integrating differential equations over a time interval, starting from a set of initial conditions related to compartments, as reported in Table A1. The time evolution of each compartment is governed by a series of hyper-parameters whose values have been estimated based on the available data and information about the first pandemic wave. The full set of hyper-parameters is reported in Table A2. BSM is based on the SEIR model that was first introduced by Tang and colleagues [15,16] and then adapted to the Italian case by Reno and colleagues [17]. BSM has been designed by redefining some differential equations and the related hyper-parameters to achieve a better modelling of the hospitalizations trend in Lombardy. In particular, the following contact function proposed by Reno and colleagues [17], where: has been entirely reconstructed, as shown in Equations (2)-(7) in the manuscript. As discussed in the manuscript, the contact function c(t) has been designed as a linear composition of two terms c w 1 (t) and c w 2 (t) as shown by (2). Specifically, c w 1 (t) models the trend of personal contacts during the first pandemic wave and is expressed as a function of two time-varying terms c 1 (t) and c 2 (t) as reported in (3). The first term c 1 (t), defined by (4), has been designed to model the exponentially decreasing number of contacts upon the adoption of containment measures for the first pandemic wave (9 March-18 May 2020). The second term c 2 (t), has been used to model the increased number of contacts after gradually lifting containment measures at the end of the first wave lockdown. Moreover, c w 2 (t) models the trend of personal contacts during the second pandemic wave and is expressed as a function of the time-varying term c 3 (t) as reported in (6). The term c 3 (t), defined by (7), has been added in the contact function to model the decreasing trend of contacts that was observed after the new containment measures put in place to control the pandemic spread during the second wave (16 October 2020 up to 6 January 2021). When formulating the contact function used by Reno and colleagues [17], both c 1 (t) and c 2 (t) terms can be left out through step function θ(t), which allows an instant modelling of the trend of contacts during the transition phase and also upon entering and exiting lockdown. On the contrary, to achieve a more realistic modelling of the trend of contacts following the adoption of containment measures both in the first and in the second lockdown, the present analysis has selected appropriate sigmoid functions that enable the representation of smoothed transitions. It should also be noted that in our formulation, transitions entering and exiting lockdown are only managed by sigmoid functions, thus leaving info about the duration of the first and the second lockdown out of contact terms c 2 (t) and c 3 (t), respectively. Finally, the c e parameter in the c 3 (t) term of the contact function allows modelling the asymptotic presumed number of personal contacts after the second wave lockdown imposed in Lombardy on 6 November 2020. Hence, c e allows estimating different hospitalization scenarios based on the pandemic evolution in Lombardy. The contact function is also parameterized by the c 0 , c b and c f rates. In the original model by Tang and colleagues [16], c 0 has been defined as equal to 14.781 referring to Wuhan region in China. In our use case, this value should have been rescaled to c 0 = 13.469, according to (A5), to fit the Lombardy context with a population of about 10.097 million people [41], where Wuhan population has been estimated at about 11.081 million individuals [44]. However, the present analysis considered a lower value of c 0 = 10 to include the preliminary containment measures that had already been introduced in Lombardy on 9 March 2020. This value has also produced the best predictions with respect to the observed hospitalizations. The parameter c b represents the lowest contact rate achieved during the lockdown period, that has been set to three according to the Italian Scientific Technical Committee report [42]. This document reports the average number of contacts per age group, considering a situation in which no containment measures are enforced on the population. As mentioned above, the lockdown was imposed on the entire Italian population from 9 March to 4 May 2020. It is safe to consider that the contact rate was only limited to home contacts in this period. Therefore, the value for c b has been obtained by calculating an average of the values for each age group. A broader discussion is needed for the c f parameter used to model the increasing trend of personal contacts between the two waves of the pandemic. Information about this parameter was not available at the time of writing, hence an intuitive approach has been applied using mobility data provided by Apple [35]. In Figure A1, the Apple mobility trends about "driving", "walking" and "transit" are shown with respect to a baseline representing a reference mobility value registered on 13 January 2020 then before the adoption of any containment measures. After the end of the first wave lockdown (4 May 2020), an increasing trend can be observed with a peak higher than 100% (in the worst case) with respect to the baseline. Hence, the authors of the analysis assumed a similar trend for personal contacts, imposing a value above 100% with respect to the reference value c 0 that was previously estimated to model contacts before the first lockdown. Moreover, the recovery rates of hospitalized people (γ H = 1/ recovery period H ) and asymptomatic infected individuals (γ A = 1/ recovery period A ) have been redefined to adapt the model to the Italian case. The recovery period for hospitalized patients and asymptomatic infected individuals lasts 13 and 10 days, as reported in [34,45] respectively. Therefore, γ H = 0.0769 and γ A = 0.1 have been considered. Furthermore, δ I (t) (the diagnosis rate), which strongly depends on the resources available for detecting new cases, has been defined as in the work by Tang and colleagues, [16]: where, as reported in Table A2, δ I0 is the diagnosis rate at the initial time (t = 0) with δ I (0) = δ I0 , δ I f is the fastest diagnose rate with lim t → ∞ δ I (t) = δ I f , and r 2 is the exponential decreasing rate. Figure A1. Apple mobility trends and the corresponding presumed contacts. Mobility trends provided by Apple, expressed as percentage variation with respect to the baseline in Italy between 13 January and 9 November 2020, are reported by straight lines with values on the y-axis, left side. On the y-axis, right side, the presumed contacts have been reported according to the mobility trend.

Evaluation of the Effective Reproduction Number
The effective reproduction number ( ) is one of the most important epidemiological indices to control the behavior of an epidemic disease. It corresponds to the average number of new infections per infectious case in a population made up of both susceptible and unsusceptible hosts in a generic time instant t. If 1, the number of cases will increase, such as at the start of an epidemic [46]. The present study used BSM to estimate hospitalization in Lombardy from 9 March to 15 October 2020.
( ) has been evaluated in each simulated day, leading to the definition of the effective daily reproduction number ( ). ( ) has been evaluated by Tang and colleagues [15,16], using the next generation matrix algorithm [39], and its expression is reported in (A7).
However, following the approach proposed by Reno et al. in [17], BSM adds a new compartment related to the isolated infected individuals (L) with respect to the SEIR model used by Tang and colleagues [15,16]. Thus, a new closed formula for ( ) has been derived by applying the next generation matrix algorithm to BSM. Each step of the algorithm is presented as follows.
Let be the vector of infected compartments, such as exposed (E), infected (I), asymptomatic (A), hospitalized (H), quarantined exposed (Eq), isolated infected (L). Let be the vector of uninfected compartments, such as susceptible (S), quarantined susceptible (Sq), and recovered (R). Let ℱ be the vector of new infection rates (transitions from Y to X) and the vector of all other rates (not new infections). Figure A1. Apple mobility trends and the corresponding presumed contacts. Mobility trends provided by Apple, expressed as percentage variation with respect to the baseline in Italy between 13 January and 9 November 2020, are reported by straight lines with values on the y-axis, left side. On the y-axis, right side, the presumed contacts have been reported according to the mobility trend.

Evaluation of the Effective Reproduction Number
The effective reproduction number R e (t) is one of the most important epidemiological indices to control the behavior of an epidemic disease. It corresponds to the average number of new infections per infectious case in a population made up of both susceptible and unsusceptible hosts in a generic time instant t. If R e > 1, the number of cases will increase, such as at the start of an epidemic [46].
The present study used BSM to estimate hospitalization in Lombardy from 9 March to 15 October 2020. R e (t) has been evaluated in each simulated day, leading to the definition of the effective daily reproduction number R d (t). R d (t) has been evaluated by Tang and colleagues [15,16], using the next generation matrix algorithm [39], and its expression is reported in (A7).
However, following the approach proposed by Reno et al. in [17], BSM adds a new compartment related to the isolated infected individuals (L) with respect to the SEIR model used by Tang and colleagues [15,16]. Thus, a new closed formula for R d (t) has been derived by applying the next generation matrix algorithm to BSM. Each step of the algorithm is presented as follows.
Let X be the vector of infected compartments, such as exposed (E), infected (I), asymptomatic (A), hospitalized (H), quarantined exposed (E q ), isolated infected (L). Let Y be the vector of uninfected compartments, such as susceptible (S), quarantined susceptible (S q ), and recovered (R). Let F be the vector of new infection rates (transitions from Y to X) and V the vector of all other rates (not new infections).
Considering the following hypothesis regarding the Disease-Free Equilibrium (DFE), which is defined as the point at which no disease is present in the population, [47]: DFE = (s * , 0, · · · , 0) it follows that S = dS dt = −(βc(t) + c(t)q(1 − β))S(I + θ A) + λS q |DFE = 0 ⇒ −(βc(t) + c(t)q(1 − β))s * (0) + λ0 = 0 Then, deriving vectors F and V with respect to X and evaluating them at the DFE point, we have: Computing the inverse of the matrix V and multiplying by F, the next-generation matrix FV −1 can be obtained as in the following: Finally, by computing the dominant eigenvalue of the next generation matrix it is possible to obtain the following closed formulation for the daily reproduction number: R d (t) = c(t)β(1 − q)(−ρδ I (t)θ − ρ I θ − ργ I θ + ργ A + δ I (t)θ + I θ + γ I θ) γ A (δ I (t) + I + γ I ) It should be noted that (A8), obtained for BSM, only differs from (A7) by Tang with regard to the rate of home isolation for infected individuals ( I ) which replaces the disease induced death rate (α).