Designing social distancing policies for the COVID-19 pandemic: A probabilistic model predictive control approach

1.   Introduction

The COVID-19 pandemic emerged in China in early December of 2019, two years ago sweeping the globe in a short period of time. During the first wave, by May 20, 2020, a total of more than 4,900,000 cases in 188 countries were reported while the death toll was 320,000 ^[1]. Today, the death toll worldwide has passed the 5.5 million people and with the appearance of the omicron variant the total infected cases have exceeded 343 million. After China, Italy became the global epicenter of the pandemic from March to mid-April, listing hundreds of deaths each day and surpassing soon the fatalities of China. By May 19, the number of confirmed cases exceeded 225,000, while the death toll passed 32,000. Of all afflicted regions, Lombardy-the nation's economic center-suffered the most; more than half of the deaths counted in Italy originated from Lombardy. The first confirmed case in the region was reported on February 21, 2020, and on March 8, the government announced a country-wide lockdown which was expanded on March 21 to include even stricter isolation measures such as the prohibition of mobility and closure of all non-essential business activities. Italy had indeed transformed into a 'red zone'. An abrupt increase in the deaths in Lombardy was mostly due to the saturation of ICUs ^[2]. The situation in Lombardy was dramatic with hundreds of deaths every day. Regarding mobility, according to Google Community Mobility Reports, the population mobility in the region during various activities was reduced as much as 80% compared to the pre-pandemic baseline activity ^[3]. As reported in a recent study ^[4] traffic flows decreased up to 55%, during the period March 1–May 3, 2020 with a peak of -87% during the Easter holiday period while distances of motorized travels decreased by one third, and in general transit demand decreased by 50%. By the end of March 2020, as the disease continued to spread rapidly, the fear of imminent collapse of the healthcare systems around the world led more than 100 countries to implement complete or partial lockdowns. Since the end of April 2020, due to the severe isolation measures, the number of daily new cases and deaths declined in many European and Asia-Pacific countries and governments began to introduce a phased reopening of the activities.

While in most of the countries there was a decline in the confirmed cases during the summer period of 2020, a second even more intense wave of the pandemic emerged during Autumn of 2020-Winter 2021. The number of dead people due to COVID-19 at the second wave surpassed 900 per day in France in mid November, 500 in Spain at the end of November and 700 in early February, 1200 in Germany at the end of December 2020, 1500 in UK in mid of January 2021, and 730 in Italy in early December 2020 ^[1]. In early November 2020, in a race against time there were also the first announcements of vaccines that have been tested successfully in the first interim analysis from phase three. However, one year later, due to the heterogeneity and the generally still low rates of vaccinations worldwide, there are concerns that the COVID-19 pandemic will become seasonal as the influenza. Thus, until a global heard immunity is achieved, the systematic scheduling of social distancing measures will still be on the table of public health policies due to the relatively high mortality ratio of the disease among the elderly and vulnerable population.

At the beginning of the pandemic efforts were focused on the development of mathematical models in an attempt to assess the dynamics for forecasting purposes, as well as to estimate epidemiological parameters that cannot be computed directly based on clinical data. These ranged from simple exponential growth models (see e.g. ^[5] for the case of China, and ^[6] for the case of Italy) and compartmental SIR-like models (see e.g. ^[7] for Wuhan, China, ^[8] for Hubei, China, ^[9] for various cities in China and ^[10] for Lombardy, Italy) to metapopulation models ^[11] including the modelling of the influence of travel restrictions and other control measures in reducing the spread (^[11]). Today, mathematical modelling efforts are focused more on the use of such mathematical models to design in a systematic way efficient coronavirus spreading control policies that account for both the effect of vaccinations and the impact of mobility policies on the COVID-19 transmission. Simple control approaches to COVID-19 involve the use of feedback control and on/off control to simulate the implementation of mobility policies ^[12]. More advanced control approaches such as Model predictive control (MPC) have also been applied, using specific epidemiological models as constraints. For example, Kohler et al. proposed a robust MPC using a compartmental model applied to Germany ^[13]; similar approaches have been used to investigate regions in China ^[14], Hungary ^[15], Italy ^[16] and the United States ^[17]. Finally, the effect of regional interventions seen as an interacting network between regions has also been proposed ^[18].

In this work, we present a probabilistic model predictive control (PMPC) approach constrained to a compartmental model extending beyond the conventional SEIR approach to account for hospitalization and quarantine and importantly for the uncertainty in the level of asymptomatic cases in the population. In particular, the model proposed by ^[10] is further developed, extending it to include hospitalized and quarantined cases. In that work, a model was presented that includes a compartment for the asymptomatic cases and proposed an optimization algorithm to assess the uncertainty that characterizes the actual number of infected cases in the total population, which is mainly due to the large percentage of asymptomatic or mildly symptomatic cases and the uncertainty regarding DAY-ZERO of the outbreak. The knowledge of both these quantities is crucial to assess the stage and dynamics of the epidemic, especially during the first growth period. By doing so, the proposed model was able to predict quite well the evolution of the dynamics of the pandemic in Lombardy two months ahead of time (see ^[10]). Here, we extend that model by incorporating two more compartments, namely the hospitalized and the quarantined cases as well as the effect of social mixing, i.e. the effect of various social activities in the potential transmission of infectious diseases as has been reported for Italy and other countries in ^[19]. The parametrization of the hospitalized compartment was based on recently-published data for the hospital-fatality risk and lengths of stay in the hospitals of Lombardy in the first wave of the pandemic ^[20]. Based on this model, we formulate and solve a PMPC problem that exploits the use of the mobility data provided by Google Community Reports for the Lombardy region to study (retrospectively) what if scenarios of the optimal control action policies, and comparing them with the control actions taken. We note that such retrospective analyses are very important as also pointed in a recently published study of the Italian exit strategy from COVID-19 lockdown ^[21], since they fill gaps in our knowledge of the various aspects of the pandemic and can help design better risk mitigation policies for combating (more) efficiently analogous situations in the future. For example, Foechesato et al. ^[22] performed a retrospective analysis for the estimation of the undetected viral circulation and the day-zero of the COVID-19 pandemic in Italy, while Osborn et al. ^[23] performed a retrospective analysis to investigate the effect of different facemask strategies during the first wave of the pandemic in New York City.

The paper is structured as follows. In section 2.1, we present the compartmental dynamical model and we justify the selection of the values of the model parameters based on publicly available epidemiological data. Based on the model, in section 2.2, we derive estimations of the basic and effective reproduction numbers. In section 2.3 we perform an analysis of the Google mobility reports for Lombardy, thus providing a model for the effect of the social distancing measures in all activities to the residential activity which is also an important factor in the transmission of infectious diseases. In section 2.4, we present the proposed PMPC scheme constrained to the impact of the social mixing and activity patterns in the potential transmission. In section 3, we report simulation results, and in section 4 we conclude with a discussion of the findings. Finally, in the Supplementary Information section, we present the sources of the data used in the present study, competing formulations to the designed PMPC one and an analysis of the effect of PMPC parameters on the resulting policies.

2.   Materials and methods

2.1.   The modelling approach

As in the PMPC scheme we considered as constraint the maximum number of persons that are hospitalized due to COVID19, we extended the model proposed in Russo et al. ^[10] to include hospitalized cases and confirmed infected cases with generally moderate symptoms that are isolated at home. Thus, the model includes two compartments for the infected cases: (1) a compartment with the cases that are symptomatic and ask for medical help, hospitalized or die; confirmed cases may also be ^[10], (2) a compartment with the asymptomatic cases ^[10] which in general are not reported and recover relatively soon enough.

Considering a mixed population with $N$ people, the state of the population is described by the following compartments: $S(t)$ represents the average number of susceptible, $E(t)$ is the average number of exposed, $I_A(t)$ is the average number of asymptomatic cases; $R(t)$ is the average number of recovered asymptomatic persons. $I_S(t)$ denotes the average number of infected cases exhibiting more severe symptoms; a part of them is hospitalized, which in our model is denoted by $H(t)$, a part of them dies (denoted by $D(t)$) and a part of them recovers (denoted by $R_H(t)$). The other part of the confirmed infected cases from the $I_S(t)$ compartment in our model is represented by the variable $Q(t)$; it is assumed to experience mild symptoms and is quarantined at home. For this compartment there are no other available data except from its specific number. The confirmed recovered (in our model denoted by $R_H(t)$) are the recovered cases dismissed from hospitals. Furthermore, there are no reported data for the number of deaths from the infected people that are isolated at home. As said this compartment represents the part of cases that experience mild symptoms and thus they aren't expected to die from the disease.

Note that a wide testing policy may also result in the identification of asymptomatic cases belonging to the compartment $I_A$ that would then be assigned to compartment $I_S$. However, as a generally reported rule in Italy, tests were conducted only for those who asked for medical care with symptoms like fever and coughing. Thus, people who did not seek medical attention were tested on rare occasions ^[24]. Hence, for all practical purposes the compartment $I_S$ reflects the cases with mild to more severe symptoms.

We should also note that in our model, first, the number of mild to severe symptomatic cases as represented by the compartment $I_S$ does not coincide with the reported/confirmed number of infected cases since, (i) not all symptomatic cases are reported as confirmed cases, (ii) the confirmed reported total positive cases is the sum of the hospitalized and guaranteed cases. Second, the number of confirmed recovered cases does not coincide with the total number of recovered cases dismissed from the hospital $R_H$ since the reported cases include also a part of the confirmed cases that were isolated at home and after a second test were identified as recovered; thus, we expect $R_H$ to be less that the reported one. Third, there are no data for the number of recovered from the compartment $Q$ and since these cases experience mild symptoms, we assume there are no deaths from this compartment. Hence, in our model only the sum of the $H$ and $Q$ compartments reflect the number of actual confirmed cases.

For our analysis, due to the short period we analyze, we assume that the total number of the population remains constant. Based on demographic data, the total population of Lombardy is ten million, i.e. $N = 10^6$. The rate at which a susceptible ($S$) becomes exposed ($E$) to the virus is proportional to the density of asymptomatic infectious persons $I_A(t)$ in the total population and the symptomatic ones $I_S(t)$ the period preceding their hospitalization or home quarantine. The proportionality constant is the "effective" disease transmission rate, $\beta: = \bar{c}\, p$, where $\bar{c}$ is the average number of contacts per day and $p$ is the probability of infection upon a contact between a susceptible and an infected.

The assumption of the model is that only a fraction, $\epsilon$, of the actual number of exposed cases $E(t)$ will experience mild to more severe symptoms and ask for medical treatment; this fraction is represented by the compartment $I_S(t)$. Thus, our discrete mean field compartmental SEASQHRD model reads:

The discrete time points $t = 1, 2, \ldots$ denote days after day-zero, with the corresponding initial condition at the very start of the pandemic ($t_0$ = 0 = DAY-ZERO) being $S(0) = N-1$, $I_{A}(0) = 1$, $E(0) = I_{C}(0) = Q(0) = H(0) = R_{A}(0) = R_{H}(0) = R_{Q}(0) = D(0) = 0$.

The values of the epidemiological parameters were fixed in the proposed model at values reported by clinical studies as follows:

● $\beta (d^{-1})$ is the "effective" transmission rate of the disease; recently published studies have reported that this rate is for all practical purposes the same for both asymptomatic and symptomatic cases (see e.g. ^[25]). The transmission rate $\beta$ cannot be obtained by clinical studies, but only by mathematical models. Here we set $\beta = 0.68$ based on a previous work for the region of Lombardy^[10].

● $\sigma (d^{-1})$ is the average per-day "effective" rate at which an exposed person becomes infectious. In many studies, it is set equal to the inverse of the mean incubation period (i.e. the time from exposure to the development of symptoms). A study in China ^[26] suggests that it may range from 2 to 14 days, with a median of 5.2 days. Another study using data from 1,099 patients with laboratory-confirmed 2019-nCoV ARD from 552 hospitals in 31 provinces/provincial municipalities of China suggested that the median incubation period is 4 days (interquartile range, 2 to 7). However, the incubation period does not generally coincide with the time from exposure to the time that someone starts to be infectious. Regarding COVID-19, it has been suggested that an exposed person can be infectious well before the development of symptoms ^[27]. Indeed, it has been suggested that for a mean incubation period of 5.2 days infectiousness starts from 2.3 days (95% CI, 0.8–3.0 days) before symptom onset ^[27]. In our model, as explained above, ${}^{1}\!\!\diagup\!\!{}_{\sigma }$ represents the period from exposure to the onset of the contagious period and accordingly has been set as $\sigma = {}^{1}\!\!\diagup\!\!{}_{3 }$.

● $h (d^{-1})$ is the average per-day rate at which an infected person is hospitalized. The inverse of the rate of hospitalization $h$, i.e. ${}^{1}\!\!\diagup\!\!{}_{h }$ is the mean time from the onset of symptoms to hospitalization for those infected cases that need hospitalization. For Lombardy, the median value of this parameter has been reported to be 7 days ^[28]. Thus, considering also the fact that the infectiousness period starts 2-3 days before the onset of symptoms during the incubation period as explained above, we set ${}^{1}\!\!\diagup\!\!{}_{h } = {4+3}$ and thus $h = {}^{1}\!\!\diagup\!\!{}_{7 }$. This coincides with the corresponding number reported by the Istituto Superiore di Sanitá, Italy ^[28].

● $\delta_A (d^{-1})$ is the average per-day "effective" recovery rate within the group of asymptomatic cases in the total population. Here, it was set equal to $\approxeq{}^{1}\!\!\diagup\!\!{}_{6.6 }$ based on studies reporting that infectiousness declines quickly within $\sim$7 days ^[27]; this coincides with the mean value of the serial interval reported for Lombardy (6.6 days) ^[29] and that for China (mean serial interval of 6.3 days (95% CI 5·2–7·6))^[30].

● $\delta_H (d^{-1})$ is the average per-day "effective" recovery rate within the subset of hospitalized cases ($H$) that finally recover. As reported recently^[20], the overall median length of stay (LoS) in hospitals, which for all practical purposes coincides with the recovery period, decreased steadily from 21.4 (20.5–22.8) days in February to 5.2 (4.7–5.8) days in June. Based on the above, we set the recovery period for the hospitalized cases as $\delta_H = {}^{1}\!\!\diagup\!\!{}_{24 }$ for the period until March 8, $\delta_H = {}^{1}\!\!\diagup\!\!{}_{18 }$ for the period March 9-March 20, $\delta_H = {}^{1}\!\!\diagup\!\!{}_{14 }$ for the period March 21-April 10, and $\delta_H = {}^{1}\!\!\diagup\!\!{}_{10 }$ for the period April 11–May 4.

● $\delta_Q (d^{-1})$ is the average per-day "effective" recovery rate within the subset of quarantined cases ($Q$). In a study that is based on 55,924 laboratory-confirmed cases, the WHO-China Joint Mission has reported a median time of 2 weeks from onset to clinical recovery for mild cases ^[31]. Hence, based on these reports, we set $\delta_Q = {}^{1}\!\!\diagup\!\!{}_{14 }$ for the quarantined cases.

● $\gamma (d^{-1})$ is the average per-day "effective" mortality rate within the subset of $H(t)$ hospitalized cases that finally die. The median time from hospitalization to death for Lombardy during the early phase of the first wave of the pandemic (until April) has been reported to be seven days through the first wave, with an inter-quartile range (3-13) ^[32] for all Italy including Lombardy and ^[20] for Lombardy. We also note that the time from hospitalization to death subsequently increased as medical personnel gained experience. Thus based on the above the time from hospitalization to death was set as $\gamma = {}^{1}\!\!\diagup\!\!{}_{5 }$ for the period until March 20, $\gamma = {}^{1}\!\!\diagup\!\!{}_{7 }$ for the period March 21-April 10, and $\delta_H = {}^{1}\!\!\diagup\!\!{}_{10 }$ for the period April 11-May 4.

● $f_1$ is the ratio of all medium to severe infected cases (represented by the $I_S$ compartment) that gets hospitalized. For the early phase of the pandemic (until March 8th), a 50% over a total 5,830 confirmed cases reported in Lombardy were hospitalized ^[29], this percentage is also confirmed by other studies (see e.g. ^[33]). Regarding the next period (March 20$^{\rm th}$ to May 4) that was characterized by a saturation of the hospitals beds, it is reasonable to expect that this ratio would fall significantly in addition due to the fear of a high-risk nosocomial infection. Indeed, for the period March 20 and on this ratio was about 20% as reported also by other studies (see e.g. ^[33]). Based on these studies and through the analysis of the provided data, this ratio was set to $f_1 = 0.65$ for the period until March 8, $f_1 = 0.60$ for the period March 9 to March 20, $f_1 = 0.5$ for the period March 21 to April 10, and $f_1 = 0.2$ for the period April 11 to May 4.

● $f_2$ is the case fatality ratio of all hospitalized cases. Based on a study considering 38,715 hospitalized COVID-19 patients in Lombardy between 21 February and 21 April 2020, this ratio was relatively high ^[34]) mainly due to the saturation of ICU beds. In particular, as reported in a recent published study ^[20], for the region of Lombardy iy is estimated to have steadily decreased from 34.6% (32.5–36.6%) in February to 29.8% (29.3–30.3%) in March, 22.0% (21.1–23.0%) in April, 14.7% (13.3–15.9%) in May and 7.6% (6.3–10.6%) in June. Based on these studies, the fatality ratio was set $f_2 = 0.27$ for the period February 24-March 8, $f_2 = 0.23$ for the period March 9-April 10, and $f_2 = 0.20$ for the period April 11 -May 4.

● $\epsilon (d^{-1})$ is the fraction of the exposed cases that enter compartment $I_S$, i.e. those that will experience mild to severe symptoms and ask for medical care. For the case of Italy and in particular that of Lombardy during the first wave of the pandemic, the ratio of the confirmed to unreported cases was estimated to be around 0.07 (interquartile range: 0.023 to 0.1), ^[10]. In our model, we set $\epsilon = 0.12$ for the period until March 8, $\epsilon = 0.10$ for the period March 9 to March 20, $\epsilon = 0.09$ for the period March 21 to April 10, and $\epsilon = 0.05$ for the period April 11 to May 4.

Finally, DAY-ZERO in Lombardy, i.e. the initial date of the introduction of the disease in the region, was set to January 15, as estimated based on a methodological framework based on genetic algorithms together with the level of asymptomatic cases in the total population ^[10]. Another study ^[35] based on genomic and phylogenetic data analysis, reports the same time period, from the second half of January 2020 to early February 2020, as the time when the novel coronavirus SARS-CoV-2 entered northern Italy.

2.2.   Estimation of the basic $R_0$ and effective $R_e$ reproduction numbers from the model

There are three infected compartments, namely $E, I_A, I_S$ that determine the fate of the pandemic. Thus, considering the corresponding equations (2.2), (2.3), (2.4), and that at the very first days of the epidemic $S\approx N$, the Jacobian of the system as evaluated at the disease-free state is:

The eigenvalues of the Jacobian matrix indicate the stability of the equilibrium dynamics. In particular, if all the eigenvalues are within the unit circle, the disease-free state is stable; otherwise the infectious disease will in an outbreak. Here, we used Jury's stability criterion ^[36]. Thus, the characteristic polynomial is given by:

where

The necessary conditions for stability read:

Thus, sufficient conditions for stability are given by the following two inequalities:

where,

It can be shown, that the second necessary condition of (2.15) and the first sufficient condition of (2.16) are always satisfied for the range of values of the epidemiological parameters considered here. The first inequality (2.14) results in the necessary condition:

It can also be shown that, for the range of the parameters considered here, the second sufficient condition of (2.17) is satisfied if the necessary condition of (2.19) is satisfied. Thus, the necessary condition of (2.19) is also a sufficient condition for stability. Hence, the disease-free state is stable, if and only if, condition of (2.19) is satisfied.

Note that in this necessary and sufficient condition of (2.19), the first term reflects the contribution of compartment $I_A$ and the second term the contribution of compartment $I_S$ in the spread of the disease. Thus, the above expression reflects the basic reproduction number $R_0$ which is qualitatively defined by $R_0 = {}^{\beta}\!\!\diagup\!\!{}_{infection \; time }$. Hence, our model results in the following expression for the basic reproduction number:

For $\epsilon$ = 0 there are no asymptomatic cases in the population; thus the $R_0$ coincides to the one obtained from a simple SIR model.

Based on the values that are used in our model, the above expression results to $R_0 = 4.38$. This estimate is close to the one reported in Russo et al. ^[10] ($R_0 = 4.56$) and it is greater than the ones reported in other studies. For example, Allieta et al. ^[37] reported $R_0 = 3.88$, while in another study regarding Italy, D'Arienzo and Coniglio ^[38] used a SIR model to fit the reported data in nine Italian cities and found that $R_0$ ranged from 2.43 to 3.10. Finally, ^[39] reported an $R_0$ around 3.4. However, we note that the above studies that provide estimates of $R_o$ are based on the reported cases which may lead to an underestimation of the basic reproduction number (see also the discussion in ^[40]), while our model takes into account the compartment of asymptomatic cases which transmit the disease.

Regarding the estimation of the effective reproduction number, $R_e$ corresponds to the average number of secondary infections from a single infectious individual during an epidemic. For its calculation, we derive the next generation matrix $G$ with elements $g_{ij}$ formed by the average number of secondary infections of type $i$ from an infected individual of type $j$, given by (see also ^[10]):

where, ${\boldsymbol{F}}$ is the vector containing the transmission rates, and ${\boldsymbol{V}}$ is the vector containing the transition rates between the infected compartments. The effective reproduction number $R_e$ is the spectral radius, i.e. the dominant eigenvalue of $G$. In our model this is:

2.3.   Google mobility report analysis

In the spring of 2020, many nations enacted policies restricting resident movement in hopes of slowing the spread of COVID-19. Although the individual compliance to these policies may vary, there is a strong link between government mobility policies and the subsequent mobility behavior in a population. For example Armstrong et al. ^[41] assessed the effect of the restriction policies in the motility of 75 canadian and american cities, and Buonomo and Della Marca investigated the effect of the lockdown during the first wave of the COVID-19 pandemic in Italy^[42].. Furthermore, these mobility patterns are correlated with changes in COVID-19 cases ^[43].

There are several dates that demarcate major changes in Lombardy mobility policies. On March 4, 2020 all schools and universities were closed. On March 8, 2020, the first major lockdown was introduced, with the encouragement of working from home and closure of recreational venues. On March 21, 2020, measures were furthered; a stay-at-home order was issued along with the banning of pedestrian activity. These measures continued until May 4, when policies were gradually relaxed; non-essential industrial activities reopened, and other pedestrian activities was allowed. These key dates can be used to divide the initial wave of COVID-19 into four approximately time periods. The first period -February 15 to March 8- started with the availability of mobility data in Lombardy and ended with the introduction of lockdown. The second period stretches from March 9 to March 20: the period of initial lockdown. The third period represents the first phase of stricter measures, ranging from March 21 to April 19, while the fourth and final period begins April 20 and continues until policy relaxation on May 4.

To describe mobility patterns in Lombardy, we used data provided by the Google Community Mobility Reports ^[3]. The reports show movement trends within a population for a given region. The mobility data is organized into six categories (i.e. activities): retail & recreation (RR), grocery/pharmacy (G), parks (P), transit (T), workplaces (W), and residential (R). We note that within this period schools & universities (SU) were closed. The mobility values are expressed as percent deviation from the baseline value; the baseline was set as the median value from January 3 to February 6, 2020. All activities except for the residential one determine the percent deviation from the change in total visitors. The residential category is determined by the change in time spent in places of residence. In order to deduce a policy that is relevant, the day to day adherence of the population to each activity should be quantified. Table 1 first presents the average activity reduction and the standard deviation for the periods with a specific social policy. Based on these calculations, the last line in the table presents the expected maximal mobility reduction and the daily adherence of the Lombardy community during the first wave. The trajectory in time of the expected mobility reduction was computed using a centered nine-point moving average of the raw Goggle Mobility Report data for each activity. Subsequently, adherence was computed as the standard deviation of the distribution of the distances of the observed mobility from their expected values. The maximal reduction was identified as the maximum of the trajectory for the expected reduction, not the maximum observed reduction in the raw Goggle Mobility Report's data.

The relative impact of social mixing on disease transmission was analyzed to inform parameter $\beta$ in the model, since the amount of personal contacts made in a given location can be used to estimate the likelihood of disease transmission in that location. Here, we used the information reported by Mossong et. al. ^[19] who analysed data including age, sex, location, duration, frequency, and occurrence of physical contacts from 7,290 participants and 97,904 contacts (thus implying 13.4 contacts per day per person) from eight European countries (Belgium, Germany, Finland, Great Britain, Italy, Luxembourg, The Netherlands, and Poland). On average, 23% of the contacts were made at home, 21% at work, 14% at school, 3% while traveling, and 16% during leisure activities. For Italy, home, work, and leisure activities each account for about 20% of all physical personal contacts reported, transport accounts for the least amount of contacts at about 5% and a 20% is accounted to multiple activities. Despite travel being responsible for the least amount of physical and verbal contacts, diminished travel has still been strongly correlated with a decrease in the reproduction number ^[44]. Therefore, based on this information, each reported mobility category likely plays a significant role in the transmission of COVID-19.

To account for the role of each social activity in the transmission of COVID-19, an aggregate value was found by appropriately weighting the six categories for a given day. Furthermore, we took into account the caution that people exhibit in their daily interactions such as physical distancing and other preventive measures (e.g. hand hygiene, mask wearing). Thus, the transmission rate is impacted by the change in mobility such that $\beta = (1-u)(1-\theta_C)\beta_0$, with $u$ a measure of the reduction in the mobility and social activities. In this formulation, when $u = 0$, there are no mobility restrictions, and when $u = 1$, there is assumed to be complete lack of mobility and social interaction in the population beyond the residential. The caution that people exhibit in their daily interactions is captured by the term $\theta_C$, with $\theta_C = 0$ being no caution. Following the first curtails in mobility on March 8$^{\rm th}$, the caution was quantified to the value $\theta_C = 0.2$. This is in agreement with what has been reported in other studies about the effect of the cautiousness in the reduction of the disease transmission rate ^[45]. Analyzing the available data, the effect of policy changes for different activities in society can thus be quantified as

Obviously, the mobility activity in the category "Residential" is affected by the other activities, including transit which in turn also includes traveling to other regions for vacation or work. In order to model and assess this dependence, a linear regression model of the following form was fit:

The parameters of the model were estimated using the available mobility data from Feb 15, to May 4 2020. In Table 2, we report the values of the parameters and their 95% confidence intervals together with the fit statistics, namely the standard error (SE) of the coefficients, the t-statistic to test the null hypothesis that the corresponding coefficient is zero against the alternative that it is different from zero, given the other predictors in the model, the corresponding p-value for the t-statistic of the hypothesis test that the corresponding coefficient is equal to zero or not. The Root Mean Squared Error was 0.0194, the R-squared 0.979, and the adjusted R-squared was 0.978. Thus, at a statistical significance threshold of $\alpha = 0.05$, the variables $RR$, $G$, $P$, $W$ and $S$ are significant for describing the residential mobility, while the transit activity is not statistically significant. Note also that based on the statistical analysis, the hypothesis that $c_0$ is different than zero cannot be rejected at the significance level of 0.05, highlighting that the residential activity is purely dependent on the rest of the activities and can be discounted in the following section. Combining (2.23) & (2.24), the resulting cumulative effect of activity curtails on the social distancing variable is presented in Fig. S1. Employing the statistical analysis results summarized at Table 1, the resulting standard deviation (i.e. adherence) for social distancing, $u(t)$, was identified to be $\sigma_u = 0.0282$ and furthermore it could be considered constant during the period Feb. 24$^{\rm th}$ to May 4$^{ \rm th}$. Note that during this period there were specific governmental policies assigning mobility curtails to each activity. From this analysis the adherence by Lombardy to the policies can be estimated, which must be accounted for.

Employing the SEASQHRD model with the identified parameters and considering the variance in social distancing due to adherence (Table 1), we simulated the progression of the disease under the social distancing policy set by the government for 200 adherence scenarios. Figure 2 presents the simulations and the reported data for Lombardy. As the results show, the predictions of the model fits fairly well the reported data. We note that the lower limit of $R_e$ on May 4 is close enough to 1. As discussed above, the difference between our estimates with the reported $R_e$ in other studies (see e.g. ^[39] who reported a value of $R_e$ around 0.6 just before the end of the lockdown on May 4) is due to the fact our model also takes into account that the asymptomatic cases also transmit the disease.

2.4.   Probabilistic model predictive Control

The developed SEASQHRD model of §2.1 was employed to propose government policy as the pandemic progresses. A systematic mechanism for the proposed policy hinges on the solution of an optimization problem that balances the impact on quality of life due to prolonged social distancing on the one hand and the impact of the pandemic on life expectancy on the other. Considerations include the question of societal adherence to government advisories as well as other unaccounted factors which affect the outcome of predictions. We thus assume that government policies need to be continually evaluated and adapted to reality. A system's approach towards this end involves the design of a probabilistic model predictive control (PMPC) structure. The design formulates the policy question as an optimization problem with a finite number of policy decisions being made and their effect being predicted by the SEASQHRD simulator for a finite time horizon forward, called the prediction horizon, $T_P$, given the current state of the epidemic. Each policy decision is effective for a finite period, called the Action Horizon, $T_A$. The number of decisions, $N_d$, multiplied by $T_A$ is called the control horizon, $T_C: = N_d\, T_A$. Once a policy has been identified, the first policy decision is enacted and by the end of period $T_A$ the state of the epidemic is measured. A new optimization problem for the period $T_A$ to $T_A+T_P$ is then formulated and solved. This recursion constitutes the MPC, where one of the advantages is that it can adapt to uncertainty in the epidemic's behavior. Note, that the state of the compartments is a prerequisite for the SEASQHRD simulator, however only three compartments are directly measured, leaving seven compartment states unknown. In the previous section the analysis of the activities resulted in a distribution of values which impacted the epidemic's evolution. Recognizing the uncertainty in the announced policy due to adherence and the limitations of the derived model, the underlying optimization problem is defined as a probabilistic one

where $\alpha(t)$ is the governmental distancing policy, consisting of curtails on activities $\alpha(t) = [RR(t)\; G(t)\; P(t)\; T(t)\; W(t)\; SU(t)]^T$; each activity in $\alpha$ is normalized to map $(0, T_P] \mapsto [0, 1]$, where $0$ implies no restrictions on the corresponding activity and $1$ implies complete social curtail on it. Function ${\bf H}(\cdot)$ denotes the standard Heaviside function. Equation (2.26) denotes that during action period $j$ policies in $\alpha_j$ are active. Following standard MPC practice the final policies $\alpha_{N_d}$ are maintained till the end of the prediction horizon; one reason for this practice is that a final time penalty term can then be neglected. The state of infection, $X$, in (2.33) is defined as $X = [S\, E\, I_A\, I_S\, H\, Q\, R_A\, R_H\, R_Q\, D]$. Term ${\bf 0}$ denotes the zero vector of appropriate size. In the specific PMPC formulation a number of epidemic progression scenarios (defined as $N_{cases}$) is investigated; the adherence of Lombardy's population to the government curtails is the uncertain parameter which is sampled based on the analysis results of Table 1. The trajectory statistics are then computed to obtain the distributions used in the formulation.

The cost function of (2.25) balances the social cost of the epidemic to the social distancing, via the relative weights $\omega_I$ and $\Omega_A$. Function ${\rm\bf E}[\cdot]$ denotes the expectation of the distribution for argument $\cdot$, while ${\rm\bf P}[\cdot]$ denotes the probability argument $\cdot$ is true. The cost function is constructed to be convex to facilitate convergence to an optimum. Constraint (2.27) reflects the need for mobility of essential workers and for essential activities to take place in order for society to keep functioning. These upper bounds where computed from the Google mobility report data (presented in Table 1), by estimating the maximal expected curtail on mobility during the period Feb-24 to May 4, during which society was in essence completely closed.

The stress on society due to mobility restrictions can be significant. Constraint (2.28) is introduced to prevent a rapid increase in the restriction of each activity. In the original formulation, we impose a constraint of 25% increase in each activity restriction with the exception of schools & universities ($SU$) which is left unconstrained (i.e., $\delta\alpha_c = [0.25\; 0.25\; 0.25\; 0.25\; 0.25\; 1.0]$). Note that there is no constraint imposed for lifting the restrictions. The societal effect of the desired policy is also captured by the use of two constraints, (2.31) limiting the overflow of patients to hospitals and (2.32) limiting the effective reproduction number, defined in (2.22) (with the basic reproduction number, $R_0$ calculated from (2.20)) so that new infection cases dynamic becomes a contraction. In Figure 2(f), we observe that during the period February 24-May 4, $R_e$ remains above one (even though it is greatly reduced) which hints at the increasing rate of asymptomatic cases that transmit the disease during the same period. The use of probabilities in constraints (2.31) & (2.32) captures the degree of risk aversion in the desired policy, by imposing a limit to the risk of violating them via tunable parameters $P^c_{bed}\in[0, 1]$ and $P^c_{rep}\in[0, 1]$. A value of $0.5$ converts this constraint to an expectation constraint. For the specific formulation, as the values of $P^c_{bed}$ and $P^c_{rep}$ increase the PMPC formulation converges to the robust MPC formulation ^[46].

3.   Simulation results

The model predictive control simulation begins on February 24$^{\rm th}$, three days after the first case was reported in Lombardy. Within this 70 day period until May 4, when the restrictions were relaxed, the horizon lengths determine the duration and frequency of policy investigation and implementation. The prediction horizon, or how far into the future the simulation investigates at a given time point, was examined at lengths of 1 to 4 weeks. As the prediction horizon increases from 1 to 4 weeks, the action during the $1^{\rm st}$ period converges to a set value. A prediction horizon of 3 weeks was chosen to balance this convergence and the uncertainty in model predictions. The control horizon, which determines the number of policies during a given prediction, was chosen as 2 weeks to balance problem complexity and solution convergence. Figure S2 of the Supporting Information showcases the effect of the different horizon choices. The action horizon, that is the actual duration of a given policy, was chosen to be 1 week. This choice of horizon allows for frequent policy changes if necessitated by the model predictions, while not being so frequent as to be unrealistic from a policy implementation standpoint.

The choice of the objective function determines the optimization calculations performed by the controller. A linear objective function of (11.19) could be used, but the quadratic objective function of (2.25) resulted in stable control pattern and was used in this work. The probabilistic form of this optimization function was used to account for inherent uncertainty in policy implementation; it is less aggressive than a robust formulation of (11.10)-(11.18), which assumes the worst-case scenario during optimization. Consequently, the choice of a PMPC formulation allows for less heavy-handed policies while still considering the uncertainty associated with population adherence.

Another consideration is the relative cost assigned to mobility restrictions, scaled by the relation between parameters $\Omega_A$ and $\omega_I$ & $\omega_R$. Focusing first on the relative cost of activity restrictions as reflected by the Google mobility report data during the period Feb-24 to May-4, was used to elucidate the preference given in the restrictions. Specifically, the need for groceries is central for a society where activities such as restaurants are curtailed. Similarly, transit is needed in order for basic social activities to take place. On the other hand, parks are not explicitly needed for society to function during epidemic conditions. Using similar arguments a convenient diagonal matrix $\Omega_A$ was defined, with the diagonal elements expressing the relative resistance to curtails being $\omega_A = [0.2\; 1.0\; 0.5\; 0.5\; 1.0\; 0.5]$. For a prediction horizon of 3 weeks, control horizon of 2 weeks, and action horizon of 1 week, a selection of $\|\Omega_A\|_F\rightarrow 0$ results in government policies of aggressive lockdown. A selection of $\|\Omega_A\|_F\approx \omega_R$ results in moderate mobility restrictions. To equally weight the impact on hospitalization and impact on disease progression we employed $\omega_I \approx 10 \omega_R$.

Figure 3 shows the simulation results of disease progression for the period Feb-24 to May-4 under the optimal policy as identified by the MPC formulation of (2.25)-(2.33) with $\omega_R = 10$, $\omega_I = 10$, $\Omega_A$ as discussed previously, $T_p = 21$, $T_c = 14$, $T_a = 7$ days. The uncertainty in the simulation is captured via presenting 200 population adherence scenarios(i.e., $N_{cases} = 200$), informed from Table 1. These parameters ultimately prescribe a government policy that corresponds to approximately a 60% reduction in mobility. The difference in outcome between the implementation of the prescribed restrictions with the historical data from Lombardy for reference by comparing Figure 3 with Figure 2. The simulation of optimal restrictions leads to less infections, hospitalizations, and deaths than the historical data. A major contributor to the difference in deaths is the date of introduction of major curtails which in our simulation starts on Feb-24 which also has the effect that population caution becomes $\theta_C = 20\%$ on that day. The historical data in contrast show the effect of casual curtails starting on Feb-21 with major (and aggressive) curtails starting on March 8; population caution becomes $\theta_C = 20\%$ on March 8.

4.   Discussion

One of the most challenging aspects of addressing the COVID-19 crisis is the scale of response, promptness, and fast adaptation to a number of factors. In a fluid situation in which the inherent uncertainly that characterizes the reported data and generally our knowledge at least at the first stages of an pandemic, mathematical modeling and control theory may play an important role for the investigation of "what if" scenarios for the implementation of efficient social distancing and vaccination policies.

Here, using a probabilistic model predictive control approach we attempted to retrospectively assess the impact of social distancing measures for various activities on the evolution of the pandemic in Lombardy during the first wave of the pandemic (Feb 21, 2020 to May 4, 2020) by quantifying the consequence on the fatalities and the pressure on hospital infrastructure. Thus, an important issue for the efficient modelling of the COVID-19 pandemic concerns the estimation of the unobserved part of infected cases that can transmit the virus ^[47]. Towards this aim, a compartmental modelling approach was presented that quantifies the unobserved number of asymptomatic cases ^[10] in the population and also the effect of social contact during different activities to the spread of infectious diseases. In particular, the proposed model extends the one proposed by Russo et al. ^[10] for the region of Lombardy, and now considers hospitalized and guaranteed populations. Furthermore, the effect of the social contact during various activities is modelled based on the Google mobility reports for the region of Lombardy as well as by metanalysis of the social contacts and their effect on the spread of infectious diseases as reported for Italy by Mossong et al. ^[19]. To increase the resolution of the PMPC formulation we also quantify the effect of the reported activities, namely, workplaces, grocery and pharmacy, retail and entertainment, parks and transit to the residential change of mobility. This was achieved by fitting a linear regression model based on the corresponding Google mobility reports.

We show that the proposed modeling framework approximates quite well the actual evolution of the pandemic in terms of reported infected cases, deaths and hospitalized cases. By activating the PMPC scheme from February 21, 2020, i.e. the day that the first official case was reported in Lombardy, we show that the obtained number of deaths would be significantly less (around 4,000 at May 4 instead of the actual number of 14,000), however the spread of the disease would not have been eradicated by May 4. Indeed our simulations show that due to the presence of a significant portion of asymptomatic cases which transmit the disease but are rarely reported, due to the mild or no symptoms, the effective reproduction number is still above 1 on May 4. However, the pressure to the hospital infrastructure would be significantly smaller, with a maximum of 3500 hospitalized cases peaking around mid April, instead of the actual number of cases that exceeded the 13000 bed capacity of the hospital system in the same period. Of course, as this is a retrospective analysis study, it should be regarded as indicative of what would happen if we knew a-priori the scale and severity of the threat. Thus, our analysis should not be considered by any means as a "ground truth" or an accurate prediction of what would happen if such measures would have been adopted. However, the results highlight the importance of very prompt responses and alertness for analogous emergent and re-emergent infectious diseases that will certainly occur into the future. Importantly, our study points to the urgent need for strengthening public preventative mechanisms including the epidemiological surveillance infrastructure that will allow for faster responses, and the enhancement of the (fragmented in most countries) primary care facilities networks for the standardized and accurate exchange of meta-information.

The proposed framework can be extended in various directions. One issue at the forefront of current policy making is vaccination, which this model does not account for. To accurately forecast policy decisions today, the model should be revised to capture vaccine-induced immunity to COVID-19. Other considerations include the probability of reinfection due to the emergence of the virus variants, which was a less significant factor in this study due to the short time period examined. Another challenge posed by the COVID-19 pandemic was the accurate state estimation. In this work, we assumed that the reported values were in principle correct. For example, we assumed that the reported cases of hospitalized and quarantined accurately presented the severity of the symptoms, and there was no bias due to the strain on the hospital network as the number of cases approached maximum capacity. This assumption can be challenged during the period when 13000 hospitalized cases were reported. Future formulations can include moving horizon estimation (MHE) that has the task of reducing the effect of measurement bias and uncertainty to estimate the state of the epidemic.

Furthermore, the PMPC formulation could be extended to quantify the uncertainty in the evolving dynamics and adapt its predictions, taking into account the reported, from clinical studies, distributions of the epidemiological parameters rather than their expected values. Combined with an MHE structure, previous observations can also be used to evaluate expected disease parameters' values. Alternatively, a further detailed compartmental model that takes also the age distribution in the population (see e.g. ^[48] can be used to. A critical point that is connected with the above is that with such a small number of infectious cases at the initial stage of the simulations, a stochastic model or a hybrid stochastic model could be more realistic in which uncertainty could be modelled in the form of realistic perturbations (see for example ^[49]).

Linked to the above, as COVID is fast becoming our new "norm", government policy makers need to examine longer periods of time. Considerations like the probability of reinfection due to the emergence of the virus variants now become dominant, as well as the slow change of parameter values as seasons change. Formulations that include adaptive MHE structures would be central in an extended PMPC formulation to accurately capture the evolution of the pandemic and assess its uncertainty. Finally, one could extend the proposed framework spatially, by extending the model to a graph of regions and their region specific interventions (for early formulations see ^[18]. These are some important aspects of the pandemic that we aim to tackle in future work.

Acknowledgments

The work was supported by the Ministero dell'Università e della Ricerca through the program of Fondo Integrativo Speciale per la Ricerca (FISR). No of the Program: FISR2020IP-02893.

Conflict of interest

The authors declare there is no conflict of interest.

Supporting Information

S.1.   Epidemiological data

All the relevant data used in this paper are publicly available and accessible at https://lab.gedidigital.it/gedi-visual/2020/coronavirus-i-contagi-in-italia/. The reported cumulative numbers of cases from February 21 to March 19 are listed in Table S1. The data from February 21 to March 8 have been used for the calibration of the model parameters and the data from March 9 to March 19 have been used for the validation of the model.

S.2.   Google mobility data

All the mobility data used in the present study are publicly available at https://www.google.com/covid19/mobility/. The mobility data for the Lombardy region between February-15 and July-30 are graphically presented in Figure S1. These data were employed to estimate the standard deviation and bias between the government guidelines and population behavior, discussed in §2.3.

S.3.   Other MPC formulations

Beyond the PMPC formulation of (2.25)-(2.33), deterministic and robust MPC problems, ((S.3.1)-(S.3.9) and (S.3.10)-(S.3.18) respectively, were also formulated and their performance was evaluated.

The deterministic problem does not consider population adherence to government guidelines, assuming perfect response. This leads to a simpler formulation at the cost of a potentially optimistic view of the progression of the epidemic that can lead to the violation of constraints (S.3.9) and importantly of (S.3.8).

At the other end of the spectrum, the robust formulation considers bounded uncertainties in the state evolution and the enforcement of the control action. It attempts to identify a policy that will not violate the constraints for all investigated scenarios. It is important to note that basic robust MPCs are designed to address problems where the uncertainty is bounded; elaborate RMPC formulations have been designed using Lyapunov arguments that guarantees the stability ^[50] and robustness ^[46] of the closed–loop system. In the presented work, this is addressed by creating a set of "scenarios"; each scenario contains values of the adherence of the population to each activity curtail and the associated expected initial state $X_0$ that will be used by the SEASQHRD model. A robust optimization problem was also formulated where a non-convex set was defined containing value ranges for each activity and unmeasurable initial state. As the results from these two robust MPC formulations were almost identical, we only present the first RMPC formulation here.

Finally, recognizing that $H$ and $\alpha$ always attain positive values, a linear objective function can be considered, compared to the widely used quadratic one:

S.4.   Effect of MPC parameters on policy design

Numerous simulations were performed to investigate the effect MPC parameters have on proposed action. Figure S2 showcases the effect of the MPC horizon parameters on the calculated optimal action. We chose horizons $T_p = 21$, $T_c = 14$, $T_a = 7$ days that balanced problem complexity and accuracy well.

The effect of the choice of weights of $\Omega_A$ can be visualized for the case when $\omega_I = 0$ & $\omega_R = 1$, i.e. they are relatively small compared to $\|\Omega_A\|_F = \approx 1.6$. We observe that when $\Omega_A = 0.6{\rm \bf I}_{6\times 6}$ stringent curtails are imposed on important actives such as groceries and work (shown in Figure S4) compared to when the "variable" weight is employed (shown in Figure S3), with the aggregate curtail being $u(t)\approx 0.57$. Considering the societal effects of activities the nominal weight for the optimization used in formulation of (2.25)-(2.33) was chosen as