Control with uncertain data of socially structured compartmental epidemic models.

The adoption of containment measures to reduce the amplitude of the epidemic peak is a key aspect in tackling the rapid spread of an epidemic. Classical compartmental models must be modified and studied to correctly describe the effects of forced external actions to reduce the impact of the disease. The importance of social structure, such as the age dependence that proved essential in the recent COVID-19 pandemic, must be considered, and in addition, the available data are often incomplete and heterogeneous, so a high degree of uncertainty must be incorporated into the model from the beginning. In this work we address these aspects, through an optimal control formulation of a socially structured epidemic model in presence of uncertain data. After the introduction of the optimal control problem, we formulate an instantaneous approximation of the control that allows us to derive new feedback controlled compartmental models capable of describing the epidemic peak reduction. The need for long-term interventions shows that alternative actions based on the social structure of the system can be as effective as the more expensive global strategy. The timing and intensity of interventions, however, is particularly relevant in the case of uncertain parameters on the actual number of infected people. Simulations related to data from the first wave of the recent COVID-19 outbreak in Italy are presented and discussed.

Introduction

From the digital tracking systems of the Koreans to UK’s initial choice of not wanting to do anything to counter the spread of the virus, passing through the militarized quarantines of the Chinese and those less authoritarian and more involved of the Italians, the reaction of different countries to the COVID-19 outbreak has shown a series of very different approaches that can also be explained considering the different cultural and political attitudes of the countries concerned.

In all cases, after an initial phase, even those governments that were less restrictive in the face of the pandemic’s inexorable progress had to take strong containment measures. There’s a graph that has become the symbol of the COVID-19 pandemic most of all. It shows in a simple and intuitive way the importance of slowing down the spread of an epidemic as much as possible (“flattening the curve”), so that the healthcare system can take care of all the sick without collapsing. Its success has helped to save many lives, raising awareness of good practices to slow down an epidemic: stay as much as possible at home, reduce social interactions and wash your hands often and well.

These ”non-pharmaceutical” intervention measures, however, entail significant social and economic costs and thus policy makers may not be able to maintain them for more than a short period of time. Therefore, a modeling approach based on a limited time horizon that takes into account the social structure of the population is necessary in order to optimize containment strategies. Most current research, however, has focused on control procedures aimed at optimizing the use of vaccinations and medical treatments (Barro et al. 2018; Bolzoni et al. 2017; Lee et al. 2010; Colombo and Garavello 2019) and only recently the problem has been tackled from the perspective of non-pharmaceutical interventions (Lin et al. 2010; Morris et al. 2020). In addition, data collected by governments are often incomplete and heterogeneous, so a high degree of uncertainty needs to be incorporated into predictive models (Capaldi et al. 2012; Capistrán et al. 2012; Chowell 2017; Jagodnik et al. 2020; Mizumoto et al. 2020; Roberts 2013). This is the case of the spreading of COVID-19 worldwide, which have been often mistakenly underestimated due to a combination of factors, including deficiencies in surveillance and diagnostic capacity, and the large number of infectious but asymptomatic individuals (Jagodnik et al. 2020; Mizumoto et al. 2020; Zhang et al. 2020).

For almost a hundred years, mathematical models have been used to describe the spread of epidemics (Kermack and McKendrick 1927). The models currently used largely originate from the model proposed by Kermack and McKendrick at the beginning of last century. Even if the model contains strong simplification assumptions, the concepts introduced through this model are essential to provide a first intuition on the dynamics of epidemics, an intuition that remains confirmed in more complex models, albeit with numerous modifications (see for example Hethcote 2000; Capasso and Serio 1978). The model provides for the division of the population into compartments, the susceptible, healthy individuals who may be infected, the infectious, who have already contracted the disease and can transmit it, and the removed, compartment that includes those who are healed and immune.

The hypothesis made by Kermack and McKendrick is that of the homogeneous “mixing”; that is, it is assumed that each individual has the same probability of contacting any other individual in the population. One understands how this hypothesis is unrealistic: we are often in contact with people from our family, our workplace, school class, group of friends and very rarely with those who live in a different place, have different ages and professions. In recent years, therefore, computational models have been developed that try to take into account additional social characteristics of individuals in order to arrive at more accurate predictions by keeping, however, the simplicity of compartmental models (Castillo-Chavez et al. 1989; Glasser et al. 2012; Lee et al. 2012; Hethcote 1996; Iannelli et al. 1992; Franceschetti and Pugliese 2008).

In this paper starting from a general compartmental model with social structure, typically the age dependence, we consider the external action of a policy maker that aims at reducing the spread of the epidemics by applying non pharmaceutical intervention measures, such as social distancing and quarantine. The mathematical problem is formulated as an optimal control problem characterized by a functional cost whose objective is to minimize the number of infectious people in a given time horizon. Through an instantaneous control strategy we compute an explicit feedback control that allows us to derive new SIR-type compartmental models capable of describing the epidemic peak reduction. Previously, this type of approach has been used successfully in the case of social models of consensus (Albi and Pareschi 2018; Albi et al. 2014, 2015a, b; Bongini et al. 2015; Caponigro et al. 2013; Düring et al. 2018).

