On the optimal control of kinetic epidemic models with uncertain social features

It is recognized that social heterogeneities in terms of the contact distribution have a strong influence on the spread of infectious diseases. Nevertheless, few data are available on the group composition of social contacts, and their statistical description does not possess universal patterns and may vary spatially and temporally. It is therefore essential to design robust control strategies, mimicking the effects of non-pharmaceutical interventions, to limit efficiently the number of infected cases. In this work, starting from a recently introduced kinetic model for epidemiological dynamics that takes into account the impact of social contacts of individuals, we consider an uncertain contact formation dynamics leading to slim-tailed as well as fat-tailed distributions of contacts. Hence, we analyse the effects of an optimally robust control strategy of the system of agents. Thanks to classical methods of kinetic theory, we couple uncertainty quantification methods with the introduced mathematical model to assess the effects of social limitations. Finally, using the proposed modelling approach and starting from available data, we show the effectiveness of the proposed selective measures to dampen uncertainties together with the epidemic trends.


Introduction
In recent years extensive research efforts have been devoted to design effective non-pharmaceutical interventions (NPIs) to mitigate the impact of the COVID-19 pandemics [4,7,23,27,32,39].In particular, several works in mathematical epidemiology shed light on the importance of the inner heterogeneity in the social structure of a population, see [5,17,19,50].In this direction, among the main factors shaping the evolution of the epidemic, the contact structure of a population has been deeply studied especially in relation to the age distribution of a population.Special attention was recently paid by the scientific community to the role and the estimate of the distribution of contacts between individuals as also a relevant cause of the potential pathogen transmission [6,9,25].Nevertheless, we have often limited information on the real social features of a population, whose characteristics are structurally uncertain and may frequently change due to exogenous processes that are also influenced by psychological factors, determining different responses in terms of individuals' protective behavior, see e.g.[20,28].
Starting from the above considerations, recent works proposed kinetic-type models to connect the distribution of social contacts with the spreading of a disease in multi-agent systems [15,17,35,49].The result is obtained by integrating a compartmental modeling approach for epidemiological dynamics with a thermalization process determining the formation of social contacts.We highlight how the advantages of kinetic modeling approaches for epidemiological dynamics rely on a clear connection between the scales of the transmission of the infection, linking agent-based dynamics with the macroscopic observable ones.Within this research framework, we mention [13,31] where epidemiological relevant states are characterized by agent-based viral load dynamics.
In this paper, we concentrate on a classical SEIR compartmentalization of the population whose contact distribution is uncertain.In particular, we introduce an interaction scheme describing the evolution in the number of social contacts of individuals.The microscopic model is based on a simple transition operator whose parameters are assumed to be uncertain.At the kinetic level, the aforementioned model is capable to identify a variety of equilibrium distributions, ranging from slim-tailed Gamma-type distributions to power-law-type distributions depending on the introduced uncertainties.In the introduced setting, the analysis of the emerging distribution is essential to define the evolution of the main moments of the system of kinetic equations via a closure approach determining the evolution of macroscopic quantities.In particular, we will consider stationary states that depend on uncertain quantities thus, the derived system of equations embeds an incomplete knowledge on the real distribution of contacts.
Therefore, the definition of effective NPIs, generally based on a generalized reduction of the number of contacts, should take into account the uncertain contact structure of a population.In particular, we aim at giving a deeper understanding of the mitigation effects due to the reduction of social interactions among individuals.To this end, we develop an approach sufficiently robust in terms of the introduced uncertainties.This is done through a combination of a kinetic epidemiological model and a control strategy whose target is to point the population towards a given target number of contacts.The development of control protocols for kinetic and mean-field equations has been deeply investigated in recent years, without pretending to review the huge literature we mention [1-3, 24, 40] and the references therein.In detail, we concentrate on modeling the lockdown policies through a selective optimal control approach.In particular, we show how the form of the implemented control may result in very different mitigation effects, that deeply depend on the heterogeneity in the contact distribution of the population.In the last part, starting from the calibrated model at our disposal, we focus on the numerical study of the proposed approach and we exploit accurate methods for the uncertainty quantification of kinetic equations.
The rest of the paper is organized as follows.In Section 2 we introduce a system of kinetic equations with SEIR compartmentalization combining the dynamics of social contacts with the spread of an infectious disease in a multi-agent system.The main features of the solution of a surrogate Fokker-Planck model are studied in Section 2.3.In Section 3 a control strategy is introduced at the kinetic level and in Section 4 we observe the effects of the control on the corresponding second-order macroscopic model.Finally, in Section 5 we investigate numerically the relationship between the kinetic epidemic model with uncertainties and its macroscopic limit.A second part is dedicated to the interface between the introduced modeling approach and available data.

