key: cord-121428-79wyxedn
authors: Dimarco, G.; Perthame, B.; Toscani, G.; Zanella, M.
title: Social contacts and the spread of infectious diseases
date: 2020-09-02
journal: nan
DOI: nan
sha: 
doc_id: 121428
cord_uid: 79wyxedn

Motivated by the COVID-19 pandemic, we introduce a mathematical description of the impact of sociality in the spread of infectious diseases by integrating an epidemiological dynamic with a kinetic modeling of population-based contacts. The kinetic description leads to study the evolution over time of Boltzmann type equations describing the number densities of social contacts of susceptible, infected and recovered individuals, whose proportions are driven by a classical compartmental model in epidemiology. Explicit calculations show that the spread of the disease is closely related to the mean number of contacts, thus justifying the lockdown strategies assumed by governments to prevent them. Furthermore, the kinetic model allows to clarify how a selective control can be assumed to achieve a minimal lockdown strategy by only reducing individuals undergoing a very large number of daily contacts. This, in turns, could permit to maintain at best the economic activities which would seriously suffer from a total closure policy. Numerical simulations confirm the ability of the model to describe different phenomena characteristic of the rapid spread of an epidemic. A last part is dedicated to fit numerical solutions of the proposed model with experimental data coming from different European countries.

1 Introduction a model which well fits with the experimental data.

An easy way to understand epidemiology models is that, given a population composed of agents, they prescribe movements of individuals between different states based on some matching functions or laws of motion. According to the classical SIR models [29] , agents in the system are split into three classes, the susceptible, who can contract the disease, the infected, who have already contracted the disease and can transmit it, and the recovered who are healed, immune or isolated. Inspired by the model considered in [14] for describing a social attitude and making use of the SIR dynamics, we present here a model composed by a system of three kinetic equations, each one describing the time evolution of the distribution of the number of contacts for the subpopulation belonging to a given epidemiological class. These three equations are further coupled by taking into account the movements of agents from one class to the other as a consequence of the infection. The interactions which describe the social contacts of individuals are based on few rules and can be easily adapted to take into account the behavior of agents in the different classes in presence of the infectious disease. Our joint framework is consequently based on two models which can be considered classical in the respective fields. From the epidemic side, the SIR and its related models are widely used in most applications (cf. [29] and the references therein). From the side of multi-agent systems in statistical mechanics, the linear kinetic model introduced in [14, 46] has been shown to be flexible and able to describe, with opportune modifications, different problems in which human behavior plays an essential role, like the formation of social contacts.

The study of the effects of the intensity of social contacts in the epidemic diffusion, among other consequences, allows to obtain in a rigorous way a variety of non-linear incidence rates of the infectious disease as for instance recently considered in [33] . It is also interesting to remark that the presence of non-linearity in the incidence rate function, and in particular, the concavity condition with respect to the number of infected has been considered in [33] as a consequence of psychological effects. Namely, the authors observed that in the presence of a very large number of infected, the probability for an infective to transmit the virus further may decrease because the population tend to naturally reduce the number of contacts. This fact will be embedded in the kinetic model producing a change in the mean value of the number of daily contacts. The importance of reducing at best the social contacts to countering the advance of a pandemic is a well-known and studied phenomenon [17] . While in normal life activity, it is commonly assumed that a large part of agents behaves in a similar way, in presence of an extraordinary situation like the one due to a pandemic, it is highly reasonable to conjecture that the social behavior of individuals is strictly affected by their personal feeling in terms of safeness. In this work, we focus on the assumption that it is the degree of diffusion of the disease that changes people's behavior in terms of social contacts, in view both of the personal perception and/or of the external government intervention. More generally, the model can be clearly extended to consider more realistic dependencies between an epidemic disease and the social contacts of individuals. However, this does not change the essential conclusions of our analysis, namely that there is a close interplay between the spread of the disease and the distribution of contacts, that the kinetic description is able to quantify. The encouraging results described in the rest of the article, finally suggest that a similar analysis can be carried out, at the price of an increasing difficulty in computations, in more complex epidemiological models like the SIDARTHE model [25, 26] , recently used to simulate the COVID-19 epidemic in Italy to validate and improve the eventual partial lockdown strategies of the government and to suggest future measures.

The rest of the paper is organized as follows. Section 2 introduces the system of three SIRtype kinetic equations combining the dynamics of social contacts with the spread of infectious disease in a system of interacting individuals. Then, in Section 3 we show that in a suitable asymptotic procedure the solution to the kinetic system tends towards the solution of a system of three SIR-type Fokker-Planck type equations with local equilibria of Gamma-type [3] . Once the system of Fokker-Planck type equations has been derived, in Section 4 we close the SIR-type system of kinetic equations around the Gamma-type equilibria to obtain a SIR model in which the presence and consequently the evolution of the social contacts leads to a non-linear incidence rate of the infectious disease satisfying the compatibility conditions introduced in [33] . Last, in Section 5, we investigate at a numerical level the relationships between the solutions of the kinetic system of Boltzmann type, its Fokker-Planck asymptotics and the SIR model. These simulations confirm the ability of the model to describe different phenomena characteristic of the trend of social contacts in situations compromised by the rapid spread of an epidemic, and the consequences of various lockdown action in its evolution. A last part is sacred to a fitting of the model with the experimental data by first estimating the parameters of the epidemic through data and successively by using these parameters in the kinetic model.