The feedback controlled models are subsequently extended to take into account the presence of uncertain infection parameters and data. In fact, to have reliable forecasts it is of paramount importance to consider the presence of uncertain quantities as a structural feature of the epidemic dynamics. This aspect is of paramount importance in the case of pandemic COVID-19, in which undetected infectious individuals play a key role in the spread of the disease. In this regard, it is worth noting that our methodology can be easily extended to the case of more complex compartmental models. The decision to limit ourselves to a simple SIR-type compartmentalization was related on the one hand to the increased complexity given by the dependence on the social structure, which proved to be crucial in the case of the COVID-19 pandemic, and on the other hand to the introduction of a systematic uncertainty in the number of infected to avoid a complex sub-compartmentalization of the infectious population and the consequent difficulties due to parameter identification and the inability to follow a data-driven approach (Giordano et al. 2020; Flaxman et al. 2020). From a mathematical point of view, we can rely on the methods of uncertainty quantification (UQ) to obtain efficient and accurate solutions based on stochastic orthogonal polynomials for the differential model with random inputs (Xiu 2010).

Few results are actually available regarding methods of UQ in epidemic systems, we mention in this direction (Capistrán et al. 2012; Chowell 2017; Capaldi et al. 2012; Roberts 2013). The main idea is to increase the dimensionality of the problem adding the possible sources of uncertainty from the very beginning of the modeling. Hence, we extrapolate statistics by looking at the so-called quantities of interest, i.e. statistical quantities that can be obtained from the solution and that give some global information with respect to the input parameter like expected solution of the problem or higher order moments. Several techniques can be adopted for the approximation of the quantities of interest, in this paper we adopt stochastic Galerkin methods that allow to reduce the problem to a set of deterministic equations for the numerical evaluation of the solution in presence of uncertainties. Compared to conventional Monte Carlo methods, based on stochastic sampling, these methods guarantee an exponential convergence in the case of smooth uncertainty distributions and allow a much more accurate and efficient estimation of random parameters. We refer the interested reader to recent surveys and monographs on the topic (Dimarco et al. 2017; Jin and Pareschi 2017; Pareschi 2021; Xiu 2010).

In particular, we consider the case in which the policy maker applies his control based on several possible estimators on the actual number of infected people. The need for long-term interventions shows that alternative actions based on the social structure of the system can be as effective as the more expensive optimal strategy. The importance of the timing and intensity of interventions is particularly relevant in the case of uncertain parameters on the actual number of infected people.

The rest of the manuscript is organized as follows. In Sect. 2 we introduce the structured social SIR model and formulate the mathematical approach for containment measures to reduce the spread of the disease. Next, a feedback controlled model used in the subsequent analysis is derived within a short time horizon approximation. In Sect. 3 we generalize the feedback controlled model to take into account the presence of uncertainties. Section 4 is dedicated to the presentation of some numerical examples including applications to the first wave of the COVID-19 epidemic in Italy. In separate Appendices we provide details on the generalizations of the present approach to more realistic epidemic models for COVID-19 including additional compartmentalizations, on the stochastic Galerkin method employed to efficiently address the uncertainties, and on the social interaction matrices characterizing the contact rates.

Control of epidemic dynamics

The starting model in our discussion is a SIR-type compartmental model with a social structure. The presence of a social structure is in fact essential in deriving appropriate sustainable control techniques from the population for a protracted period, as in the case of the recent COVID-19 epidemic. We will discuss in Sect. 3 how to modify the model through the introduction of a stochastic parameter that takes into account the dependence on uncertain data, and thus implicitly introduce the role of undetected infectious in the dynamics.

Compartmental models with social structure

The heterogeneity of the social structure, which impacts the diffusion of the infective disease, is characterized by the vector characterizing its social state and whose components summarize, for example, the age of the individual, its number of social connections or its economic status (Hethcote 1996, 2000). We denote by , and , the distributions at time of susceptible, infectious and recovered individuals, respectively in relation to specific social characteristics. We assume that the rapid spread of the disease and the low mortality rate allows to ignore changes in the social structure, such as the aging process, births and deaths.

Consequently, for a given population of total number N, we have thatwhere is the total distribution of the social features defined by the vector . Hence, we recover the total fraction of the population which belong to the susceptible, infected and recovered as followsIn a situation where changes in the social features act on a slower scale with respect to the spread of the disease, the socially structured compartmental model follows the dynamicswhere the function represents the uncertain interaction rate among individuals with different social features and the recovery rate which may depend on the social feature.

Often, in socially structured models the interaction rate between people is assumed to be separable, and proportionate to the activity level of the social feature (Hethcote 1996, 2000), as followswith the average number of people contacted by a person with social feature per unit time. Alternative approaches are based on preferential mixing (Glasser et al. 2012; Castillo-Chavez et al. 1989). Specific examples of age-dependent social interaction matrices are reported in Appendix C.

We introduce the usual normalization scalingand observe that the quantities S(t), I(t) and R(t) satisfy the conventional SIR dynamicswhere the fraction of recovered is obtained from . We refer to Hethcote (1996), Hethcote (2000) for analytical results concerning model (2) and (4). In the following we will adopt the simple compartmental model (2) to derive our feedback controlled formulation in presence of uncertainty. The extension to more realistic compartmental models in epidemiology, designed specifically for the COVID-19 pandemic, can be carried out in a similar way (Gatto et al. 2020).

