Testing and isolation to prevent overloaded healthcare facilities and reduce death rates in the SARS-CoV-2 pandemic in Italy.

Background: During the first wave of COVID-19, hospital and intensive care unit beds got overwhelmed in Italy leading to an increased death burden. Based on data from Italian regions, we disentangled the impact of various factors contributing to the bottleneck situation of healthcare facilities, not well addressed in classical SEIR-like models. A particular emphasis was set on the undetected fraction (dark figure), on the dynamically changing hospital capacity, and on different testing, contact tracing, quarantine strategies.
Methods: We first estimated the dark figure for different Italian regions. Using parameter estimates from literature and, alternatively, with parameters derived from a fit to the initial phase of COVID-19 spread, the model was optimized to fit data (infected, hospitalized, ICU, dead) published by the Italian Civil Protection.
Results: We show that testing influenced the infection dynamics by isolation of newly detected cases and subsequent interruption of infection chains. The time-varying reproduction number (R t) in high testing regions decreased to <1 earlier compared to the low testing regions. While an early test and isolate (TI) scenario resulted in up to ~31% peak reduction of hospital occupancy, the late TI scenario resulted in an overwhelmed healthcare system. Conclusions: An early TI strategy would have decreased the overall hospital usage drastically and, hence, death toll (∼34% reduction in Lombardia) and could have mitigated the lack of healthcare facilities in the course of the pandemic, but it would not have kept the hospitalization amount within the pre-pandemic hospital limit.

Introduction

The COVID-19 outbreak created a worldwide pandemic causing more than 4,000,000 deaths and over 190 million total cases worldwide as of July 2021[1]. Many countries implemented non-pharmaceutical interventions (NPIs), which were effective in reducing virus spreading. This was further supported by social distancing, mask duty, and hygiene measures. Though different models of NPIs and their implementation methods have been proposed, their impact and effectiveness on disease dynamics are under scrutiny and remain a matter of global discussion[2-6]. Singapore and Hong Kong were able to contain the virus by aggressive testing[7], while South Korea adopted a trace, test, treatment strategy[8]. In a different but similarly effective approach, Japan averted the risk of contagion by isolating the whole contact clusters and by heavily relying on the self-awareness and discipline of the population[9].

The COVID-19 outbreak originated in Wuhan, People’s Republic of China, in early December 2019. Within two months, it erupted and unfolded with tremendous speed in Italy, which became the European epicenter of disease spreading, forcing the government to impose a lockdown on March 9th, 2020. On March 19th, 3405 people had already died in Italy, thereby surpassing China, while 41035 people were diagnosed as COVID-19 positive. This induced Italy to shut down all non-essential businesses on March 21st. Despite the strict measures applied, in Lombardia alone a total of 28545 symptomatic people were infected by April 8th, accounting for 12976 hospital admissions, followed by Emilia Romagna (4130), Piemonte (3196), and Veneto (1839)[10]. These large numbers led to the complete collapse of the healthcare system within a few weeks of the first detection of COVID-19 cases, most notably in Lombardia where even funeral homes had been overwhelmed and were incapable of responding in a reasonable time[11]. Even though the state expanded the hospital and intensive care unit (ICU) capacities, it could not prevent the bottleneck situation of the healthcare system and presumably caused a large number of deaths for a prolonged period.

Many factors aggravated the COVID-19 situation in Italy, among which the distinct demographic structure of Italy with nearly 23% of the population of age 65 years or older[12], larger household size, and the prevalence of three-generation households compared to Germany[13] as well as limited hospital and ICU capacities. At the beginning of the pandemic, Italy focused on testing symptomatic patients only, which resulted in a large proportion of positive tests and high case fatality rates (CFR) compared to other countries[14]. A large proportion of cases remained undetected, which became a major driver of new infections. A different study estimates that in Italy the actual number of total infections was around 30-fold higher than reported, while for Germany it was less than ten-fold[15] (data up to March 17th 2020).