Our first goal is to build a kinetic system which suitably describes the spreading of an infectious disease under a dependence of the contagiousness parameters on the individual number of social contacts of the agents. As in classical SIR models [29] , the entire population is divided into three classes: susceptible, infected and recovered individuals.

Aiming to understand social contacts effects on the dynamics, we will not consider in the sequel the role of heterogeneity in the disease parameters (such as the personal susceptibility to a given disease), which could be derived from the classical SIR model, suitably adjusted to account for new information [12, 38, 48] . Consequently, agents in the system are considered indistinguishable [40] . This means that the state of a person in each class at any instant of time t ≥ 0 is completely characterized by the number of contacts x ≥ 0, measured in some unit.

While x is a natural positive number at the individual level, without loss of generality we will consider x in the rest of the paper to be a nonnegative real number, x ∈ R + , at the population level. We denote then by f S (x, t), f I (x, t) and f R (x, t), the distributions at time t > 0 of the number of social contacts of the population of susceptible, infected and recovered individuals, respectively. The distribution of contacts of the whole population is then recovered as the sum of the three distributions

As outlined in the Introduction, we do not consider for simplicity of presentation disease related mortality as well as the presence of asymptomatic individuals. Therefore, we can fix the total distribution of social contacts to be a probability density for all times t ≥ 0

As a consequence, the quantities

denote the fractions of susceptible, infected and recovered respectively. For a given constant α > 0, we also denote with m α (t) the moment of order α of the distribution of the number of contacts

and the moments of order α for the distributions of the number of contacts for each class

Unambiguously, we will indicate the mean values, corresponding to α = 1, by m(t) and, respectively,

In what follows, we assume that the evolution of the densities is built according to the classical SIR model [29] , and that the various classes in the model could act differently in the social process constituting the contact dynamics. The kinetic model then follows combining the epidemic process with the contact dynamics, as modeled in [14, 46] . This gives the system

where γ is the constant recovery rate while the transmission of the infection is governed by the function K(f S , f I ), the local incidence rate, expressed by

In (3) the contact function κ(x, y) is a nonnegative function growing with respect to the number of contacts y of the population of infected, such that κ(x, 0) = 0. Then, the spreading of the epidemic depends heavily on the function κ(·, ·) used to quantify the rate of possible contagion in terms of the number of social contacts between the classes of susceptible and infected. The evolution of the mass fractions J(t), J ∈ {S, I, R} then obeys to the classical SIR model by choosing κ(x, y) ≡ β > 0 [29] , thus considering the spreading of the disease independent of the intensity of social contacts. A leading example is obtained by choosing a rank 1 matrix

where α, β are positive constants, that is by taking the incidence rate directly proportional to the product of the number of contacts of susceptible and infected people. When α = 1, for example, the incidence rate takes the simpler form

In equations (2) the operators Q J , J ∈ {S, I, R} characterize the thermalization of the distribution of social contacts in the various classes in terms of repeated interactions [14, 46] . The Q J , J ∈ {S, I, R} are integral operators that modify the distribution of contacts f J (x, t), J ∈ {S, I, R}. Their action on observable quantities is given by

where B(x) measures the interaction frequency, · denotes mathematical expectation with respect to a random quantity, and x * J denote the post-interaction value of the number x of social contacts of the J-th population. Last, the constant τ in front of the interaction integral measures the time scale of the thermalization of the distribution of social contacts. Remark 1. The relaxation constant τ plays an important role in the time evolution of the model (2), since it relates the time scale of the epidemic evolution with that of the statistical formation of social contacts. Small values of the constant τ will correspond to a fast adaptation of people to a steady situation. Hence, since it is reasonable to assume that adaptation of people is faster that the dynamics of epidemic, in what follows we will assume τ 1.

The process of formation of the social contacts distribution is obtained by taking into account the typical aspects of human being, in particular the search, in absence of epidemics, of opportunities for socialization. In addition to that, social contacts are due to the common use of public transportations to reach schools, offices and, in general, places of work and to basic needs of interactions due to work duties. As shown in [3] , this leads individuals to stabilize on a characteristic number of daily contacts depending on the social habits of the different countries, (represented by a suitable valuex M of the variable x, which can be considered as the mean number of contacts relative to the population under investigation).