Kinetic epidemic models with uncertain contact distribution
In this section, we introduce a compartmental model describing the spreading of an infectious disease coupled with a kinetic-type description of the contact evolution of a system of individuals [17,18,35,49].In addition, we will also take into account uncertainties collecting the missing information on the contact distribution.
In more details, we consider a system of agents that can be subdivided into the following relevant epidemiological states [8,14,30]: susceptible (S) agents are the ones that can contract the disease, infectious agents (I) are responsible for the spread of the disease, exposed (E) agents have been in contact with infectious ones but still may or may not become contagious; finally, removed (R) agents cannot spread the disease.
To incorporate the impact of contact distribution in the infectious dynamics, we denote by f J = f J (z, x, t) the distribution of the number of contacts x ∈ R + at time t ≥ 0 of agents in compartment J, where J ∈ C := {S, E, I, R}.The random vector z ∈ I z ⊆ R dz , with d z ∈ N, collects all the uncertainties of the system and we suppose to know its distribution p(z) such that We define the total contact distribution of a society as while the mass fractions of the population in each compartment and their moment of order r > 0 are given by In the following, to simplify notations we will indicate with m J (z, t), J ∈ C, the mean values corresponding to r = 1.Hence, we assume that the introduced compartments in the model could act differently at the level of the social process constituting the contact dynamics. of the functions f J (z, x, t) follows by combining the epidemic process with the contact dynamics.This gives the system where the operators Q J (f J ) characterizes the emergence of the distribution of social contacts in the compartment J ∈ C. The transmission of the infection is governed by the local incidence rate defined as where κ(x, x * ) is a nonnegative contact function measuring the impact of contact rates among different compartments.A leading example for κ(x, x * ) is obtained by choosing with β > 0 and α > 0. In the following, we will stick to the case α = 1 for simplicity so that This choice formalizes an incidence rate that is proportional on the product of the number of contacts of susceptible and infected people.The other epidemiological parameters characterizing the spread of the disease are ζ > 0, the transition rate of exposed individuals to the infected class and γ > 0, the recovery rate.The introduced parameters have been summarized in Table 1.Finally, the relaxation parameter 0 < τ ≪ 1 represents the frequency at which the agents modify their contact distribution in response to the epidemic dynamics.As we will see, we are assuming that the social dynamics is much faster than the epidemic dynamics [50].

Contact formation dynamics
The total number of contacts can be viewed as a result of the superimposition of repeated updates and possible deviations due to aleatoric uncertainty, see [26,37].In particular, similarly to [17,18] we consider the following microscopic scheme where x ′ J − x is the elementary variation of the number of contacts and Φ δ ε defines the transition function with ε > 0. In (5) we introduced a constant µ > 0 linked to the maximum variability of the function and the centered random variable η ε such that η 2 ε = εσ 2 , being ⟨•⟩ the expectation with respect to the introduced random variable.The constant ε > 0 tunes the strength of interactions.We remark that the microscopic model (4) depends on a parametric uncertainty and δ = δ(z), such that δ(z) ∈ [−1, 1], for any z ∈ R dz .The transition function ( 5) is defined such that it is simpler to reach a high number of daily contacts while it is very unlikely to go under a certain threshold.This type of asymmetry is typical of human and biological phenomena as shown e.g. in [17,18,29,33,41,43].In the regime ε ≪ 1 we have Note also that the function Φ δ ε is such that 1] and ε > 0. Clearly, the choice µ < 1 implies that, in absence of randomness, the value x ′ J remains positive if x is positive.It is interesting to observe that Φ δ ε is asymmetric around that value x/m J with respect to different distributions of δ.In particular, Φ δ ε is increasing and convex for any x/m J ≤ 1 if δ > 0 whereas, if δ < 0, the transition function becomes concave in an interval [0, x], x/m J < 1, and then convex.
Once the microscopic process (4) is given, the time evolution of the distribution of the number of social contacts f follows by resorting to kinetic collision-like approaches, see [11,37], that quantify the variation of the density of the contact variable in terms of an interaction operator, for any time t ≥ 0. The time evolution of f is given by the following kinetic equation written in weak form where we indicated with φ : R + → R, φ(x) ∈ C ∞ (R + ) an observable quantity.In the following, we will consider an uncertain interaction kernel expressing a multiagent system in which the frequency of changes in the number of social contacts depends on x through the following law being in particular We observe that the kernel (8) mimics the fact that a priori information on the frequency of interaction of a system of agents is missing, see [34].
Remark 2.1.If we consider φ(x) = 1 in (7) we easily get the conservation of the mass.Furthermore, if Therefore, if we exploit the form of the interaction kernel (8) we have that m J (z, t) is a conserved quantity of (7) if δ is a discrete random variable such that δ(z) ∈ {−1, 1} for all z ∈ R dz .A possible example that we will study in the following is given by δ(z) = 1 − 2z, where z ∼ Bernoulli(p).

Fokker-Planck scaling and steady states
In general, it is difficult to compute analytically the equilibrium state of the kinetic model (7).A possible approach has its roots in the so-called grazing collision limit of the classical Boltzmann equation [11,45].
In this direction, a deeper insight on the steady states can be obtained through a quasi-invariant technique [26,37,42].The goal is to derive a simplified Fokker-Planck model from the introduced Boltzmann-type dynamics.For such surrogate model, the study of asymptotic properties is much easier.The idea is to scale simultaneously interactions and interaction frequency.Hence, the equilibrium in contact distribution is reached faster than the time scale of the epidemic dynamics.In details, assuming φ ∈ C ∞ we may observe that for ε ≪ 1 the difference x ′ J − x, J ∈ C, is small and we can perform a Taylor expansion Plugging the above expansion in the interaction operator Q J (f J )(z, x, t) in (7) and thanks to the scaling (6) we get where we have defined the remainder term Assuming |η 3 ε | < +∞ we can prove that, in the limit ε → 0 + , the remainder vanishes thanks to the smoothness of the function φ proceeding as in [12].Hence, in the quasi-invariant scaling regime, we can show that the solution to model (9) converges to Integrating back by parts (11) we obtain the Fokker-Planck model complemented by no-flux boundary conditions We can observe now that the steady state of equation ( 12) depends on the parametric uncertainty of the model and is given by corresponding to generalized Gamma density with C δ,σ 2 ,µ,m J > 0 normalization constant.In particular, we can observe that in the limit δ → 0 we get where again C (0) which is a Gamma distribution.On the other hand, if δ(z) ≡ −1 from ( 14) we get corresponding to an inverse Gamma distribution.More generally, we may observe that the distribution ( 14) exhibits different behaviors depending on the uncertain parameter δ(z).In particular, for each realization of the random variable δ(z) such that δ < 0 the equilibrium density exhibits fat tails with a polynomial decrease for x → +∞.On the other hand, for each realization of the random variable δ(z) such that δ ≥ 0, the equilibrium density is characterized by slim tails.From the modelling point of view, a fat-tailed distribution of contacts defines a society where a non-negligible portion of agents has a high number of contacts.Therefore, the fact that the parameter δ characterizing the tails of the distributions is uncertain means that we take into account the lack of knowledge on the behaviour of the society.
Remark 2.2.In the present context we have neglected effects related to opinion-type dynamics that may influence the process of contact formation.Recent experimental results have shown that social norm changes are often triggered by opinion alignment phenomena.In particular, the perceived adherence of individuals' social network has a strong impact on the effective support of protective behaviour.Therefore, the individual responses to threat are a core question to set up effective measures in the presence of cases escalation.

Uniqueness of the solution
In this subsection, we prove some properties of the solutions of the Cauchy problem (1) for any z ∈ R dz .Let us first concentrate on the Cauchy problem defined by the Fokker-Planck-type problem (12) with given initial condition f J (z, x, 0) ≥ 0. We may apply the arguments of [10,22] to show the positivity of the solution of ( 12).Proposition 2.3.Let f J be a solution of the Cauchy problem where Proof.Let us consider a positive constant ε > 0. We introduce an increasing approximation of the sign function sign ε (f J )(z, x), z ∈ R dz , x ∈ R + , with J ∈ {S, E, I, R}, and define the approximation . Hence, we write the Fokker-Planck equation in weak form where we consider the smooth function φ = sign ε (f J )(z, x) to obtain where we recall that δ = δ(z).Since the boundary terms sign x=0 vanish in view of the boundary conditions, we have Next we observe that for all z ∈ R dz Therefore we have Hence, integrating by parts the first term of the above equation we obtain that in the limit ε → 0 + such term vanishes and for all z ∈ R dz and for all J ∈ C. Therefore, for all t ≥ 0, if we take another solution g J (x, t) of the Cauchy problem (1) with initial condition g 0 J = g J (z, x, 0), we have Corollary 2.4.Let f J be a solution of the Cauchy problem (18) , for all t ≥ 0 and z ∈ R dz .
Proof.The result follows from a similar proof presented in [10].
Now, we concentrate on the epidemic dynamics proving the positivity of the solution of the SEIR-type compartmental system in absence of the collision operators Q J , J ∈ C, see [21].
with the initial data f J (z, x, 0) ≥ 0 for all x ≥ 0 and z ∈ R dz , and K(f S , f I ) defined as Proof.We proceed by contradiction.Let us suppose that there exists a time instant t 0 > 0 such that there exists a point x 0 > 0 such that and for all x ∈ R + , z ∈ R dz .Then, f E (z, x, t) ≥ 0 for all t ∈ [0, t 0 ) and x ≥ 0. If not, there must be a time t 1 ∈ [0, t 0 ) such that there exists a value x 1 > 0 for which and for all x ∈ R + , z ∈ R dz .Hence, integrating the third equation of ( 25) we get Then we have that is not coherent with the hypothesis.As a consequence, it holds f E (z, x, t) ≥ 0 for all t ∈ [0, t 0 ), all x ≥ 0 and all z ∈ R dz .Furthermore, we also have that f I (z, x, t) ≥ 0 for all t ∈ [0, t 0 ) and x ≥ 0. If not, there should be a time t 2 ∈ [0, t 0 ) such that there exists a position x 2 > 0 for which for all x ≥ 0 and z ∈ R dz .Proceeding as before we get that is not coherent with the hypothesis.It follows that f I (z, x, t) ≥ 0 for all t ∈ [0, t 0 ) and x ≥ 0. In view of the results on f E and f I , we get f R (x, t) ≥ 0 for all t ∈ [0, t 0 ) and x ≥ 0.
To conclude, we observe that which is the desired contradiction.Therefore, f S (z, x, t) ≥ 0 for all t ∈ [0, t 0 ) and x ≥ 0.
Once proved the positivity of the contact formation model and of the epidemiological dynamics, we can conclude that the solution of the general Cauchy problem (1) with a non-negative initial data f J is positive a.e. for all t ≥ 0 and z ∈ R dz .
In the following we concentrate on the uniqueness of the solution of the introduced model.
Theorem 2.6 (Uniqueness of the solution).Let f J , g J , with J ∈ {S, E, I, R}, be two solutions of the Cauchy problem where we take K(f S , f I ) as in Proposition 2.5 and constant, positive epidemiological parameters β, ζ, γ > 0. Furthermore, we assume the existence of a positive constant κ > 0 such that ∥κ(x, Proof.In the following, we drop the dependence on x ∈ R + , t ≥ 0 and z ∈ R dz for brevity.We first observe that the difference between two solutions is itself solution of the system From the proof of Proposition 2.3 we get Now, we can rewrite K(f S , f I ) − K(g S , g I ) as follows from which we have with c > 0. This allows us to write Then there exists which, by Gronwall's inequality, gives the claim.

Selective control of the kinetic epidemic model
In Section 2 we introduced and discussed a variety of kinetic models to describe the contact formation dynamics in a society.The main brick of the construction relies on the choice of the transition functions (5) embedding uncertainties in the elementary updates (4), and characterizing the growth in terms of an uncertain parameter δ = δ(z).In particular, it was shown that, for negative values of the parameter δ, the resulting equilibrium contact distribution is given by a distribution with polynomial tails (17).On the other hand, slim tailed distributions can be obtained for positive values of δ, see ( 15)- (16).
In this section, we will investigate the possibility to control the dynamics of contact formation mimicking the action of non-pharmaceutical interventions which should then mitigate the risk factors linked to the transmission of the infection.The new kinetic description allows to enlighten the effects of interventions of the policy maker by acting on the contact distribution of the society of which partial information is available.It is worth to mention that the control of multiagent systems has been recently investigated as a natural follow-up issue in the description and modeling of their self-organization ability, see e.g.[1-3, 24, 33] and the references therein.

The controlled model
To mimic the action of non-pharmaceutical interventions, we add to the microscopic evolution of the social contacts a second update dynamics, implementing an additive control term u, to limit selectively the social activities, see [2,41].Hence, the contact formation is influenced by the uncertain dynamics defined in (4) and, in parallel, by the elementary interaction under control where x ′′ J − x is variation of social contacts in the presence of the control u and S(x) ≥ 0 is a selective function which depends on the number of contacts.
The small parameters τ and ε represent, respectively, the speed at which the contact dynamics equilibrium is reached and the limit from the Boltzmann dynamics to the Fokker-Planck one [18].Two different speed values need to be considered since it is reasonable to assume that such interventions share the time scale with the epidemics, wich is much faster than the contact formation process.We remark also that this second interaction scheme is independent by the uncertain parameter z ∈ R dz and it is linked to a new additive Boltzmann collisional operator which scales with the epidemic dynamics.
The optimal control u * is such that under the constraint (28), where U is the set of the admissible controls, i.e., the set of controls such that x ′′ J ≥ 0. We define the cost J J as follows being κ > 0 a penalization coefficient and x T,J > 0 the desired target number of contacts to reach in each compartment.We remark that the introduced penalization can depend by the compartment of the agent and that the control obtained from (29) subject to (28) is independent on z ∈ R dz .Typical choices for the cost function J J are obtained for p = 1, 2 and a clear analytical understanding is generally difficult for general convex functions and suitable numerical method should be developed.
Let us consider the simple case p = 2. Hence, the minimization of ( 29) can be achieved via a Lagrangian multiplier approach.We define the Lagrangian where θ ∈ R is the multiplier associated to the constraint (28).Then we compute which yields the optimal control Thus, plugging u * into (29) we obtain the controlled update which is a non negative quantity, as required.The kinetic equation expressing the introduced control strategy in the presence of elementary transitions of the contact formation dynamics is a sum of collision operators where the first term on the rhs has been defined in (7) and the second operator describes the impact of non-pharmaceutical interventions on the formation of social contacts.In (32) we have introduced also a second kernel B(z, x), in principle different from B(z, x), describing the frequency of interactions of the agents under the action of the control.
Similarly to what we have done in the uncontrolled scenario, under the grazing limit ε → 0 and scaling the penalization as κ = τ ν, ν > 0, we get a surrogate Fokker-Planck model accounting for an additional drift term quantifying the impact of the control where see [18], whose steady state is given by corresponding to a generalized Gamma density.
Remark 3.1.It is interesting to observe that if B ≡ 1 we can easily determine a S(x) to force a slim tailed equilibrium even in the case δ < 0 for any z ∈ R dz .In particular, we have that any selection function S(x) with superlogarithmic growth is sufficient to ensure that f ∞ J (z, x) is slim-tailed.

Damping effects on the model uncertainties
It is of interest to quantify the effects of the introduced controls on the uncertainties of the kinetic model.
Under suitable hypothesis, it has been observed how the lack of information of system of agents can be dampened for small penalizations, see e.g.[33,44].In the following, we concentrate on the damping effects of the control in terms of the introduced uncertainties by choosing a Maxwellian kernel for the control operator, i.e.B(z, x) ≡ 1, and considering two possible selective functions.We consider the uniform control case S(x) ≡ 1 and the possible selective control that is increasing with x ∈ R + , S(x) = √ x.Let us consider the model (32) and we introduce the time scale ξ = εt.We restrict our analysis to the case in which δ(z) is a discrete random variable such that δ(z) ∈ {−1, 1}.We recall that the mean is conserved in time as observed in Remark 2.1.

By indicating m
Hence, by considering the scaled penalization κ = ντ we get in the limit ε → 0 whose large time behavior is We have since |Φ δ | ≤ µ.In (36) we used the notation m ∞ r,J (z) to indicate the moment of order r > 0 of compartment J ∈ C at the equilibrium, i.e.
Let us consider two cases: • If we consider S(x) ≡ 1, then we get where the quantity m ∞ 1−α(δ),J (z) is finite under the assumption δ(z) ∈ {−1, 1}.Therefore, from bound (37), we have that a vanishing penalization ν leads to a relaxation of the mean to the target x T,J .
Therefore, looking at the variance with respect to the uncertainties z ∈ R dz , we have for all • If we consider now S(x) = √ x, from Jensen's inequality we have Therefore, we obtain the estimate which again, for vanishing penalization ν, implies that the mean reaches the target.Considering the variance with respect to the random variables z ∈ R dz , we obtain Hence, we argue that the introduced controls are capable of damping the variability due to the presence of uncertainties in the distribution of social contacts.Furthermore in the case of zero diffusion case σ 2 = 0 we have In the limit ε → 0 and t → +∞ and with the scaled penalization κ = ντ we obtain .
• Considering S(x) ≡ 1, we get observing that in the limit ν → 0 we have m ∞ (z) → x T,J .
where again, in the limit ν → 0 + , we have m ∞ (z) → x T,J .
Therefore, we can observe that the introduced controls push the energy E ∞ (z) towards the square of the mean number of contacts m ∞ (z).In other words, the steady state converges to a Dirac delta distribution centered at x = x T,J .

Controlled kinetic epidemic model
Once defined the control of the social dynamics, we can define a new kinetic epidemic model embedding the presence of non-pharmaceutical interventions.Following the discussions of Section 3.1, we combine the epidemic process with the controlled contact dynamics as As discussed in Section 2, the transmission of the infection is governed by the local incidence rate K(f S , f I ) defined in (2), the thermalization of the distribution of social contacts in each compartment is given by Q J (f J ) together with the operators C J (f J ) defined in (33).
It is interesting to observe how, under the introduced scaling, the definition of non-pharmaceutical interventions acts at the same time scale of the epidemic dynamics.Hence, the equilibrium states of the dynamics of social contacts result unaltered by the introduction of the control.This fact will be essential in the subsequent section to derive second order macroscopic models describing the evolution of the conserved moments of (39).

Observable effects of non-pharmaceutical interventions
Epidemiological data are typically macroscopic quantities characterizing the evolution of a subset of the introduced compartments.In the following, we derive a macroscopic model which is consistent with the introduced kinetic epidemic model.
We recall here that in [17,18,49] one of the underlying assumptions was that the contact distribution of the population could be fruitfully estimated as an experimentally consistent Gamma distribution [6].In this work, we put uncertainty precisely on the nature of the tail of the contact distribution, which in principle changes the characteristic of the related macroscopic system, thus changing also the efficacy of the containment strategies.

Derivation of the macroscopic model
Recalling that the operators Q J and C J , coupled with no-flux boundary conditions, are mass-preserving, let us integrate system (39) with respect to x to obtain under the assumption on the local incidence rate (3).In (40) we obtained a system for the evolution of the mass fractions.However, we can observe that the system is not closed like in the ones in the classical compartmental framework, since the evolution of ρ J (z, t) depends on the evolution of the first order moment of the distribution functions f J (z, x, t).The evolution of the momentum reads where from (34) we get The hierarchical coupling of moments is a well-known problem in kinetic theory.The closure can, however, be obtained formally by resorting to a limit procedure.Indeed, assuming that the time scale involved in the process of contact formation is τ ≪ 1, we obtain a fast thermalization of the contact distribution of agents with respect to the evolution of the epidemics.Therefore, for τ ≪ 1 the distribution function f J (z, x, t) reaches fast the steady state equilibrium, which is a generalized Gamma distribution with mass fractions ρ ∞ J and local mean values m ∞ J .As observed in Remark 2.1, the case in which δ(z) is a discrete random variable such that δ(z) ∈ {−1, 1} is particularly interesting in the present modeling approach since the mean is conserved.In the following, we stick to this choice and we assume also such that Under this assumption, we can express the second order moment of the generalized Gamma distributions in terms of the mean , where we recall that we fixed λ = µ/σ 2 .Therefore, at the macroscopic level, we obtain the following system of equations for the time evolution of the first order moments in each compartment In (42) the terms G J (f J ), J ∈ C, embed the action of the control at the level of the mean number of social contacts and read We observe now that ( 40) and ( 42) describe in closed form the time evolution of an epidemic where the transition between compartments depend on the mean number of social contacts in the population.
In particular, in the cases S 2 (x) ≡ 1 and S 2 (x) = x, we have For small penalization of the control ν → 0 + , the mean number of connections stabilizes towards the values Therefore, a selective strategy may outperform the uniform one depending on the value of Λ δ (z).We observe that, for vanishing penalizations, the expected number of connections in the compartment , indeed exploiting the information in (41) we get

Numerical examples
In this section, we present several numerical results.We first construct an implicit structure preserving (SP) method [38,48] with a stochastic-Galerkin approach [16,46,51] for system (39).This kind of methods are spectrally accurate in the space of the random parameters under suitable regularity assumptions.For a survey on available methods for the uncertainty quantification of kinetic models we mention [36] and the references therein.In particular, we study the influences of the uncertainties in the spreading of an epidemics and the capability of the designed control strategies in reducing both the peak of the epidemics and the variability of the results given by the random parameters.Furthermore, we consider the macroscopic system of ODEs ( 40)-( 42) and we estimate relevant parameters characterizing non-pharmaceutical interventions based on real epidemiological data.We first estimate the relevant epidemiological parameters thanks to the dataset of the John Hopkins University1 .Hence, we evaluate the impact of different control strategies during the first wave of infection in Italy.

Stochastic Galerkin methods
In order to solve numerically system (39), let us rewrite it in vector form where f = {f J } J , Q = {Q J } J , C = {C J } J , J = {S, E, I, R}, and P is the vector whose components are the transitions rates between the compartments.).Initial conditions given by (49).
Stochastic Galerkin (sG) methods are based on the approximation of the solution f (z, x, t) on a set of polynomials {Ψ h (z)} M h=0 of degree less or equal to M ∈ N, orthonormal with respect to the distribution of the random parameters, such that The polynomials are chosen following the so-called Wiener-Askey scheme [46,47].In the previous relation, we denote by fh (x, t) = { fh,J (x, t)} J the projections of the solution along the linear space generated by the polynomial of degree h where we denote by Ω ⊆ R dz the space of the random parameters.
We discretize the time domain [0, T ] with a time step of size ∆t > 0 and we denote by f n (x) an approximation of f (x, t n ) with t n = n∆t.The first order time splitting method reads: Contact & control dynamics: Epidemic exchange: We plug f M into (45)-( 46) and we project against Ψ h (z)p(z)dz on Ω for each h = 0, . . ., M .Hence, we end with two systems of M + 1 vector equations for the coefficients of the expansion.
The sG reformulation of the contact embedding the control dynamics reads We discretize (47) with a central finite differences approach and we apply a fully-implicit-in-time scheme following the construction presented in [16,48].The epidemic exchange system is where System ( 48) is then integrated through a first order Euler method.
To show the spectral convergence property of the designed sG method, we consider a contact dynamics in the uncontrolled scenario, i.e. with S(x) = 0, of a generic compartment J, that is, we take a single component of (47).We compute a reference solution with N = 25001 grid points of size ∆x = 0.02 in the x-domain [0, 500], ∆t = 0.1, τ = 10 −5 and sG expansion up to order M = 40.We fix the parameters as µ = 0.5, σ 2 = 0.1, being λ = µ/σ 2 , and we consider a one-dimensional uncertainty in a way that δ(z) = z with z ∼ U([−1, 1]).Since the distribution of z is uniform, we consider Legendre polynomials.The initial distribution is a deterministic Gamma with m J = 10.Then, we compute the L 2 error on the first order moment of the distribution at fixed time T = 1 for increasing M .
In Figure 1, we may observe the decay of the numerical error in the space of the random parameter as the order of accuracy increases.We observe that we reach essentially the machine precision within a finite order M .

Test 1: Uncontrolled model
In this section, we focus on the uncontrolled scenario, i.e., system (44) with S(x) = 0. We fix the parameters as β = 0.0025, γ = 0.1, ζ = 0.3, µ = 0.5, σ 2 = 0.1, with λ = µ/σ 2 , we consider a onedimensional uncertainty in a way that δ(z) = z and we investigate the behavior of the model for a uniform random variable z with different support.The x-domain is [0, 500], discretized with N = 25001  grid points of size ∆x = 0.02, the time domain [0, 150] is discretized with the time step ∆t = 0.1; the scale parameter is τ = 10 −5 .We fix the sG expansion up to order M = 5 in all the simulations.The initial conditions for the f 0 J (x) are deterministic Gamma distributions with ρ 0 S = 0.97, ρ 0 E = ρ 0 I = ρ 0 R = 0.01 and m 0 J = 10 for every compartment J.In Figure 2 we show the time evolution of the masses of the compartments for different choices of the random parameter, namely: We observe that the choice c) is associated to a contact equilibrium with fat tails, indicating that there exists a higher probability that agents possess a great number of contacts.Indeed we observe that this choice generates at the equilibrium the smallest number of Susceptible and the highest number of Removed with respect to the other ones, indicating that the epidemics has spread more.Moreover, note also how the peaks of the Infected and Exposed are above the others.The choice b), associated to contact equilibrium with both fat and slim tails, exhibits an intermediate behavior with respect to a), which is associated to slim tails, and c), as expected.

Test 2: Consistency of the macroscopic limit
We consider the coupled system ( 40)- (42) with the underlying assumption that the random variable z follows a Bernoulli distribution of parameter p as in (41).In the following we will fix p = 1/2.We numerically check the consistency of the derived macroscopic closure of the kinetic model which leads to the system ( 40)-( 42) in the limit τ → 0 + .We solve the coupled ODEs with a fourth-order Runge-Kutta method with ∆t = 0.05, the kinetic system ( 44) is solved with the same discretization described in Section 5.2, with the initial conditions (50).The epidemiological parameters are summarized in Table 2.In Figure 3 we observe that smaller values of the time scale τ corresponds to better time-by-time accordance between the kinetic equations and the macroscopic model.

Test 3: Controlled model and uncertainty damping
Let us consider now the controlled model.We concentrate first on the contact dynamics of a single generic compartment J without epidemic exchange, i.e. a component of (47).We are indeed interested in evaluating the effectiveness of the designed control in reducing the tails of the distributions and damping the uncertainties of the system.The parameters, the space and time discretization and the initial conditions are chosen as in Section 5.1.We fix M = 5 and x T,J = 5.
We choose two different selective functions, in the first case we assume S(x) ≡ 1, corresponding to a control that is uniform over the population being independent from the number of contacts.We consider then the selective case with S(x) = √ x, the resulting control has a stronger impact on agents with a higher numbers of contacts.To quantify the effectiveness in reducing uncertainty of the adopted control strategy, we define an index that measures the distance from the target x T and the variability at a given time On the top row of Figure 4, we show the expectation of G ν (z) versus the penalization coefficient ν, for the chosen selective functions.We observe that the control S(x) = √ x is more efficient than the uniform selection, in the sense that reduces more both the variability and the distance from the target for a fixed penalization ν, as discussed in Section 3.2.
On the bottom row of Figure 4, we display in semilogarithmic scale the expectation of the uncontrolled distribution f ∞ (z, x) at the equilibrium, together with the expectations of the numerical solution of (47) at the fixed time T f = 1, for penalizations ν = 1, 10.Note how the introduced control is capable to change the behavior of the tails of the distribution.
Then, we consider the full model ( 44) with epidemic exchange.In particular, we are interested in understanding whether the control on the contact dynamics is able to reduce the spreading of the epidemics and the variability due to the uncertain parameter.To this end, we consider the computational setting of Section 5.2 with z ∼ U([−0.5, 0.5]), T f = 150 and the selective functions S(x) = 1, √ x.In Figure 5 we compare the time evolution of the expectations of the masses ρ J (z, t) in the uncontrolled scenario (black) and under the action of the control with ν = 10 3 , 10 2 (blue and red), for all the compartments.We observe that the control is able to increase the fraction of Susceptible (first row) at the equilibrium and to reduce the Removed (fourth row), but also to dampen the peaks of Exposed (second row) and Infected (third row), meaning that the epidemics has spread less.As expected, with a fixed penalization,  In both cases, we fix p = 1/2 and the epidemiological parameters as in Table 2. Kinetic equations solved with ∆x = 0.02 in the interval [0, 500], and ∆t = 0.1 with T = 20.Initial distribution as in (50).the selective control S(x) = √ x is more efficient than the uniform one, and it is also capable of reducing the uncertainties on the results, as we can notice from the right column, red lines, of Figure 5 5

.5 Test 4: A data-oriented approach
As remarked at the beginning of the section, and following the approach proposed in [18], we will focus on the first wave of the SARS-CoV-2 epidemic during the first half of 2020, particularly in the case of Italy.There, the first detected case was on January, 30th, while the first containment measures were applied on March, 9th.

Test 4a: Calibration of the model
The first step of the calibration is to estimate the unknown epidemiological parameters in the unconstrained regime, assuming that no restriction on the number of contacts was having place, which translates into having G J (f ∞ J )(z, t) ≡ 0. We fixed the known clinical parameters in agreement with the available literature of the field (see, e.g., [18,27] and references therein).In all subsequent figures, we highlighted the evolution of system ( 40)-( 42) obtained in the deterministic cases δ ≡ −1 or δ ≡ 1.Also, we choose the case p = 1/2, to show the performance in an intermediate case.
As done in [18], we solved a least square problem to minimize the relative L 2 norm of the difference between the reported number of infected ρ I and recovered ρ R , and the theoretical evolution of the model ρ I (t) and ρ R (t), with t varying in the timespan [t 0 , t L ] preceding the lockdown regime.For what concerns the initial data, we assume that ρ E (t 0 ) = ρ I (t 0 ) = ρ R (t 0 ) = 1, i.e., t 0 marks nearly the start of the epidemics, while for the average initial number of contacts we set m S (t 0 ) = m E (t 0 ) = m R (t 0 ) = 10, in agreement with the experimentally observed mean number of contacts in a Western country before the pandemic [6].In order to take into account illness and quarantine periods for infected individuals, we fixed their mean number of contacts to be m I (t) ≡ 3 throughout their infection, which corresponds to the average number of family contacts.Thus, the constrained minimization problem is the following: where θ ∈ [0, 1], ∥ • ∥ L 2 ([t0,t L ]) is the relative norm over the time horizon [t 0 , t L ], while we constrained β to belong to the interval [0, 0.01] and λ to satisfy 3 < λ ≤ 10.
In Table 2 we report the parameters obtained by solving problem (52) for different choices of the parameter p, where we fixed the norm coefficient θ = 10 −3 , to better observe the trend with respect to the infectious individuals.