Compartmental models have been widely used to describe the dynamics of epidemics, for example, SIR models[16] that consider three compartments, namely susceptible, infected, and recovered, or more complex SEIR models[17,18] that take susceptible, exposed, infectious, and recovered compartments into account. Typically, these models either exclude the undetected index cases[4,17,18], or ignore their dynamic nature[19], and structurally these models are not developed to address the load on the healthcare system. Besides these epidemic models, simple algorithms exist in the literature for estimating the time-varying reproduction number and have been widely used in the context of many infectious diseases (e.g., measles, H1N1 swine flu, polio, etc.)[20]. Several studies[21,22] estimated the undetected case number in Italy, but its dynamics in the context of different testing strategies and implications on the healthcare system were not considered. Even though the general compartmental SIR and SEIR type models are useful in inferring epidemic spread and public health interventions, we needed to introduce additional compartments to investigate how the pandemic is shaped by several influential factors (e.g., dark numbers, regional testing strategies, hospital beds); for instance, we included a specific compartment for infected undetected cases (IX) to analyse the impact of the region-wise undetected cases upon the evolution of the . Similarly, hospital (HU and HR) and ICU (UD and UR) compartments were introduced to monitor the load on the healthcare system.

Additionally, the absence of reliable symptom onset data and heterogeneity in the actual infectious period among the asymptomatic, pre-symptomatic, and symptomatic individuals require a more complex model that not only can accurately portray the dynamics of COVID-19 spread but also can disentangle the impact of intertwined factors like the variation of undetected components, limited and changing hospital and ICU availability.

Existing modeling studies that analyze the COVID-19 situation in Italy[19,23] or other regions[4,24-26], in general, did not address some fundamental aspects of the ongoing pandemic like temporal dynamics of undetected infections, the benefits of a high testing and isolation strategy, or the impact of a limited and dynamically changing healthcare capacity on the lives lost. In this study, we address the bottleneck situation of the healthcare facility, the benefits of extending hospital infrastructure, and the impact of an early testing and isolation strategy on the healthcare system with a COVID-19-specific mathematical model. To evaluate the COVID-19 situation in Italy in a realistic framework, we first estimated the undetected fraction (dark figure) of infections across different regions of Italy. We used this information to determine the parameters of the model and showed that our model is structurally identifiable. We studied the influence of the dark figure and implemented NPIs on the time-dependent reproduction number, . With data about regional hospital and ICU bed capacities, we estimated that an extra 25% of people died in Lombardia due to the overwhelmed healthcare system. We investigated the impact of early testing strategy and, alternatively, of contact tracing combined with quarantine (~10 fold more isolation of infected) policy in the setting of elevated hospital capacity as it currently stands. This strategy would have reduced the death toll by 20% to 50%.

Methods

SECIRD-model

To understand the impact of potential aggravating factors, namely infections from undetected index cases, early vs late testing strategy, and limited healthcare facilities on disease progression, we developed a COVID-19-specific SECIRD-model parametrized for Italy. The SECIRD-model distinguishes healthy individuals without immune memory of COVID-19 (susceptible, S), infected individuals without symptoms but not yet infectious (exposed, E), and infected individuals without symptoms who are infectious (carrier, CI, CR). The carriers are distinguished into a fraction α of asymptomatic (CR) and (1 − α) of pre-symptomatic infected (CI). The latter are categorized into a fraction μ of detected symptomatic (IH and IR) and (1 − μ) of undetected mild-symptomatic (IX). Out of the CI, a fraction ρ gets hospitalized (IH), and (1 − ρ) become symptomatic but recover without hospitalization (IR). Further, compartments for hospitalization (H) and intensive care units (U) were introduced to monitor the load on the healthcare system. A fraction ϑ of H requires treatment in ICU (HU) while a fraction (1 − ϑ) recovers from hospital without ICU treatment (HR). δ and (1 − δ) represent the fraction of patients in ICU who subsequently die (UD) or recover (UR). The compartment (R) consists of patients recovered from different infection states. The Reference Model (Fig. 1; equations are in the Supplementary Methods 1.3) was solved with parameters in Table 1.

Fig. 1

Model schemes.

a The Reference Model distinguishes healthy individuals with no immune memory of COVID-19 (susceptible, S), infected individuals without symptoms but not yet infectious (exposed, E), infected individuals without symptoms who are infectious (carrier, CR,I, asymptomatic and pre-symptomatic, respectively), infected (IX,H,R), hospitalized (HU,R) and Intensive Care Units (ICU) (UD,R) patients, dead (D) and recovered (RX,Z), who are assumed immune against reinfection. This scheme also applies to the Asymptomatic Model. b The Capacity Model is a modified branch of the Reference Model to investigate the impact of limited hospital and ICU access onto the death toll. and are steep exponential functions diverting the flux from IH and HU to D, respectively, when hospital and ICU occupancy reached their respective current capacities and . (c) The Testing Model is a modified branch of the Reference Model used to evaluate the impact of increasing case detection and isolation onto infection dynamics; IXD and IX describe newly detected and undetected cases, respectively. Rx with x ∈ [2,…,9] are per day transition rates between different states. Behavioral parameters (ρ, ϑ, δ, R1 and R10) are subject to contingent factors, like Non-Pharmaceutical Interventions (NPIs), self-awareness, availability of hospital beds, etc., and, hence, are functions of time.