Then, following a consolidated path in the kinetic theory of social phenomena, one aims in obtaining the characterization of the distribution of social contacts in a multi-agent system (the macroscopic behavior), starting from some universal assumption about the personal behavior of agents (the microscopic behavior). Indeed, as in many other human activities, the daily amount of social contacts is the result of a repeated upgrading based on well-established rules. To this extent, it is enough to recall that the daily life of each person is based on a certain number of activities, and each activity carries a certain number of contacts. Moreover, for any given activity, the number of social contacts varies in consequence of the personal choice. The typical example is the use or not of public transportations to reach the place of work or the social attitudes which scales with the age. Clearly, independently of the personal choices or needs, the number of daily social contacts contains a certain amount of randomness, that has to be taken into account. Also, while it is very easy to reach a high number of social contacts attending crowded places for need or will, going below a given threshold is very difficult, since various contacts are forced by working or school activities, among others causes. This asymmetry between growth and decrease, as exhaustively discussed in [14, 27, 28] , can be suitably modeled by resorting to the value function description of the elementary variation of the x variable.

It is important to outline that, in the presence of an epidemic, the characteristic mean number of daily contactsx M reasonably changes in time, even in absence of an external lockdown intervention, in reason of the perception of danger linked to social contacts that people manifest. Consequently, even if not always explicitly indicated, we will assumex M =x M (t). Also, this time dependent value can be different depending on the class to which agents belong.

A further non secondary aspect of the formation of the number of social contacts is that their frequency is not uniform with respect to the values of x. Indeed, it is reasonable to assume that the frequency of interactions is inversely proportional to the number of contacts x. This relationship takes into account that it is highly probable to have at least some contacts, and the rare situation in which one reaches a very high values of contacts x.

The choice of a variable interaction frequency has been fruitfully applied in a different context in a different situation [22] . This was done to better describe the evolution in time of the wealth distribution in a western society. There, the frequency of the economic transactions has been proportionally related to the values of the wealth involved, to take into account the low interest of trading agents in transactions with small values of the traded wealth. As discussed in [22] , the introduction of a variable kernel into the kinetic equation does not modify the shape of the equilibrium density, but it allows a better physical description of the phenomenon under study, including an exponential rate of relaxation to equilibrium for the underlying Fokker-Planck type equation.

Following [14, 27, 28, 46] , we will now illustrate the mathematical formulation of the previously discussed behavior. In full generality, we assume for the moment that individuals in different classes can have a different mean number of contacts. The microscopic updates of social contacts of individuals in the class J, J ∈ {S, I, R} will be taken of the form

In a single update (interaction), the number x of contacts can be modified for two reasons, expressed by two terms, both proportional to the value x. In the first one, the coefficient Φ (·), which can assume both positive and negative values, characterizes the typical and predictable variation of the social contacts of agents, namely the personal social behavior of agents. The second term takes into account a certain amount of unpredictability in the process. The usual choice is to assume that the random variable η is of zero mean and bounded variance, expressed by η = 0, η 2 = λ, with λ > 0. Small random variations of the interaction (5) will be expressed simply by multiplying η by a small positive constant √ , with 1, which produces the new (small) variance λ. Notice that this formulation makes an essential use of the choice x ∈ R + . This scaled variable is the variable η appearing in (5) . The function Φ plays the role of the value function in the prospect theory of Kahneman and Twersky [30, 31] , and contains the mathematical details of the expected human behavior in the phenomenon under consideration. In particular, the hypothesis on which it is built is that it is normally easier to increase the value of x than to decrease it, in relationship with the mean valuē x J , J ∈ {S, I, R}. In terms of the variable s = x/x J we consider as in [14] the class of value functions given by

where 0 < δ ≤ 1 and 0 < µ < 1 are suitable constants characterizing the intensity of the individual behavior, while the constant > 0 is related to the intensity of the interaction. Hence, the choice 1 corresponds to small variations of the mean difference x * J − x and making it common to both effects, randomness and adaptation, permits to equilibriate their effects as can be seen in Section 3. In (7), the value µ denotes the maximal amount of variation of x that agents will be able to obtain in a single interaction. Note indeed that the value function Φ δ (s) is such that

Clearly, the choice µ < 1 implies that, in absence of randomness, the value of x * J remains positive if x is positive. As proven in [14] , the value function (6) satisfies

These properties are in agreement with the expected behavior of agents, since deviations from the reference point (s = 1 in our case), are bigger below it than above. Letting δ → 0 in (6) allows to recover the value function Φ 0 (s) = µ s − 1 s + 1 , s ≥ 0.

introduced in [27, 28] to describe phenomena characterized by the lognormal distribution [2, 35] . Once the elementary interaction (5) is given, for any choice of the value function, the study of the time-evolution of the distribution of the number x of social contacts follows by resorting to kinetic collision-like models [8, 40] , that quantify at any given time the variation of the density of the contact variable x in terms of the interaction operators. For a given density f J (x, t), J ∈ {S, I, R}, the action on the density of the interaction operators Q J (f )(x, t) in equations (2) is fruitfully written in weak form. The weak form corresponds to say that for all smooth functions ϕ(x) (the observable quantities) the action of the operator Q J on ϕ is given by

