Italy was the first European country experiencing a substantial outbreak of the coronavirus disease 2019 (COVID-19) and
consequently applying lockdown emergency measures to contain the spread of the epidemic, similarly to the Chinese strategy. In particular, as the number of infections was growing exponentially [1] , the Italian authorities decided to close all schools on March 5th 2020 and introduced community containment lockdown measures starting on March 9th. Spain, France, Germany, United Kingdom (UK) and United States (US) are the countries with the highest number of infections and deaths. All these countries adopted public health measures and, most of them, implemented a complete lockdown, but such measures were introduced at different time points following local outbreaks and with different strengths [1] .
The aim of our work is to model the evolution of the COVID-19 pandemic in different regions, taking into account major confounders such as undetected cases, related to the different case fatality ratio measured and taking into account background death data, and different mobility changes as effect of public health measures.
Assuming that isolation is successfully applied to the positive detected cases, undetected and asymptomatic cases will represent the main source of infections in the near future, potentially determining new outbreaks. This poses a crucial challenge for public health decision makers, namely to determine when and how to relax the control measures implemented so far, taking into account the whole number of estimated infectious individuals and not only the tested cases.
Our model provides an accurate estimation of the undetected and asymptomatic cases, therefore allowing to predict future epidemic peaks and design the most effective way to relax the control measures, in terms of timing and strength in term of mobility restrictions and adoption of personal protective behavior. Moreover, it can be expanded and tailored when more and new region-specific data will become available.
We model the spread of the COVID-19 pandemic in different geographical regions using a modified SEIR compartmental model [2] , fitted using publicly available data from the 2019-nCoV Data Repository by Johns Hopkins CSSE and from the Italian official repository updated as of 24th April [3] In a classical SEIR model the population is divided into four classes: Susceptible, Exposed, Infected, and Recovered. In our specific case, illustrated in Supplementary Methods, we added the classes of Positive tested, Undetected infected and, both for detected and undetected ones, the two possible ending point, Healing and Death ( Figure 1A ). Disease-specific parameters were derived from the recent literature, including a mean incubation time ( −1 ) of 5. We compared available Italian death records data from years 2015-2020 (considering only February and March) and assumed that the increase in deaths registered in 2020 until 28th March in each Italian region is entirely related to the COVID-19 pandemic, even though only a fraction was confirmed as positive. In particular, we estimated that undetected COVID-19 cases in Italy represent 53% of the registered COVID-19 deaths (Table 1) , which is consistent with other scenarios. For instance, on 17th
April 2020 the Chinese government announced that in Hubei province there are 1290 COVID-19-related deaths (occurred from the beginning of the pandemic), which were not previously registered as COVID-19 deaths, and that represent about 40% of the officially registered deaths [9] . This fact supports the use of the Italian estimation for all analyzed countries. Therefore, by combining the background death records data with the measured Naïve Case fatality Rates (nCFRs) up to 10th March, ranging from 2.26% in Germany to 18.27% in Lombardy (Italy), to obtain an aCFR of 1.4% (CI95% 1.2-1.7%) [8] , it is possible to give a reasonable estimate of the proportion of undetected cases in different regions.
We introduced in the model the public health measures, taking into account the effects of mobility changes as well as the potential effects of adopting a personal protective behavior, like frequent hands washing and wearing a mask. The former is modelled with a decreasing logistic function, where the target plateau is the average mobility change reported by Google COVID-19
Community Mobility Report in each geographical area over the last 7 available observations (11th-17th April) and the shape is calibrated knowing the date of the implementation of public health measures (social distancing advisement or school closures) [1] and date of the maximum reduction of the mobility ( Figure 1B-D) .
A B C D The effect of the adoption of a personal protective behavior as frequent handwashing or masks wearing, recommended worldwide [10] , was derived from a metanalysis of studies performed in the context of SARS epidemic [11] , estimating an Odds Ratio (OR) of the use of masks of 0.32 (CI95% 0.25-0.40) and an OR of frequent handwashing of 0.45 (CI95% 0.36-0.57) for the risk of infection. The models are fitted until 17th April leaving the last 7 observed days for the accuracy analysis of models' predictions through Mean Absolute Percentage Error (MAPE) analysis [12] . The model is deterministic, thus the reasonable confidence intervals of fitted estimates were derived considering the CI95% of the aCFR and the CI95% of personal protective measures
ORs derived from literature data [8, 11] Mathematical details of the model, including the equations, parameters calibration and fitting procedure are illustrated in the Supplementary Methods.
In the Italian scenario the undetected proportion of cases ranges from 76% in Lazio to 96% in Sardinia, being 93% for the whole Italy (Table 1, Figure 2A ). Epidemic curves obtained using the calibrated SEIR models are reported in Figure 2B ; the accuracy analysis of forecasting showed a median MAPE of 9.2% [Interquartile Range (IQR) 3.6-14.5%; range 0.9-34.8%; Table 1 ].
The forecasting procedure permits an estimate on 24th April of the whole fraction of infected population, including undetected cases, ranging from 0.4% (CI 0.3-0.5%) for Campania to 13.6% (CI 13.4-13.9%) for Lombardy, being ≈ 5 to 40-fold higher than positive tested cases, as shown in Table 1 and in heatmaps in Figure 2A . This is coherent with recent literature results estimating the COVID-19 prevalence, until 7th April, ranging from 0.35% to 11.2% across different Italian regions [13] .
In the international context the number of supposed undetected cases ranges from 60% in Germany to 92% in UK. The application
of the SEIR models, tested in the whole Italy, France, Germany, UK, Spain, and US was satisfactory, permitting also to show the trajectories of overall cases, including undetected ones ( Figure 3A ). Analyzing the forecasting accuracy in the international (0), transmission rate at time 0;
(0) , basic
reproductive number imputable to positive tested cases; (0) , basic reproductive number imputable to undetected cases;
(0), overall basic reproductive number;
1∕ 2 , fitted recovery notification time; MAPE, mean absolute percentage error; CI, reasonable confidence interval;
∕ , fraction of detected cases derived from undetected ones until 24/4/2020; ∕ , fraction of the overall new cases derived from undetected ones until 24/4/2020; * since 13th April 2020 public data regarding reported recovery cases are not anymore available. April to 24th April are used for judging models' accuracy. The shadow regions represent the reasonable confidence intervals.
The best fit in each region of the recovery notification time (1∕ 2 ), has median value of 63 days (IQR 38-99; range 15-10'000 days), possibly representing different efficiency of recovery reporting in each geographical scenario. The estimation of the basic reproductive number ( 0 ) was derived from the model through the Next Generation Matrix Method [14] . The contribution of positive detected cases or undetected ones for the generation of secondary cases is weighted to obtain the overall reproductive number ( 0 ), representing the overall number of newly cases generated along the infectious time of an affected case ( Table 1 ). The mean 0 related to detected cases ( 0 ) is 7.5 ± 2.1, the one related to undetected ones ( 0 )15.7±4.5, determining an overall weighted 0 value of 15.1± 4.3. Figure 3B ). According to our model, there were 2502 (CI 2251-3205) infectious subjects on 7th April, most of those being undetected (96.4%). We show how the infection dynamics depend on the level of target mobility that is restored, compared to the baseline, and on the adoption of personal protective measures. These two factors can make the difference between a huge and rapid secondary epidemic peak or a more delayed and manageable one ( Figure 3B ).
Since its first discovery as a new human virus on January 8th in Wuhan, China, the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread worldwide, determining the COVID-19 pandemic that up to 6th May 2020 has caused over 3 588 773 confirmed cases and 245 503 death [16] . The main public health issue regarding COVID-19 disease is that 5%
to 10% of the patients need intensive care recovery with artificial ventilation for onset of pneumonia complicated with acute respiratory distress syndrome (ARDS); as a consequence, the number of intensive care beds, which is often limited, can become insufficient [17] .
Despite the proven effectiveness of lockdown measures, as shown in the Chinese epidemic, the challenge is now to assess their efficacy and to plan strategies to progressively relax them [18, 19] . Mathematical models can help to determine the trajectories of the pandemic, considering the confirmed cases and estimating the undetected ones. In particular, the SIR or SEIR compartmental models have the advantage of providing useful insights for health care decision makers and of being easily manageable from a mathematical point of view [2, 20] . Moreover, if adequately calibrated with good literature-based parameters the simplicity of these models can avoid overfitting results [2] . The supposed proportion of undetected cases is reported to be highly variable. Its estimation from case fatality rates was already proposed [1] , but its integration with background death records data can further improve this estimate.
In the available published literature, searching PubMed between Dec 1st 2019 and April 19th 2020, we found several papers modelling the COVID-19 pandemic using SEIR models, most of them studying the evolution of the epidemic in Wuhan or
Hubei or in single geographical context, and proposing different scenarios depending on different transmission rates. However, only few of them take into account a measurable effect of public health intervention (e.g., mobility changes) [21] . The work by Gatto et. al. [21] is of major interest, as they apply a modified SEIR model (SEPIA) using data of geo-located smartphone users to estimate the effect of lockdown intervention measures on the Italian epidemic evolution [22] . Our results in the Italian case-study report the cumulative number of overall affected cases as of March 25th and March 28th equal to 1.46*106 (CI To the best of our knowledge our is the first SEIR compartmental model taking into account the mobility changes measured by Google in order to track the effect of lockdown measures. The decreasing logistic function permitted to solve the system analytically and to integrate the effect of the adoption of personal protective measures (hand washing an masks wearing) derived from literature estimates [11] .
Fitting the recovery reporting time after the isolation of positive tested cases permitted to depict wide differences among countries (IQR 36-95 days) but also in-between Italian regions (IQR 38-99 days); of note the fitted value of 10?000 days for UK is the consequence of the interruption of disclosing publicly recovery data since the beginning of April 2020. Recovery?s reporting time heterogeneity is of main interest as the removal time of recovered individuals is expected to influence the shape of the epidemic curves of active cases before and after relaxing the lockdown measures. Furthermore, as shown in the Hubei case-study ( Figure 3B ), the outputs of our model in term infectious groups (undetected or not already confirmed cases) can be Even though specific and quantitative predictions from a mathematical model need to be taken with caution, some key messages can be derived, as the estimation and the evolution of complex phenomena that are hardly imaginable by the common sense, for instance the hidden dynamic of undetected cases parallel to the dynamic of detected ones [23] . Despite having some limitations, not considering spatial-temporal dynamics and stochastic features, the main advantage of our system is its easy application, that can be even improved by further updating mobility information, nowadays publicly available, and acquiring region-specific background death records data.
None.
None declared.
Because this study involved the use of publicly available aggregate data, approval by a research ethics board was not required.
All data used in this manuscript are publicly available. COVID-19 epidemiological data for Italy are available at https://github.com/pcm-dpc/COVID-19, and for countries at https://github.com/CSSEGISandData/COVID-19. Italian death records data are available at https://www.istat.it/it/archivio/240401#Decessidel2020-0, (downloaded on 11th April 2020).
Google COVID-19 Community Mobility Report are available at https://www.google.com/covid19/mobility/ (downloaded on 23rd April 2020).
Herein, we introduce the mathematical details of the epidemiological model presented in this study. It is a modified version of the well-known SEIR (susceptible-exposed-infected-recovered) compartmental model, accounting both for specific features of the COVID-19 infection and for the impact of public health protective measures on the dynamic of the transmission rate. The population is divided into nine classes representing the possible states of the individual in relation to the diseases:
• Susceptible (S), it is the class containing individuals that can potentially contract the disease;
• Exposed (E), includes individuals that have been exposed to disease (infected) but are not yet capable of transmission (not yet infectious);
• Infected (I), includes individuals that are both infected and infectious;
• Positive Tested (P), contains individuals recognized as positive by the health national system and assumed to be well isolated (either hospitalized or quarantined) and not infectious anymore;
• Recovered (R), includes subjects healed and whose recovery is notified by the health system;
• Dead (D);
• Hidden Infected (̃ ), includes subjects affected by the disease (either asymptomatic or symptomatic) but not officially tested and registered as infected by the health national system;
• Hidden Recovered (̃ ), includes individuals recovered form the disease, but previously not registered as positive by the health national system;
• Hidden Dead (̃ ) includes individuals not officially registered by the health national system, who were infected and died;
The mathematical model is described by the following set of nine ordinary differential equations (ODEs), representing the time rates of change of the introduced variables in a population of size :
In the equations above, the number of Susceptible individuals ( ) decreases over time due to contacts with the infectious population ( +̃ ). The force of the infection is driven by the transmission rate parameter defined as the probability of disease transmission in a single contact times the average number of contacts per unit time per person.All the other parameters represent the rates at which individuals pass from a compartment to the subsequent one; they can be seen as the inverse of the number of days on average spent in the former class. Specifically, for an Exposed individual, is the rate of change from the class to either the Infected ( ) or the hidden Infected (̃ ) ones. 1 is the rate of change from the Infected class to the class. The parameter 2 and̃ represent, respectively, the recovery rate of a hidden infected or the recovery's notification rate by the health system of a positive tested subject. Similarly, and̃ are the mortality rates.
As said before, the number of days spent in each class is given by the inverse of the corresponding rate: the incubation period = −1 , the infectious period of visible Infected 1 = −1 1 , the average diagnosis time for an Exposed individual to being Positive tested 2 = −1 2 and the average dead period of both Positive tested = −1 and hidden Infectious ̃ =̃ −1 . The recovery periods of the hidden infectious is given by ̃ =̃ −1 . We stress the fact that the recovery period of the visible infected individual is given by = 1 + 2 −1 .
As regards the other parameters, and̃ represents the probability for an individual respectively in or̃ to end in the dead class or̃ , whereas the parameter ℎ is the fraction of that becomes Hidden Infectious̃ . In Figure 1A of the main text we report the compartmental diagram illustrating the model, as well as the introduced parameters. The ODEs system is supplemented by the initial conditions: Note that the variables ,̃ , and̃ are cumulative and that the system has the mass conservation property: since
the total population is constant over time: In our model, the average number of contacts is described by considering an initial constant value 0 reduced over time by a decreasing function ( ) with (0) = 1:
Similarly, for the probability of transmission we have:
The modelling of the mobility time dependent part ( ) is given by a decreasing logistic function with a plateau defined as:
where the normalization constant granting (0) = 1 is:
In Eq.(10), the parameter is the abscissa of the inflection point and the parameter is related to the first derivative calculated in that point, that is the steepness of the function. The parameter is the plateau of the function:
The parameter is calibrated using the Google COVID-19 Community Mobility Report data by computing the average mobility changes for the period 11th-17th April for the four categories: retail and recreation, grocery and pharmacy, transit stations and workplaces mobility. As regards the second effect of personal protective measures, we do not assume an explicit behaviors of the ( ) function. Instead, we incorporated its effect as a reduction of the plateau . The consequent transformation of the mobility decreasing function can be described by considering the aggregate product of both the time dependent terms as:
Note that, despite this transformation, the second derivative is invariant and the abscissa of the inflection point is still the same of ( ). In order to fit the function with the actual data, we note that the derivative in the inflection point is equal to The initial date of the transition can be written as a function of and :
and similarly, the final date of the transition is:
In Figure 1B 
The parameter 0 is defined as the âĂŸaverage number of secondary cases arising from a typical primary case in an entirely susceptible population. We exploit the next generation matrix method (Van den Driessche 