Table 1

Parameter ranges used in the Reference Model: determination of the boundaries for literature-based parameter set was based on the interpretation of the values given in the references[27,28,32].

Parameter	Comments/References	Description	Parameter ranges from literature
			Min	Max
R1	Time-dependent	transmission probability of COVID-19 per each contact made with an infectious person (CI, CR,IX, IR, and IH in the model)
R2	[54,64–66]	the inverse of R2 represents latent, non-infectious period following the transmission of COVID-19. 1/R2 = 5.2 − 1/R3; median incubation period is 5.2 days
R3	[54,64–66]	the inverse of R3 represents the pre-symptomatic infectious period	14.2	25.2
R4	[67–69]	the inverse of R4 represents the infectious period for the mild symptomatic cases without requiring hospitalization (including the undetected symptomatic people (IX))	114	17
R5	[27,28,70]	the inverse of R5 represents the duration for which the hospitalized cases stay in general hospital care before discharge without requiring further intensive care	116	15
R6	[28,71]	the inverse of R6 represents the duration a patient stays at home before hospitalization	17	0.9
R7	[28,70,71]	the inverse of R7 represents the time spent in general hospital care before admission to ICU	13.5	1
R8	[27,72]	the inverse of R8 represents the time spent in ICU before recovery	116	13
R9		the inverse of R9 represents the duration for which the asymptomatic cases remained infectious following their latent non-infectious period	1R_9=1R_3+0.5×1R_4
R10	Time-dependent[28,71,73]	the inverse of R10 represents the time spent in ICU before dying	110	0.9
α	fixed,[36–38]	undocumented asymptomatic fraction	0.4	0.4
β	Assumed	the risk of infection from the registered and quarantined (IH+IR) patients	0.05	0.25
ρ	Time-dependent	the fraction of documented infections that require hospitalization	0.01	0.9
ϑ	Time-dependent[74–76]	the fraction of hospitalized patients that require further intensive treatment	0.01	0.7
δ	Time-dependent[74–76]	the fraction of ICU patients that have fatal outcome	0.3	0.9
μ¯	this fraction represent the total undocumented infection including the asymptomatic cases, estimated through MLE method of the Bayesian framework
μ		documented symptomatic fraction	μ=1−μ¯1−α

Initial condition

Italian regions started documenting epidemiological data at different dates, at the earliest February 24th, 2020. As we considered a fixed incubation period of 5.2 days in our model, we assumed that at minimum, the first entry in the dataset (number of total cases) was the exposed number 5.2 days earlier. In addition to the documented infection, we calculated undocumented cases based on the estimated region-wise dark number by a Bayesian MCMC framework (Supplementary Methods 1.2). The sum of documented and undocumented cases was the initial exposed population. We used regional population as the initial susceptible population. For all other compartments, we started the simulation with zero. The simulation began from −5.2 days with the aforementioned initial conditions, and undetected cases were split between asymptomatic and symptomatic undetected cases according to the parameters used for those compartments.

Parameterization

We distinguished physiological and behavioral model parameters. Physiological parameters depend on the nature of the virus (Rx, x = 2, ..., 9) and remain unchanged throughout the analysis of the pandemic. We first determined the range of physiological values for each parameter from the literature[27,28] (Table 1) and then estimated parameters’ values from a fit to the exponential growth of case numbers (infection, hospitalized, ICUs, and death cases) during the first two weeks of the pandemic. In total, 56 data points (14 daily data for those four observables) were used to estimate the physiological parameters. This initial phase was not yet affected by NPIs, public awareness, or an overwhelmed healthcare system and, thus, reflects viral properties. Some of the physiological parameters may be internally linked. For instance, hospitalization and ICU cases increase with infection cases, and disentangling those internal relations is difficult with limited data availability. It is likely that the best combination of parameters contains those internal relations. We assumed that the virus variant remained the same during the investigation period, therefore we kept physiological parameters constant throughout. However, environmental factors, testing policy, interventions, public behavior, self-isolation, hospitalization, etc. might have distorted such relations (e.g., the rate of increase in hospital and ICU cases as infection increases differ in different phases of the pandemic) and substantially altered the disease dynamics by impacting the transmission probability, dark number, hospitalizations, and death rate. These contingent factors affect the behavioral parameters (ρ,ϑ,δ, R1, R10). We estimated the behavioral model parameters by minimizing the sum of squared differences between the observed data (active infections, hospitalized, ICU patients, and death numbers (Italy Data on Coronavirus 2020[29])) and model simulations using Matlab’s nonlinear least-squares optimizer. This procedure was repeated separately for each region in Italy in moving time windows of 7 days to account for local specifics and temporal changes in disease transmission. This moving-window technique with the size of a calendar week reduces periodic fluctuations that are an artifact of the unequal distribution of tests among the weekdays.