Here expectation · takes into account the presence of the random parameter η in the microscopic interaction (5) . The function B(x) measures the interaction frequency. The right-hand side of equation (8) quantifies the variation in density, at a given time t > 0, of individuals in the class J ∈ {S, I, R} that modify their value from x to x * J (loss term with negative sign) and agents that change their value from x * J to x (gain term with positive sign). In [14] , the interaction kernel B(x) has been assumed constant. This simplifying hypothesis, which is not always well justified from a modeling point of view, is the common assumption in the Boltzmann-type description of socio-economic phenomena [21, 40] . In particular, the role of a non constant collision kernel B(x) has been analyzed in its critical aspects in [22] .

Following the approach in [22, 46] , we express the mathematical form of the kernel B(x) by assuming that the frequency of changes which leads to increase the number x of social contacts is inversely proportional to x. Hence, we consider collision kernels in the form

for some constants a > 0 and b > 0. This kernel assigns a low probability to interactions in which individuals have already a large number of contacts and assigns a high probability to interactions in which the value of the variable x is small.

The constants a and b can be suitably chosen by resorting to the following argument [46] . For small values of the x variable, the rate of variation of the value function (6) is given by

Hence, for small values of x, the mean individual rate predicted by the value function is proportional to x δ−1 . Then, the choice b = δ would correspond to a collective rate of variation of the system independent of the parameter δ which instead characterizes the individual rate of variation of the value function.

A second important fact is that the individual rate of variation (9) depends linearly on the positive constant , and it is such that the intensity of the variation decreases as decreases. Then, the choice a = 1 , is such that the collective rate of variation remains bounded even in presence of very small values of the constant . With these assumptions, the weak form of the interaction integrals in (8), is given by

Note that, in consequence of the choice made on the interaction kernel B(·), the evolution of the density f (x, t) is tuned by the parameter , which characterizes both the intensity of interactions and the interaction frequency. This choice implies

Consequently, the actions of both the value function and the random part of the elementary interaction in (5) survive in the limit → 0.

Remark 3. The approach just described can be easily adapted to other compartmental models in epidemiology like SEIR, MSEIR [6, 13, 29] and/or SIDHARTE [25, 26] . For all these cited models, the fundamental aspects of the interaction between social contacts and the spread of the infectious disease we expect not to change in a substantial way.

The limit procedure induced by (11) corresponds to the situation in which elementary interactions (5) which produce an extremely small modification of the number of social contacts (quasi-invariant interactions) are prevalent in the dynamics. At the same time, the frequency of these interactions is suitably increased to still permitting to see their effect.

In kinetic theory, this is a well-known procedure with the name of grazing limit. We point the interested reader to [11, 20, 40] for further details. Expanding the difference ϕ(x * J ) − ϕ(x) in Taylor series, and then passing to the limit into the right-hand side of (10), one obtains, for J ∈ {S, I, R}

If we impose at x = 0 the boundary conditions

and a suitable rapid decay of the densities f J (x, t) at x = +∞ [21] , the limit operators Q δ J , J ∈ {S, I, R}, coincide with the Fokker-Planck type operators

characterized by a variable diffusion coefficient and a time-dependent drift term. Indeed, in view of Remark 2 we havex J =x J (t). Thus, system (2), in the limit → 0 corresponds to the following Fokker-Planck system in strong form

The Fokker-Planck system (14) is complemented with the boundary conditions (12) at x = 0. Clearly, the steady distributions satisfy the ordinary differential equations corresponding to equations (14) with time derivatives set equal to zero. It is interesting to remark that the explicit form of the equilibria is easily obtained when both β and γ are set equal to zero, and the contact function κ(x, y) = 0, so that the incidence rate K(f S , f I ) vanishes. In this simple case, by assuming that the mean valuesx J , J ∈ {S, I, R} are constant, and by setting ν = µ/λ, the equilibria are given by the functions

J ∈ {S, I, R}. By fixing the mass of the steady state (15) equal to one, in agreement with [3] , the consequent probability densities are generalized Gamma f ∞ (x; θ, χ, δ) defined by [34, 45] 

characterized in terms of the shape χ > 0, the scale parameter θ > 0, and the exponent δ > 0 that in the present situation are given by

It has to be remarked that the shape χ is positive, only if the constant ν = µ/λ satisfies the bound

Note that condition (18) holds, independently of δ, when µ ≥ 4λ, namely when the variance of the random variation in (5) is small with respect to the maximal variation of the value function. Note moreover that for all values δ > 0 the moments are expressed in terms of the parameters denoting respectively the meanx J , J ∈ {S, I, R}, the variance λ of the random effects and the values δ and µ characterizing the value function φ δ defined in (6) . The standard Gamma and Weibull distributions are included in (15) , and are obtained by choosing δ = 1 and, respectively δ = χ. In both cases, the shape χ = ν, and no conditions are required for its positivity.

The Fokker-Planck system (14) contains all the information on the spreading of the epidemic in terms of the distribution of social contacts. Indeed, the knowledge of the densities f J (x, t), J ∈ {S, I, R}, allows to evaluate by integrations all moments of interest. However, in reason of the presence of the incidence rate K(f S , f I ), as given by (3), the time evolution of the moments of the distribution functions is not explicitly computable, since the evolution of a moment of a certain order depends on the knowledge of higher order moments, thus producing a hierarchy of equations, like in classical kinetic theory of rarefied gases [8] . We will come back to this question later on.

Before discussing the closure, we establish some choices in the model in order to obtain results which are in agreement with the ones reported in [3] . We fix the value δ = 1, and the incidence rate as in (4) . In particular, the choice δ = 1 gives, for J ∈ {S, I, R}

In this case (17) imply χ = ν and θ =x J /ν, and the steady states of unit mass, for J ∈ {S, I, R}, are the Gamma densities

With this particular choice, the mean value and the variance of the densities (19), J ∈ {S, I, R}, are given by

We now go back to the problem of modeling the time evolution of the densities f J (x, t). We assume that due to the presence of the epidemic, the population tends to reduce the typical average number of contacts which exhibits in standard situations. This can happen due to two main reasons: on voluntary basis for preventing contagion or by authorities decision through lockdown measures. This average reduction effect can be introduced by assuming that the mean number of daily contacts x J (t), J ∈ {S, I} depends on the proportion of infected

where the function H(r) is decreasing with respect to r, 0 ≤ r ≤ 1, starting from H(0) = 1. In particular, in our model, we do not consider a time-dependent mean number of social contacts for the class of recovered. This is due to the fact that the behavior of epidemic does not depend on it. For J ∈ {S, I, R}, andx J (t) as in (20), we now definē

With these notations, system (14) with δ = 1 reads

Integrating both sides of equations in (21) with respect to x, and recalling that the Fokker-Planck type operators are mass-preserving, we obtain the system for the evolution of the proportions J(t) defined in (1), J ∈ {S, I, R}

As anticipated, unlike the classical SIR model, system (22) is not closed, since the evolution of the mass fractions J(t), J ∈ {S, I, R} depends on the mean values m J (t).

Similarly to the derivation of macroscopic equations in kinetic theory, the closure of system (22) can be obtained by resorting, at least formally, on a limit procedure. In fact, as outlined in the Introduction, the typical time scale involved in the social contact dynamic is τ 1 which identifies a faster adaptation of individuals to social contacts with respect to the evolution time of the epidemic disease. The choice of the value τ 1 pushes the distribution function f S (x, t) towards the Gamma equilibrium density with a mass fraction S(t) and momentum

as it can be easily verified from the differential expression of the interaction operatorQ S . Indeed, if τ 1 is sufficiently small, one can easily argue from the exponential convergence of the solution of the following expression

towards the equilibrium f ∞ S (x; θ, ν) (see [47] for details), that the solution f S (x, t) remains sufficiently close to the Gamma density with mass S(t) and momentum given by (23) for all times. This equilibrium distribution f ∞ S (x; θ, ν) can be plugged into the first equation of (21) and then one can integrate it with respect to the variable x. This procedure gives formally the first equation of (22). Analogous analysis and formal limit procedure can be done with the second equation in system (21) , which leads to relaxation towards a Gamma density with mass fraction I(t) and momentum given by m I (t) =x I I(t)H(I(t)).

Consequently, we obtain the closure ∂S(t) ∂t = −β S(t) I(t) H 2 (I(t)),

In system (26) we definedβ = βx IxS . This identifies the classical transmission parameter of the SIR model, where however now the difference is that this quantities are not postulated but instead derived starting from microscopic considerations. In the following, we refer to this model as the social contact-SIR model (S-SIR).

Remark 4. The outlined closure strategy can be obtained also by resorting to the splitting method, a very popular numerical approach for the Boltzmann equation [23, 39] . If at each time step (t, t + ∆t) we consider sequentially the population-based interaction and relaxation operators in the first equation in (14) , during this short time interval we recover the evolution of the density from the joint action of relaxation ∂f S (x, t) ∂t

and SIR interaction

We recall once again that the typical time to consider is τ 1, which identifies a faster adaptation of individuals to social contacts with respect to the evolution time of the epidemic disease. Since the operatorQ S is mass preserving, in the considered time interval the relaxation (27) with a value τ 1 pushes the solution of equation (24) towards the Gamma equilibrium density with the same mass fraction S(t) of the initial datum, and momentum m S (t + ∆t) =x S S(t)H(I(t)), (29) as it can be easily verified from the differential expression of the interaction operatorQ S . Indeed, if τ 1 is sufficiently small, one can easily argue from the exponential convergence of the solution to (24) towards the equilibrium [47] , that the solution f S (x, t + ∆t) is sufficiently close to the Gamma density with mass S(t) and momentum given by (29) , and this density can be used into the SIR step (28) to close the splitted system (24)- (28) . Analogous procedure can be done with the second equation in system (21) , which leads to relaxation towards a Gamma density with mass fraction I(t) and momentum given by

Consequently, substituting into system (22) we obtain the closed system (26) .

It remains to quantify the action of the social contacts on the evolution of the epidemic. For this reason, let us introduce for a given constant N 1 the decreasing function [33] 

which describes a possible way in which, in presence of the spread of the disease, the susceptible and the infected population tend to reduce the mean number of daily social contactsx J , J ∈ {S, I}. This choice produces a SIR model with global incidence rate

that fulfils all the properties required by the non-linear incidence rates considered in [33] . 

for all S, I > 0. In addition to the form (31), we can also consider the following function

This function satisfies the same properties than (31) in terms of the incidence rate requests detailed in [33] and takes into account memory effects on the population's behavior. In fact, people may adapt their life style in terms of possible daily contacts to answer to the actual pandemic situation as To conclude, let us introduce the basic reproduction number R 0 of this model, i.e. the average number of secondary cases produced by a single infected agent introduced into an entirely susceptible population. This is given by

where γ is the recovery rate of infected. According to the analysis of [33] , an autonomous compartmental epidemiological model with the non-linear incidence rate (32) under the constant population size assumption is stable. Such system has either a unique and stable endemic equilibrium state or no endemic equilibrium state at all. Since the incidence rate D(S, I) satisfies the conditions (33)- (34) , if R 0 > 1, then the endemic equilibrium state Q * = (S * , I * ) of system (26) is asymptotically stable. If R 0 ≤ 1, then there is no endemic equilibrium state, and the infection-free equilibrium state is asymptotically stable.

In addition to analytic expressions, numerical experiments allow us to visualize and quantify the effects of social contacts on the SIR dynamics used to describe the time evolution of the epidemic. More precisely, starting from a given equilibrium distribution detailing in a probability setting the daily number of contacts of the population, we show how the coupling between social behaviors and number of infected people may modify the epidemic by slowing down the number of encounters leading to infection. In a second part, we discuss how some external forcing, mimicking political choices and acting on restrictions on the mobility, may additionally improve the reduction of the epidemic trend avoiding concentration in time of people affected by the disease and then decreasing the probability of hospitalization peaks. In a third part, we focus on experimental data for the specific case of COVID-19 in different European countries and we extrapolate from them the main features characterizing the contact function H(·).

The starting point is represented by a population composed of 99% of susceptibles and 0.01% of infected. The distribution of the number of contact is described by (15) with ν = 4, δ = 8 and x J = 15 while the epidemic parameters are k(x) = 0.8

x J x and γ = 0.8. The kinetic model (14) is solved by a splitting strategy together with a Monte Carlo approach where the number of samples used to described the population is fixed to M = 10 6 . The time step is fixed to ∆t = 10 −2 and the scaling parameter is τ = 10 −2 . These choices are enough to observe the convergence of the Boltzmann dynamics to the Fokker-Planck one as shown in Figure 1 where the analytical equilibrium distribution is plotted together with the results of the Boltzmann dynamics. We considered also uniform initial distribution

being χ(·) the indicator function. In the introduced setting, we then compare two distinct cases: in the first one we suppose that the social contacts do not affect the solution, meaning H(I(t)) = 1, while the second includes the effects of the function H(I(t)) given in (31) with N = 10. The results are depicted in Figure 2 . The top right images show the time evolution of the distribution of the number of contacts for the two distinct cases, while the middle images report the corresponding evolution of the epidemic. For this the second case, the function H(I(t)) as well as the distribution of contacts for respectively the susceptible and the recovered are shown at the bottom of the same figure. We clearly observe a reduction of the peak of infected in the case in which the dynamics depends on the number of contacts with H(I(t)) given by (31) . We now repeat the same simulation by only changing the considered epidemic parameters. In particular, we consider a lower infection and recovery rate given by k(x) = 0.2/x J x and γ = 1, respectively. In Figure 3 , we show the evolution of two epidemic profiles in time for the case in which social contacts do not affect the solution and for the case in which H(I(t)) is a function of the number of infected as for the previous situation. Results show the same qualitative behavior: peak reduction and spread of the number of infected over time is observed when sociality is taken into account. The total number of infected is also reduced in the second case.

Next, we compare the effects on the spread of the disease when the population adapts its habits with a time delay with respect to the onset of the epidemic. This kind of dynamics corresponds to a modeling of a possible lockdown strategy whose effects are to reduce the mobility of the population and, correspondingly, to reduce the number of daily contacts in the population.

The setting is similar to the one introduced in Section 5.1 and we consider a switch between H = 1 to H(I(t)) = 1/ 1 + N I(t) when the number of infected increases. The social parameters are ν = 4, δ = 8 and x J = 15, as before, while the epidemic parameters are k(x) = 0.4/x J x and γ = 0.6, the final time is fixed to T = 15. The initial distribution of contacts is also assumed to be of the form (37) .

We consider three different settings, in the first one H = 1 up to t < T /4, in the second one up to t < T /2 while in the third one we prescribe a lockdown for a limited amount of time (1 < t < 2) and then we relax back to H = 1. The results are shown in Figure 4 for both the distribution of daily contacts over time and the SIR evolution. We can identify three scenarios. The first on the top gives a slightly change to the standard epidemic dynamics. Indeed, we can observe around t = 7.5 a reduction of the speed of the infection. For the second, we observe an inversion in the trend of the epidemic around t = 4, while for the third case we first observe inversion and then the resurgence of the number of infected when the lockdown measures are relaxed.

In this part, we consider data about the dynamics of COVID-19 in three European countries: France, Italy and Spain. For these three countries, the evolution of the disease, in terms of reported cases, evolved in rather different ways. The estimation of epidemiological parameters of compartmental models is an inverse problem of generally difficult solution for which different approaches can be considered. We mention in this direction a very recent comparison study [36] . It is also worth to mention that often the data are partial and heterogeneous with respect to their assimilation, see for instance discussions in [1, 7, 9, 44] . This makes the fitting problem challenging and the results naturally affected by uncertainty.

The data concerning the actual number of infected, recovered and deaths of COVID-19 are publicly available from the John Hopkins University GitHub repository [16] . For the specific case of Italy, we considered instead the GitHub repository of the Italian Civil Protection Department [41] . In the following, we present the adopted approach which is based on a strategy with two optimisation horizons (pre-lockdown and lockdown time spans) depending on the different strategies enacted by the governments of the considered European countries.

In details, we considered first the time interval t ∈ [t 0 , t ], being t the day in which lockdown started in each country (Spain, Italy and France) and t 0 the day in which the reported cases hit 200 units. The lower bound t 0 has been imposed to reduce the effects of fluctuations caused by the way in which data are measured which have a stronger impact when the number of infected is low. Once the time span has been fixed, we then considered a least square problem based on the minimization of a cost functional J which takes into account the relative L 2 norm of the difference between the reported number of infected and the reported total casesÎ(t),Î(t) +R(t) and the evolution of I(t) and I(t) + R(t) prescribed by system (26) with H ≡ 1. In practice, we solved the following constrained optimisation problem min β,γ J (Î,R, I, R) (38) where the cost functional J is a convex combination of the mentioned norms and assumes the following form

We choose p = 0.1 and we look for a minimum under the constraints 0 ≤ β ≤ 0.6, 0.04 ≤ γ ≤ 0.06.

In Table 1 we report the results of the performed parameter estimation together with the resulting attach rate R 0 defined in (36) .

Once that the contagion parameters have been estimated in the pre-lockdown time span, we successively proceeded with the estimation of the shape of the function H from the data.To estimate The second optimisation problem has been solved up to last available data for each country with daily time stepping h = 1 and over a time window of three days. This has been done with the scope of regularizing possible errors due to late reported infected and smoothing the shape of H. Both optimisation problems (38)- (39) have been tackled using the Matlab functions fmincon in combination with a RK4 integration method of the system of ODEs.

In Figure 5 , we present the result of such fitting procedure between the model (26) and the experimental data. The evolution of the estimated H(t), t ∈ [t , T ], is instead presented in the left column of Figure 6 . From this figure, it can be observed in the case of Italy how, even if the daily number of infected decreases after May 1st, the estimated H remains quite stable after this day. This behavior cannot be reproduced by using a function H(I(t)) as the one given in (31) . Instead, a function H(t, I(t)) of the form (35) , which takes into account both the instantaneous number of confirmed infections N I(t) and the total number of infected in the population N s<t I(s) up to time t, is able to better better fit the observed results. In fact, it appears that taking into account possible memory effects on the contact function is a good path to follow. In the right column of Figure 5 , the two discussed contact functions have been denoted with H 1 = H 1 (N I(t)) and (39) . The corresponding R 2 coefficient is also reported.

. Both are assumed of the form

A least square minimization approach has been finally developed to determine a > 0, b ∈ R that best fit the estimated curve whose results are presented in Table 2 . We can easily observe in the right column of Figure 6 how the function H 2 is capable to better explain the estimated values of H especially after the epidemic peak. To evaluate the goodness of fit we can finally use the so-called coefficient R 2 , as reported in Table 2 . Results show that the function H 2 appears more suitable in terms of the fitting results for all tested situations. This fact may indicate that people are rather fast to apply social distancing, and therefore to reduce their average number of contacts, whereas they tend to restore the pre-pandemic average contact rate more slowly, possibly due to further memory effects.

In this last part, we discuss the results of the S-SIR model when the contact function has the shape extrapolated in the previous paragraph. In particular, we devoted it to show whether the obtained extrapolated function H, which depends on the product between the current number of infected and the total number of infected, produces qualitative trends that are in agreement with the data. As discussed in Section 5.3, it appears that the curve of infection may be better explained in the case in which the contact function is both a function of the instantaneous number of infected and of the total number of people that contracted the disease up to time t. In order to compare qualitatively the observed curve of infected and the theoretical one, we consider the following setting for the three countries under study: ν = 4, δ = 8 and x J = 15, ∆t = 0.01, τ = 0.01, M = 10 5 . Moreover, we suppose S(t = 0) and I(t = 0) to match the relative number of susceptible and infected of each country at the time in which we start our comparison. Finally, we consider 

where we use the parameters of Table 2 that we recall here for the seek of clarity: (a, b) = (1.41, 0.15) in the case of France, and (a, b) = (0.017, 0.25) in the case of Italy. The case of Spain will be discussed later. In Figure 7 we show the profiles of the infected over time together with the shape of the function H again over time. The results show that with the choices done for the contact function, it is possible to reproduce at least qualitatively the shape of the trend of infected during the pandemic observed in Italy and in France.

It is worth to remark that the considered social parameters have been estimated only the in the case of France, see [3] , and we assumed that the initial contact distribution is the same for the Italian case. We now consider the case of Spain. For this country, according to Figure 5 , the trend of infected undergoes a deceleration during the lockdown period. This can be also clearly observed in Figure 6 where the extrapolated shape of the contact function H is shown. Let also observe that while the global behavior of this function is captured by the fitting procedure, we however lose the minimum which takes place around end of April. This minimum is responsible of the deceleration in the number of infected and can be brought back to a strong external intervention in the lifestyle of Spain country with the scope of reducing the hospitalizations. This effect can be reproduced by our model by imposing the same behavior in the function H. To that aim, the Figure 8 reports finally the profile of the infected over time together with the shape of the function H again over time for this last case. The results show that also in this case, the S-SIR model is capable to qualitatively reproduce the data. 

The development of strategies for mitigating the spreading of a pandemic is an important public health priority. The recent case of COVID-19 pandemic has seen as main strategy restrictive measures on the social contacts of the population, obtained by household quarantine, school or workplace closure, restrictions on travels, and, ultimately, a total lockdown. Mathematical models represent powerful tools for a better understanding of this complex landscape of intervention strategies and for a precise quantification of the relationships among potential costs and benefits of different options [17] . In this direction, we introduced a system of kinetic equations coupling the distribution of social contacts with the spreading of a pandemic driven by the rules of the SIR model, aiming to explicitly quantify the mitigation of the pandemic in terms of the reduction of the number of social contacts of individuals. The kinetic modeling of the statistical distribution of social contacts has been developed according to the recent results in [3] , which present an exhaustive description of contacts in the France population, divided by categories. The numerical experiments then show that the kinetic system is able to capture most of the phenomena related to the effects of partial lockdown strategies, and, eventually to maintain pandemic under control.

Control with uncertain data of socially structured compartmental models

The Log-normal Distribution

The French Connection: The First Large Population-Based Contact Survey in France Relevant for the Spread of Infectious Diseases

The theory of the nonlinear, spatially uniform Boltzmann equation for Maxwellian molecules

Economic and social consequences of human mobility restrictions under COVID-19

Mathematical models in epidemiology. With a foreword by Simon Levin

Parameter estimation and uncertainty quantification for an epidemic model

The Boltzmann equation and its applications

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecast

X: Interaction of maturation delay and nonlinear birth in population and epidemic models

On a Kinetic Model for a Simple Market Economy

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

Mathematical epidemiology of infectious diseases: model building, analysis and interpretation

Kinetic modeling of alcohol consumption

Wealth distribution under the spread of infectious diseases

An interactive web-based dashboard to track COVID-19 in real time. The Lancet Infectious Diseases

Strategies for mitigating an influenza pandemic

Estimating the number of infections and the impact of non-pharmaceutical interventions on COVID-19 in 11 European countries

Inferring the Structure of Social Contacts from Demographic Data in the Analysis of Infectious Diseases Spread

Fokker-Planck equations in the modelling of socio-economic phenomena

Fokker-Planck equations in the modelling of socio-economic phenomena

Non-Maxwellian kinetic equations modeling the evolution of wealth distribution

Relaxation schemes for nonlinear kinetic equations

A simple SIR model with a large set of asymptomatic infectives

Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Call center service times are lognormal. A Fokker-Planck description

Human behavior and lognormal distribution. A kinetic description

The Mathematics of Infectious Diseases

Prospect theory: an analysis of decision under risk

Choices, values, and frames

Determining the best population-level alcohol consumption model and its impact on estimates of alcoholattributable harms

Non-linear incidence and stability of infectious disease models

A physical basis for the generalized Gamma distribution

Log-normal distributions across the sciences: keys and clues

The reproductive number of COVID-19 is higher compared to SARS coronavirus

Social contacts and mixing pat-terns relevant to the spread of infectious diseases

On the spread of epidemics in a closed heterogeneous population

Time Relaxed Monte Carlo Methods for the Boltzmann Equation

Interacting multiagent systems: kinetic equations and Monte Carlo methods

Dipartimento della Protezione Civile. GitHub: COVID-19 Italia -Monitoraggio Situazione

Gmel, Gerhard; Statistical modeling of volume of alcohol exposure for epidemiological studies of population health: the US example

Transmission dynamics of the etiological agent of SARS in Hong Kong: Impact of public health interventions

Epidemic models with uncertainty in the reproduction

A generalization of the Gamma distribution

Trails in Kinetic Theory: foundational aspects and numerical methods

Entropy-type inequalities for generalized Gamma densities

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission