During the health crisis due to the COVID-19 pandemic, the development of reliable mathematical models, supported by the availability and analysis of complete and accurate data, prove to be fundamental tools for the interpretation and understanding of the epidemic outbreak, as well as for providing support to digital health [11] . If fed by accurate input data, epidemiological mathematical models can enable forecasting the trend of the outbreak and allow a description of the dynamics of the infection, in particular for those figures related to number of infected, hospitalized, healed and deceased individuals. Even when available, data may be difficult to polish, organize, and analyse in meaningful ways that allow for a coherent and immediate interpretation. This can be very useful, from one hand, to inform the general public and, from another hand, to provide institutional bodies and authorities with factual quantitative information that may enable the implementation of relevant decisions. This paper has a twofold goal. The first one is to present a newly developed mathematical dashboard (accessible online at the address https://www.epimox. polimi.it) that gathers all the most relevant data concerning infected, hospitalized, healed and deceased individuals, as well as to characterize the evolution of the epidemic outbreak through suitable fittings of the former figures, together with their first and second rates of variation. This analysis, which focuses on the situation in Italy at both national and regional scales -for all the 20 Italian regions -will enhance the interpretation and transparency of available data thanks to an immediate appreciation of the past epidemic development, together with its current trends. The second goal is to enable a comparison between the first epidemic wave, the one that initiated on February 21st, 2020 and exhausted in early June 2020, and the new (second) wave that took off in early October and is currently violently hitting Italy this Fall. The two waves feature indeed several analogies as well as significant differences that we aim to highlight and analyze. Moreover, still on the ground of this comparison, a preview on the expected trend of the epidemic outbreak for the near future will be discussed, based on a data extrapolation after curve registration from one side, and of a model-based epidemiological model from another side.
Our analysis highlights several important features. The more relevant are:
(i) In the early phase of the exponential outbreak, the timeline for the implementation of containment measures is crucial. It is more efficient to operate early with moderate restrictions than later with stricter ones; (ii) The second epidemic wave has a slower pace but a much larger diffusion that will eventually yield far worse figures (in terms of fatalities, hosted in intensive care units, etc.); (iii) The epidemiological mathematical model allows the investigation of various scenarios that conform to the different operational modes devised by the Italian government. A comparison among them allows to identify those that are potentially more effective to contrast the near future epidemic development.
An outline of this paper is as follows: in Section 1 a description of the COVID-19 epidemic data time series is supplied together with a description of the data processing tools (filtering, derivative, scaling, normalization) used to ease a straightforward identification of the main characteristics of each time series. The first epidemic wave is analysed and the results are compared with the corresponding analysis of the second epidemic wave. In Section 2, a data-based strategy to forecast the evolution of the epidemic is proposed. A model-based alternative for the forecast is finally introduced in Section 3, where it is used to compare different near future scenarios. Conclusions follow.
The data used in this paper are those made available by the Italian Dipartimento della Protezione Civile through the open data repository https://github.com/ pcm-dpc/COVID-19. Data are communicated on a daily basis and comprise: the number of individuals who are currently positive, isolated at home, hospitalized, hosted in ICUs (Intensive Care Units), the daily new positive cases, the cumulative number of deaths and recovered since the beginning of the pandemic, and the number of swabs performed. All these data are supplied at the regional and national levels, while the available data at a finer scale (provinces) are unfortunately limited to the count of total positive cases since the beginning of the epidemic.
From now on we will refer to these data as raw data. These raw data are then smoothed by resorting to a local polynomial regression based on the Savitzky-Golay convolution filter [14] : at day we attribute the value of the polynomial of degree that approximates, in the least squares sense, the 2 +1 values of the raw data centered on day n, i.e. in the range [ − , + ], where ≤ 2 . Here, we consider a cubic least squares polynomial and a window size of 21 days (that is, = 3 and = 10). A standard approach based on a weekly moving average could be obtained by taking = 0 and = 3. Finally, a backward moving average can also be considered where the average for day n is computed on the time window [ − 6, ] .
From now on, when not differently specified, the time series that will be presented will be the smoothed curves obtained using the Savitzky-Golay convolution filter.
The second step is to calculate the first and second rates of variation of the different compartments. First rate of variation describes how fast a trend is increasing or decreasing: change of sign from positive to negative for the first rate of variation indicates switching from increasing to decreasing in the corresponding curves, whereas a change on the second rate of variations denotes a change of convexity (a point of inflection).
The time series of the raw data (dashed lines) and filtered data (solid lines) for some relevant indicators at the national Italian level, namely the daily new positive cases, the daily deaths, the number of individuals that are hospitalized with symptoms and those hosted in ICUs, are reported in Figure 1 .
By analysing the time series of the total positive cases using a log scale (see Figure 2 ) along the full evolution of the epidemic in Italy, we can observe, for the first epidemic wave, an exponential growth ( ) = 1 (which is linear in log scale with slope 1 ) located in the first two weeks of March with a doubling time of approximately 3 days. During the second half of August, very likely because of the relaxation of social distances during holidays and newly imported cases from abroad, a few days of exponential growth can be observed, although featuring a much lower growth factor 2 , yielding a doubling time of approximately 15 days. Finally, what we will refer to as the second wave had its exponential growth during October, with growth factor 3 and a doubling time of about 8 days. We believe this is likely connected to the increased number of contacts associated with the school opening in September and related commuting, the restart of recreational activities in closed ambiences, and the full recovery of the working activities that were dramatically reduced during the spring lockdown, as well as (although at a minor extent) during the summer period.
We first focus on the first epidemic wave that occurred in Italy across the Spring 2020. In our dashboard, a specific time range can be selected and each time series can be normalised with respect to its maximum attained in the prescribed time range. This allows for an immediate identification of the day the different peaks have occurred and of their relative positions (relative delays). As displayed in Figure 3 , it can be noticed that the first compartment that reaches a peak is that of the number of daily new positive cases (on March 24th), followed 5 days later on March 29th by the peak of daily deaths. The peaks for the number of individuals hosted in ICUs and those hospitalized occur on April 1st and 3rd, respectively, that is after 7 and 9 days since the peak of daily new positive cases. Even if, at a first glance, the fact that the peak of daily deaths occurs before those of the hospitalized and hosted in ICUs may appear surprising, this is instead reasonable since the latter data do not refer to new daily entries but to the total number of individuals who are hosted in hospitals or ICUs at a specific date. Unfortunately, information (raw data) on the daily admissions in hospitals and ICUs are not available. Finally, the peak of total positive cases was reached on April 21st, that is almost 1 month past the peak of daily new positive cases. A similar analysis carried out on the rate of change (first derivative) of the different indicators may be used to determine when the initial exponential phase for each indicator is over. In Figure 4 , the normalized first derivatives of the same indicators discussed above are presented. In particular, we notice that the growth rate of the daily new positive cases reached its maximum on March 15th, while the growth rate of the number of hosted in ICUs is attained just 4 days later (March 19th) . This corresponds to an inflection point in the time series and indicates the time at which the growth rate starts decreasing.
It is worth noticing that, even if a maximum on the rate of change indicates that the evolution of the epidemic has overcome the initial exponential phase, relaxing too early the containment measures might turn not being wise. For instance, on March 15th the number of cumulative deaths at the Italian national level was 1, 809, and it raised up to 10, 779 on March 29th when the peak of daily new positive cases was 
The second wave of the COVID-19 epidemic in Italy is still in its rising phase at the current time (November 19th). At the date of the current analysis the only compartment that seems to have reached its peak is the daily new positive cases. However, since the raw data on this compartment is very noisy (see Figure 6) , due to the different number of cases tested across the weekend, it may be possible that the actual peak is not yet reached. A different indicator that can be used to assess whether the peak of new infections has been attained is the ratio between the daily new positive cases and the daily number of swabs (see Figure 7 ). As expected, this quantity has weekly oscillations with lower amplitudes and it also display a peak on November 11th.
All the other compartments, including that of the patients hosted in ICUs, that of the hospitalized and that of the daily deaths, are still rising (the associated normalized curves are displayed in Figure 8 ).
The first hint of a possible reduction of the growing rate has been observed for some compartments (Figure 9 ), such as the daily new positive cases, with a peak on the growth rate on October 24th, and on the the hosted in ICUs, with a peak attained on November 5th, that is 12 days later. We observe that this delay is larger (12 days against 4 days) than the one observed in the first wave. This may suggest a possible slower dynamics of the second wave, a fact that ought to be considered when trying to perform any forecasting of the near future evolution of the epidemic.
In particular, a possible estimate of the peak dates for the different compartments may be guessed. In particular, if the peak of daily new positive cases can be assumed on November 12th, then we may expect the peak of patients hosted in ICUs 20-25 days later (December 2nd-7th). Although this estimate is based on a merely speculative argument, nonetheless it will be confirmed by an approach based on curve registration that will be introduced in Section 2. 12 days Fig. 9 : Rate of change of daily new cases and hosted in ICUs highlighting the time shifts between the inflection points during the first epidemic wave
The same analyses that have been proposed at the national Italian level can be carried out at the regional scale, that is for each of the 20 Italian Regions. In Figure 10 , the normalized time series of daily new positive cases, hospitalized, hosted in ICUs and daily deaths are presented for 8 Italian regions, among which Lombardia, Emilia Romagna, Piemonte, Veneto that were severely hit by the first epidemic outbreak, and Lazio, Puglia, Campania, Sicilia that are far more evidently affected by the latter epidemic wave than by the former. The aim is to help highlighting some peculiar dynamics and trends that occur at the local, regional level. For instance, in Figure 11 the normalized time series restricted to the second wave are displayed. It can be noticed that in most regions the curves of patients hosted in ICUs and those hospitalized follow a similar trend, while in Campania the evolution of the hosted in ICUs has a different behaviour at the beginning of November: the trend may seem to indicate a possible saturation effect on the availability of ICUs (assuming that the raw data that have been communicated were correct). 
Peaks dates estimates based solely on similarities between the two waves may not be rigorously justified, especially if we consider that the two waves occurred under very different conditions, e.g. different individual behaviours (masks wearing, social distancing habits, ...), different social distancing rules and very different initial conditions. Still, they are most likely yielding better estimates than those that could be obtained by a naive polynomial extrapolation of the most recent observed data. As a matter of example, in Figure 12 , we display (in dashed lines) the forecast on the next 50 days for hosted in ICUs and daily deaths, based on a polynomial fit of degree 4 on the last recorded 20 days. The abrupt reduction that is obtained drives to zero these quantities in just a few weeks, which is totally unrealistic having in mind how typically epidemic waves decrease. It is clear that this forecast can only be considered when looking at a short time horizon, namely few days, while it becomes unmeaning for predictions on a longer time -a well known numerical misbehaviour of polynomial extrapolation (see [12] ).
A better forecasting strategy can be obtained by exploiting the similarities between the first and second waves, not only limited to the location of the curve peaks (as done at the end of the previous section), but considering the full curve shape. This can be accomplished by curve registration (see [13] for an overview on this subject). The registration procedure is sketched in Figure 13 and is performed by first computing the Exponentially Modified Gaussian (EMG) curve (blue line) that best fits the first wave. We denote this curve as ( ), 0 ≤ ≤ 1 , with 0 the first day of the recorded data (February 24th) and 1 equal to August 1st. Then a second minimization problem is solved to compute the time shift and scaling factors to apply to the computed EMG curve (red line) to best fit the rising portion of the second wave in the time range [ , ], with coinciding with October 15th and with the last recorder date. Namely, we look for the optimal time shifth, and the scaling factors¯1 and¯2 such that
where is the value of the considered data series at day .
For each data series, the fitted EMG curve and the optimal values for shift and scaling factors are computed and, in this way the shape of the first wave can be used to complete the second wave for the different compartments.
The extrapolation based on this curve registration procedure is implemented in the dashboard and, in Figure 14 , it has been used for the forecast of two relevant compartments such as hosted in ICUs and daily deaths. 
Epidemic forecasting is a difficult task because of the intrinsic variability and uncertainty of the Covid-19 pandemic. Uncertain, incomplete or inaccurate data (with regards to both initial conditions and time series of the different compartments) represent a serious threat, another being due to the partial knowledge of the behaviour of the specific infecting agent, not to say about the dynamic evolution of environmental and social conditions.
Since the seminal work [7] , where the first SIR compartmental model based on a system of nonlinear ordinary differential equations was proposed, a large variety of models have appeared in the literature -see, for instance, [6, 2, 8, 5, 4 , 1] -each trying to cope with specific aspects of the problem. If in the original SIR model the epidemic evolution is described by the number of individuals belonging to the susceptible (S), infected (I) and removed (R) compartments, several models consider increasing the number of compartments to include, for instance: • the exposed individuals (those who have already been exposed to the infecting agent but are not yet infectious); • possible splitting of the infected compartment into different classes according to the actual level of severity; • a distinction of the removed compartment in recovered and dead.
Although these models are typically working under many simplifying assumptions intrinsically related to their compartmental nature, they enable forecasting analysis which go beyond the simple extrapolation of a measured data. Indeed, once a model has been calibrated to be able to reproduce the past, observed time evolution until the present time, simulation regarding the future dynamics may include the investigation of different scenarios where the change of conditions (associated, for instance, to changing social distancing rules) are included by suitably modeling how the parameters change. One of the most critical aspects in the development of complex compartmental models is indeed related to their calibration based on available data. This critical aspect is twofold: on the one hand, data related to the different compartments may not be available (or they may not even be collected); on the other hand, the resulting data assimilation problem may suffer from limited identifiability of the parameters, as discussed in the recent paper [10] focusing on SIR-like models for COVID-19. A new compartmental model named SUIHTER has been recently developed by these authors (and will be presented in a forthcoming paper) with the goal of facing the first of these two issues, i.e. defining a model best suited for the data actually available. In the context of the COVID-19 epidemic outbreak, the time series that are daily collected and made available (which have already been described discussing the dashboard for the data analysis) lead us to consider for the new model the following compartment:
• S (susceptibles); • U (undetected), infectious individuals who have not yet been identified; • I (isolated), infectious individuals that have been quarantined at home; • H (hospitalized), infectious individuals that have been hospitalized with symptoms; • T (threatened): infectious individuals that are hosted in ICUs;
Based on the SUIHTER model, a preliminary analysis has been carried out with the objective of investigating the possible impact of containment measures, with various levels of restriction. The considered scenarios range between the current state (at November 19th) to a strict lockdown (similar to the one imposed during the Spring 2020), in particular we identify the following:
• Scenario 0: the current situation (at November 19th) in which some limitations are adopted at national level (distance learning from 9th grade, bar and restaurant with limited activity during the evening, . . . ), while additional stricter limitations (confinement within municipality limits, distance learning from 7th grade and universities, . . . ) are adopted only in some regions (denoted as "red" regions) ; • Scenario 1: the stricter limitations are imposed in the whole national territory (all regions becoming "red") starting on November 20th; • Scenario 2: a complete lockdown is imposed starting on November 20th.
Three compartments which are considered relevant for monitoring the evolution of the epidemic are displayed in Figure 15 , namely the isolated (representing the 95% of the total detected cases), the hosted in ICUs (which should not exceed the available ICUs places) and the cumulative deaths (that can be regarded as a dramatic synthetic indicator of the epidemic effect).
The model has been calibrated using a Monte-Carlo Markov Chain (MCMC) [3] approach implemented in the Python library pymcmcstat [9] looking for a robust estimate of the parameters minimizing the least square error with respect to the data. For each parameter considered in the calibration (in our case the transmission rates during the different phases of the epidemic), the MCMC calibration starts from a prior, in our case a uniform distribution around an initial guess, and provides its posterior probability density function.
The solid curves displayed in Figure 15 represent the time evolution of the different compartments in the 3 considered scenarios, while the shaded region surrounding each curve represents the 95% confidence interval associated with the probabilistic characterization of the model parameters. As expected, when more severe restrictions are imposed the epidemic curve evolves faster towards its flattening, even if these results seem to indicate that, considering the timing for this stricter restrictions, it may be by now too late to obtain a significant drop of the number of fatalities. Indeed, the cumulative deaths that amounts today at 47, 870 will overcome the impressive threshold of approximately 90, 000 on February 1st for Scenario 0, while Scenarios 1 and 2 would allow to lower this estimate down to 86, 617 and 84, 076, respectively.
A second set of simulations has been performed using the calibrated model and imposing that the same restrictions would be kept for a shorter time range and in fact relaxed on December 15th. (This strategy would be pursued to comply with an increasing demand for reopening most commercial and social activities in view of Christmas holidays.)
A detailed description of the limitation introduced have been published on the Gazzetta Ufficiale on November 4th and can be found at https://www.gazzettaufficiale.it/eli/gu/2020/ 11/04/275/so/41/sg/pdf This restriction relaxation is considered in three new scenarios (named 3, 4 and 5), which evolve before December 15th as scenarios 0, 1 and 2, respectively. The results, which are presented in Figure 16 , show that interventions with a limited time range of application may impair the effect of the restrictions. Public authorities should carefully weight whether a lighter restriction over a longer time frame could be more effective than a stricter one that however ends too early.
To better appreciate the comparison between the 6 scenarios, in Figure 17 we report the differences in terms of cumulative deaths with respect to Scenario 0. The positive contribution given by the additional restrictions in Scenarios 1 and 2 are clearly highlighted, as well as the negative effect of the relaxation at December 15th for all scenarios which partially impairs the results obtained by the previous stricter restrictions.
A similar comparison on the number of patients hosted in ICUs is reported in Figure 18 where the effect of the relaxation on December 15th amounts to an increase exceeding 200 units for both Scenarios 3 and 4.
In this paper we have examined the epidemic outbreak that occurred in northern Italy since early Spring 2020 and is still severely affecting the entire Country. Data were reported for several compartments (total and daily new positive cases, isolated at home, hospitalized, hosted in ICUs, cumulative deaths, recovered) that allow to provide a synthetic yet informative description of the epidemic evolution. Several analyses were carried out. A careful comparison between the two waves (that of last Spring, and the one that took off in early Fall) is made, followed by an extrapolation of the second wave along the next few months based on a registration procedure operated on the two curves. Finally, we have reported a few results obtained by using the new epidemiological differential mathematical model SUIHTER with the aim of providing a forecast on the epidemic curves for the next couple of months. 