Perturbation and parameter identifiability

To generate the standard deviation for , we perturbed the behavioral parameters (ρ, ϑ, δ, R1, R10) 10% of their optimized value and sampled uniformly within this range such that the total parameter variation, κ, defined as [30], remains within 10% of its reference value. k, and L represent the parameters of the altered system, the reference system and the total number of parameters, respectively. We generated dynamics for 100 perturbed parameter sets for the statistical analysis.

We addressed parameter identifiability in two ways: Structural identifiability based on synthetic outbreak data and practical identifiability based on real data. In the first method, we randomly sampled parameters within a range specified in Table 1 and by using random initial conditions of model state variables. Then we used the resulting dynamics of the state variables as model observables and checked for a unique solution in the parameter space. We repeated this procedure 100 times (Supplementary Fig. 1 for a typical result). In the second method, we considered nationwide Italy data for the period February 24th to May 22nd, 2020, and fixed the physiological parameters as described in the Parameterization section. As we are estimating behavioral parameters (ρ, ϑ, δ, R1, R10) by considering a moving time windows of 7 days, we checked practical identifiability of these parameters in each time window. We found that the parameters are identifiable in more than 75% of the cases (Supplementary Fig. 2 for a typical results when all parameters are identifiable; and the Supplementary Methods 1.1 for more details).

Basic reproduction number

The basic reproduction number is defined as the expected number of secondary infections produced by a single infection in a population where everyone (assuming no immune memory) is susceptible[31] and reflects the transmission potential of a disease. For COVID-19, the dynamics of the pandemic was influenced by several factors, like, the self-awareness in the community, interventions and policies implemented by the authorities and immunisation of the population. Therefore, the time-dependent reproduction number that describes the expected number of secondary cases per infected person at a given time of the epidemic, is a more practically useful quantity to understand the impact of interventions, behavioral changes, seasonal effects, etc. on the disease dynamics[20,32].

In multi-compartmental epidemic models, can be derived with the next generation matrix method[33-35], where the Jacobian matrix consists of two factors, rate of appearance of new infections into the infection compartment () and transfer of infected into other compartments (). The elements Gij of represent the expected number of secondary infections in compartment i caused by a single infected individual of compartment j. The reproduction number is given by the dominant eigenvalue of G (the derivation of is provided in the Supplementary Methods 1.4):where N0 is the total population and S0 is the susceptible population, both at the start of the pandemic (parameters are listed in Table 1). was calculated using the parameter set estimated by the initial fit that considers only the first two weeks of data points as described in the Parameterization section. An ensemble of parameter sets (as described in the Perturbation and parameter identifiability section) was used to calculate the standard deviation in the .

In order to understand the impact of awareness in the population, NPIs and policies implemented by the authority upon the development of the time-varying reproduction number [20], we fitted the model parameters to data in shifting time windows of one week. This approach has two advantages: first, the reproduction number is determined as a time-dependent variable and thus reflects the impact of NPIs on the infection dynamics; second, the moving-window dampens sudden jumps in the data because of reporting delays. In each time window, a best fit of the model parameters was found based on the cost function value (squared difference between data and simulation). In the next window, the fitting was repeated with initial conditions given by the model state in the previous time window. As described in the Perturbation and parameter identifiability section, an ensemble of perturbed parameter sets was used to calculate the standard deviation in the .

in time window k reads[32]:where ρ(t) denotes the hospitalized fraction of identified symptomatic cases in the kth time window. Data analysis of the clinical state of all infected cases (up to June 22nd) by the Istituto Superiore di Sanità (ISS) showed ~30% asymptomatic cases, with increasing tendency [36-38]. In a study performed in Vo’ Euganeo, Veneto, the percentage of asymptomatic cases was found to be in the range of 40%[39]. We set the asymptomatic fraction to α = 0.4. The fraction of undetected cases (Estimation of undetected cases section in the Supplementary Methods 1.2 and Table 2) is by definition:

Table 2

Estimation of the total number of infections, the Infection Rate (IR), the Infection fatality rate (IFR)1.

Areas	IFR in % (95% CI)	Estimated total Infections (Undetected %)	IR in % (95% CI)	CFR in %	Detected Infections
Italy	1.58 (1.04–1.84)	2627807 (93.73%)	4.37 (3.8–6.64)	13.11	165155
Emilia Romagna	1.84 (1.03–2.24)	252985 (91.69%)	5.79 (4.84–10.22)	13.26	21029
Liguria	2.08 (1.15–2.6)	85924 (93.09%)	5.63 (4.57–10.01)	13.6	5936
Lombardia	1.66 (1.03–1.9)	1390759 (95.53%)	13.83 (12.16–22.19)	18.3	62153
Marche	1.88 (0.88–2.47)	58555 (90.62%)	3.93 (3.05–8.11)	13.56	5503
Piemonte	1.73 (0.78–2.12)	258792 (92.94%)	6.1 (5.06–13.4)	11.05	18229
Toscana	1.63 (0.69–2.36)	62671 (87.77%)	1.43 (0.99–3)	7.25	7666
Valle d’Aosta	1.54 (0.73–2.34)	9785 (90.19%)	9.74 (6.4–17.94)	12.63	958
Veneto	1.3 (0.57–1.71)	141466 (89.67%)	2.77 (2.19–6.09)	6.43	14624

1Based on the data provided by ISTAT up to April 15th[42,43]. Age specific IFRs are reported in Supplementary Fig. 5.

Asymptomatic Model

In the Asymptomatic Model, all symptomatic cases are detected, i.e., μ = 1. We compared the results from this model with those from the Reference Model to understand, in an ideal situation, the implication of detecting all symptomatic cases for the pandemic development.

Testing Model

In order to understand the influence of extra testing on infection dynamics, we adopted a model where a fraction of the undetected infected cases (IX) is detected (IXD) via testing and hence, contained. The newly detected infected (IXD) contribute to new infections with a frequency reduced by a factor β (β < 1) but the infectious period remains unaltered (1/R4) (Testing Model in Figs. 1 and 2). Li et al.[40] have demonstrated, in the context of COVID-19 transmission in China, that strict control measures (travel restrictions, enhanced testing, self-quarantine, contact precautions, etc.) helped in improving the fraction of all documented infections from 14% to 65%, ~4−5 fold. For Italy, we estimated that 90% of the infections remained undetected (Table 2). We assumed that enhanced testing reduced this dark figure by daily 2% until reaching 60% (assuming similar efficiency as in China, i.e., documented infections increasing 10% to 40%). We introduced a time-dependent fraction of undetected infections, which, starting from , was decreased daily by steps of 2% down to 60%. The asymptomatic fraction was fixed as in the Reference Model (α = 0.4). The undetected portion of symptomatic is instead modified so that the fraction of undetected cases, μ1(t), and the fraction of newly detected cases, μ2(t), satisfies:with

Fig. 2

Flowchart of the study design including features and purposes of the SECIRD models.

The flowchart illustrates the steps followed to obtain the results. Each solid arrow points to the result obtained through the step from which the arrow starts, while each dotted arrow links the input to the step where that input was used.

The parameters obtained by fitting the data with the Reference Model were transferred into the Testing Model. This maintains the compartmental flow of the Reference Model and thus ensures that the result reflects the sole effect of isolating a fraction of undetected infections. Correspondingly, becomes:

Capacity Model

To estimate the impact of capacity limitations of the healthcare system, we implemented time-varying capacity constraints on the hospital () and ICU () accessibility (Capacity Model in Figs. 1 and 2), using available data for the number of hospital and ICU beds in the different regions. Table 3 reports the pre-pandemic capacity and the increased capacity, specifically allocated to COVID-19 patients, together with the date of accomplished installation. Some regions doubled their capacity and, presumably, this extension of infrastructure has been implemented in a step-wise manner. We assumed a linear increase of the hospital () and ICU () capacity from three days before exhaustion until reaching the maximum capacity on the date of accomplished installation. This new capacity was available thereafter. The exhaustion date was determined from the data and refers to the day at which the number of hospitalized and ICU patients became larger than the initial capacity. Before the pandemic, 85% of the hospital beds and 50% of the ICU beds were occupied[41]. In the Capacity Model, 15% and 50% of the pre-pandemic total capacity (Table 3) was considered as the baseline capacity of hospital and ICU beds, i.e., the starting values of and , respectively.

Table 3

Hospital bed and ICU capacity before and in the course of the pandemic[41,77]1,2.

Regions	ICU	Beds	Added ICU	Date ICU	Added beds	Date beds
Abruzzo	109	4410	67	31/03/2020	537	23/04/2020
Basilicata	49	1861	24	17/03/2020	139	17/03/2020
Calabria	153	5739	60	11/04/2020	126	11/04/2020
Campania	506	17977	104	11/04/2020	773	14/04/2020
Emilia Romagna	449	17295	259	24/03/2020	2189	24/03/2020
Friuli Venezia Giulia	127	4333	102	02/04/2020	358	08/05/2020
Lazio	557	20817	323	24/04/2020	1527	21/04/2020
Liguria	186	5690	127	07/04/2020	1241	01/04/2020
Lombardia	859	37767	939	03/04/2020	11673	12/04/2020
Marche	115	5183	132	31/03/2020	638	06/04/2020
Molise	31	1225	12	28/03/2020	31	07/04/2020
Piemonte	317	16313	500	08/03/2020	4451	16/04/2020
Puglia	302	12531	297	11/04/2020	1027	26/04/2020
Sardegna	123	5739	40	14/04/2020	92	07/04/2020
Sicilia	392	15821	312	23/04/2020	1632	04/05/2020
Toscana	377	12021	247	06/04/2020	1350	05/04/2020
Umbria	70	3259	35	25/03/2020	131	11/04/2020
Valle d’Aosta	12	481	25	03/04/2020	262	03/04/2020
Veneto	487	17512	331	17/03/2020	1910	17/03/2020
Bolzano (AP)2	40	2047	66	16/04/2020	442	03/04/2020
Trento (AP)2	32	2113	70	02/04/2020	382	07/04/2020

1ICU and normal Beds represent the pre-pandemic total beds. In the simulation we used 50% of ICU and 15% of normal beds as baseline capacity. Added ICU and Added beds represent increased allocation specifically for COVID-19 patients. Date ICU and Date beds is the date when the additional beds and ICUs were in place.

2AP autonomous province.

Upon reaching the capacity limit, the influx should be stopped until a vacant bed is available. This could be achieved by introducing a Heaviside step function or any other piecewise method, but this type of function introduces discontinuities and makes solving the ODEs computationally demanding and error-prone. Here, we introduced two functions (f and f) that behave like a step function but are continuous in nature. f and f return 1 as long as there are free hospital or ICU beds and 0 otherwise. In the Capacity Model, we introduced the compartment ID, to which the flux from IH is directed when the hospital capacity is reached. Then, in a sharp transition, patient access to the hospital or ICU is reduced by the fractions and , respectively:

Both factors increase fatal outcomes of infections when hospital and ICU capacities are reached. We assumed that inaccessibility of hospital or ICU leads to faster and more frequent death. In particular, when the ICU capacity is reached, people in the hospital compartment (HU) die after 1/R7 days which is faster than via the hospital-ICU-dead route ((1/R7 + 1/R10)). Similarly, when the hospital capacity is reached, people in the infected compartment (IH) die after 1/R6 + 1/R7 days, satisfying 1/R6 < 1/R6 + 1/R7 < 1/R6 + 1/R7 + 1/R10 (see Reference Model in the Supplementary Methods 1.3).

The MaxCap Model is defined by the same equations as the Capacity Model, but the parameters, and were set to the maximum hospital and ICU capacity, respectively, from the beginning of the simulations.

Data and code

Italy COVID-19 data of infected cases, hospitalized and ICU patients, and death numbers were provided by the Protezione Civile Italiana[29]. Demographic and mortality data used to estimate IFR, are available from the Italian Institute of Statistics (ISTAT) website[42,43]. ISTAT collects mortality data from the Italian National register office for the resident population (ANPR). An automated method was implemented, and parameter estimation was carried out in Matlab 2019b[44] with a combination of the Data2Dynamics framework[45]. The code is available at https://github.com/arnabbandyopadhyay/COVID-19-in-Italy, and has been archived on Zenodo at ref. [46]. For the Bayesian estimation of COVID-19 IFR of Italian regions, see Supplementary Methods 1.2.