A simplified math approach to predict ICU beds and mortality rate for hospital emergency planning under Covid-19 pandemic.

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.

Introduction

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 28,000 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 and Bucher, 2020; Remuzzi and 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 13,800 official deaths in the hospitals and almost 1400 ICU beds (Grasselli et al., 2020a, 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 and 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.

Dynamics of Covid-19 phenomenon

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, 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 and 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 non-pharmaceutical interventions (Bayham and 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.

Simple math models for ICU beds and deaths

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 (Molinari et al., 2020) 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 and 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 (Molinari et al., 2020).

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 h 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).

Qualitative features of ICU beds and deaths models

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 and 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 (R 2). 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 R 2 values with either two or three nines as decimals, almost equal to 1, i.e. .

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 Fig. 1 ).

Fig. 1

Doubling time of ICU patients in Italy and in Lombardy in the very first days of the pandemic. The lower the values the higher the pressure on the health system. Italy shows a time delay of few days respect to Lombardy where the pandemic manifested earlier.

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.

Quantitative features of ICU beds and deaths models

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, a, b, c are adaptive parameters and 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.

Eq. (1) in semilogarithmic coordinates assumes the form:which is the equation of a line whose slope is b, the intercept is , and .

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

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

The logistic curve Eq. (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).

Table 1

Inflection time and halfway time for the logistic and Gompertz curves.

	Logistic	Gompertz
Inflection time	t_infl=ln(b)c	t_infl=ln(b)c
Halfway time	t_halfway=ln(b)c	t_halfway=−ln(ln(2)b)c

Both the logistic and the Gompertz models have the same maximum plateau when t → ∞, .

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 Eq. (6) through the following nonlinear regression procedure (Hosmer et al., 2011):

The minimization of Eq. (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., if the desired percentage is 98%). An easier approach to solve Eq. (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).

(1), (2), (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:

Equally, the EMG model can be written as the product of three terms:where is the multiplying factor, 10 is the exponential term and is the gaussian term. Eq. (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 b, c 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 ().

The first derivative of the exponential function is:

The first derivative of the logistic function is:

The first derivative of the Gompertz function is:

The first derivative of the EMG function is:

Eqs. (10-13) are good when drawn on a graph and compared with daily fluctuations of both ICU and death figures. Indeed, they help understanding from a differential point of view the trend of the phenomenon as well as the point where the derivative changes sign (from positive, increasing trend, to negative, decreasing trend which shows the time when the maximum plateau of ICU beds occurs). Equally, the maximum and minimum values of the derivative correspond to the inflection points (i.e. maximum rate of increase and decrease of the phenomenon respectively).

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 t 0 are known as they are the plateau value and the corresponding time when the maximum of ICU beds is reached, respectively. Consequently, Eqs. (14) and (15) reduce to two adaptive parameters, b, c. Both Eqs. (14) and (15) exhibit the same final null plateau when the pandemic expires.

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

Model identification

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 Eq. (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 Eq. (18). Equally, one can evaluate the either the mean (MAE) or the median (MEDAE) absolute errors via Eqs. (19) and (20) respectively (K. Zhang et al., 2015).

Analysis of the dynamics of ICU patients and deaths in Italy and in Lombardy

Dynamics of ICU patients

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 Eq. (1). This assumption proved true after observing the linear trend of dependent variables (i.e. ICUs and deaths) in semilogarithmic coordinates Eq. (4) and evaluating the determination coefficient.

Fig. 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 R 2 that approximates 1 (i.e. ).

Fig. 2

Linear trend in semilog coordinates of ICU patients in Italy in the first 15 days of pandemic. Cyan line with dots reports real data. Orange line with dots shows extrapolated values. Dashed line is the linear trend line that evidences the exponential growth of the phenomenon.

As time progresses, the phenomenon leaves progressively the purely linear trend in semilog coordinates and moves towards a quadratic behavior (in those same semilog coordinates) that is embodied by an EMG model (see Eqs. (8) and (9)).

Fig. 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. ). 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, based on daily variations of the same variables, which exhibit significant fluctuations. A further attention point is the logarithmic ordinate axis that flattens the fluctuations of high values. The fluctuations appear rather small on the semilog diagram that flattens the natural oscillations of the pandemic.

Fig. 3

Quadratic trend in semilog coordinates of ICU patients in Italy in the day 13–36 interval. Cyan line with dots reports real data. Dashed line is the quadratic trend line i.e. EMG model of the phenomenon.

The logistic and Gompertz models, both direct Eqs. (2) and (3) and reverse Eqs. (14) and (15), are further mathematical representations that allow discriminating among models. Besides the logistic and Gompertz models (direct and reverse), Fig. 4 shows the EMG model that describes both the plateau and the downhill phases. For the sake of clarity, the downhill part of EMG practically 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).

Fig. 4

ICU patients in Lombardy. The uphill trend is described by the logistic and Gompertz models (dashed lines) from day 1 until day 32 (24-mar). The plateau and downhill trends are described by the EMG model from green line) from day 32 (24-mar). The reverse logistic and Gompertz models describe the downhill trend (cyan and blue lines) from the maximum plateau (day 41, 2-apr).