Test 4b: Assessment of different restriction strategies
Once all the epidemiological parameters are estimated, we can focus on the constrained regime, i.e., the subsequent lockdown phase.Within the framework of our model, we can interpret the lockdowns enforced during the first wave of the pandemic in Western Europe as a form of control strategy whose associated selection function S(x) is uniform with respect to the number of contacts x.With this perspective, it is Comparison between data relative to reported infected and recovered people (respectively, black crosses and black plus signs) and time evolution of the mass fractions of infectious agents ρI (t) (red solid line) and removed agents ρR(t) (blue solid line), as prescribed by system ( 40)-( 42) with target xT (t) obtained by solving problem (29).We also reported the evolution in time of the mass fraction of the exposed compartment ρE(t) (green solid line and green circles).In all cases both selection functions S(x) ≡ 1 and S(x) = √ x were employed, but we report distinct curves only for the exposed agents for better clarity, since in all cases we obtain nearly superimposable results, also with respect to different choices of p. Epidemiological parameters as in Table 2.
interesting to compute the optimal target value x T in the control term which permits to fit the data.As a simplifying assumption, we assume this value x T equal for each compartment and we study the two cases of uniform and selective restrictions, which can be obtained by fixing in the dynamics S(x) = 1 and S(x) = √ x, respectively, while we can compute G J (f ∞ J )(z, t) by equation (43).Hence, we solve an optimization problem in the lockdown timespan [t L + 1, t f ], for a sequence of time steps t n over a moving time window of one week (we tried to keep the notation consistent with the one in [18]).Again, it is a constrained least-square problem: with k L = 3, k r = 4.We report the result of such fitting in Figure 6, along with the associated estimated evolution of the exposed compartment for both selection functions S(x) ≡ 1 and S(x) = √ x.In this case, we report the results only for the value p = 1/2, since the fitting procedure gives almost indistinguishable results with respect to the choice of p ∈ [0, 1].
We observe that the estimated value for the target x T is higher when a selective lockdown is enforced, meaning that employing a non-uniform control strategy would achieve the same effects with respect to the number of infected people while allowing greater sociality, especially for the first part of the restriction period.
We also computed the total cost of such measures as the sum of the functionals J S + J E + J R , where J H is defined for H ∈ C as that is, the functional (30) can be seen as the instantaneous approximation of J H , which is obtained by considering (31) in the limit ε, τ → 0 + (see [18]).We see that the cost is strongly influenced by the

Figure 4 :
Figure 4: Test 3. Top row: expectation of Gν (z) defined in (51) with xT,J = 5, versus the penalization ν, for S(x) = 1 (left) and S(x) = √ x (right).Bottom row: details of the expected values of the distributions: the black line represents the uncontrolled distribution at the equilibrium, the blue and red lines are the controlled distribution for ν = 1, 10, at the fixed time T f = 1, for the selective functions S(x) = 1 (left) and S(x) = √ x (right).In all the simulations we choose ∆x = 0.02 in the interval [0, 500], ∆t = 0.1 with T f = 1 and τ = 10 −5 .The sG expansion is M = 5, the uncertain parameter is δ(z) = z with z ∼ U ([−1, 1]).Initial conditions given by (49).

Figure 6 :
Figure 6: Test 4a.Comparison between data relative to reported infected and recovered people (respectively,

Table 1 :
The kinetic model defining the time evolution Parameters definition in the SEIR model (1).