In order to simplify the description, we will consider the one-dimensional case and set the social dependence as the age a of the individual because of its importance in epidemic dynamics. It is clear, however, that similar containment procedures can impact also on other social features, like the wealth of the individuals (Dimarco et al. 2020). We will first formulate the feedback controlled SIR model in the deterministic case and subsequently extend our approach to the presence of uncertain parameters.

Optimal control of structured compartmental model

In order to define the action of a policy maker introducing a control over the system based on social distancing and other containment measures linked to the social structure we consider an optimal control framework. The choice of an appropriate functional is problem dependent (Lee et al. 2010).

In our setting, we account the minimization of the total number of the infected population I(t) through the an age dependent control action depending both on time and pairwise interactions among individuals with different ages.

Thus, we introduce the optimal control problemsubject towith initial condition , and .

The number of infected individual is measured by a monotone increasing function such that . This function models the policy maker’s perception of the impact of the epidemic by the number of people currently infected and in the sequel will be referred to as perception function. For example , for implies an underestimation of the actual number of infected corresponding to . The control aims to minimize this measure of the total infected population by reducing the rate of interaction between individuals. We consider a quadratic cost for its actuation.

Such control is restricted to the space of admissible controlswhich ensure the admissibility of the solution for (6). The above restriction on admissible controls can be relaxed if we consider controls that violate the previous condition locally but preserve the inequality in integral form after integration against .

The solution to problem (5–6) is computed through the optimality conditions obtained from the Euler-Lagrangian as follows where are the associated Lagrangian multipliers. By computing the variations with respect to (s, i, r) we retrieve the adjoint systemwith terminal conditions and . Note that the contribution of vanishes since the control does not act directly on population R, and the removed population is not considered in the minimization of the functional. The optimality condition readsThe optimality conditions (7–8) are first order necessary conditions for the optimal control . In order to be admissible then the control readswhere .

The approach just described, however, is generally quite complicated when there are uncertainties as it involves solving simultaneously the forward problem (5–6) and the backward problem (7–8). Moreover, the assumption that the policy maker follows an optimal strategy over a long time horizon seems rather unrealistic in the case of a rapidly spreading disease such as the COVID-19 epidemic. Let us emphasize that extending the above optimal control formulation to more complex compartmental models designed specifically for COVID-19, like SEPIAR or SIDHARTE (Gatto et al. 2020; Giordano et al. 2020), can be done by generalizing the control functional (5) to include, for example, the hospitalized compartment, or other specific indicators that can be measured from the data. For all of these models, the feedback control strategy described in the next section does not change substantially. We refer the reader to Appendix A for more details.

Feedback controlled compartmental models

In this section we consider short time horizon strategies which permits to derive suitable feedback controlled models. These strategies are suboptimal with respect the original problem (5–6) but they have proved to be very successful in several social modeling problems (Albi and Pareschi 2018; Albi et al. 2015a, b, 2014; Düring et al. 2018). To this aim, we consider a short time horizon of length and formulate a time discretize optimal control problem through the functional restricted to the interval , as followssubject toBy recalling that the macroscopic information on the infected is we can derive the minimizer of computing . We retrieve the following nonlinear equationIn order to pass to the limit we must rescale the penalization term as so that we can introduce the above instantaneous strategy directly in the discrete system (10)-(11).

The resulting controlled dynamic, corresponding to the feedback controlled continuous system (6), reads as follows In what follows we provide a sufficient conditions for the admissibility of the instantaneous control in terms of the penalization term . Indeed we want to assure that the dynamics preserve the monotonicity of the number of susceptible population s(a, t).

Proposition 1

Let , and a monotonically non decreasing function, then for all and , solutions to (13) are admissible if the penalization satisfies the following inequalitywhere and are respectively the peak reached by the infected of age and by the total population of infected.

Proof

By imposing the non-negativity of the controlled reproduction rate inside the integral we haveThis inequality has to be satisfied for every time . Next we observe that the number of susceptible is decreasing in time, therefore for all t. Moreover reaches a peak before decreasing to 0 (note that this peak can also be in ), say for the macroscopic variable and . Thus, thanks to the monotonicity of , we can restrict the previous inequality as follows

In Fig. 1 we report the phase diagram of susceptible-infected trajectories for the controlled model with homogeneous mixingwith . The dynamics are similar to the classical SIR model but with a nonlinear contact rate. In particular, the trajectories are flattened when the value of the control is such that . Note, however that this status is not an equilibrium point of the system.

Fig. 1

Phase diagram of susceptible-infected trajectories for the controlled SIR-type model with homogenous mixing and . Different choices of the penalization term are reported. Left plot the case , right plot . The line markers point out the peaks of the infected population for each choice of

To understand this, let us observe that an equilibrium state for (15) satisfies the equationsAn equilibrium point corresponds to the classical state in which we have the extinction of the disease and arbitrary and defined by the asymptotic state of the dynamics (Hethcote 2000). Now, let’s suppose that , , we can look for solutions where control is able to perfectly balance the spread of the diseaseconsequentlywhich is satisfied only for when . The stability analysis of this unique equilibrium point can be performed by standard arguments and we omit the details.

Control of epidemic dynamics with uncertainties

Since the beginning of the outbreak of new infectious diseases, the actual number of infected and recovered people is typically underestimated, causing fatal delays in the implementation of public health policies facing the propagation of epidemic fronts. This is the case of the spreading of COVID-19 worldwide, often mistakenly underestimated due to deficiencies in surveillance and diagnostic capacity (Jagodnik et al. 2020; Remuzzi and Remuzzi 2020). Health systems are struggling to adopt systematic testing to monitor actual cases. Moreover, another important epidemiological factor affecting data reliability is the large proportion of asymptomatic (Jagodnik et al. 2020; Mizumoto et al. 2020; Zhang et al. 2020).

Among the common sources of uncertainties for dynamical systems modeling epidemic outbreaks we may consider the followingIn the following we consider the effects on the dynamics of uncertain data, such as the initial conditions on the number of infected people or the interaction and recovery rates. On the numerical level we consider techniques based on stochastic Galerkin methods, for which spectral convergence on random variables is obtained under appropriate regularity assumptions (Xiu 2010). For simplicity, the details of the numerical method that allows to reduce the uncertain dynamic system to a set of deterministic equations are reported in Appendix B.

noisy and incomplete available data;

structural uncertainty due to the possible inadequacy of the mathematical model used to describe the phenomena under consideration.

Socially structured models with uncertain inputs

We introduce the random vector whose components are assumed to be independent real valued random variableswith the Borel set. We assume to know the probability density characterizing the distribution of . Here, is a random vector taking into account various possible sources of uncertainty in the model.

In presence of uncertainties we generalize the initial modeling by introducing the quantities , and representing the distributions at time of susceptible, infectious and recovered individuals. The total size of the population is a deterministic conserved quantity in time, i.e.Hence, the quantitiesdenote the uncertain fractions of the population that are susceptible, infectious and recovered respectively.

The social structured model with uncertainties readsTo illustrate the impact of uncertainties let us consider the simple following example. In the case of homogeneous mixing with uncertain contact rate , , distributed as p(z) the model readswith deterministic initial values and . The solution for the proportion of infectious during the initial exponential phase is (Roberts 2013) and its expectation where the function represents the statistical correction factor to the standard deterministic exponential phase of the disease . If z is uniformly distributed in we can explicitly compute More in general, if z has zero mean then by Jensen’s inequality we have for , so that the expected exponential phase is amplified by the uncertainty (see Roberts 2013).

In a similar way, keeping constant, but introducing a source of uncertainty in the initial data , and distributed as p(z) the solution in the exponential phase reads and then its expectation where is the mean of the variable z. Therefore, the expected initial exponential growth behaves as the one with deterministic initial data . Of course, if both sources of uncertainty are present the two effects just described sum up in the dynamics.

Remark 1

The presence of a large number of undetected infected is at the basis of the construction of numerous epidemiological models with an increasingly complex compartmental structure in which the original compartment of the infected is subdivided into further compartments with different roles in the propagation of the disease (Giordano et al. 2020; Gatto et al. 2020; Flaxman et al. 2020). To clarify the relationships to other deterministic compartmental models, let us consider, for simplicity, model (17) in absence of a social structureand with a one-dimensional random input distributed as p(z). Furthermore, for a function F(z, t) we will denote its expected value as . Now, starting from a discrete probability density functionwe have , with . Taking the expectation in (21), we can writewith , . For example, in the case n=2, by identifying and with the compartments of detected and undetected infectious individuals, respectively, we can formulate the equivalent partitioningwhich has the same structure of a SIAR compartmental model including the undetected (or the asymptomatic) class. In the following, we will not rely on discrete probability distributions, but on continuous representations that can be associated with the overall probability of having a certain number of infectious individuals (detected or undetected). The additional dependence of the epidemiological parameters from the random variable allows us to take into account changes in the corresponding dynamics of disease transmission and recovery.

The feedback controlled model with random inputs

In presence of uncertainties the optimal control problem (5–6) is modified as follows being a suitable operator taking into account the presence of the uncertainties . Examples of such operator that are of interest in epidemic modeling are the expectation with respect to uncertaintiesor relying on data which underestimate the number of infectedwhere is a given value such that , and .

In (24) the space of admissible controls is defined as or in a more relaxed form if we consider the above inequality after integration against .

The minimization problem (24) is subject to the following dynamics with initial condition , and .

The implementation of instantaneous control for dynamics in presence of uncertainties follows from the derivation presented in Sect. 2.3. We can derive the minimizer of computing from the restriction of the minimization problem (24) in the interval or equivalentlywhere we assumed , to obtain the following nonlinear identityThe above assumption on is clearly satisfied by (25) and (26), where in the case of (26) we used the notationBy introducing the scaling we can pass to the limit for to getwhich defines the feedback controlled model in presence of uncertainties.

Examples from the COVID-19 outbreak in Italy

In this section we present several numerical tests on the constrained compartmental model with uncertain data. Details of the numerical method used are given in Appendix B. In an attempt to analyze sufficiently realistic scenarios, in the following examples we will refer to values taken from Italian data on the COVID-19 epidemic (del Consiglio dei Ministri 2020). More precisely, in the first test case we illustrate the behavior of the model in a simplified setting in absence of uncertainty and social structure and without trying to reproduce scenarios closely related to current data. In the second test case, following a progressively more realistic approach, we consider the impact of the presence of uncertain data in the controlled model with homogeneous social mixing and calibrated on Italian data. The same setting is then considered in Test 3 taking into account the additional effects given by the social structure of the system, characterized by suitable social interaction functions and an age-dependent recovery rate. The final scenario, explored in Test 4, examines the influence on the spread of infectious disease induced by relaxed confinement measures related to the social structure of the system.

Test 1: Containment in homogeneous social mixing dynamics

To illustrate the effects of controls introduced that mimic containment procedures, let us first consider the case where the social structure is not present. Furthermore, to simplify further the modeling, in this first example we neglect any dependence on uncertain data.

We consider as initial small number of infected and recovered , . These normalized fractions refer specifically to the first reported values in the case of the Italian outbreak of COVID-19, even if in this simple test case we will not try to match the data in a quantitative setting but simply to illustrate the behavior of the feedback controlled model. Based on recent studies (Jagodnik et al. 2020; Zhang et al. 2020; Liu et al. 2020), the initial infection rate of COVID-19 has been estimated between 2 and 6.5. Here, to exemplify the possible evolution of the pandemic we consider a value close to the lower bound, taking and , namely a recovery rate of 10 days, so that .

In Figs. 2 and 3 we report the infected and recovered dynamics based on the activation of the control in two different time frames. In Fig. 2 the activation for , which means that after 100 days we suppose that all containment restrictions are cancelled. In Fig. 3 we consider a larger activation time frame .

Fig. 2

Test 1. Evolution of the fraction of infected (left) and recovered (right) based on the feedback constrained model (15) for , a perception function , and several penalizations . The choice corresponds to the unconstrained case. In last row the normalized case , with

Fig. 3

Test 1. Evolution of the fraction of infected (left) and recovered (right) based on the feedback constrained model (15) for , a perception function , and several penalizations . The choice corresponds to the unconstrained case. In last row the normalized case , with

With the choice of the perception function , , we can observe how the control term is able to flatten the curve even if, as expected, the case gives rise to a weaker control action. To make the two controls, and , comparable for the same penalization factor we also consider the normalized case, where , and the constant is a normalization constant such thatwith an estimate of the maximum number of infected in absence of control. In particular, in our test case, from , assuming we obtain .

Note that, if the activation time is too short the control is not able to significantly change the total number of infected (and therefore recovered). On the other hand, by enlarging the activation time in combination with a sufficiently small penalty constant, the peak infection is not only reduced, but the total number of infected people is decreased. To achieve this, the control should be kept activated for a sufficiently long time and with the right intensity in a kind of plateau regime where there is a perfect balance between the containment effect and the spread of the disease. On the contrary, if the control is too strong, the majority of the population remains susceptible and consequently the disease will start spreading again forming a second wave after the containment policy is removed. Similar conclusions (that may appear counterintuitive) have been shown also by other authors (see for example (Britton et al. 2020; Lunelli et al. 2009)).

The cost functional depends on the value of q and can be evaluated summing up contributions in (9) with explicit form of the control given by (12). In Fig. 4 the cost of the two interventions is compared. We can see how a higher cost is associated with that can be obtained with the control only for smaller penalizations. As expected, the normalized case essentially realigns the cost of interventions. Then, in Fig. 5 we compare the performance of the two controls in with a deactivation time . We consider for both and normalized. It can be observed that there is a minimum control horizon for both strategies, in order to avoid the onset of a second infection peak. A sufficiently long control horizon is therefore necessary to reduce the impact of the infection.

Fig. 4

Test 1. Evaluation of the cost functional J(u) for the dynamics in Fig. 3

Fig. 5

Test 1. Evolution of the of the fraction of infected (left) and recovered (right) based on the feedback constrained model (15) for a perception function , and normalized, for different control actions in with a deactivation time and a fixed penalization

Test 2: Impact of uncertain data on the epidemic outbreak

Next we focus on the influence of uncertain quantities on the controlled system with homogeneous mixing focusing on available data for COVID-19 outbreak in Italy, see del Consiglio dei Ministri (2020). The estimation of epidemiological parameters is a very difficult problem that can be addressed with different approaches (Capaldi et al. 2012; Chowell 2017; Roberts 2013). It is worth to mention that, in the case of COVID-19, the number of infected and recovered has been largely underestimated, especially in the early phases of the epidemic, see (Jagodnik et al. 2020; Mizumoto et al. 2020). Here, we restrict ourselves to identifying the deterministic parameters of the model through a suitable data fitting procedure, considering possible deviations due to such underestimates as part of the subsequent uncertainty quantification process.

The data fitting process

In details, we have adopted a two-level approach in estimating the parameters in absence of uncertainties. In the phase preceding the lockdown we estimated the epidemic parameters in an unconstrained regime, where we assumed no social containment procedure was activated. This estimate was then kept int the subsequent lockdown phase where we estimated as a function of time the value of the control penalty parameter. Both these two calibration steps were analyzed under the assumption of homogeneous mixing, therefore model (15) has been used in lockdown phase.

First, we estimated the parameters in the time interval solving a lest square problem based on the minimization of the relative norm of the difference between the reported number of infected and recovered , and the theoretical evolution of the unconstrained model whose solution at time is indicated with I(t) and R(t). More precisely, we considered the following minimization problemwhere and denotes the relative norm over a time horizon [t, s].

Problem (30) has been solved with the constraints and . Indeed, according to several studies the time to viral clearance during the early phase of the epidemic, corresponding to the time from the first positive test to the first negative test, can approximately span in average from 10 to 24 days, see (Chen 2020; Gatto et al. 2020; Lavezzo 2020).

At the end of the above optimization process we obtained the values , computed by averaging the optimization results with and . The choice of a small value for is due to the low reliability on the recovered data at this early stage.

Next, we estimate the penalization in time by solving (15), in the lockdown time interval and for a sequence of time steps , the corresponding least square problems in where are integers, and where we fix the values estimated in the first optimization step. In details, we solve the following minimization problemover a window of seven days corresponding to and for regularization along one week of available data. Both minimization problems (30–31) have been solved testing various numerical methods in combination with adaptive solves for the systems of ODEs. The results have been obtained using Matlab functions fmincon in combination with ode45. The available data start on February 24 2020, when moderate social restrictions were enforced by the Italian government, and since the lockdown started on March 9 2020, thus we considered (days).

The corresponding time dependent values for the expected penalization for a perception function are reported in Fig. 6. After an initial adjustment phase the penalty terms converge towards a constant value that we can assume as fixed in predictive terms for future times in a lockdown scenario. This is consistent with a situation in which society needs a certain period of time to adapt to the lockdown policy.

Fig. 6

Test 2. Estimated control penalization terms over time from reported data on number of infected and recovered in the case of COVID-19 outbreak in Italy

Remark 2

Finally, we remark that the data fitting procedure is easily generalizable to epidemic models with additional compartments (Gatto et al. 2020; Giordano et al. 2020). Let us denote with the compartments that are related to the reported data , , as the number of actual cases, hospitalized, deaths, etc. In the first step one minimizes a weighted norm in the formwhere are suitable weights such that and is the vector of the model parameters that need to be estimated. Since the problem may admit multiple minima leading to unrealistic solutions one usually perform the above optimization process under constraints on the range of values of some parameters such as recovery rate, incubation period, etc.

An important difference, compared to a simple SIR compartmentalization, is the presence of compartments that are not data-driven as exposed, pre-symptomatic, asymptomatic, etc. that prevent the realization of the data fitting since their values are unknown. A way to overcome this difficulty is to solve the differential model starting from an unknown time using as initial data the presence of a single individual in the first compartment promoting the infection (for example the exposed) and zero individuals in all other compartments. The idea is to simulate the early phase of the epidemic, by optimizing also the initial time in problem (32).

After this, in the second step one considers the corresponding feedback controlled model (see Appendix A) and solves the minimization problemwhere is the vector of penalization terms that need to be estimated. Even in this case, assumptions on the range of values of the penalization terms may be necessary to avoid unrealistic solutions. Finally, we underline that most epidemic models are not data-driven, but one can always assume that the total number of reported cases underestimates the actual number of cases and perform the above data fitting to obtain a lower bound on the evolution of the epidemic. Then including a suitable data uncertainty, as in the present work, allows the recovery of realistic scenarios for the pandemic progression.

Introducing data uncertainty

To account for the impact of uncertainties in the data and parameters we then consider a two-dimensional random variable with independent components such thatandwhere are chosen distributed as symmetric Beta functions in [0, 1], and are the initial number of reported cases and recovered taken from (GitHub: https://github.com/pcmdpc/COVID-19) on February 24, 2020. The choice of a Beta distribution for is coherent with other authors (Roberts 2013; Xiu 2010). However, different probability distribution functions may be considered if additional information on the nature of the uncertainties are available.

It should be noted that the estimated reproduction number, computed in the first optimization step, corresponds to , which is at the upper limit of the values reported in the literature for COVID-19 (Jagodnik et al. 2020; Zhang et al. 2020; Liu et al. 2020). Being aware of the limitations of the data fitting on reported data, to consider a more realistic range of values we assumed a stochastic dependence in and taking into account the faster recovery of asymptomatic individuals (Gatto et al. 2020) and the fact that asymptomatic individuals might be slightly less infectious than symptomatic cases (Sayampanathan et al. 2021). In the simulations we take , and in (35). Under this assumptions the reproduction number covers a range of values approximatively in [3.13, 6.27] with an expected value around 4.25.

The reproduction number in the feedback controlled model is estimated fromIn (36) the time is the lockdown time, in the case under study March 9th, and the indicator function. The feedback control u(t) is defined from (29) in the case of homogeneous mixing and assuming as in (26), where is the total number of infected reported at time t. This leads toin agreement with the fact that the confinement restrictions have been implemented accordingly to the reported data. Note that, otherwise the action of the uncertainty translates into the control and leads to the unrealistic effect that the largest is the number of unreported infected the largest is the action of the control in the system.

In Fig. 7 we report the expected value of together with the and confidence levels with respect to the variable . The estimated reproduction number relative to data fitting is reported with x-marked symbols and represents an upper bound for . The results show that the reproduction number, thanks to the containment actions, has been drastically reduced and its expected value fell below one between March 23rd and March 29th for both and . After March 29th the observed is stably below unity. Note that, in both controls the results are very similar, without any need of renormalization for due to the data fitting process.

Fig. 7

Test 2. Estimated reproduction number from the feedback controlled model with uncertain data (35) for a perception function , (left) and (right) together with the confidence bands. We mark with dash-dotted green lines the days in which the lower band and the expected fell below one, and with x-markers the estimated reproduction number relative to data fitting

Next we considered the evolution of the uncertain number of infected. In the following we assumed a strongly underestimated initial number of infected (including asymptomatic), taking so that the reported infected along the time horizon of the simulation represent approximately a portion of the total infected persons computed by the feedback controlled model. This is in accordance with the WHO findings that around of infected are asymptomatic1 and with the results of preliminary serological campaigns promoted in Italy2.

In Fig. 8 we represent the evolution of the expected value of the number of infected obtained by the controlled model with perception function in presence of uncertain contact and recovery rates (35) and initial uncertain data (34) assuming and . We represent the expected values of the current cases (left) and of the total cases (right) along with the confidence level with respect to the variables and . The shaded color band is relative to the variability in whereas the shaded gray band is relative to the variability in . The bars below the graph are the reported values of the number of infected on which the model has been calibrated. The results with do not highlight significant differences with respect to the case and therefore have been omitted.

Fig. 8

Test 2. Evolution of expected current cases (left) and of the expected total cases (right) and their 95% confidence bands with respect to (shaded color) and (shaded gray) for the feedback controlled model with perception function , and uncertain initial data (34–35)

Test 3: The effect of social contacts in the population

Subsequently, we analyze the effects of the inclusion of age dependence and social interactions in the above scenario. More precisely we consider the social interaction rate , recovery rate and uncertain initial number of infected. These functions were normalized using the estimated parameters and in accordance withwhere f(a) is the age distribution with , .

The age dependent social interaction rate is defined as follows,where , thus for we recover the homogeneuos mixing, whereas for we have a full social mixing behavior. The social interaction function, , accounts for the interactions due to specific activities and is defined by (55) in Appendix C. However, since after the discovery of the first case (February 21), schools, and many places of aggregation were closed in most regions of Northern Italy, we assume that is 0 from February 24 onwards, while is reduced by a factor one-half from March 9 onwards.

The choice of age dependent recovery rate involves a certain degree of arbitrariness, nevertheless it is reasonable to account such heterogeneity as observed in different studies (Faes et al. 2020; Voinsky et al. 2020; Verity 2020; Wang et al. 2020; Russell et al. 2020; Paradisi and Rinaldi 2020). In order to account fast recovery rate of young people, and slow recovery of the eldest we chose to be constant up to a specific age and then a decreasing function of the age. We express mathematically the recovery rate aswith and such that (38) holds.

We divided the computation time frame into two zones and used different models in each zone, in accordance with the policy adopted by the Italian Government. The first time interval defines the period without any form of containment from 24 February to 9 March, the second the lockdown period from 9 March. In the first zone we adopt the uncontrolled model with homogeneous mixing. Hence, in the second zone we compute the evolution of the feedback controlled age dependent model (27–29) with matching (on average) interaction and recovery rates accordingly to (38) and with the estimated value of the control penalization as reported in Fig. 6 until April 30. After April 30 the computation advances in time using as penalization term the constant asymptotic value reached by . The initial values for the age distributions of susceptible and infectious individuals are shown in Fig. 9 in agreement with reported data3.

Fig. 9

Test 3. Distribution of age in Italy (left) and distribution of infected (right) together with the corresponding continuous approximations

In Fig. 10 we report the results of the expected number of infected with the related confidence bands in case of homogeneous mixing and different levels of social mixing () for the constant recovery rate . The homogeneous mixing hypothesis leads to a lower estimate of the maximum number of infected and shows a slower decay over time in the case , whereas for the decay is comparable. In Fig. 11 we compare the case of constant and age-dependent recovery rates for the social mixing scenario. The expected number of infected are shown in the top row, bottom row depicts the total number of recovered and infected people. The heterogeneity of the recovery rate makes the epidemic more persistent and causes an increase in the total number of cases with respect to the homogeneous recovery rate. Finally, in Fig. 12 we report the expected age distribution of infectious individuals in time for . It is evident the effect of the age dependent recovery rate in the increase of cases among the oldest part of the population. Note that, this effect is partially compensated by the strength of the social mixing which reduces interaction among the elderly.

Fig. 10

Test 3. Expected number of infected in time for the perception function , (left) and (right) and a constant recovery rate together with the confidence bands for homogeneous mixing (), mild social mixing and full social mixing ()

Fig. 11

Test 3. Expected number of infected and total cases of infected and recovered in time for a perception function , (left) and (right) together with the confidence bands for the social mixing scenario with age-independent or age-dependent recovery rate

Fig. 12

Test 3. Expected age distribution of infectious individuals for a perception function with mild (left) and full (right) social mixing. In the top row and in the bottom row defined in (40)

Test 4: Reducing the epidemic through relaxed social containment

One of the major problems in the application of very strong containment strategies, such as the lockdown applied in Italy, is the difficulty in maintaining them over a long period, both for the economic impact and for the impact on the population from a social point of view. In order to analyze sustainable control strategies, therefore, it is necessary to resort to models with a social structure and control methods based on specific forms of social distancing that allow the economy to restart and the population to dedicate itself, albeit in a limited way, to its pre-pandemic activities.

In accordance with the interaction function introduced in the Appendix C, we considered the following age dependent relaxation functionwhere the interval defines the age group related to a stronger relaxation of the restrictions (in the sequel we assume ), and the parameters characterize the intensity of the heterogeneity of the relaxation process over the different age classes. Hence we relax the control parameter defined in (29) according to The evolution of the infection is considered within two different relaxation times, at May 4 and at June 3, as actuated by the Italian Government during the first wave of the pandemic. We report in Table 1 the specific choice of the parameter and associated to different periods. Note that, these values are directly related with an increase of the disease transmission rather than to a relaxation of restrictions. In fact, it is clear that under suitable safety precautions a relaxation of restrictions may not contribute to the progress of the epidemic.

Table 1

Reduction of the feedback control (42) over different time periods due to the relaxation of the lockdown processes by the choice of the parameter and of the age dependent function defined in (41)

	Until March 9	March 9–May 3	May 4–June 2	from June 3
p_s×100%	–	0%	15%	20%
p_w×100%	–	0%	5%	10%

In Fig. 13 (top row) we report the evolution of the age-controlled model with perception function , with mild social mixing and with full social mixing and homogeneous recovery rate. Each plot compare the evolution with relaxation of the containment measure and without, respectively indicated with and E[I]. It is easily observed how the relaxation process increases the number of infected persons. Although the expected value remains under control, the wider confidence bands highlight the risk of a resumption of the epidemic. Of course, a further relaxation of the values in Table 1 will lead to a higher risk of restarting the epidemic wave. In Fig. 13 (bottom row) we report also the evolution of the age-dependent expected number of infected individuals, both for mild social mixing and full social mixing. We remark that, as expected, the relaxation process focuses the increase of infections in the interval characterized by .

Fig. 13

Test 4. Expected number of infected with relaxed control (42) characterized by Table 1, perception function , and homogeneous recovery rate, using mild social mixing (left) and with full social mixing (right). In the bottom row the corresponding expected age distribution of infectious individuals is reported

As can be seen, a gradual strategy can keep the average number of infections under control and have an outcome comparable to the fully controlled model at a potentially lower social cost. The timing and intensity of interventions, however, are crucial to prevent the restart of the epidemic wave.

Conclusions

Quantifying the impact of uncertain data in the context of an epidemic emergency is essential in order to design appropriate containment measures. Such containment measures, implemented by several countries in the course of the COVID-19 epidemic, have proved effective in reducing the reproduction number to below or very close to one. These large-scale non-pharmaceutical interventions vary from country to country but include social distancing (banning large mass events, closing public places and advising people not to socialize outside their families), closing borders, closing schools, measures to isolate symptomatic individuals and their contacts, and the large-scale lock-down of populations with all but essential prohibited travel.

One of the main problems is the sustainability of these interventions, which until the introduction of a vaccine will have to be maintained in the field for long periods. However, estimating the reproductive numbers of SARS-CoV-2 is a major challenge due to the high proportion of infections not detected by health care systems and differences in test application, resulting in diverse proportions of infections detected over time and between countries. Most countries initially had only the capacity to test a small proportion of suspected cases and tests were reserved for severely ill patients or high risk groups. The available data therefore offer a systematic partial overview of trends.

In this article, starting from a SIR-type compartmental model with social structure, we developed new stochastic mathematical models describing the actions of a government agency to contain infections among the population in presence of an uncertain number of infectious individuals. More precisely, in the present model, the state of the infected is defined by a multi-dimensional time-dependent function characterized by the age, which proved essential in the description of the COVID-19 pandemic, and by a systematic uncertainty which permits to avoid additional sub-compartmentalization of the infectious population.

These assumptions allows to derive a socially structured model that contains the control action in feedback form based on the perception of the policy maker of the spread of the disease. Subsequently, the uncertainty in the model has been approximated by expanding the state variables into orthogonal polynomials in the random space, reducing the problem to a set of deterministic equations for the distribution of the solution through the course of the epidemic. The resulting controlled stochastic dynamical system is then solved using this deterministic formulation, which in the case of sufficiently regular uncertainty distributions allows efficient and accurate estimations of the random parameters.

The numerical simulations, performed using data from the recent COVID-19 outbreak in Italy, show, on the one hand, the model’s ability to characterize the presence of the asymptomatic population trough the introduction of an uncertainty in the number of infected and in their epidemic role, and, on the other hand, in presence of control to well describe the effects of non pharmaceutical interventions aimed at flattening the infection curve. In particular, identifying some scenarios in agreement with government actions, containment measures by the population based on the resumption of certain occupational activities characterized by specific age groups and social interaction matrices were studied.

Further studies in this direction will aim to consider more comprehensive epidemic models, including the effects of other clinical compartments of interest, along with generalization of the control term to different objective functions or to the case of multiple control terms for each social activity, in order to design optimal strategies to mitigate the overall epidemic impact.

Table 2

The different polynomial expansions connected to the probability distribution of the random component ,

Probability law	Expansion polynomials	Support
Gaussian	Hermite	(-∞,+∞)
Uniform	Legendre	[a, b]
Beta	Jacobi	[a, b]
Gamma	Laguerre	[0,+∞)
Poisson	Charlier	N

Table 3

The parameters defining the details of the interaction functions used in the simulations

Contact function	Parameters
β_F(a,a_∗)	β_0	σF2	λ_F,1	λ_F,2	μ	C_±μ	C_0	ω_F
	0.04	0.125	0.02	100	0.3	0.5	1	1
β_E(a,a_∗)	β_0	σE2	λ_E					ω_E
	0.04	0	0.105					0
β_P(a,a_∗)	β_0	σP2	λ_P					ω_P
	0.04	0.00125	0.4					0.5