Fig. 5 shows the daily variations of ICU patients in Lombardy in the most exacting period of ICU emergency that includes the maximum plateau. The daily variations of ICU patients highly oscillate and fluctuate, and depend on a number of factors that influence both the input and output terms (i.e. patients given access to ICU, deaths, weaned-from-ventilator patients who are transferred to step-down units, and possibly new installed beds).

Fig. 5

Daily variations of ICU patients in Lombardy. The lines show the first derivatives of the EMG, logistic, and Gompertz models.

The model first derivatives describe the continuous rate of change of model predictions as a function of time and approximate the real daily variations. Positive values of first derivatives mean that the slope of the tangent of the curve is uphill, negative values mean downhill slopes. Consequently, when the first derivative is equal to zero (as for the EMG model in Fig. 5) the maximum plateau of the corresponding model occurs. In Lombardy, the EMG model of Fig. 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 Fig. 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 Fig. 6 ).

Fig. 6

Extrapolation of ICU patients in Lombardy over a long-term horizon. The lines and colors are the same as in Fig. 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.

Fig. 6 shows on a longer time scale the predictions of reverse logistic and reverse Gompertz models on a long-term horizon. One can easily extract valuable information from that diagram by setting monitoring thresholds when the ICU phenomenon falls below a suitable value. According to the model predictions, the number of ICU beds in Lombardy should fall below 10% of the maximum measured value on day 92–96 (23–27 May), which is 50–54 days after the top of ICU maximum allocation (which topped 1381 beds). Similar simulations can be carried out to forecast when the models predict the last 10 patients in a Covid-19 ICU ward. For Lombardy this would happen on 26-Jun and 8-Jul for the reverse logistic and reverse Gompertz models respectively (i.e. 84 and 96 days after the maximum plateau (i.e. at days 126 and 138 respectively), based on 30-Apr real data).

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).

Dynamics of deaths

The description of death dynamics can ground on the same modeling approach based on the logistic, Gompertz, and EMG curves (Eqs. (2), (3), and (8), (9) respectively). For the sake of brevity, the initial exponential growth Eq. (1) will not be reported also because the aforementioned models feature in their formulation an initial part where the curve follows the exponential trend.

Fig. 7 shows the real trend of deaths in Italy since day 1 and the model predictions until the end of June. The Gompertz model is the only reliable one over such long-term horizons.

Fig. 7

Deaths in Italy. Last available real data on 30-Apr.

The consistency between real data and model predictions is quite high. This model allows determining that the inflection point (when the fatalities reach the highest daily increment) occurs at day 38 (30 March) whilst the halfway point is reached at day 44 (5 April). The Gompertz model allows also predicting (as of latest data available on 30 April) that the total expected deaths toll will be 34,037 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.

