Control with uncertain data of socially structured compartmental epidemic models
CControl with uncertain data of socially structuredcompartmental epidemic models
Giacomo Albi ∗ Lorenzo Pareschi † Mattia Zanella ‡ Abstract
The adoption of containment measures to reduce the amplitude of the epidemic peak is akey aspect in tackling the rapid spread of an epidemic. Classical compartmental models mustbe modified and studied to correctly describe the effects of forced external actions to reducethe impact of the disease. In addition, data are often incomplete and heterogeneous, so a highdegree of uncertainty must naturally be incorporated into the models. In this work we addressboth these aspects, through an optimal control formulation of the epidemiological model inpresence of uncertain data. After the introduction of the optimal control problem, we for-mulate an instantaneous approximation of the control that allows us to derive new feedbackcontrolled compartmental models capable of describing the epidemic peak reduction. Theneed for long-term interventions shows that alternative actions based on the social structureof the system can be as effective as the more expensive global strategy. The importance ofthe timing and intensity of interventions is particularly relevant in the case of uncertain para-meters on the actual number of infected people. Simulations related to data from the recentCOVID-19 outbreak in Italy are presented and discussed.
Keywords : Compartmental model, uncertainty quantification, social structure, optimal con-trol, containment policies, COVID-19.
From the digital tracking systems of the Koreans to UK’s initial choice of not wanting to doanything to counter the spread of the virus, passing through the militarized quarantines of theChinese and those less authoritarian and more involved of the Italians, the reaction of differentcountries to the COVID-19 outbreak has shown a series of very different approaches that can alsobe 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 ofthe pandemic’s inexorable progress had to take strong containment measures. There’s a graph thathas become the symbol of the COVID-19 pandemic most of all. It shows in a simple and intuitiveway the importance of slowing down the spread of an epidemic as much as possible (”flattening thecurve”), so that the healthcare system can take care of all the sick without collapsing. Its successhas helped to save many lives, raising awareness of good practices to slow down an epidemic: stayas 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 eco-nomic costs and thus policy makers may not be able to maintain them for more than a short periodof time. Therefore, a modelling approach based on a limited time horizon that takes into accountthe 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 ofvaccinations and medical treatments [5, 7, 14, 15] and only recently the problem has been tackledfrom the perspective of non-pharmaceutical interventions [29, 32]. In addition, data collected by ∗ Computer Science Department, University of Verona, Italy. ( [email protected] ) † Mathematics and Computer Science Department, University of Ferrara, Italy. ( [email protected] ) ‡ Mathematics Department, University of Pavia, Italy. ( [email protected] ) a r X i v : . [ q - b i o . P E ] A p r overnments are often incomplete and heterogeneous, so a high degree of uncertainty needs tobe incorporated into predictive models [9, 11, 16, 25, 31, 38]. This is the case of the spreading ofCOVID-19 worldwide, which have been often mistakenly underestimated due to a combinationof factors, including deficiencies in surveillance and diagnostic capacity, and the large number ofinfectious but asymptomatic individuals [25, 31, 41].For almost a hundred years, mathematical models have been used to describe the spread ofepidemics [27]. The models currently used largely originate from the model proposed by Kermackand McKendrick at the beginning of last century. Even if the model contains strong simplificationassumptions, the concepts introduced through this model are essential to provide a first intuitionon the dynamics of epidemics, an intuition that remains confirmed in more complex models,albeit with numerous modifications (see for example [10, 23]). The model provides for the divisionof 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 thepopulation. One understands how this hypothesis is unrealistic: we are often in contact with peoplefrom our family, our workplace, school class, group of friends and very rarely with those who livein a different place, have different ages and professions. In recent years, therefore, computationalmodels have been developed that try to take into account additional social characteristics ofindividuals in order to arrive at more accurate predictions by keeping, however, the simplicity ofcompartmental models [13, 20–22, 24, 28].In this paper starting from a general compartmental model with social structure we considerthe external action of a policy maker that aims at reducing the spread of the epidemics by ap-plying non pharmaceutical intervention measures, such as social distancing and quarantine. Themathematical problem is formulated as an optimal control problem characterized by a functionalcost 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 allowsus to derive new SIR-type compartmental models capable of describing the epidemic peak reduc-tion. Previously, this type of approach has been used successfully in the case of social models ofconsensus [1–4, 8, 12, 18].The feedback controlled models are subsequently extended to take into account the presenceof uncertain infection parameters and data. In fact, to have reliable forecasts it is of paramountimportance to consider the presence of uncertain quantities as a structural feature of the epidemicdynamics. From a mathematical point of view, we can rely on the methods of uncertainty quanti-fication (UQ) to obtain efficient and accurate solutions based on stochastic orthogonal polynomialsfor the differential model with random inputs [40]. Few results are actually available regardingmethods of UQ in epidemic systems, we mention in this direction [9, 11, 16, 38]. The main ideais to increase the dimensionality of the problem adding the possible sources of uncertainty fromthe very beginning of the modelling. Hence, we extrapolate statistics by looking at the so-calledquantities of interest, i.e. statistical quantities that can be obtained from the solution and thatgive some global information with respect to the input parameter like expected solution of theproblem or higher order moments. Several techniques can be adopted for the approximation ofthe quantities of interest, in this paper we adopt stochastic Galerkin methods that allow to reducethe problem to a set of deterministic equations for the numerical evaluation of the solution inpresence of uncertainties. We refer the interested reader to recent surveys and monographs on thetopic [17, 26, 34, 40].In particular, we consider the case in which the policy maker applies his control based on severalpossible estimators on the actual number of infected people. The need for long-term interventionsshows that alternative actions based on the social structure of the system can be as effective as themore expensive optimal strategy. The importance of the timing and intensity of interventions isparticularly relevant in the case of uncertain parameters on the actual number of infected people.The rest of the manuscript is organized as follows. In Section 2 we introduce the structuredsocial SIR model and formulate the mathematical approach for containment measures to reduce2he spread of the disease. Next, a feedback controlled model used in the subsequent analysisis derived within a short time horizon approximation. In Section 3 we generalize the feedbackcontrolled model to take into account the presence of uncertainties. Section 4 is dedicated tothe presentation of some numerical examples including applications to the COVID-19 epidemicin Italy. The most up-to-date information available to date has been used for these simulations.This illustrates the potential application of the model as a real-time forecasting tool. In separateAppendices we provide details of Galerkin’s stochastic method employed to efficiently address theuncertainties and the social interaction matrices characterizing the contact rates. The starting model in our discussion is a SIR-type compartmental model in epidemiology witha social structure. The presence of a social structure is in fact essential in deriving appropriatesustainable control techniques from the population for a protracted period, as in the case of therecent COVID-19 epidemic.
The heterogeneity of the social structure, which impacts the diffusion of the infective disease, ischaracterized by the vector a ∈ Λ ⊆ R d a characterizing its social state and whose componentssummarize, for example, the age of the individual, its number of social connections or its economicstatus [22,23]. We denote by s ( a , t ), i ( a , t ) and r ( a , t ), the distributions at time t > N , we have that s ( a , t ) + i ( a , t ) + r ( a , t ) = f ( a ) , (cid:90) Λ f ( a ) d a = N, where f ( a ) is the total distribution of the social features defined by the vector a . Hence, werecover the total fraction of the population which belong to the susceptible, infected and recoveredas follows S ( t ) = (cid:90) Λ s ( a , t ) d a , I ( t ) = (cid:90) Λ i ( a , t ) d a , R ( t ) = (cid:90) Λ r ( a , t ) d a . (1)In a situation where changes in the social features act on a slower scale with respect to thespread of the disease, the socially structured compartmental model follows the dynamics ddt s ( a , t ) = − s ( a , t ) 1 N (cid:90) Λ β ( a , a ∗ ) i ( a ∗ , t ) d a ∗ ddt i ( a , t ) = s ( a , t ) 1 N (cid:90) Λ β ( a , a ∗ ) i ( a ∗ , t ) d a ∗ − γ ( a ) i ( a , t ) ddt r ( a , t ) = γ ( a ) i ( a , t ) , (2)where the function β ( a , a ∗ ) ≥ γ ( a ) ≥ β ( a , a ∗ ) = b ( a ) b ( a ∗ ) (cid:82) + ∞ b ( a ) f ( a ) d a (3)with b ( a ) the average number of people contacted by a person with social feature a per unittime. Alternative approaches are based on preferential mixing [13, 21]. Specific examples of age-dependent social interaction matrices are reported in Appendix B.3e introduce the usual normalization scaling s ( a , t ) N → s ( a , t ) , i ( a , t ) N → i ( a , t ) , r ( a , t ) N → r ( a , t ) , (cid:90) Λ f ( a ) d a = 1 , and observe that the quantities S ( t ), I ( t ) and R ( t ) satisfy the conventional SIR dynamic ddt S ( t ) = − (cid:90) Λ s ( a , t ) (cid:90) Λ β ( a , a ∗ ) i ( a ∗ , t ) d a ∗ d a ddt I ( t ) = (cid:90) Λ s ( a , t ) (cid:90) Λ β ( a , a ∗ ) i ( a ∗ , t ) d a ∗ d a − (cid:90) Λ γ ( a ) i ( a , t ) d a , (4)where the fraction of recovered is obtained from R ( t ) = 1 − S ( t ) − I ( t ). We refer to [22, 23]for analytical results concerning model (2) and (4). In the following we will adopt the simplecompartmental model (2) to derive our feedback controlled formulation in presence of uncertainty.The extension to more realistic compartmental models in epidemiology such as SEIR and/orMSEIR can be carried out in a similar way.In order to simplify the description, we will consider the case d a = 1 and set the socialdependence as the age a of the individual because of its importance in epidemic dynamics. It isclear, however, that similar containment procedures can be applied also on the basis of other socialfeatures. We will first formulate the feedback controlled SIR model in the deterministic case andsubsequently extend our approach to the presence of uncertain parameters. In order to define the action of a policy maker introducing a control over the system based onsocial distancing and other containment measures linked to the social structure we consider anoptimal control framework. The choice of an appropriate functional is problem dependent [15].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 interactionsamong individuals with different ages.Thus, we introduce the optimal control problemmin u ∈U J ( u ) := (cid:90) T ψ ( I ( t )) dt + 12 (cid:90) T (cid:90) Λ × Λ ν | u ( a, a ∗ , t ) | dada ∗ dt, (5)subject to ddt s ( a, t ) = − s ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ ddt i ( a, t ) = s ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ − γ ( a ) i ( a, t ) ddt r ( a, t ) = γ ( a ) i ( a, t ) (6)with initial condition i ( a,
0) = i ( a ), s ( a,
0) = s ( a ) and r ( a,
0) = r ( a ).The number of infected individual is measured by a monotone increasing function ψ ( · ) suchthat ψ : [0 , → R + . This function models the policy maker’s perception of the impact of theepidemic by the number of people currently infected. For example ψ ( I ) = I q /q , for q > q = 1. The control aims to minimize this measure of the total infectedpopulation by reducing the rate of interaction between individuals. We consider a quadratic costfor its actuation.Such control is restricted to the space of admissible controls U ≡ (cid:8) u | ≤ u ( a, a ∗ , t ) ≤ min { M, β ( a, a ∗ ) } , ∀ ( a, a ∗ , t ) ∈ Λ × [0 , T ] , M > (cid:9) , L ( s, i, r, p s , p i , p r , u ) = J ( u )+ (cid:90) T (cid:90) Λ p s · (cid:18) ddt s + s ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ (cid:19) da dt + (cid:90) T (cid:90) Λ p i · (cid:18) ddt i − s ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ + γ ( a ) i ( a, t ) (cid:19) da dt + (cid:90) T (cid:90) Λ p r · (cid:18) ddt r − γ ( a ) i ( a, t ) (cid:19) da dt where p s ( a, t ) , p i ( a, t ) , p r ( a, t ) are the associated lagrangian multipliers. Computing the variationswith respect to ( s, i, r ) we retrieve the adjoint system ddt p s = ( p s − p i ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ ddt p i = ψ (cid:48) ( I ( t )) + (cid:90) Λ ( p s ( a ∗ , t ) − p i ( a ∗ , t ))( β ( a ∗ , a ) − u ( a ∗ , a, t )) s ( a ∗ , t ) da ∗ + γ ( a ) p i (7)with terminal conditions p s ( a, T ) = 0 , p i ( a, T ) = 0 and p r ( a, T ) = 0. Note that the contributionof p r ( a, t ) vanishes since the control does not act directly on population R , and the removedpopulation is not considered in the minimization of the functional. The optimality conditionreads ν ( a, a ∗ ) u ( a, a ∗ , t ) = ( p s − p i ) s ( a, t ) i ( a ∗ , t ) . (8)The optimality conditions (7)-(8) are first order necessary conditions for the optimal control u ∗ ( a, t ). In order to be admissible then the control reads u ( a, a ∗ , t ) = max (cid:26) , min (cid:26) p s − p i ν ( a, a ∗ ) s ( a, t ) i ( a ∗ , t ) , φ M,β ( a, a ∗ ) (cid:27)(cid:27) , where φ M,β ( a, a ∗ ) = min { β ( a, a ∗ ) , M } .The approach just described, however, is generally quite complicated when there are uncertain-ties 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 along time horizon seems rather unrealistic in the case of a rapidly spreading disease such as theCOVID-19 epidemic. In this section we consider short time horizon strategies which permits to derive suitable feedbackcontrolled models. These strategies are suboptimal with respect the original problem (5)-(6) butthey have proved to be very successful in several social modeling problems [1–4, 18]. To this aim,we consider a short time horizon of length h > J h ( u ) restricted to the interval [ t, t + h ], as followsmin u ∈U J h ( u ) = ψ ( I ( t + h )) + 12 (cid:90) Λ × Λ ν ( a, a ) | u ( a, a ∗ , t ) | dada ∗ . (9)subject to s ( a, t + h ) = s ( z, a, t ) − hs ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ (10)5 ( a, t + h ) = i ( a, t ) + hs ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ − hγ ( a ) i ( a, t ) . (11)Recalling that the macroscopic information on the infected is I ( t + h ) = I ( t ) + h (cid:90) Λ (cid:20) s ( a, t ) (cid:90) Λ ( β ( a, a ∗ ) − u ( a, a ∗ , t )) i ( a ∗ , t ) da ∗ − γ ( a ) i ( a, t ) (cid:21) da, we can derive the minimizer of J h computing D u J ( u ) ≡
0. We retrieve the following nonlinearequation ν ( a, a ∗ ) u ( a, t ) = hs ( a, t ) i ( a ∗ , t ) ψ (cid:48) ( I ( t + h )) . (12)Introducing the scaling ν ( a, a ∗ ) = hκ ( a, a ∗ ) we can pass to the limit for h →
0. Hence (12) reducesto the instantaneous control u ( a, t ) = 1 κ ( a, a ∗ ) s ( a, t ) i ( a ∗ , t ) ψ (cid:48) ( I ( t )) . (13)The resulting controlled dynamics is retrieved by applying the instantaneous strategy directly inthe continuous system (6) as follows ddt s ( a, t ) = − s ( a, t ) (cid:90) Λ (cid:16) β ( a, a ∗ ) − κ ( a, a ∗ ) s ( a, t ) i ( a ∗ , t ) ψ (cid:48) ( I ( t )) (cid:17) i ( a ∗ , t ) da ∗ (14a) ddt i ( a, t ) = s ( a, t ) (cid:90) Λ (cid:16) β ( a, a ∗ ) − κ ( a, a ∗ ) s ( a, t ) i ( a ∗ , t ) ψ (cid:48) ( I ( t )) (cid:17) i ( a ∗ , t ) da ∗ − γ ( a ) i ( a, t ) (14b) ddt r ( a, t ) = γ ( a ) i ( a, t ) (14c)In what follows we provide a sufficient conditions for the admissibility of the instantaneous controlin terms of the penalization term κ ( a, a ∗ ). Indeed we want to assure that the dynamics preservethe monotonicity of the number of susceptible population s ( a, t ). Proposition 2.1.
Let β ( a, a ∗ ) ≥ δ > , then for all ( a, a ∗ ) ∈ Λ × Λ solutions to (14) are admissibleif the penalization κ satisfies the following inequality δκ ( a, a ∗ ) ≥ s ( a )¯ i ( a ∗ ) ψ (cid:48) ( ¯ I ) , ∀ ( a, a ∗ ) ∈ Λ × Λ (15) where ¯ i ( a ) and ¯ I are respectively the peak reached by the infected of class a and by the totalpopulation.Proof. By imposing the non-negativity of the controlled reproduction rate inside the integral wehave β ( a, a ∗ ) κ ( a, a ∗ ) ≥ s ( a, t ) i ( a ∗ , t ) ψ (cid:48) ( I ( t )) . This inequality has to be satisfied for every time t ≥
0. Second we observe that the number ofsusceptible s ( a, t ) is decreasing in time therefore s ( a ) ≥ s ( a, t ) for all t . Moreover i ( a, t ) reachesa peak before decreasing to 0 (note that this peak can also be in t = 0), say ¯ I for the macroscopicvariable and ¯ i ( a ). Thus, thanks to the monotonicity of ψ ( · ), we can restrict the previous inequalityas follows δκ ( a, a ∗ ) ≥ s ( a )¯ i ( a ∗ ) ψ (cid:48) ( ¯ I ) . In Figure 1 we report the phase diagram of susceptible-infected trajectories for the controlledmodel with homogeneous mixing ddt S ( t ) = − (cid:18) β − κ S ( t ) I ( t ) ψ (cid:48) ( I ( t )) (cid:19) S ( t ) I ( t ) ddt I ( t ) = (cid:18) β − κ S ( t ) I ( t ) ψ (cid:48) ( I ( t )) (cid:19) S ( t ) I ( t ) − γI ( t ) , (16)6 Figure 1: Phase diagram of susceptible-infected trajectories for the controlled SIR-type model withhomogenous mixing and ψ ( I ) = I q /q . Different choices of the penalization term κ are reported.Left plot the case q = 1, right plot q = 2.with ψ ( I ) = I q /q . The dynamic is similar to the classical SIR model but with a nonlinearcontact rate. In particular, the trajectories are flattened when the value of the control is suchthat κβ ≈ S ( t ) I ( t ) ψ (cid:48) ( I ( t )) and the reproductive number is close to zero. Note, however that thisstatus is not an equilibrium point of the system.To understand this, let us observe that an equilibrium state ( S ∗ , I ∗ ) for (16) satisfies theequations ( κβ − S ∗ I ∗ ψ (cid:48) ( I ∗ )) S ∗ I ∗ = 0(( κβ − S ∗ I ∗ ψ (cid:48) ( I ∗ )) S ∗ − κγ ) I ∗ = 0 . An equilibrium point corresponds to the classical state in which we have the extinction of thedisease I ∗ = 0 and S ∗ arbitrary and defined by the asymptotic state of the dynamics [23]. Now,let’s suppose that I ∗ (cid:54) = 0, S ∗ (cid:54) = 0, we can look for solutions where control is able to perfectlybalance the spread of the disease κβ = S ∗ I ∗ ψ (cid:48) ( I ∗ ) κβS ∗ = ( S ∗ ) I ∗ ψ (cid:48) ( I ∗ ) + κγ, consequently κβS ∗ = ( S ∗ ) κβS ∗ + κγ = 0 , which is satisfied only for γ = 0 when κ (cid:54) = 0. Since the beginning of the outbreak of new infectious diseases, the actual number of infectedand recovered people is typically underestimated, causing fatal delays in the implementation ofpublic health policies facing the propagation of epidemic fronts. This is the case of the spreadingof COVID-19 worldwide, often mistakenly underestimated due to deficiencies in surveillance anddiagnostic capacity [25, 37]. Health systems are struggling to adopt systematic testing to monitoractual cases. Moreover, another important epidemiological factor that pollutes the available datais the proportion of asymptomatic [25, 31, 41]. 7mong the common sources of uncertainties for dynamical systems modeling epidemic out-breaks we may consider the following • noisy and incomplete available data • structural uncertainty due to the possible inadequacy of the mathematical model used todescribe the phenomena under consideration.In the following we consider the effects on the dynamics of uncertain data, such as the initialconditions on the number of infected people or the interaction and recovery rates. On the numericallevel we consider techniques based on stochastic Galerkin methods, for which spectral convergenceon random variables is obtained under appropriate regularity assumptions [40]. For simplicity, thedetails of the numerical method that allows to reduce the uncertain dynamic system to a set ofdeterministic equations are reported in the Appendix A. We introduce the random vector z = ( z , . . . , z d z ) whose components are assumed to be independ-ent real valued random variables z k : (Ω , F ) → ( R , B R ) , k = 1 , . . . , d z with B R the Borel set. We assume to know the probability density p ( z ) : R d z → R d z + characterizingthe distribution of z . Here, z ∈ R d z is a random vector taking into account various possible sourcesof uncertainty in the model.In presence of uncertainties we extend the initial modelling by introducing the quantities s ( z , a, t ), i ( z , a, t ) and r ( z , a, t ) representing the distributions at time t ≥ s ( z , a, t ) + i ( z , a, t ) + r ( z , a, t ) = f ( a ) , (cid:90) Λ f ( a ) da = N. Hence, the quantities S ( z , t ) = (cid:90) Λ s ( z , a, t ) da, I ( z , t ) = (cid:90) Λ i ( z , a, t ) da, R ( z , t ) = (cid:90) Λ r ( z , a, t ) da, (17)denote the uncertain fractions of the population that are susceptible, infectious and recoveredrespectively.The social structured model with uncertainties reads ddt s ( z , a, t ) = − s ( z , a, t ) (cid:90) Λ β ( z , a, a ∗ ) i ( z , a ∗ , t ) N da ∗ ddt i ( z , a, t ) = s ( z , a, t ) (cid:90) Λ β ( z , a, a ∗ ) i ( z , a ∗ , t ) N da ∗ − γ ( z , a ) i ( z , a, t ) ddt r ( z , a, t ) = γ ( z , a ) i ( z , a, t ) (18)To illustrate the impact of uncertainties let us consider the simple following example. In the caseof homogeneous mixing with uncertain contact rate β ( z ) = β + αz , α > z ∈ R distributed as p ( z ) the model reads ddt S ( z, t ) = − ( β + αz ) S ( z, t ) I ( z, t ) ddt I ( z, t ) = ( β + αz ) S ( z, t ) I ( z, t ) − γI ( z, t ) , (19)8ith deterministic initial values I ( z,
0) = I and S ( z,
0) = S . The solution for the proportion ofinfectious during the initial exponential phase is [38] I ( z, t ) = I e γ [( β + αz ) S t − t ] , and its expectation E [ I ( · , t )] = I (cid:90) R e γ [( β + αz ) S t − t ] p ( z ) dz = I e γ ( βS t − t ) W ( t ) , (20)where the function W ( t ) = (cid:90) R e γ ( αzS t − t ) p ( z ) dz represents the statistical correction factor to the standard deterministic exponential phase of thedisease I e γ ( βS t − t ) . If z is uniformly distributed in [ − ,
1] we can explicitly compute W ( t ) = γ sinh (( αS t ) /γ ) αS t > , t > . More in general, if z has zero mean then by Jensen’s inequality we have W ( t ) > t >
0, sothat the expected exponential phase is amplified by the uncertainty (see [38]).In a similar way, keeping β constant, but introducing a source of uncertainty in the initial data I ( z,
0) = I + µz , µ > z ∈ R distributed as p ( z ) the solution in the exponential phase reads I ( z, t ) = ( I + µz ) e γ [ βS t − t ] , and then its expectation E [ I ( · , t )] = (cid:90) R ( I + µz ) e γ [ βS t − t ] dz = ( I + µ ¯ z ) e γ ( βS t − t ) , (21)where ¯ z is the mean of the variable z . Therefore, the expected initial exponential growth behavesas the one with deterministic initial data I + µ ¯ z . Of course, if both sources of uncertainty arepresent the two effects just described sum up in the dynamics. In presence of uncertainties the optimal control problem (5)-(6) is modified as followsmin u ∈U J ( u ) := (cid:90) T R [ ψ ( I ( · , t ))] dt + 12 (cid:90) T (cid:90) Λ × Λ ν | u ( a, a ∗ , t ) | dada ∗ dt, (22)being R [ ψ ( I ( · , t ))] a suitable operator taking into account the presence of the uncertainties z .Examples of such operator that are of interest in epidemic modelling are the expectation withrespect to uncertainties E [ ψ ( I ( · , t ))] = (cid:90) R dz ψ ( I ( z , t )) p ( z ) d z (23)or relying on data which underestimate the number of infectedmin z ∈ R dz ψ ( I ( z , t )) . (24)In (22) U the space of admissible controls is defined as U = (cid:110) u | ≤ u ( a, a ∗ , t ) ≤ min { M, min z β ( z , a, a ∗ ) } , ∀ ( a, a ∗ , t ) ∈ Λ × [0 , T ] , M > (cid:111) . ddt s ( z , a, t ) = − s ( z , a, t ) (cid:90) Λ ( β ( z , a, a ∗ ) − u ( a, a ∗ , t )) i ( z , a ∗ , t ) da ∗ ddt i ( z , a, t ) = s ( z , a, t ) (cid:90) Λ ( β ( z , a, a ∗ ) − u ( a, a ∗ , t )) i ( z , a ∗ , t ) da ∗ − γ ( z , a ) i ( z , a, t ) ddt r ( z , a, t ) = γ ( z , a ) i ( z , a, t ) (25)with initial condition i ( z , a,
0) = i ( z , a ), s ( z , a,
0) = s ( z , a ) and r ( z , a,
0) = r ( z , a ).The implementation of instantaneous control for dynamics in presence of uncertainties followsfrom the derivation presented in Section 2.3 and we omit the details. The resulting feedbackcontrol u ( a, a ∗ , t ) reads u ( a, a ∗ , t ) = 1 κ ( a, a ∗ ) R [ s ( · , a, t ) i ( · , a ∗ , t ) ψ (cid:48) ( I ( · , t ))] , (26)and defines the feedback controlled model in presence of uncertainties. In this section we present several numerical tests on the constrained compartmental model withuncertain data. Details of the numerical method used are given in Appendix A. In an attempt toanalyse sufficiently realistic scenarios, in the following examples we will refer to values taken fromItalian data on the COVID-19 epidemic [36]. More precisely, in the first test case we illustrate thebehavior of the model in a simplified approach in the absence of uncertainty and social structureand without trying to reproduce scenarios closely related to current data. In the second testcase, following a progressively more realistic approach, we consider the impact of the presence ofuncertain data in the controlled model with homogeneous social mixing and calibrated on Italiandata. The same setting is then considered in Test 3 taking into account the additional effect on thespreading of the infectious disease given by the social structure of the system described by suitablesocial interaction functions. The final scenario, explored in Test 4, examines the possibility ofreducing the amplitude of the epidemic peak with the application of relaxed confinement measuresrelated to the social structure of the system.
To illustrate the effects of controls introduced that mimic containment procedures, let us firstconsider the case where the social structure is not present. Furthermore, to simplify further themodeling, in this first example we neglect any dependence on uncertain data.We consider as initial small number of infected and recovered i (0) = 3 . × − , r (0) =8 . × − . These normalized fractions refer specifically to the first reported values in the caseof the Italian outbreak of COVID-19, even if in this simple test case we will not try to match thedata in a predictive setting but simply to illustrate the behavior of the feedback controlled model.Based on recent studies [25, 41], the initial infection rate of COVID-19 has been estimated ofan R = β/γ between 2 and 6 .
5. Here, to exemplify the possible evolution of the pandemic weconsider a value close to the lower bound, taking β = 25 and γ = 10, so that R = 2 . t ∈ [0 . , t > t ∈ [0 . , ψ ( I ) = I q /q , q ≥
1, we can observe how the control term is able to flattenthe curve even if, as expected, the case q = 2 gives rise to a weaker control action. Note, however,that if the activation time is too short the control is not able to significantly change the total10igure 2: Test 1 . Evolution of the fraction of infected (left) and recovered (right) based on theconstrained model in the interval t ∈ [0 . ,
1] for ψ ( I ) = I q /q , q = 1 , κ = 10 − , − , − . The choice κ = + ∞ corresponds to the unconstrained case.number of infected (and therefore recovered). On the other hand, by enlarging the activation timein 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 keptactivated for a sufficiently long time and with the right intensity in a kind of plateau regimewhere there is a perfect balance between the containment effect and the spread of the disease. Onthe contrary, if the control is too strong, the majority of the population remains susceptible andconsequently the disease will start spreading again forming a second wave after the containmentpolicy is removed.The cost functional depends on the value of q and can be evaluated summing up contributionsin (9) with explicit form of the control given by (13). In Figure 4 the cost of the two interventionsis compared. We can see how a higher cost is associated with q = 1 that can be obtained with thecontrol q = 2 for weaker penalizations. Then, in Figure 5 we compare the performance of the twocontrols with a deactivation time in the range t ∈ [1 , κ = 10 − for q = 1 and κ = 10 − for q = 2. It can be observed that there is a minimumcontrol horizon for both strategies, in order to avoid the onset of a second infection peak. Asufficiently long control horizon is therefore necessary to reduce the impact of the infection. Next we focus on the influence of uncertain quantities on the controlled system with homogeneousmixing focusing on available data for COVID-19 outbreak in Italy, see [36]. According to recentresults on the diffusion of COVID-19 in many countries the number of infected, and thereforerecovered, is largely underestimated on the official reports, see e.g. [25, 31].The estimation of epidemiological parameters is a very difficult problem that can be addressedwith many different approaches [9, 16, 38]. In our case, we have limited ourselves to identifyingthe deterministic parameters of the model through a standard fitting procedure, considering thepossible uncertainties due to such an estimation as part of the subsequent uncertainty quantific-11igure 3:
Test 1 . Evolution of the fraction of infected (left) and recovered (right) based on theconstrained model in the interval t ∈ [0 . ,
2] for ψ ( I ) = I q /q , q = 1 , κ = 10 − , − , − . The choice κ = + ∞ corresponds to the unconstrained case.ation process. It has been reported, in fact, that deterministic methods based on compartmentalmodels overestimate the effective reproduction number [30].In order to calibrate the models with the reported quantities we solved two separate constrainedoptimization problems in absence of uncertainties. First we estimated β and γ by considering onthe time interval [ t , t u ] the following least square problemmin β,γ (cid:18)(cid:90) t u t | I ( t ) − ˆ I ( t ) | dt (cid:19) / , (27)namely the L norm of the difference between the observed number of infected ˆ I ( t ) and thetheoretical evolution of the unconstrained model ( κ = + ∞ ). Problem (27) has to be solved underthe constraints β ≥ γ ≥
0. This problem is solved over the time span [ t , t u ] where we assumedno social containment procedure was activated.Next, we estimate the penalization κ = κ ( t ) in time by solving over a sequence of time steps t i of size h in the controlled time interval t ∈ [ t u , t c ] the problemmin κ ( t i ) (cid:32)(cid:90) t i + k r ht i − k l h | I ( t ) − ˆ I ( t ) | dt (cid:33) / , (28)under the constraint κ ( t i ) > k l , k r ≥ β and γ estimated in the first optimization step (27). In details, since the available datastart on February 24 2020, when no social restrictions were enforced by the Italian government,and since the lockdown started on March 9 we considered t u = 14 (days). The second fittingprocedure has been activate up to last available data with daily time stepping ( h = 1) and awindow of eight days ( k l = 3, k r = 5) for regularization along one week of available data. To solveproblem (27) we tested different optimization methods in combination with adaptive solvers forthe system of ODEs. Once calibrated the model we estimated β ≈ β e = 31 . γ ≈ γ e = 4 . .5 1 1.5 210 -3 -1 -3 -1 Figure 4:
Test 1 . Evaluation of the costs functional J ( u ) for ψ ( I ) = I q /q , q = 1 (left) and q = 2 (right) on the activation frame t ∈ [0 . ,
2] for the dynamic in Figure 3. We considered threepenalizations κ = 10 − , − , − . (a) Infected (b) Recovered(c) Infected (d) Recovered Figure 5:
Test 1 . Evolution of the constrained model for ψ ( I ) = I q /q and q = 1 (left column) and q = 2 (right column) for a deactivation time of the control in [1 , κ = 10 − for q = 1 and κ = 10 − for q = 2 in order to align the costs of interventions.13 ar 15 Apr 1 Apr 1510 -9 -8 -7 -6 -5 -4 Figure 6:
Test 2 . Estimated control penalization terms over time from reported data on numberof infected people in the case of COVID-19 outbreak in Italy.Figure 7:
Test 2 . Evolution of expected number of infected and their confidence bands for thecalibrated control model with ψ ( I ) = I q /q , q = 1 , µ = 10and uncertain reproduction number (30) with α = 1 in the case of COVID-19 outbreak in Italy.which corresponds to an initial attack rate R ≈ .
25 which agrees with other observations [30].The corresponding time dependent values for the control in the case ψ ( I ) = I q /q are reportedin Figure 6. After an initial adjustment phase the penalty terms converge towards a constantvalue 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 thelockdown policy.In order to have an insight on global impact of uncertain parameters we consider a two-dimensional uncertainty z = ( z , z ) with independent components such that i ( z ,
0) = i (1 + µz ) , z ∼ U ([0 , , µ > , (29)and β ( z ) = β e + αz , z ∼ U ([ − , , α > , (30)where i is the same as in Test 1 and taken from [36] on February 24, 2020 and β e and γ e areestimated from (27). Of course, other probability distributions p ( z ) = p ( z ) p ( z ) defined onlyfor z , z ≥ Test 2 . Estimated reproduction number R from the controlled model for ψ ( I ) = I q , q = 1 and q = 2 together with the confidence bands in the case of COVID-19 outbreak in Italy.parameter for each random variable, so we limit ourselves to the simple situation of uniformuncertainty.In all the considered evolutions we adopted a stochastic Galerkin approach, see Appendix Awith M = 10 and a fourth order Runge-Kutta method for the time integration. The feedbackcontrolled model has been computed using as estimation of the total number of infected reported,namely β ( z ) − k min z ∈ R dz { S ( z , t ) I ( z , t ) q } , (31)in agreement with the fitting procedure obtained from the lower bounds of the uncertain initialdata.In Figure 7 we represent the evolution of the expected value of the number of infected obtainedby the controlled model in the presence of initial random data (29) and uncertain contact frequency(30). The value µ = 10 have been chosen accordingly to the WHO suggestions that around 80%are asymptomatic . The uncertainty in the contact rate has been modeled taking α = 1.We represented the expected values of the number of infected along with the confidence bandsobtained from the overall variance in z . The range of true infected as seen increases dramaticallyin both cases q = 1 and q = 2 with a strong uncertainty on the effective numbers. The barsbelow the graph are the reported values of the number of infected on which the model has beencalibrated.From (31) we can estimate the effect of the lockdown policies on the value of R in time. Theeffect of the uncertain initial estimated value R is translated in the confidence bands from z reported in the left plot of Figure 8. The results show that the R reproduction number, thanksto drastic containment actions, has been drastically reduced and its expected value is now justbelow one. This shows the importance to continue with containment policies to avoid a restart ofthe spread of the epidemic. We analyze the effects of the inclusion of age dependence and social interactions in the abovescenario. More precisely we consider the social interaction functions reported in Appendix Band uncertain initial number of infected. These functions were normalized using the previouslyestimated parameters β e and γ e in accordance with β e = 1 a (cid:88) j ∈A (cid:90) Λ × Λ β j ( a, a ∗ ) da da ∗ , γ e = 1 a max (cid:90) Λ γ ( a ) da, (32) Q&A: Similarities and differences COVID-19 and influenza.
20 40 60 80 10000.0050.010.0150.02 0 20 40 60 80 10000.0050.010.0150.02
Figure 9:
Test 3 . Distribution of age in Italy (left) and distribution of infected (right) togetherwith the corresponding continuous approximations .Figure 10: Test 3 . Expected number of infected in time for ψ ( I ) = I q , q = 1 (left) and q = 2(right) together with the confidence bands in the case of COVID-19 outbreak in Italy for homo-geneous mixing or social mixing (top row) and, in the social mixing scenario, age-independent orage-dependent recovery rate (bottom row). The abscissae are measured in days starting from thebeginning of the epidemic. 16igure 11: Test 3 . Expected number of recovered in time for ψ ( I ) = I q , q = 1 (left) and q = 2(right) together with the confidence bands in the case of COVID-19 outbreak in Italy for the socialmixing scenario, age-independent or age-dependent recovery rate. The abscissae are measured indays starting from the beginning of the epidemic.where Λ = [0 , a max ], a max = 100 and where γ ( a ) accordingly to recent studies on age-relatedrecovery rates [39] has been chosen as a decreasing function of the age as follows γ ( a ) = γ e + Ce − ra , (33)with r = 5 and C ∈ R such that (32) holds. Clearly, this choice involves a certain degree ofarbitrariness since there are not yet sufficient studies on the subject, nevertheless, as we will seein the simulations, it is able to reproduce more realistic scenarios in terms of age distribution ofthe infected without altering the behaviour relative to the total number of infected.We divided the computation time frame into two zones and used different models in eachzone, in accordance with the policy adopted by the Italian Government. The first time intervaldefines the period without any form of containment from 24 February to 9 March, the secondthe lockdown period from 9 March. In the first zone we adopted the uncontrolled model withhomogeneous mixing for the estimation of epidemiological parameters. Hence, in the second zonewe compute the evolution of the feedback controlled age dependent model (25)-(26) with matching(on average) interaction and recovery rates (32) and with the estimated control penalization κ ( t )as reported in Figure 6 until April 23rd. After April 23rd the computation advances in time usingas penalization term the constant asymptotic value ¯ κ reached by k ( t ). The initial values for theage distributions of susceptible and infectious individuals are shown in Figure 9 in agreement withreported data .In Figure 10 we report in the top line the results of the expected number of infected withthe related confidence bands in case of homogeneous mixing and social mixing for the constantrecovery rate γ e . At the bottom we compare the case of constant and age-dependent recoveryrates for the social mixing scenario. In general, the homogeneous mixing hypothesis leads to anunderestimation of the maximum number of infected and shows a slower decay over time, thelatter effect also found in the presence of an age-dependent recovery rate. The expected totalnumber of recovered people is shown in Figure 11. Finally in Figure 12 we report the expected agedistribution of infectious individuals in time. We remind that all the simulations presented in thissection have been obtained assuming to maintain the lockdown measures introduced on March 9for all subsequent times. One of the major problems in the application of very strong containment strategies, such as thelockdown applied in Italy, is the difficulty in maintaining them over a long period, both for the Source ISTAT ( ) and Istituto Superiore di Sanit`a ( ) Test 3 . Expected age distribution of infectious individuals over a long time horizon of200 days in the social mixing scenario with γ = γ e (left) and age dependent recovery γ ( a ) definedin (33).economic impact and for the impact on the population from a social point of view. In order toanalyze sustainable control strategies, therefore, it is necessary to resort to models with a socialstructure and control methods based on specific forms of social distancing that allow the economyto restart and the population to dedicate itself, albeit in a limited way, to its pre-pandemicactivities.In accordance with the interaction function characterizing the productive activities introducedin the Appendix B, see also Figure 14, we considered the following age dependent penalization κ ( a, a ∗ ) = (cid:26) κ w a, a ∗ ∈ Λ P κ s elsewherewhere Λ P defines the typical age group related to productive activities, for example Λ P = [20 , /κ s > /κ w characterize two control actions related to a strong and a weaker containmentof social distance. We will assume 1 /κ s = p s / ¯ κ where ¯ κ is the asymptotic value calculated for thepenalization parameter in the lockdown period, while κ w has been relaxed compared to κ s , i.e.1 /κ w = p w / ¯ κ , with 1 ≥ p s ≥ p w .In Figure 13 we report the evolution of the age-controlled model in the case q = 1. We haveconsidered two possible scenarios that reproduce the evolution of the infected in the hypothesis ofrelaxation of the lockdown measures at May 4 as foreseen by the Italian Government.In the first scenario, on the left, we consider two intervals that simulate the two differentcontrolled periods measured in days from the beginning of the epidemic, T = [14 ,
69] with p s = p w = 1, and T = [70 , p w = 0 . p s = 0 .
9. In Figure 13 (left) we compare this choice with a situation where the relaxation of theso-called phase 2 control measures in the interval T are further loosened by chosing p w = 0 . p s = 0 .
6. The respective numbers of infected in these two situations are denoted by T u , I u . Itis easily observed how we systematically increase the number of infected by successive relaxationsof the initial strict social-distancing measures. In Figure 13 (bottom, left) we report also theevolution of age-dependent expected number of infected individuals.The second scenario assumes that after an initial opening in which few productive activitiesare allowed, the government will gradually loosen the containment procedures with a gradualapproach after a certain period of time. In this case we assumed three control zones: T = [14 . p s = p w = 1 , the interval T = [70 .
98] with p s =0.9, p w =0.6 and finally a third period T = [99 . p s =0.8, p w =0.5. This last period will correspond to a further relaxation ofconfinement measures on June 1st. The solution related to the number of infected is indicated with I u and is shown in Figure 13. (right). For comparison we also report the full lockdown solutiondenoted with I u . In the lower part of the figure on the right the corresponding distributions ofinfected individuals by age over time are also reported. As can be seen, a gradual strategy allows18igure 13: Test 4 . Top row: evolution of the total number of infected in a two-zone scenarioidentified by κ and κ (left) and in a three zones scenario identified by κ and κ . Bottom row:age distribution of infected in the two-zones (left) and three-zones (right) scenario.to contain the peak of contagions and to have a result comparable with the fully controlled modelat a much lower social cost. The timing and intensity of the interventions, however, are crucial toprevent a restart of the epidemic wave. Quantifying the impact of uncertain data in the context of an epidemic emergency is essentialin order to design appropriate containment measures. Such containment measures, implementedby several countries in the course of the COVID-19 epidemic, have proved effective in reducingthe R reproduction number to below or very close to one. These large-scale non-pharmaceuticalinterventions vary from country to country but include social distance (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-scalelock-down of populations with all but essential prohibited travel.One of the main problems is the sustainability of these interventions, which until the introduc-tion of a vaccine will have to be maintained in the field for long periods. However, estimating thereproductive numbers of SARS-CoV-2 is a major challenge due to the high proportion of infec-tions not detected by health care systems and differences in test application, resulting in diverseproportions of infections detected over time and between countries. Most countries currently haveonly the capacity to test a small proportion of suspected cases and tests are reserved for severely19ll patients or high risk groups. The available data therefore offer a systematic partial overview oftrends.In this article, starting from a SIR-type compartmental model with social structure, we havedeveloped new mathematical models describing the actions of a government agency to control thepopulation in order to reduce the estimated number of infected people in the presence of uncertaindata. This hypothesis allows to derive a model that contains the control action in feedback formbased on the perception of the policy maker of the spread of the disease. Subsequently the modelhas been effectively solved in the presence of uncertain data, expanding the state variables intoorthogonal polynomials in the uncertainty space, reducing the problem to a set of deterministicequations for the distribution of the solution through the course of the epidemic. The resultingcontrolled dynamic system then results in a deterministic stochastic solution that enables efficientestimation of uncertain parameters.The numerical simulations, carried out using data from the recent COVID-19 outbreak in Italy,show, on the one hand, the model’s ability to well describe lockdown scenarios aimed at flatteningthe infection curve and, on the other hand, how the high uncertainty of the data on the number ofinfected people makes it very difficult to provide long-term quantitative forecasts. By identifyingsome plausible scenarios in accordance with the literature, sustainable containment measures bythe population based on the resumption of certain occupational activities characterized by spe-cific age groups and social interaction matrices have been studied. The results demonstrate theeffectiveness of such approaches based on the social structure of the system capable of achievingresults similar to much more restrictive policies in reducing the risk that the virus may return tospread when the restrictions are lifted but at a significantly lower socio-economical cost.Further studies in this direction will be aimed at considering more realistic epidemic modelsthan the one analyzed in this work together with multiple control terms specific for each socialactivity in order to design optimal strategies to mitigate the overall epidemic impact. Acknowledgements
G. A. and L. P. acknowledge the support of PRIN Project 2017, No. 2017KKJP4X Innovativenumerical methods for evolutionary partial differential equations and applications. This researchwas partially supported by the Italian Ministry of Education, University and Research (MIUR)through the “Dipartimenti di Eccellenza” Programme (2018-2022) – Department of Mathematics“F. Casorati”, University of Pavia. M.Z. is member of GNFM (Gruppo Nazionale per la FisicaMatematica) of INdAM (Istituto Nazionale di Alta Matematica), Italy.
A Stochastic Galerkin approximation
In this Appendix we give the details of the stochastic Galerkin (sG) method used to solve thefeedback controlled system (25)-(26) with uncertainties. To this aim, we consider a random vec-tor z = ( z , . . . , z d ) with independent components and whose distribution is p ( z ) : R d z → R d z + .The stochastic Galerkin approximation of the differential model (25)-(26) is based on stochasticorthogonal polynomials and provides a spectrally accurate solution under suitable regularity as-sumptions, see [40]. We consider the linear space P M of polynomials of degree up to M generatedby a family of polynomials { Φ h ( z ) } M | h | =0 that are orthonormal in the space L (Ω) such that E [Φ h ( · )Φ k ( · )] = δ kh , ≤ | k | , | h | ≤ M being k = ( k , . . . , k d ) a multi-index, | k | = k + · · · + k d z with δ kh the Kronecker delta function, and E [ · ] the expectation with respect to p ( z ). The construction of the polynomial basis { Φ h ( z ) } M | h | =0 depends on the distribution of the uncertainties and must be chosen in agreement with the Askeyscheme [40]. We summarise in Table 1 several polynomials bases in connection with the law of arandom component of z . 20robability law Expansion polynomials SupportGaussian Hermite ( −∞ , + ∞ )Uniform Legendre [ a, b ]Beta Jacobi [ a, b ]Gamma Laguerre [0 , + ∞ )Poisson Charlier N Table 1: The different polynomial expansions connected to the probability distribution of therandom component z k , k = 1 , . . . , d z .Assuming now s ( z , a, t ), i ( z , a, t ) and r ( z , a, t ) in L (Ω) we may approximate these termsthrough a generalized polynomial chaos expansion in the random space as follows s ( z , a, t ) ≈ s M ( z , a, t ) = M (cid:88) | k | =0 ˆ s k ( a, t )Φ k ( z ) i ( z , a, t ) ≈ i M ( z , a, t ) = M (cid:88) | k | =0 ˆ i k ( a, t )Φ k ( z ) r ( z , a, t ) ≈ r M ( z , a, t ) = M (cid:88) | k | =0 ˆ r k ( a, t )Φ k ( z ) , (34)where the quantities ˆ s k , ˆ i k , ˆ r k are projections in the polynomial spaceˆ s k = E [ s ( · , a, t )Φ k ( · )] , ˆ i k = E [ i ( · , a, t )Φ k ( · )] , ˆ r k = E [ r ( · , a, t )Φ k ( · )] . The sG formulation of system (25) is obtained first by replacing the solutions with theirstochastic polynomial expansions ddt s M ( z , a, t ) = − s M ( z , a, t ) (cid:90) Λ (cid:0) β ( z , a, a ∗ ) − u M ( a, a ∗ ) (cid:1) i M ( z , a ∗ , t ) da ∗ ddt i M ( z , a, t ) = s M ( z , a, t ) (cid:90) Λ (cid:0) β ( z , a, a ∗ ) − u M ( a, a ∗ ) (cid:1) i M ( z , a ∗ , t ) da ∗ − γ ( z , a ) i M ( z , a, t ) ddt r M ( z , a, t ) = γ ( z , a ) i M ( z , a, t ) , (35)with u M ( a, a ∗ ) = 1 κ ( a, a ∗ ) R (cid:2) s M ( z , a, t ) i M ( z , a ∗ , t ) ψ (cid:48) ( I M ( z , t )) (cid:3) , (36) I M ( z , t ) = M (cid:88) | k | =0 ˆ I k ( t )Φ k ( z ) , ˆ I k ( t ) = (cid:90) Λ i M ( z , a, t ) da, and where s M , i M , r M are defined by (34). Then, thanks to the orthonormality of the polynomialbasis of P M , multiplying by Φ m , for all | m | ≤ M , and taking the expectation with respect to p ( z )we obtain the following coupled system of M + 1 deterministic equations for the evolution of each21rojection ddt ˆ s k ( a, t ) = − M (cid:88) | m | =0 B km ( a, t )ˆ s m ( a, t ) ddt ˆ i k ( a, t ) = M (cid:88) | m | =0 B km ( a, t )ˆ s m ( a, t ) − M (cid:88) | m | =0 G km ˆ i m ( a, t ) ddt ˆ r k ( a, t ) = M (cid:88) | m | =0 G km ˆ i m ( a, t ) (37)where B km = M (cid:88) | l | =0 (cid:90) R dz (cid:18)(cid:90) Λ (cid:0) β ( z , a, a ∗ ) − u M ( a, a ∗ ) (cid:1) ˆ i l ( a ∗ , t ) da ∗ (cid:19) Φ k ( z )Φ m ( z )Φ l ( z ) p ( z ) d zG km = (cid:90) R dz γ ( z , a )Φ k ( z )Φ m ( z ) p ( z ) d z . (38)The above system is then integrated in time directly in the space of projections. We remark thatstatistical quantities of interest, such as expectation and variance of infected, can be recovered as E [ i M ( · , a, t )] = (cid:90) R dz i M ( z , a, t ) p ( z ) d z = M (cid:88) | k | =0 ˆ i k ( a, t ) E [Φ k ( · )] = ˆ i ( a, t ) , whereas for the variance we getVar( i M ( · , a, t )) = E [ i M ( · , a, t ) ] − E [ i M ( · , a, t )] = M (cid:88) | k | =0 M (cid:88) | h | =0 ˆ i k ˆ i h E [Φ k ( · )Φ h ( · )] − ˆ i ( a, t ) = M (cid:88) | k | =0 ˆ i k ( a, t ) − ˆ i ( a, t ) . B Social mixing functions
In this appendix we report the details of the social interaction functions that characterize thedynamics of social mixing. These characteristics are in fact crucial for a correct prediction ofoutcomes, especially in diseases transmitted by close contacts. Several large-scale studies havebeen designed in the last decade to determine relevant age-based models in social mixing. Withoutattempting to review the vast literature on this topic, we mention [6, 33, 35] and the referencestherein.The number of contacts per person generally shows considerable variability according to age,occupation, country and even day of the week, in relation to the social habits of the population.Nevertheless, some universal behaviours can be extracted which emerge as a function of specificsocial activities. Social mixing is highly age-related, which means that people usually tend tointeract with other people of a similar age. Young people have a high rate of contact with adultsaged around 30-39 and older people over 60, i.e. their parents and grandparents respectively.Contact rates are indeed very high at home and at school. On the other hand, professional mixingis weakly assortative by age and tends to be determined by uniform interactions, approximatelybetween people from 25 and 60 years old.For these reasons we consider an interaction function determined by three main sub-functionsthat characterize the family, the school and the professional mixing. Therefore, a stylized functionapproximating a realistic contact matrix can be written as follows β ( a, a ∗ ) = (cid:88) j ∈A β j ( a, a ∗ ) , Figure 14: From left to right, contour plot of the social interaction functions β F , β E and β P taking into account the different contact rates of people with ages a and a ∗ in relation to specificsocial activities. The function β F characterizes the family contacts, β E the education and schoolcontacts, and β P the professional contacts.where the functions β j ( a, a ∗ ) take into account the different contact rates of people with ages a and a ∗ in relation to specific social activities of the type A = { F, E, P } , where we identifyfamily contacts with F , education and school contacts with E and professional contacts with P .The particular structure of these social interaction matrices was determined empirically in [6, 35].Here, according to these observations, we propose suitable mathematical functions that can becalibrated to reproduce empirical observations.In details, familiar contacts tend to concentrate on a three-band matrix with a peak aroundyounger ages. This can be reproduced considering the function˜ β F ( a, a ∗ ) = β + λ F, λ F, + ( a + a ∗ ) · e − λ F, ( a + a ∗ ) (1 + ( a − a ∗ ) ) λ F, , λ F, , λ F, , λ F, , λ F, > . Hence, we define the family interactions as β F ( a, a ∗ ) = C F (cid:32) ˜ β F ( a, a ∗ ) + α (cid:88) µ ,µ = ± µ ˜ β F ( a + µ , a ∗ + µ ) (cid:33) , being µ > β E ( a, a ∗ ) = C E (cid:18) β + exp (cid:26) − σ E (cid:0) ( a − λ E ) + ( a ∗ − λ E ) (cid:1)(cid:27)(cid:19) β P ( a, a ∗ ) = C P (cid:18) β + exp (cid:26) σ P (cid:0) ( a − λ P ) + ( a ∗ − λ P ) (cid:1)(cid:27)(cid:19) , being C F , C E , C P > (cid:82) Λ β ( a, a ∗ ) da da ∗ = 1 , λ E > λ P > β ( a, a ∗ ). The details of the parameters used in the simulations are reported in Table 2. References [1] G. Albi, L. Pareschi. Selective model-predictive control for flocking systems.
Commun. Appl.Ind. Math. , (2): 4–21, 2018.[2] G. Albi, L. Pareschi, M. Zanella. Uncertainty quantification in control problems for flockingmodels. Math. Probl. Eng. , : Art. ID 850124, 2015.23igure 15: The social interaction function β = β F + β E + β P .Contact function Parameters β F ( a, a ∗ ) β λ F, λ F, λ F, λ F, µ α / /
100 8 100 3 /
10 1 / β E ( a, a ∗ ) β σ E λ E / /
20 21 / β P ( a, a ∗ ) β σ P λ P / /
40 2 / Commun. Math. Sci. , (6): 1407–1429, 2015.[4] G. Albi, L. Pareschi, M. Zanella. Boltzmann-type control of opinion consensus through lead-ers. Philos. Trans. R. Soc. Lond. Ser. A: Math. Phys. Eng. Sci. , (2028): 20140138, 2014.[5] M. Barro, A. Guiro, D. Ouedraogo. Optimal control of a SIR epidemic model with generalincidence function and a time delays. Cubo , (2): 53–66, 2018.[6] G. B´eraud, S. Kazmercziak, P. Beutels, D. Levy-Bruhl, X. Lenne, N. Mielcarek, Y. Yazdan-panah, P.Y. Bo¨elle, N. Hens, B. Dervaux. The French connection: The first large population-based contact survey in France relevant for the spread of infectious diseases. PLoS ONE , (7): e0133203, 2015.[7] L. Bolzoni, E. Bonacini, C. Soresina, M. Groppi, Time-optimal control strategies in SIRepidemic models, Math. Biosci. , : 86–96, 2017.[8] M. Bongini, M. Fornasier, D. Kalise. (Un)conditional consensus emergence under perturbedand decentralized feedback controls. Discrete Contin. Dyn. Syst. , (9): 4071-4094, 2015.[9] A. Capaldi, S. Behrend, B. Berman, J. Smith, J. Wright, A. L. Lloyd. Parameter estimationand uncertainty quantification for an epidemic model. Math. Biosci. Eng. , (3): 553–576,2012. 2410] V. Capasso, G. Serio. A generalization of the Kermack-McKendrick deterministic epidemicmodel. Math. Biosci. , (1): 43–61, 1978.[11] M. A. Capistr´an, J. A. Christen, J. X. Velasco-Hern´andez. Towards uncertainty quantificationand inference in the stochastic SIR epidemic model, Math. Biosci. , (2): 250–259, 2012.[12] M. Caponigro, M. Fornasier, B. Piccoli, E. Tr´elat. Sparse stabilization and optimal control ofthe Cucker-Smale model. Math. Control Relat. Fields , (4): 447-466, 2013.[13] C. Castillo-Chavez, H. W. Hethcote, V. Andreasen, S. A. Levin, W. M. Liu. Epidemiologicalmodels with age structure, proportionate mixing, and cross-immunity. J. Math. Biol. , (3):233–258, 1989.[14] R. M. Colombo, M. Garavello. Optimizing vaccination strategies in an age structured SIRmodel. Math. Bios. Eng. , (2): 1074–1089, 2019.[15] S. Lee, G. Chowell, C. Castillo-Chvez. Optimal control for pandemic influenza: the role oflimited antiviral treatment and isolation. J. Theor. Biol. , (2): 136–150, 2010.[16] G. Chowell. Fitting dynamic models to epidemic outbreaks with quantified uncertainty: Aprimer for parameter uncertainty, identifiability, and forecasts. Infect. Dis. Model. , (3): 379–398, 2017.[17] G. Dimarco, L. Pareschi, M. Zanella. Uncertainty quantification for kinetic models in socio-economic and life sciences. In Uncertainty Quantification for Hyperbolic and Kinetic Equa-tions , Editors S. Jin, and L. Pareschi, SEMA SIMAI Springer Series, vol. 14, pp 151–191,2017.[18] B. D¨uring, L. Pareschi, G. Toscani. Kinetic models for optimal control of wealth inequalities.
Eur. Phys. J. B , , Article No. 265, 2018.[19] S. Flaxman et al. Estimating the number of infections and the impact of non-pharmaceuticalinterventions on COVID-19 in 11 European countries, Report 13. Imperial College COVID-19Response Team , 30 March 2020.[20] A. Franceschetti, A. Pugliese. Threshold behaviour of a SIR epidemic model with age structureand immigration.
J. Math. Biol. , (1): 1–27, 2008..[21] J. Glasser, Z. Feng, A. Moylan, S. Del Valle, C. Castillo-Chavez. Mixing in age-structuredpopulation models of infectious diseases. Math. Bios. , (1): 1–7, 2012.[22] H. W. Hethcote, Modeling heterogeneous mixing in infectious disease dynamics, in Modelsfor Infectious Human Diseases.
V. Isham and G. F. H. Medley, eds., Cambridge UniversityPress, Cambridge, UK, (1996): pp. 215–238.[23] H. W. Hethcote, The mathematics of infectious diseases.
SIAM Rev. , (4): 599–653, 2000.[24] M. Iannelli, F. A. Milner, A. Pugliese. Analytical and numerical results for the age-structuredS-I-S epidemic model with mixed inter-intracohort transmission. SIAM J. Math. Anal. , (3)662–688, 1992.[25] K. Jagodnik, F. Ray, F. M. Giorgi, A. Lachmann. Correcting under-reportedCOVID-19 case numbers: estimating the true scale of the pandemic. Preprint medRvix:2020.03.14.20036178 .[26] S. Jin, L. Pareschi. Uncertainty Quantification for Hyperbolic and Kinetic Equations , SEMA-SIMAI Springer Series, 2017.[27] W.O. Kermack, A.G. McKendrick. A Contribution to the Mathematical Theory of Epidemics.
Proc. Roy. Soc. Lond. A , : 700–721, 1927.2528] S. Lee, M. Golinski, G. Chowell. Modeling optimal age-specific vaccination strategies againstpandemic influenza. Bull. Math. Biol. , (4): 958-980, 2012.[29] F. Lin, K. Muthuraman, M. Lawley. An optimal control theory approach to non-pharmaceutical interventions. BMC Infect. Dis. , (1), p. 32, 2010.[30] Y. Liu, A. A. Gayle, A. Wilder-Smith, J. Rockl¨ov. The reproductive number of COVID-19 ishigher compared to SARS coronavirus. J. Travel Med. , (2): 1–4, 2020.[31] K. Mizumoto, K. Kagaya, A. Zarebski, G. Chowell. Estimating the asymptomatic proportionof coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship,Yokohama, Japan, 2020. Euro. Surveill. , (10): 2000180, 2020.[32] D. H. Morris, F. W. Rossine, J. B. Plotkin, S. A. Levin. Optimal, near-optimal, and robustepidemic control. Preprint arXiv:2004.02209, 2020.[33] J. Mossong, N. Hens, M. Jit, P. Beutels, K. Auranen, R. Mikolajczyk, M. Massati, S. Salmaso,G. Scalia Tomba, J. Wallinga, J. Heijne, M. Sadkowska-Todys, M. Rosinska, W. J. Edmunds.Social contacts and mixing patterns relevant to the spread of infectious diseases. PLoS Med. , (3): e75, 2008.[34] L. Pareschi. An introduction to uncertainty quantification for kinetic equations and relatedproblems. In Trails in Kinetic Theory: Foundational Aspects and Numerical Methods , SEMA-SIMAI Springer Series, Eds. G. Albi, S. Merino-Aceituno, A. Nota, M. Zanella, to appear.[35] K. Prem, A. R. Cook, M. Jit. Projecting social contact matrices in 152 countries using contactsurveys and demographic data.
PLoS ONE , (9): e1005697, 2017.[36] Presidenza del Consiglio dei Ministri, Dipartimento della Protezione Civile. GitHub: COVID-19 Italia - Monitoraggio situazione, https://github.com/pcm-dpc/COVID-19 , 2020.[37] A. Remuzzi, G. Remuzzi. COVID-19 and Italy: what next? Lancet , :1225–1228, 2020.[38] M. G. Roberts. Epidemic models with uncertainty in the reproduction. J. Math. Biol. , :1463–1474, 2013.[39] S. Wang , F. Zhong , W. Bao , Y. Li , L. Liu , H. Wang , Y. He. Age-dependent risks ofIncidence and Mortality of COVID- 19 in Hubei Province and Other Parts of China Hongdou. medRxiv preprint doi: https://doi.org/10.1101/2020.02.25.20027672 , 2020.[40] D. Xiu. Numerical Methods for Stochastic Computations: A Spectral Methods Approach , Prin-ceton University Press, 2010.[41] S. Zhang, M. Diao, W. Yu, L. Pei, Z. Lin, D. Chen. Estimation of the reproductive numberof novel coronavirus (COVID-19) and the probable outbreak size on the Diamond Princesscruise ship: A data-driven analysis.
International Journal of Infectious Diseases ,93