key: cord-321735-c40m2o5l
authors: Manca, Davide; Caldiroli, Dario; Storti, Enrico
title: A simplified math approach to predict ICU beds and mortality rate for hospital emergency planning under Covid-19 pandemic
date: 2020-06-04
journal: Comput Chem Eng
DOI: 10.1016/j.compchemeng.2020.106945
sha: 
doc_id: 321735
cord_uid: c40m2o5l

The different stages of Covid-19 pandemic can be described by two key-variables: ICU patients and deaths in hospitals. We propose simple models that can be used by medical doctors and decision makers to predict the trends on both short-term and long-term horizons. Daily updates of the models with real data allow forecasting some key indicators for decision-making (an Excel file in the Supplemental material allows computing them). These are beds allocation, residence time, doubling time, rate of renewal, maximum daily rate of change (positive/negative), halfway points, maximum plateaus, asymptotic conditions, and dates and time intervals when some key thresholds are overtaken. Doubling time of ICU beds for Covid-19 emergency can be as low as 2-3 days at the outbreak of the pandemic. The models allow identifying the possible departure of the phenomenon from the predicted trend and thus can play the role of early warning systems and describe further outbreaks.

Covid-19 is the most exacting pandemic since the Spanish flu of more than a century ago. The fast outbreak of Covid-19 and the wide spread all over the world transformed a local disease, initially located in China, into a global problem; thus the name: pandemic (Fauci et al., 2020) . Italy (60.5 million inhabitants) is one of the most plagued nations by this pandemic with almost 28000 official deaths in the hospitals (as of 30-apr-2020) and over 4000 ICU beds to treat patients at the peak of Covid-19 emergency (Livingston & Bucher, 2020; Remuzzi & Remuzzi, 2020) . Likewise, Lombardy (10 million inhabitants) is the most crowded region of Italy and the most hit and afflicted region of Europe with over 13800 official deaths in the hospitals and almost 1400 ICU beds (Grasselli et al., 2020a; Grasselli et al., 2020b) .

Before the pandemic, the number of ICU beds for any treatments in Lombardy was about 700 (Grasselli et al., 2020a; Manca, 2020a) and in winter months those beds are usually 85-90% occupied (Grasselli et al., 2020a) . The pandemic called for doubling that number at an incredible pace to cope with the repeated tsunami waves of very complicated patients affected by DIP (i.e. diffuse interstitial pneumonia). Those patients required ever-increasing treatments ranging from oxygen masks, to helmet C-PAP, to NIV, and eventually tracheal intubation (Desai & Aronoff, 2020) . In some hospitals, such as the Lodi hospital in Lombardy (the one where the first Covid-19 italian patient was diagnosed and sheltered in ICU, a 38 year-old sporty male with no comorbidities who remained in the ICU ward for 18 days before moving to PIP, i.e. post-intensive care, for other 14 days) before the pandemic there were 6 ICU beds (Cutuli, 2020) . In a matter of few weeks, the beds became 18 and eventually 27 at the peak of the emergency (i.e. 4.5 times more). This called for huge efforts in terms of burden on medical doctors, nurses, and managing team that literally transformed the original wards, operating rooms, and recovery rooms to set up and operate that high number of new ICU beds. The mathematical models describing the phenomenon dynamics, which are reported in this paper, have allowed understanding the fast evolution and preparing daily for the continuously increasing burden. Besides the predicted numbers, those models allowed also forecasting the different phases of the pandemic and quantifying some basic indicators about the daily variations, the key times, the key figures, the expected decrease, the progressive reach of a maximum plateau before facing with the decrease of ICU beds for Covid-19 which we are measuring right now. Together with the ICU beds, which are reported daily by the Italian civil protection department every evening at 6 PM CET (Dipartimento della Protezione Civile, 2020), we also monitored the dynamics of official deaths.

Official deaths are the ones that occur to patients sheltered in hospitals after they result positive to a nasal swab and possibly to a CT scan. Official deaths are important as they measure one of the two possible outcomes of hospital treatment, either success or failure. Physicians have to know and somehow forecast not only the daily number of deaths but also the asymptotic value predicted at the end of the pandemic. Such a number is a big burden for the psychology of those who struggle for life (e.g., medical doctors, nurses) but it can prepare them to deal with that fatal outcome and somehow find an upper bound of the phenomenon, as shown in the following. As detailed in Section 2, the mathematical models can be incomplete in the sense that further models might describe with the same or even higher efficacy the dynamics of real data. Nonetheless, we will show the reliability, robustness, and quality of the suggested models that can describe the phenomenon evolution with a high degree of precision.

The object of this paper consists in making the reader aware of the main features of a pandemic in terms of its dynamic evolution based on key points, key intervals, key dates, and key numbers. The mathematical models can be used to get a background of past events, understand the current situation, and forecast the pandemic evolution over either short-term or long-term horizons. The qualitative and quantitative assessment of the pandemic dynamics allows medical decision-makers to plan for the emergency, prepare for severe times, and finally relax the safety measures and revert to standard elective medicine when the pandemic deflates. The models can outline possible deviations of the pandemic from the expected evolutions and therefore play the role of early warning systems in case of new outbreaks.

Mathematical modeling can be of real help in describing, understanding, and eventually forecasting how a virus diffuses within a domain (e.g., population, province, region, country, and continent). Mathematical models feature a set of mathematical equations that include a number of adaptive parameters that can be determined numerically grounding on available real data (Panovska-Griffiths, 2020) . Once these models are refined and the adaptive parameter computed, one can use these models to understand what happened in the past (e.g., lessons learnt) and even more important to forecast what might happen in the future (e.g., emergency planning, resource allocation) (He et al., 2020; Poston et al., 2020; Steinberg et al., 2020) . This Section focuses on describing the time evolution (i.e. dynamics) of Covid-19 pandemic centering on only two variables that are the most reliable and therefore robust data, namely the number of ICU patients (aka ICU beds) and official deaths in hospitals. These data are reported daily by the Italian civil protection department (Dipartimento della Protezione Civile, 2020). Day 1 is when the first patient was sheltered in the intensive care unit (22-Feb-2020 in Italy).

We chose to describe and follow, but also to predict, the dynamics of those real values with the simplest mathematical models available. Such models were carefully chosen by grounding on a keep-it-simple approach. The models are regression-based, which means that they minimize the distance between model predictions and real data. We selected the models that best describe the dynamics of the phenomenon and feature the lowest number of adaptive parameters. Indeed, we only used either two or three adaptive parameters to keep the approach simple, avoid overparameterization, and preserve numerical robustness and stability (Manca, 2020a (Manca, , 2020b .

This approach allows implementing those models on any computers without the necessity of relaying on dedicated programming tools. Intentionally, we implemented a set of math models that can run on Excel as it is widely available on any operating systems (e.g., Windows, Mac, Linux) and works on several platforms (e.g., computers, tablets, mobile phones). In addition, Excel (even though it is not a real programming software) requires a lower learning curve respect to other programming tools such as Matlab or Mathematica or programming languages such as Fortran, C/C++, Python (just to cite a few).

The supplemental material to this manuscript includes an Excel file with the proposed models that can be used to extract important information about past, present and future situation, and to draw the diagrams that describe the phenomenon dynamics and allow understanding its qualitative and quantitative evolution.

Once the models are identified (i.e. regressed respect to real data), they show new values of the adaptive parameters which play a role in quantifying some important indicators about the pandemic.

We intentionally did not use any epidemiologic models such as SIR, SIRD, SEIR, SEIRD models that are based on a set of few differential equations with initial conditions and a number of adaptive parameters as well as strong assumptions and simplifications (C.H. Li et al., 2014; G.H. Li & Zhang, 2017; Magal et al., 2016; Xia et al., 2009) . The family of SIR models shows a higher detail of description of the phenomenon. However, according to our opinion, their use is good for running parametric predictive scenarios based on a number of assumptions and hypotheses that are quite sensitive to the selection of proper values of the adaptive parameters and functional description.

Reliable values of those parameters will be available only at the end of the pandemic and depend significantly on the political decisions endorsed at different times and intensity by each country (or even each region locally) on social-distancing measures and more in general on nonpharmaceutical interventions (Bayham & Fenichel, 2020; Cowling et al., 2020) .

Conversely, the math models proposed in this paper can be used daily (by updating the ICU and deaths data respectively) and automatically adapt to the evolution of those values.

The phenomena behind ICU beds and deaths follow two different evolutionary curves in case of a pandemic. ICU beds start from zero at the very beginning of the pandemic, they keep on increasing once the pandemic deploys its intensity up to a maximum value. After the plateau, the ICU beds start decreasing and reach a final null value when the pandemic expires.

Conversely, the deaths number continues increasing monotonically and reaches a final value, which is the fatalities toll the country/region has to pay to the pandemic plague. Actually, the deaths number is a cumulated value that increases with the daily fatalities.

As far as big regions and countries are concerned, the number of ICU beds increases monotonically up to the maximum plateau. This is what both Lombardy and Italy (Onder et al., 2020) dealt with. Conversely, small regions in terms of population and/or positive cases may experience non-monotonically increasing trends. For instance, this is the case of Umbria a central region in Italy with 882,000 inhabitants living in small municipalities distributed on its territory, which experimented a trend of ICU beds quite different from the most involved Italian regions, the northern ones. Indeed, Umbria saw often the number of ICU beds remaining constant for a few days. Now and then, that number decreased and then increased again with rather small fluctuations (i.e. few units per day) due to the tiny number of infected people. For this kind of regions/countries, the models and trends discussed in this paper are not recommended and should be avoided.

Back to Lombardy, after the fast explosion of ICU patients following an exponential growth (Manca, 2020a; Remuzzi & Remuzzi, 2020) , it experienced a saturation condition protracted over several days (Manca, 2020b) before reaching the maximum plateau and turning into the descent trajectory. That saturation condition was not constant and the healthcare system was flexible enough to increase at a lower extent (than the necessary one) the number of ICU beds. The "saturation" term means that the capacity of daily increasing the number of beds, to cope with the tsunami wave of new patients requiring an ICU treatment, was lower than the compulsory number of new ICU beds. It is worth observing that the continuous creation of new ICU beds took to a subtle drift in terms of treatment quality as new beds were created under huge pressure with an above-limited capacity in already existing wards and intensive/step down areas reconfigured in a matter of few days. Coupling the saturation condition with the quality of new ICU beds had an impact on the treatment quality of ICU patients and consequently on the fatalities toll .

Indeed, there was also a limiting factor played by the number of physicians and nurses subject to over-intensive activities, with shifts of 15-16 hours per day, day over day for more than one month, sometimes being themselves hit by Covid-19 (more than 150 medical doctors left their lives due to Covid-19) (ADNkronos, 2020).

For the sake of simplicity, let us focus first on the qualitative trend of ICU patients. The first days of the pandemic see an exponential increase (Remuzzi & Remuzzi, 2020 ). An exponential curve is characterized by a doubling period that remains constant if no external measures or disturbances occur. It is rather easy to show that a phenomenon is exponential when it fits a straight line in a semilogarithmic diagram (i.e. a diagram where x-axis is linear (time, number of days) and y-axis reports the logarithm of the real values (ICU beds). It is possible to measure the quality of the fitting curve (i.e. straight-line) respect to real data by means of the determination coefficient ( 2 R ).

This is a nondimensional coefficient ranging from 0 to 1 (Hahs-Vaughn, 2016). The higher the value, the better the consistency of the model with real data. In case of both Italy and Lombardy, the first epidemic days saw 2 R values with either two or three nines as decimals, almost equal to 1, i.e. 2 0.99, ,0.999 R  .

The exponential growth remains stable for about 15-16 days. For the sake of correctness, every day the slope of the straight line decreases slightly and the intercept increases slightly (i.e. the point where the line cuts the y-axis). This is intrinsic to the epidemic dynamics and shows that the phenomenon starts reducing its momentum, although it remains exponential (see also Figure 1 ).

There is a time, days 16-22, when new values of ICU beds start departing from the straight line and begin to follow a parabola on the semilog plane. This is the case of the so-called exponentially modified Gaussian (EMG) trend (Golubev, 2010) . On a standard Cartesian plane, the phenomenon is well described by two distinct regions. The first one with upward concavity and the second one with downward concavity. The upward concavity reflects the exponential trend that becomes progressively linear at the inflection point where the concavity changes from upward to downward. At the inflection point, the maximum positive velocity of change occurs. In math terms, we say that at the inflection point the first derivative (i.e. maximum daily change of the observed variable, ICU beds) has the highest positive slope (i.e. maximum increment). These are the toughest days to stand the tsunami wave as day over day the increment of ICU beds is the highest one. In addition, the bigger the region/country with its number of infected patients, the higher its inertia. After the inflection point, the phenomenon continues to increase but at a slower pace. In this phase, the increase is monotonic, which means that the number of ICU beds continues to climb, but the daily growth reduces progressively up to a maximum plateau (around days 39-43). Again, the inertia of the system is usually high and instead of having a train at the top of a rollercoaster that quickly starts falling, the plateau remains stable for a few days. That maximum plateau is critical for the whole system as the ICU stay is rather long (Murthy et al., 2020) . Usually, patients remain in ICU wards at least fifteen days (with twenty-day stay the standard value) (Cutuli, 2020) and, respect to Covid-19 emergency, this quite a long time allows describing the whole ICU beds inflation period with curves such as the logistic (Hosmer et al., 2013) or the Gompertz (Panik, 2014) ones. At the maximum plateau, small oscillations may occur but finally the system starts decreasing (days 42-45). The reverse trend of the first ascending period occurs in this second descending phase up to a final null value when the pandemic is out. When the descent starts, the concavity is downward and becomes upward after a new inflection point (with negative slope, i.e. downhill). Before the inflection point (days 60-70), the phenomenon increases the velocity of reduction of ICU beds.

After the inflection, the reduction pace slows down and approaches progressively the final null plateau (days 105-120). For the sake of simplicity, the two ascending and descending tracts of the overall ICU beds phenomenon can be described with two separate segments of suitable models such as (i) exponentially modified Gaussian (EMG), (ii) logistic, (iii) Gompertz curves. All these models ground on either physical or biological foundations. The logistic model (Hosmer et al., 2013) was originally proposed in 1838 by Pierre François Verhulst to describe the growth of populations where the rate of reproduction depends on both the existing population and the amount of available resources. The Gompertz model (Panik, 2014) is similar to the logistic one and was designed by Benjamin Gompertz in 1825 to describe the law of human mortality. The EMG model (Golubev, 2010) fits well processes involving normally-distributed inputs and exponentially distributed outputs. EMG characterizes the transition probability of cellular cycles and embraces both the deterministic and probabilistic contributions.

The same qualitative discussion made to model the inflation period of ICU beds may be adapted to model the regional/national fatalities. Contrary to ICU beds, the deaths curve is only monotonically increasing and finally approaches the maximum plateau once the pandemic expires.

The deployment of that curve is longer and shows a lag time respect to the ICU curve of about 12-14 days in proximity of the first ascending inflection point. It is worth observing that the fatalities curve is a pure cumulative curve (i.e. integral curve) that sums up the single death tolls experienced daily.

Quantitatively, the inflection point, when the maximum daily fatalities occur, arrives at days 34-38

and it takes a total of 100-120 days to reach 98% of the final maximum plateau. Mathematically, the logistic and Gompertz models predict the practical extinction of deaths after days 150-200.

The hiatus between the logistic and Gompertz models increases as the pandemic deploys its dynamics. The predicted numbers should be taken with a grain of salt as (i) these are just forecasts when this manuscript was written; (ii) the fatalities phenomenon depends heavily on the social distancing measures enacted by regions and countries together with the very uncertain outcomes that the so-called phase 2 will produce when people progressively reduce social-distancing and start again to live and work as before the pandemic although with a higher awareness about safety and health related risks. Only time will tell if these forecasts approach the real phenomenon or further modeling issues should be accounted for.

The mathematical description of the exponential model is:

The mathematical description of the logistic model is:

The mathematical description of the Gompertz model is:

where t is time, ,, abc are adaptive parameters and 2.718281828459045235... e  is the Euler's number. y is the dependent variable that is predicted by the model (i.e. ICU patients or deaths). The exponential model features two parameters, whilst the logistic and Gompertz ones feature three parameters. As discussed in Section 2.2 the approach to models parameterization is the minimal one and in line with the keep-it-simple philosophy. As far as the independent variable is concerned, t can be considered either continuous or more correctly discontinuous as it corresponds to the number of days since the start of the pandemic in a specific region/country. Equation (1) in semilogarithmic coordinates assumes the form:

which is the equation of a line whose slope is b , the intercept is

It is straightforward to determine the doubling time of the exponential growth from Equation (4):

The higher the slope of the line (Equation 4), the shorter the doubling time of the exponential phenomenon.

The logistic curve (Equation 2) is symmetric respect to the inflection point and the sigmoid function is a special case of the logistic function (Hosmer et al., 2013) . The Gompertz curve is similar qualitatively to the logistic function but it is not symmetric and takes a longer time to reach the final asymptotic plateau. In addition, the Gompertz plateau is higher than the logistic one (once the parameters of both curves are identified to minimize the distance from the same set of real data). Table 1 shows the analytic formulae that allow computing the inflection and halfway points. The inflection point identifies the time when the rate of change is highest, whilst the halfway point determines the time when the phenomenon reaches a value that is half than the maximum plateau. It is worth observing that the time when the halfway condition occurs does not mean that after an equivalent amount of time the whole phenomenon completes. Indeed, the halfway condition refers to the dependent variable ( y ) and not to the independent one ( t ).

Both the logistic and the Gompertz models have the same maximum plateau when t , ya  . For the sake of clarity, each curve has its own a value that is computed by minimizing the sum of the squared distances between the real data (whose cardinality is NP ) and the model predictions (Equation 6) through the following nonlinear regression procedure (Hosmer et al., 2011) :

The minimization of Equation (6) is multidimensional and unconstrained even though the optimizing procedure may be eased by specifying that the degrees of freedom (i.e. the adaptive parameters) are positive. Once the model is identified, it is possible to answer the following questions: how much time does the phenomenon take to reach a given percentage respect to the maximum plateau? For instance, at what time the death phenomenon will reach 98% or 99% of the asymptotic condition (i.e. the maximum final plateau when the pandemic expires)?

Mathematically one has to solve the following nonlinear algebraic equation:

Where p is the desired fraction (e.g., 0.98 p  if the desired percentage is 98%). An easier approach to solve Equation (7) consists in reporting in a table the model predictions (e.g., an Excel spreadsheet) as a function of time and to search for the first time when the phenomenon overtakes that percentage (as the asymptotic value a is known). Equations (1-3) are good at describing the monotonically increasing phase of ICU patients and the whole fatalities phenomenon.

Once the ICU phenomenon touches its maximum value it starts decreasing as discussed in Section 2.2. The increase-plateau-decrease region can be described by an exponentially modified Gaussian (EMG) model whose mathematical description is: 

where 10 c C  is the multiplying factor, 10 bt is the exponential term and 2 10 at is the gaussian term. Equation (8) is even more flexible as it can be applied to the uphill part of the ICU phenomenon (besides the downhill one). For both the uphill and downhill parts, the EMG model always showed positive values for , bc and negative values for a which is a mandatory condition for the gaussian contribution. The weak point of EMG is that it is rather precise in predicting values over short-time intervals, whilst it is less reliable over longer periods when applied to the uphill portion of the phenomenon. However, across the maximum plateau and the descending section, the EMG predictive performance is rather good.

The evaluation of the daily rate of change of the phenomenon can be carried out either in a discrete way by computing the difference between two consecutive model predictions (for instance in the Excel spreadsheet) or by analytic differentiation ( dy y dt

).

The first derivative of the exponential function is:   ln 10 10 bt y a b  

The first derivative of the logistic function is:

The first derivative of the Gompertz function is:

The first derivative of the EMG function is:

Equations ( Once the logistic and Gompertz curves reach the maximum plateau, they conclude their scope as far as ICU beds are concerned. However, their intrinsic monotonic nature can be suitably exploited to describe the descent towards the end of pandemic when ICU beds are null. Both logistic and Gompertz models can be rewritten according to a reverse formulation respect to the original one.

It is sufficient to translate them, change their sign, and move the initial condition to the maximum plateau as reported in the following.

The mathematical description of the reverse logistic model is:

The mathematical description of the reverse Gompertz model is:

For the sake of clarity, parameters a and 0 t are known as they are the plateau value and the corresponding time when the maximum of ICU beds is reached, respectively. Consequently, Equations (14-15) reduce to two adaptive parameters, , bc. Both Equations (14-15) exhibit the same final null plateau when the pandemic expires.

Finally, the first derivatives of the reverse logistic and Gompertz models are respectively: 

The models reported in Section 2.3 can be selected and used critically according to their performance and consistency with real data published daily in most countries. It is necessary to determine the values of the adaptive parameters by means of an identification procedure that is based on either a linear or a nonlinear regression of real data depending on the mathematical nature of the model (H.H. Zhang et al., 2011) . The good news is that Excel implements a multidimensional optimizer, which is capable of solving the minimization problem of Equation (6).

Any other programming environments can solve as well that same problem provided a multidimensional unconstrained optimization routine is available. The interested reader can refer to the Supplemental material to identify the models and their parameters.

One of the performance indicators to discriminate the models and find the most representative is RMSE , i.e. the root mean square error (Equation 18). Equally, one can evaluate the either the mean ( MAE ) or the median ( MEDAE ) absolute errors via Equations (19) and (20) 

This section shows how to choose, use, and extract important information from the proposed models. The case study is applied to both Italy and Lombardy in terms of ICU patients and deaths.

In the very first days of the pandemic, the phenomenon is purely exponential (Equation 1 ). This assumption proved true after observing the linear trend of dependent variables (i.e. ICUs and deaths) in semilogarithmic coordinates (Equation 4) and evaluating the determination coefficient. Figure 2 shows how, in the first days of the pandemic, the ICU beds grow exponentially as the real data lay almost perfectly on a straight line. The linear trend is even stronger after day 10. It is worth remarking that the y-axis is logarithmic which means that the integer values are powers of ten in linear coordinates. For the sake of clarity, 1 means 10 ICU patients, 2 means 100 ICUs and 3 means 1000 ICUs. The linear trend is confirmed by the determination coefficient 2 R that approximates 1 (i.e. 2 = 0.9917229 R ).

As time progresses, the phenomenon leaves progressively the purely linear trend in semilog coordinates and moves towards a quadratic behaviour (in those same semilog coordinates) that is embodied by an EMG model (see Equations 8-9). Figure 3 shows the high affinity of real data with the EMG model, which is also confirmed by the very high value of the determination coefficient (i.e. 2 = 0.9991630 R ). It is not common to have real data (i.e. experimental data) that are so consistent with a model. Somehow and to a first sight, data appear tamed or even worst manipulated to make them smoother, although this is absolutely not true. Actually, this smoothing effect is the result of large numbers and cumulated curves. In the following, we will illustrate some diagrams, 

coincides with the logistic predictions and, as a result, it is almost invisible in the diagram. The trend of real data in Lombardy fluctuates a bit more than the corresponding Italian trend due to the extreme pressure and saturation condition exerted by the Covid-19 emergency on that region.

When Lombardy exceeded the threshold of 900-1000 patients in intensive care units, it was a daily struggle with Covid-19 pandemic to create further intensive care beds to shelter the continuously increasing wave of patients requiring ICU treatment and tracheal intubation. That was the challenge that absorbed the most from the medical doctors and nurses of intensive care wards (Cutuli, 2020) . Consequently, when the first derivative is equal to zero (as for the EMG model in Figure 5 ) the maximum plateau of the corresponding model occurs. In Lombardy, the EMG model of Figure 5 identifies the maximum plateau of the ICU patients' trend at day 41 (2-apr) and the same happens for Italy. Equally, the minimum values of the first derivatives for the three models of Figure 5 occur for Lombardy around day 58-61 (19-22 April) and for Italy around day 56-57 (17-18 April). These are the days when the highest daily decreases of ICU patients are expected. Afterward, the daily decrease continues but its intensity lowers and finally becomes negligible in proximity of the null plateau (see also Figure 6 ). Figure 4 . The red horizontal bottom line shows the 10% threshold respect to the maximum measured value and identifies the times when most of the ICU wards should be empty. Last available real data on 30-Apr.

Equally, the reverse logistic and reverse Gompertz models predict that for Italy the descent below the 10% threshold of ICU patients will occur on 25-May and 28-May respectively (i.e. 52 and 55 days after the maximum plateau, which is also at days 94 and 97 respectively). These values are in line with the predictions for Lombardy. Equally, the last remaining 10 ICU patients would turn out on 18-Jul and 31-Jul for the reverse logistic and reverse Gompertz models respectively (i.e. 106

and 119 days after the maximum plateau, which is also at days 148 and 161 respectively, based on 30-Apr real data availability).

The description of death dynamics can ground on the same modeling approach based on the logistic, Gompertz, and EMG curves (Equations 2, 3, and 8-9 respectively). For the sake of brevity, the initial exponential growth (Equation 1) will not be reported also because the aforementioned models feature in their formulation an initial part where the curve follows the exponential trend. allows also predicting (as of latest data available on 30 April) that the total expected deaths toll will be 34037 and that 98% of that value will be reached on day 111 (11 June). It is worth observing that the Gompertz curve is not symmetric and that after the inflection point the change of concavity makes it approach the final maximum plateau slower than other curves such as the logistic and EMG ones. At the end of June (day 130), the Gompertz model predicts that 99.17% of the whole phenomenon is manifested and that 282 fatalities remain before the pandemic is out. It is straightforward to remark that these are predictions based on a mathematical model that is subject to the availability of reliable data. The more the real data the better. Indeed, every day the national and regional reports release new data and one can carry out forecasts that are more reliable. In addition, the model forecasts ground on the assumption that future values will be registered and made available according to the same boundary conditions in terms of social distancing measures and data collection conditions. It is worth observing that the behaviour of Italy in terms of pandemic dynamics all over its territory is not uniform. It would be a huge mistake to assume a uniform distribution of ICU patients and fatalities in its regions as the death toll (as of day 69) sees more than 49% 

The models of Section 2.3 applied to the case study of Lombardy and Italy proved their efficiency in reproducing real data and were used to forecast the evolution of key parameters as the number of ICU patients and deaths on both short and long-time horizons. The same models can be applied to different countries and regions (Hopman et al., 2020) if reliable and timely data are collected and shared by the public bodies of civil protection and health. The same models can also be applied to describe the dynamics of other variables provided they are reliable and collected according to the same standards and methodologies. For instance, at least in Italy, the swabs done to test possibly infected people adopted different approaches in different regions and they are not representative of the expected number of infected people (Grasselli et al., 2020b; Molinari et al., 2020) . Indeed, a patient initially positive to a swab has to receive two negative in a row results before being declared no more contagious (e.g., a medical doctor of a Lombardy hospital was tested with 9 consecutive swabs before being declared negative. He had to wait 46 days from the first positive swab). Therefore, the total number of swabs is not representative of the number of people tested. In addition, a large number of probably infected people have not been tested with a swab and nonetheless have to undergo a quarantine of at least 14 days. This happens to a large number of people living together at home with just released patients from hospitals and to asymptomatic or paucisymptomatic individuals (Kimball et al., 2020; Nicastri et al., 2020) .

The three main models proposed in this paper (i.e. EMG, logistic, and Gompertz), either direct or reverse, can be used for qualitative and quantitative purposes in the Covid-19 emergency. They play two separate roles. First, they track real data and allow discriminating among models to find the most reliable one(s) but also they allow understanding if any unexpected trends start occurring. Being analytic and fully developed in terms of functional dependency, they are continuous and feature analytic derivatives of any order. This means that not only the cumulated values but also the daily variations can be observed and monitored to understand if any drifts of the phenomenon are occurring. Second, these models can be used to predict the evolution of the phenomenon over either short or long-term horizons. During the very first tough days, the most important forecasts are required over short-time intervals to understand the doubling times of ICU beds and allocate/prepare suitable resources to cope with the repeated tsunami waves of patients with acute respiratory distress syndrome. When the pressure is at least partially released, the models can be used to predict when some safer thresholds will be reached. In those cases, it will be possible to decrease the number of Covid-19 beds, to close dedicated wards, to reallocate human resources, to restart elective hospital activities that were intentionally shutdown or reduced to the minimum to focus on the demanding emergency. Also, the reported models have been exploited for decision making (see Figure 9 ) to understand if some extreme decisions in terms of resources allocation would timely meet the dynamics of phenomenon (e.g., reallocation of wards, building new hospitals, moving ICU patients to other hospitals and regions/countries).

Models can be compared for their consistency and precision in describing the phenomenon by evaluating suitable key performance indicators (KPIs) as the ones reported in Equations (18-20) (K. Zhang et al., 2015) . Table 2 and Table 3 There is not a clear best-performer among the models used to describe the dynamics of ICU patients. However, the Gompertz model behaves on average better than the logistic one both uphill and downhill. The EMG model can predict the transit across the maximum plateau before leaving the stage to the reverse logistic and Gompertz models. The performance of the EMG model is better in Lombardy than in Italy and it is comparable with the uphill and downhill trends of direct and reverse logistic and Gompertz models.

It is worth remarking that each KPI is characterized by a specific functional dependency and therefore the values of a KPI should be considered just to compare different models belonging to the same region/country (K. Zhang et al., 2015) . For instance, it is reasonable that the KPIs of Italy are higher than those of Lombardy are as the involved numbers are higher for the whole nation respect to that single region. Indeed, the maximum plateau counted 4068 patients in Italy and 1381 patients in Lombardy. Since the unit of measure of all the KPIs reported in Table 2 and Table   3 is the number of ICU patients, it is possible to observe that proportionally the same models were a bit less accurate in the predictions of Lombardy respect to those of Italy. However, the prediction capability of all the proposed models is rather accurate as they describe the allocation of thousands of ICU beds and can cope with the intrinsic oscillations of biological and massive organisms (e.g., regions and nations) when they are hit by a pandemic that exerts paramount pressure and burden on their vital resources (Livingston & Bucher, 2020) .

The comparison of models for the prediction of fatalities is easier as the phenomenon is somehow simpler since it is monotonically increasing and once it passes the inflection point one has to wait only for the final plateau. In case of fatalities, the Gompertz model is the clear winner in terms of precision and reliability over the whole horizon of available real data (from the very beginning of pandemic until day 69, 30 April 2020 for Italy and Lombardy). Conversely, the logistic model had to be discarded after day 55 as its predictions were too optimistic and it calculated that final plateau would occur earlier and with a too low fatalities toll. Since the number of deaths in Lombardy and Italy were rather high at day 69 (27967 and 13772 respectively), the KPI values reported in Table 4 are pretty low respect to the dimension of the phenomenon (where the expected final number of fatalities was 34037 for Italy and 15519 for Lombardy according to the projections of day 69, 30 April).

The paper presented and discussed a few regression models to predict two of the most important variables in a pandemic from the point of view of decision-making and emergency planning. These are the number of ICU patients who must be sheltered in dedicated Covid-19 wards and the number of fatalities. These models can be applied to different regions and countries as the pandemic phenomenon has the same qualitative features. The two, maximum three, adaptive parameters allow describing quantitatively the dynamics of those different regions and countries.

Indeed, every region and every country are characterized by different features bound to their territory nature, population distribution (in terms of age, density, life-styles, human interactions, family habits), and political decisions (in terms of progressive/immediate, relaxed/inflexible lockdowns, social distancing, and other non-pharmaceutical interventions). The mathematics behind the proposed models is rather simple and can be implemented in an Excel file that can be used by most decision makers and medical doctors. Based on real data that most countries/regions produce daily, these regression models can be identified in terms of their adaptive parameters and used to forecast the trends of ICU and deaths on either short or longterm horizons. These same models can also be adapted to track other variables as far as the reliability of those variables is good enough to preserve their consistency. Once identified, the models can determine (i) the doubling time of the phenomenon, (ii) the inflection point where the daily increment of the phenomenon is either maximum positive (uphill) or maximum negative (downhill), and (iii) the maximum plateau. In addition, these models can evaluate the time when some conditions occur, such as the achievement of a fraction of the maximum plateau. By monitoring some suitable KPIs, the user can assess the performance of the models and understand when they should be quit or if the reverse version of the model (in case of the logistic and Gompertz models) would fit better the new real data. The periodic update of these models and the important details that one can extract by observing the diagrams that compare the prediction capabilities respect to real data allow detecting possible drifts of the phenomenon and can play the role of early-warning systems if the pandemic derails from the expected evolution.

Coronavirus, 150 medici morti in Italia

Impact of school closures for COVID-19 on the US health-care workforce and net mortality: a modelling study. The Lancet

Impact assessment of nonpharmaceutical interventions against coronavirus disease 2019 and influenza in Hong Kong: an observational study

COVID-19: the Lodi model (in Italian)

Masks and Coronavirus Disease 2019 (COVID-19)

COVID-19 Italia-Monitoraggio situazione

Covid-19 -Navigating the uncharted

Exponentially modified Gaussian (EMG) relevance to distributions related to cell proliferation and differentiation

Critical Care Utilization for the COVID-19 Outbreak in Lombardy, Italy: Early Experience and Forecast during an Emergency Response

Baseline Characteristics and Outcomes of 1591 Patients Infected with SARS-CoV-2 Admitted to ICUs of the Lombardy Region

Applied multivariate statistical concepts

How to transform a general hospital into an "infectious disease hospital" during the epidemic of COVID-19

Managing COVID-19 in Low-and Middle-Income Countries

Applied Survival Analysis: Regression Modeling of Time to Event Data

Applied Logistic Regression: Third Edition

Asymptomatic and presymptomatic SARS-COV-2 infections in residents of a long-term care skilled nursing facility -King County

Analysis of epidemic spreading of an SIRS model in complex heterogeneous networks

Dynamic behaviors of a modified SIR model in epidemic diseases using nonlinear incidence and recovery rates

Coronavirus Disease 2019 (COVID-19) in Italy

Final size of an epidemic for a two-group SIR model

Analysis of the number growth of ICU patients with Covid-19 in Italy and Lombardy

Dynamics of ICU patients and deaths in Italy and Lombardy due to Covid-19

SARS-CoV-2: The Lombardy scenario in numbers

Care for Critically Ill Patients with COVID-19

Coronavirus disease (COVID-19) in a paucisymptomatic patient: Epidemiological and clinical challenge in settings

Case-Fatality Rate and Characteristics of Patients Dying in Relation to COVID-19 in Italy

Growth curve modeling: Theory and applications

Can mathematical modelling solve the current Covid-19 crisis?

Management of Critically Ill Adults with COVID-19

COVID-19 and Italy: what next? The Lancet

Calculated decisions: COVID-19 calculators during extreme resource-limited situations

Epidemics of sirs model with nonuniform transmission on scale-free networks

Linear or nonlinear? Automatic structure discovery for partially linear models

A comparison and evaluation of key performance indicator-based multivariate statistics process monitoring approaches

The authors acknowledge the valuable discussions with MD Piergiorgio Villani (Lodi hospital), MD Giovanni Mistraletti (San Paolo di Milano hospital), and MD Francesco Trotta (Lodi hospital).