Fig. 8 shows on a semilogarithmic diagram the trends of both cumulated fatalities and daily changes. After the first days, the consistency of cumulated real values with model predictions is decidedly high. It is worth observing that the discrepancy observed until 7–8 March (days 15–16) is not much significant as it is of few units and can be observed only because the ordinates are on a logarithmic scale. Conversely, the oscillations of daily changes of deaths are manifest and show the complex nature of the pandemic. The maximum value of the first derivative of the Gompertz model (in Fig. 8) corresponds to the inflection point (in Fig. 7) and occurs at day 37 (29 March).

Fig. 8

Deaths in Italy on a semilogarithmic diagram. Comparison of cumulated fatalities and daily variations. Last available real data on 30-Apr.

Very similar inferences on fatalities can be carried out for Lombardy based on what presented for Italy. The trends reported in Fig. 7 and Fig. 8 for Italy are pretty similar to those of Lombardy. Only the absolute values change. For the sake of detail, grounding on the real data available at day 69 (30 April), the Gompertz model identifies that the inflection point occurred at day 35 (27 March), the halfway point at day 40 (1 April) and that 98% of the total fatalities toll (maximum plateau of 15,519 deaths) will be reached at day 98 (29 May).

It is worth observing that the behavior 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% of Italian deaths concentrated in Lombardy although the inhabitants of Lombardy account for less than 17% of the whole national population. Indeed, some southern regions of Italy suffered a minor if not negligible fatalities toll. For instance, in Campania (where Naples is the first city) the fatalities were 359 (as of day 69) with a population of 5.8 million inhabitants. Conversely, Lombardy the most scourged region (with Milan as first city) suffered 13,772 fatalities with 10 million inhabitants, which is a 22-fold higher value if scaled to population.

Discussion

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 Fig. 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).

Fig. 9

Allocation of ICU beds to Covid-19 wards in Lombardy as a function of the number of patients admitted to ICUs in the day 57–66 period. On 27-apr the total number of treated ICU patients was 4269 (1864 dismissed, 1661 deaths, 744 sheltered. Age min 0, max 92, mean 62 (std dev 11), median 63, 25% percentile 56, 75% percentile 69 [y]).

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 (18), (19), (20) (K. Zhang et al., 2015). Table 2 and Table 3 report the KPIs of the ICU models for both Italy and Lombardy.

Table 2

Key performance indicators of different models of the ICU beds dynamics in Italy. RMSE, MAE, and MEDAE are defined in (18), (19), (20).

ITALY	Logistic	Gompertz	EMG	Reverse logistic	Reverse Gompertz
Time interval	Days 1–32	Days 1–32	Days 32–53	Days 41–69	Days 41–69
RMSE	32.1540	38.0560	49.9381	31.9778	29.6996
MAE	28.9102	34.0622	42.0433	25.3397	24.2865
MEDAE	28.9426	29.1963	42.0536	21.6841	18.5969

Table 3

Key performance indicators of different models of the ICU beds dynamics in Lombardy. RMSE, MAE, and MEDAE are defined in (18), (19), (20).

LOMBARDY	Logistic	Gompertz	EMG	Reverse logistic	Reverse Gompertz
Time interval	Days 8–32	Days 8–32	Days 32–69	Days 41–69	Days 41–69
RMSE	18.8038	14.1331	17.0645	16.7249	16.2730
MAE	15.7513	11.4652	13.2487	14.4411	14.1539
MEDAE	14.6162	8.5640	11.3843	15.3614	13.4041

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 and 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 (27,967 and 13,772 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 34,037 for Italy and 15,519 for Lombardy according to the projections of day 69, 30 April).

Table 4

Key performance indicators of different models of the fatalities both in Italy and in Lombardy. RMSE, MAE, and MEDAE are defined in (18), (19), (20).

	ITALY		LOMBARDY
	Logistic	Gompertz	Logistic	Gompertz
Time interval	Days 3–55	Days 3–69	Days 3–55	Days 3–69
RMSE	325.8346	158.9923	154.1105	91.1673
MAE	284.2970	120.8251	130.5953	74.1202
MEDAE	286.4472	68.9499	111.0064	63.8503

Conclusions

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 long-term 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.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: