The spread of COVID-19 has been thwarted in most countries through non-pharmaceutical interventions. In particular, the most effective measures in this direction have been the stay-at-home and closure strategies of businesses and schools. However, population-wide lockdowns are far from being optimal, carrying heavy economic consequences. Therefore, there is nowadays a strong interest in designing more efficient restrictions. In this work, starting from a recent kinetic-type model which takes into account the heterogeneity described by the social contact of individuals, we analyse the effects of introducing an optimal control strategy into the system, to limit selectively the mean number of contacts and reduce consequently the number of infected cases. Thanks to a data-driven approach, we show that this new mathematical model permits us to assess the effects of the social limitations. Finally, using the model introduced here and starting from the available data, we show the effectiveness of the proposed selective measures to dampen the epidemic trends.

This article is part of the theme issue ‘Kinetic exchange models of societies and economies’.

1. Introduction

Non-pharmaceutical interventions (NPIs) have had a strong impact in limiting the spread of SARS-CoV-2 [1]. The introduced NPIs span from basic personal hygiene precautions to more severe actions such as the mandatory stay-at-home and business and schools closures.

These measures were intended to diminish transmission rates and, to this end, aimed at reducing person-to-person contacts via social distancing. In fact, even if the precise role played by those interventions has still to be carefully understood in view of a deeper comprehension of the transmission mechanisms of the disease [2,3], it has been recognized how the social attitudes of a population are a crucial feature for the spreading of airborne diseases [4,–8]. Thus, it appears that understanding the interplay between social heterogeneity and clinical characteristics of a disease is an essential feature to design efficient and robust countermeasures in the near future.

The classical epidemiological models [9] are usually based on a compartmentalization of the population, denoted in the following by $C$ , in which the agents composing the system are split into epidemiologically relevant states, like susceptible, who can contract the disease, and several other classes of agents who may transmit it. Therefore, to model the transmission of a pathogen through a population, if $N_{i} (t)$ represents the number of individuals in the compartment $i \in C$ with $| C | = n > 1$ , one generally considers a system of differential equations that can be written in compact form as

\frac{d N (t)}{d t} = P (N (t)),

1.1

where

N (t)

is the

n

-dimensional vector with components

N_{i} (t)

, and

P

is the

n

-dimensional vector whose components

P_{i} (N)

i \in C

, are the transition rates between compartments, namely the variation of the number of individuals in the

i

th compartment due to interactions with individuals in other compartments.

Several variations of system (1.1) have been introduced to forecast the penetration of a virus and to prevent its spread, see e.g. [9,–12]. Generalizations of the mentioned compartmental models include additional factors characterizing the transmission dynamics, such as the introduction of age structure, contacts heterogeneity and mobility, see e.g. [–4,13,15].

In this work, we propose a kinetic approach to derive consistent transition rates in a population characterized by a heterogeneous number of contacts together with selective NPIs. Indeed, the number of social contacts has been recently identified by the scientific community as one of the main relevant causes of the potential pathogen transmission [6,–8,16]. In this respect, the NPIs imposed by the different governments were mainly based on restrictions whose scope was to reduce as much as possible the number of contacts among individuals.

In more detail, we aim at giving a deeper understanding of the mitigation effects due to the reduction of the social interactions among individuals in the dynamics of SARS-CoV-2. This is done through a combination of a multiscale mathematical model which takes into account both the statistical information on the distribution of social contacts in a society as well as a control strategy whose target is to point the population towards a given benchmark number of contacts. The underlying theoretical framework that we consider is the one coming from the kinetic theory of collective social phenomena [17]. The diffusion of epidemics is here described as a result of the interactions between a large number of individuals, each one having a different social propensity. In this respect, as shown in Dimarco et al. [18] and Zanella et al. [19], the fundamental tools of statistical physics are capable of providing a useful insight into epidemiological dynamics linking classical compartmentalization approaches with a statistical understanding of the social aspects. In fact, the essential multiscale nature of kinetic theory permits us to determine the resulting macroscopic (or aggregate) and measurable characteristics of the evolution of the disease [19]. One of the key directions we pursue in this research is given by the combination of consolidated theoretical and modelling tools [17] with data-driven techniques. More precisely, starting from the data at our disposal, we calibrate our mathematical model passing through the determination of the relevant epidemiological parameters of the epidemic from the data. Subsequently, we concentrate on modelling the lockdown policies through an optimal control approach. In particular, we show how different kinds of 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 dedicate ourselves to the study of the effects induced by alternative control strategies.

We stress that, even if in the rest of the paper we mainly refer our results to a compartmental model with exposed and asymptomatic, i.e. the so-called SEIAR compartmentalization, the present approach is quite general, and applies both to simple [9,–12] as well as to more complex epidemiological models [20,–22]. For this reason, we will keep the presentation as general as possible.

2. Modelling epidemic dynamics in the presence of contact heterogeneity

We discuss in this part a possible improvement of the classical compartmental description of epidemic spreading which takes into account the statistical description of the social behaviour of individuals. This permits us to determine the distribution of the average number of connections in a system of $N$ interacting individuals, see [18,19], and to relate this quantity with the infection dynamics. In this new approach, agents which are supposed to belong to a given compartment in term of response to a given pathogen are additionally characterized by their daily number of contacts $x \geq 0$ , whose statistical description is obtained through the distribution functions $f_{i} (x, t)$ , $t \geq 0$ , where $i \in C$ denotes the $i$ th compartment (susceptible, infected, asymptomatic, etc.), such that

\sum_{i \in C} f_{i} (x, t) = f (x, t), \int_{R_{+}} f (x, t) d x = N (t) .

In this setting, the time evolution of the functions

f_{i} (x, t)

i \in C

is obtained by integrating the compartmentalization epidemiological description with the simultaneous characterization of the formation of social contacts in a society [18]. This merging is expressed by the system

\frac{\partial f (x, t)}{\partial t} = P (x, f (x, t)) + \frac{1}{τ} Q (f (x, t)) .

2.1

As in (1.1),

P

is the

n

-dimensional vector whose components

P_{i} (x, f)

i \in C

, are the transition rates between compartments (now depending both on the contact variable

x

and on the vector

f (x, t)

), while

Q

is the

n

-dimensional vector whose components

Q_{i} (f_{i})

i \in C

are suitable relaxation operators which describe the formation of equilibrium distribution of contacts of the underlying compartments. Details about this emerging equilibrium distribution, together with its main properties, will be provided in the following. Lastly,

τ > 0

measures the time scale at which these equilibrium distributions of social contacts are reached by relaxation.

In general, the description of the time evolution of epidemic spreading in terms of the distributions $f_{i} (x, t)$ , as given by (2.1), is richer than the classic one obeying (1.1). This is due to the fact that, in (2.1), the system depends on an additional variable allowing a more realistic characterization of the transition rates between compartments. In particular, by assuming that the pathogen transmission does not depend upon the contact among individuals, the proposed description (2.1) collapses to the standard approach (1.1) as shown in [18].

In our model, the transition rates can now be assumed to be dependent on the intensity of social contacts in the compartments. As detailed later, in fact, they are obtained by taking into account, in a very natural way, the heterogeneity of the population in the form of the variance of the distributions $f_{i} (x, t)$ [18,23]. This point of view provides a better understanding of the role of the different individual behaviours in the spreading of a disease [5].

In general, the knowledge of the densities $f_{i} (x, t)$ allows us to evaluate various macroscopic quantities, independent from the number of contacts, the so-called moments of the distributions, obtained by integration with respect to the variable $x \in R_{+}$ . In this respect, the zero-order moment gives the number of individuals $N_{i} (t)$ in each compartment at time $t \geq 0$

N_{i} (t) = \int_{R^{+}} f_{i} (x, t) d x, i \in C .

2.2

In the same way, one can quantify the mean numbers

{\bar{x}}_{i} (t)

of contacts of individuals in each compartment at time

t \geq 0

evaluating

{\bar{x}}_{i} (t) = \frac{1}{N_{i} (t)} \int_{R^{+}} x f_{i} (x, t) d x, i \in C .

2.3

Let now $f_{S}$ , $S \in C$ , indicate the contact distribution of susceptible people (people at risk of becoming infected), and let the subset $I \subset C$ , $| I | < | C |$ include the compartments of people that may transmit the infection. Since the transition rate $P_{S} (x, f (x, t))$ contains the details about the spreading of the infection among susceptible individuals, the new formulation (2.1) allows us in particular to characterize this transmission in terms of the intensity of social contacts between the compartments of susceptible and infectious individuals. This intensity is measured at time $t \geq 0$ by relating the densities $f_{S} (x, t)$ and $f_{j} (x, t)$ , $j \in I$ . Following [18], the transition rate $P_{S}$ is then expressed by introducing the function $K (f_{S}, {f_{j}}_{j \in I})$ , the so-called local incidence rate, given by

K (f_{S}, {f_{j}}_{j \in I}) (x, t) = f_{S} (x, t) \int_{R^{+}} \sum_{j \in I} κ_{j} (x, y) f_{j} (y, t) d y .

2.4

In (2.4), the contact functions

κ_{j} (x, y)

are non-negative functions growing with respect to the number of contacts

x

and

y

of the populations of susceptible and infectious agents in the compartment

j \in I

, and such that

κ_{j} (x, 0) = 0

. In particular, the choice

κ_{j} (x, y) = β_{j} x^{α} y^{α}, j \in I

with constants

α, β_{j} > 0

corresponds to consider an incidence rate dependent on the product of the number of contacts of susceptible and infectious people. When

α = 1

, for example, the incidence rate takes the simpler form

K (f_{S}, {f_{j}}_{j \in I}) (x, t) = x f_{S} (x, t) \sum_{j \in I} β_{j} {\bar{x}}_{j} (t) N_{j} (t) .

2.5

We discuss now the details of the dynamics of formation of the social contacts which are modelled by the operators $Q_{i} (f_{i})$ in (2.1). These operators are of Fokker–Planck-type with variable coefficient of diffusion and linear drift

Q_{i} (f_{i} (x, t)) = \frac{1}{2} \frac{\partial}{\partial x} [σ^{2} \frac{\partial (x f_{i} (x, t))}{\partial x} + μ (\frac{x}{{\bar{x}}_{i} (t)} - 1) f_{i} (x, t)], i \in C,

2.6

coupled with no-flux boundary conditions such that the Fokker–Planck-type equation

\frac{\partial}{\partial t} f_{i} (x, t) = Q_{i} (f_{i}) (x, t)

is mass and momentum preserving. In (2.6),

σ

and

μ

are positive constants which characterize the intensities of the diffusion and drift terms. The above description through Fokker–Planck operators is obtained thanks to an upscaling of an underlying microscopic dynamics which takes into account the typical aspects of human behaviour for what concerns the social attitudes [18]. In particular, the microscopic description takes into account the interactions due to the basic daily needs such as work duties, schools and the search, in the absence of epidemics, for opportunities for socialization. These aspects are embedded in the microscopic model by relying on the so-called prospect theory of Kahneman & Tversky [24] as detailed later.

The kernels of the operators $Q_{i}$ furnish the class of equilibrium densities, that are obtained by solving the first-order differential equations

σ \frac{\partial (x f_{i}^{\infty} (x))}{\partial x} + μ (\frac{x}{{\bar{x}}_{i}} - 1) f_{i}^{\infty} (x) = 0, i \in C, \sum_{i \in C} N_{i} (t) = N (t) .

It is immediate to verify that, in this setting, the equilibrium distributions

f_{i}^{\infty} (x)

of given mass

N_{i}

and local mean

{\bar{x}}_{i}

are the Gamma distributions

f_{i}^{\infty} (x) = f_{i}^{\infty} (x, {\bar{x}}_{i}, λ) = N_{i} {(\frac{λ}{{\bar{x}}_{i}})}^{λ} \frac{1}{Γ (λ)} x^{λ - 1} \exp {- \frac{λ}{{\bar{x}}_{i}} x},

2.7

characterized by the shape parameter

λ = μ / σ

. It is important to stress the role of the positive constant

λ

which is a measure of the heterogeneity of the population, statistically distributed according to (2.7). It is also remarkable that the above equilibrium dynamics, satisfying the condition

λ > 1

, is in agreement with the experimental results in [6]. This study represents one of the first attempts in which a large scale experimental measure of the relation between contact and spread of infectious diseases has been done. To better characterize the equilibrium distribution, we also note that the condition on

λ

identifies the Gamma densities for which

f_{i}^{\infty} (x = 0) = 0

. Moreover, higher-order moments of this class of distributions depend only on

N_{i}

{\bar{x}}_{i}

and

λ

. In particular, the local variance is expressed by

V_{i} = \frac{1}{N_{i}} \int_{R^{+}} {(x - {\bar{x}}_{i})}^{2} f_{i}^{\infty} (x, {\bar{x}}_{i}, λ) d x = \frac{{\bar{x}}_{i}^{2}}{λ} .

2.8

Hence, expression (2.8) permits us to relate the large values of

λ

with small heterogeneity of the population as claimed before.

Before introducing in the model the aspects related to the control policies, we briefly discuss how the Fokker–Planck-type operators (2.6) can be derived from collisional kinetic equations describing the evolution of social connections [25,26]. These kinetic equations are obtained by looking to a possible characterization of the personal behaviour of a single agent in terms of daily contacts starting from some well-established assumptions. Indeed, the total amount of contacts can be simply viewed as the result of a repeated upgrading: the daily life of each person is based on a certain number of activities each carrying a certain number of contacts. Moreover, independently of the personal choices or needs, this number contains a certain amount of randomness, that has to be taken into account. These assumptions give, letting $x^{'}$ be the microscopic number of social contacts resulting from the interaction,

x^{'} = x - Φ_{ϵ} (\frac{x}{{\bar{x}}_{i}}) x + η_{ϵ} x,

2.9

with

η_{ϵ}

a centred random variable such that

E [η_{ϵ}^{2}] = ϵ σ^{2}

with bounded moments up to order three, and

Φ_{ϵ} (\cdot)

a transition function of the form

Φ_{ϵ} (s) = μ \frac{e^{ϵ (s - 1)} - 1}{e^{ϵ (s - 1)} + 1}, s \geq 0.

2.10

The transition function has been constructed by taking into account that while is very easy to reach a high number of social contacts for need or will, going below a given threshold is very difficult, since contacts are forced by basic activities which cannot be avoided in general. This asymmetry between growth and decrease, as exhaustively discussed in [18,25,28], can be suitably modelled by (2.10) which plays the role of the so-called value function in the prospect theory of Kahneman & Tversky [24] and implements the advocated asymmetry in the contact distribution formation.

Hence, with the above discussed choices, at the aggregate scale the evolution of a contact distribution $f_{i, ϵ}$ is given by the Boltzmann-type equation in weak form

\frac{\partial}{\partial t} \int_{R_{+}} φ (x) f_{i, ϵ} (x) d x = \frac{1}{ϵ} E [\int_{R_{+}} B (x) (φ (x^{'}) - φ (x)) f_{i, ϵ} (x, t) d x],

2.11

where

B (x) = 1 / x

is the interaction kernel assigning a low probability to interactions for individuals exhibiting a large number of contacts and assigning a high probability to interactions when

x

is small. Letting finally

ϵ \to 0^{+}

(the limit of grazing interactions) one shows that (2.11) converges, up to extraction of a subsequence, to the Fokker–Planck equation (2.6), see e.g. [18,29]. This result is particularly useful since it permits us to characterize the equilibrium state of the system, i.e. it permits us to find the equilibrium Gamma distribution (2.7) as discussed in [18] and experimentally found in [6]. Moreover, it allows us to analytically consider and study control problems at the kinetic level, i.e. to define effective actions at the microscopic level the scope of which is to modify the behaviour of the single individuals. Subsequently, thanks to well established kinetic theory tools [17], it is possible to upscale the microscopic actions in order to compute its reflection at the aggregate, measurable, scale of description and of interest.

(a) Observable effects of selective social restrictions

In this part, we model the lockdown measures intended to reduce the social interactions by introducing an additive term $u$ at the level of the microscopic formation of the contacts (2.9) to limit selectively the social activities. This gives the controlled elementary interaction

x^{'} = x - Φ_{ϵ} (\frac{x}{{\bar{x}}_{i}}) x + \sqrt{ϵ τ} W (x) u + η_{ϵ} x, W \geq 0.

2.12

The function

W (x)

is such that if

W \equiv 1

the control variable

u

is independent from the number

x \geq 0

of the social contacts, while if

W (x) \propto x

the control is considered heavier for people with a higher number of contacts. Note that in (2.12), the size of the controlled variable is tuned by the small parameters

ϵ

and

τ

, which ensures that the control acts at the right scale with respect to both the grazing variable

ϵ

characterizing the passage from a Boltzmann to a Fokker–Planck dynamics and the relaxation variable

τ

characterizing the speed at which an equilibrium state for the distribution of contact is reached, compared to the speed at which the epidemic dynamics travels.

Then, the so-called optimal control $u^{*}$ is determined as the minimizer of a given cost functional

u^{*} = arg min_{u \in U} J (x^{'}, u),

2.13

subject to the constraint (2.12). In (2.13), the minimum is taken on the set

U

of all admissible controls. A classical choice for

J

is a quadratic cost functional of the form

J_{i} (x^{'}, u) = \frac{1}{2} E [{(x^{'} - x_{T, i})}^{2} + ν u_{i}^{2}] i \in C,

2.14

ν > 0

being a penalization coefficient measuring the economic and social cost of the control and

x_{T, i} > 0

the desired target number of social contacts for the compartment

i \in C

. According to (2.14), the cost of the control increases quadratically with the distance to the desired target of average number of contacts. The presence in (2.14) of the mean operator is related to the fact that, by hypothesis, a control cannot act on the random fluctuations introduced to mimic the unanticipated daily interactions among agents. Let us observe that other convex cost functionals may be considered as well. However, they often lead to problems whose analytical solution cannot be obtained explicitly. Proceeding now as in [30], the control introduced in (2.12) can be solved explicitly giving

u = - \frac{\sqrt{ϵ τ} W (x)}{ν + ϵ τ W^{2} (x)} (x - x_{T, i} - Φ_{ϵ} (\frac{x}{{\bar{x}}_{i}}) x),

2.15

which generates the optimal constrained law of growth of the social contacts

x_{*}^{'} = x - \frac{ν}{ν + ϵ τ W^{2} (x)} Φ_{ϵ} (\frac{x}{{\bar{x}}_{i}}) x - \frac{ϵ τ}{ν + ϵ τ W^{2} (x)} (x - x_{T, i}) W^{2} (x) + η_{ϵ} x .

2.16

Now, in the presence of the above described microscopic controlled interactions, one can consider as before the construction of a kinetic collisional equation for what concerns the description of the formation of the social structure

\frac{\partial}{\partial t} \int_{R_{+}} φ (x) f_{i, ϵ} (x) d x = \frac{1}{ϵ} E [\int_{R_{+}} \frac{(φ (x_{*}^{'}) - φ (x))}{x} f_{i, ϵ} (x, t) d x],

2.17

where we recall that the frequency of interaction has been chosen such that

B (x) = 1 / x

. Then, thanks to the passage to the grazing limit one can show that the controlled kinetic model converges to a Fokker–Planck type equation which contains an additional drift term, that quantifies the role of the control mimicking lockdown measures. In fact, since by hypothesis

ϵ ≪ 1

, we can write

φ (x_{*}^{'}) - φ (x) = (x_{*}^{'} - x) \partial_{x} φ (x) + \frac{1}{2} {(x_{*}^{'} - x)}^{2} \partial_{x}^{2} φ (x) + \frac{1}{6} {(x_{*}^{'} - x)}^{3} \partial_{x}^{3} φ (\hat{x}),

with

\hat{x} \in (min (x, x_{*}^{'}), max (x, x_{*}^{'}))

. Then using (2.16) one gets

\begin{aligned} \frac{\partial}{\partial t} \int_{R_{+}} φ (x) f_{i, ϵ} (x) d x & = - \frac{1}{ϵ} \int_{R_{+}} \frac{ν}{ν + ϵ τ W^{2} (x)} Φ_{ϵ} (\frac{x}{{\bar{x}}_{i}}) \partial_{x} φ (x) f_{i, ϵ} (x, t) d x \\ - \int_{R_{+}} \frac{τ}{ν + ϵ τ W^{2} (x)} \frac{(x - x_{T, i})}{x} W^{2} (x) \partial_{x} φ (x) f_{i, ϵ} (x, t) d x \\ + \frac{σ^{2}}{2} \int_{R_{+}} \partial_{x}^{2} φ (x) x f_{i, ϵ} (x, t) d x + R_{Φ_{ϵ}} (f_{i, ϵ}) (x, t), \end{aligned}

where the remainder

R_{Φ_{ϵ}} (f_{i, ϵ}) (x, t)

can be shown to go to zero in the limit

ϵ \to 0^{+}

since by assumption

φ

and its derivatives are bounded in

R_{+}

and the stochastic variable

η

is assumed to have bounded moments at least of order 3 as previously stated. We point the reader to [17,18] for a detailed discussion on this aspect. In the limit

ϵ \to 0^{+}

, one gets

\begin{aligned} \frac{\partial}{\partial t} \int_{R_{+}} φ (x) f_{i} (x) d x & = - \int_{R_{+}} Φ (\frac{x}{{\bar{x}}_{i}}) \partial_{x} φ (x) f_{i} (x, t) d x \\ - \int_{R_{+}} \frac{τ}{ν} \frac{(x - x_{T, i})}{x} W^{2} (x) \partial_{x} φ (x) f_{i} (x, t) d x + \frac{σ^{2}}{2} \int_{R_{+}} \partial_{x}^{2} φ (x) x f_{i} (x, t) d x \end{aligned}

2.18

since for

ϵ ≪ 1

we have

Φ_{ϵ} (\frac{x}{{\bar{x}}_{i}}) \approx ϵ \frac{μ}{2} (\frac{x}{{\bar{x}}_{i}} - 1) = ϵ Φ (\frac{x}{{\bar{x}}_{i}}) .

Integrating back by part (2.18), and considering the no flux boundary condition at

x = 0^{+}

to preserve mass and momentum, the relaxation equations take the form

\frac{\partial f_{i} (x, t)}{\partial t} = \frac{1}{τ} Q_{i} (f_{i} (x, t)) + C_{i} (f_{i} (x, t)), i \in C,

2.19

where the operator

C_{i} (f_{i})

in (2.19) is given by

C_{i} (f_{i} (x, t)) = \frac{1}{ν} \frac{\partial}{\partial x} [\frac{x - x_{T, i}}{x} W^{2} (x) f_{i} (x, t)],

2.20

and

Q_{i} (f_{i} (x, t))

remains the same as in (2.6). Note in particular that in (2.20) the control term does not depend on the time scale

τ

of the Fokker–Planck operator.

Therefore, in the presence of social restrictions modelled by the control term introduced above, system (2.1) is modified as

\frac{\partial f (x, t)}{\partial t} = P (x, f (x, t)) + C (f (x, t)) + \frac{1}{τ} Q (f (x, t)),

2.21

where the

n

-dimensional vector

C

quantifies the effect of the control on the epidemic spreading. According to the previous analysis, the components of the vector control

C

are given by the operators of type (2.20), where the values of the penalization coefficients and of the target values may depend on the specific compartment

i \in C

In the next part, we focus on a specific compartmental model for which we give the details of the controlled strategy for the reduction of the social contacts.

Remark 2.1.

We observe that the cost functional (2.14) introduced at the level of microscopic interactions can be seen as an instantaneous approximation of a macroscopic functional obtained by considering (2.15) in the limit $ϵ, τ \to 0^{+}$ whose form is

J_{i} = \frac{1}{2} \int_{R_{+}} (1 + \frac{W^{2} (x)}{ν}) {(x - x_{T, i})}^{2} f_{i}^{\infty} (x) d x, i \in C .

In the above functional, we consider the contributions at the final time

T

and the steady distribution

f_{i}^{\infty}

of contacts defined in (2.7) instead of the instantaneous single interaction among agents.

(b) A SEIAR-type compartmentalization

We consider the case $C = {S, E, I, A, R}$ , i.e. we concentrate on a model with exposed and asymptomatic compartments and we discuss the observable effects of the introduced control at the level of the mean number of connections in two leading examples.

The choice $C = {S, E, I, A, R}$ gives, according to (2.21), the following system of kinetic equations relating the epidemiological dynamics with the number of daily contacts expressed in terms of a probability distribution $f_{i} (x, t)$

\begin{aligned} \frac{\partial}{\partial t} f_{S} (x, t) & = - x f_{S} (x, t) \sum_{j \in {I, A}} β_{j} {\bar{x}}_{j} N_{j} (t) + C_{S} (f_{S}) (x, t) + \frac{1}{τ} Q_{S} (f_{S}) (x, t), \\ \frac{\partial}{\partial t} f_{E} (x, t) & = x f_{S} (x, t) \sum_{j \in {I, A}} β_{j} {\bar{x}}_{j} N_{j} (t) - γ_{E} f_{E} (x, t) + C_{E} (f_{E}) (x, t) + \frac{1}{τ} Q_{E} (f_{E}) (x, t), \\ \frac{\partial}{\partial t} f_{I} (x, t) & = ξ γ_{E} f_{E} (x, t) - γ_{I} f_{I} (x, t) + C_{I} (f_{I}) (x, t) + \frac{1}{τ} Q_{I} (f_{I}) (x, t), \\ \frac{\partial}{\partial t} f_{A} (x, t) & = (1 - ξ) γ_{E} f_{E} (x, t) - γ_{A} f_{A} (x, t) + C_{A} (f_{A}) (x, t) + \frac{1}{τ} Q_{A} (f_{A}) (x, t), \\ and \frac{\partial}{\partial t} f_{R} (x, t) & = γ_{I} f_{I} (x, t) + γ_{A} f_{A} (x, t) + C_{R} (f_{R}) (x, t) + \frac{1}{τ} Q_{R} (f_{R}) (x, t), \end{aligned}}

2.22

where, to have a completely definite system, it remains to clarify the role of the parameter

ξ

which gives the fraction of exposed individuals becoming successively infected and its counterpart

(1 - ξ)

which gives the fraction of asymptomatic individuals after the exposure to the virus. In addition, the functions

γ_{I}

and

γ_{A}

represent the recovery rates of infected and asymptomatic individuals respectively, which give the frequency at which individuals move to the compartment of recovered individuals. In the present setting, the choice of the constants

β_{I}, β_{A} > 0

takes into account that the population of the infected compartment is considered recognized by the health system through a swab test and therefore subject to isolation. We collect in the asymptomatic compartment unrecognized infections that are not subject to isolation. We highlight that we preferred to neglect the compartment of deceased people since they do not influence the spreading of the disease. Anyway, the inclusion of this new compartment to assess the risks related to the infection can easily be considered.

Since the $Q_{i}$ and $C_{i}$ components in (2.22), coupled with no-flux boundary conditions, are mass preserving, we can integrate the equations of the system with respect to the contact variable $x$ to obtain a macroscopic system of equations describing the evolution of the fraction of individuals in each compartment independently of the number of daily social contacts. This gives the system

\begin{aligned} \frac{d N_{S}}{d t} & = - {\bar{x}}_{S} N_{S} \sum_{j \in {I, A}} β_{j} {\bar{x}}_{j} N_{j}, \\ \frac{d N_{E}}{d t} & = {\bar{x}}_{S} N_{S} \sum_{j \in {I, A}} β_{j} {\bar{x}}_{j} N_{j} - γ_{E} N_{E}, \\ \frac{d N_{I}}{d t} & = ξ γ_{E} N_{E} - γ_{I} N_{I}, \\ \frac{d N_{A}}{d t} & = (1 - ξ) γ_{E} N_{E} - γ_{A} N_{A}, \\ and \frac{d N_{R}}{d t} & = γ_{I} N_{I} + γ_{A} N_{A}, \end{aligned}}

2.23

which is not closed since the evolution in time of the vector solution

N (t)

depends on the mean number of contacts

{\bar{x}}_{i}

. Since the

Q_{i}

components, under the same hypothesis on the boundary fluxes, are also momentum preserving, by multiplying the distributions

f_{i} (x, t)

x

and integrating with respect to the contact variable, it is possible to recover a further macroscopic set of equations describing the constrained dynamics of the mean number of social contacts in each compartment. However, one immediately realizes that this second set of macroscopic observables depends upon higher order moments of the distributions

f_{i} (x, t)

. In other words, in general, the system of macroscopic equations which can be obtained by integration of (2.22) over the number of social contacts is not closed. In fact, system (2.23) depends on the time evolution of

{\bar{x}}_{i} (t), i = S, E, I, A, R

while the system giving the time evolution of

{\bar{x}}_{i} (t), i = S, E, I, A, R

depend upon second-order moments of the

f_{i} (x, t), i = S, E, I, A, R

. This represents nothing else than the classical so-called closure problem often encountered in kinetic theory [17].

However, letting $τ \to 0$ forces the distributions $f_{i} (x, t)$ to converge towards the Gamma vector distribution of components (2.7) with mass fractions $N_{i} (t)$ , and local mean values ${\bar{x}}_{i} (t)$ , $i \in C = {S, E, I, A, R}$ . This amounts to assuming that social dynamics are much faster than epidemic dynamics. In this situation, it is possible to obtain a closed set of macroscopic equations. In particular, the second order moments of the distributions $f_{i} (x, t)$ can be explicitly computed in terms of the mass fractions and mean values resorting to (2.8). This gives the system

\begin{aligned} \frac{d {\bar{x}}_{S}}{d t} & = - \frac{1}{λ} {\bar{x}}_{S}^{2} \sum_{j \in {I, A}} β_{j} {\bar{x}}_{j} N_{j} + C (f_{S}^{\infty}), \\ \frac{d {\bar{x}}_{E}}{d t} & = - {\bar{x}}_{S} \frac{N_{S}}{N_{E}} (\frac{λ + 1}{λ} {\bar{x}}_{S} - {\bar{x}}_{E}) \sum_{j \in {I, A}} β_{j} {\bar{x}}_{j} N_{j} + C (f_{E}^{\infty}), \\ \frac{d {\bar{x}}_{I}}{d t} & = ξ γ_{E} ({\bar{x}}_{E} - {\bar{x}}_{I}) \frac{N_{E}}{N_{I}} + C (f_{I}^{\infty}), \\ \frac{d {\bar{x}}_{A}}{d t} & = (1 - ξ) γ_{E} ({\bar{x}}_{E} - {\bar{x}}_{A}) \frac{N_{E}}{N_{A}} + C (f_{A}^{\infty}), \\ and \frac{d {\bar{x}}_{A}}{d t} & = γ_{I} ({\bar{x}}_{I} - {\bar{x}}_{R}) \frac{N_{I}}{N_{R}} + γ_{A} ({\bar{x}}_{A} - {\bar{x}}_{R}) \frac{N_{A}}{N_{R}} + C (f_{R}^{\infty}), \end{aligned}}

2.24

where

C (f_{i}^{\infty}) (t) = \frac{1}{ν N_{i} (t)} \int_{R_{+}} (\frac{x_{T, i}}{x} - 1) W^{2} (x) f_{i}^{\infty} (x, t) d x .

2.25

Equations (2.24) coupled with the system (2.23) describe in closed form the time evolution of an epidemic where the transition functions from one compartment to another depend upon the mean level of social contacts in the population. In addition, these latter depend on the heterogeneity of the population through the parameter

λ = μ / σ

characterizing the distributions

f_{i} (x, t)

, while the term (2.25) models the effects of the lockdown measures.

In the rest of the paper, we focus on two specific benchmark cases for which it is possible to compute explicitly the effect of the social restrictions on the system. In particular, we first concentrate on the case of uniform restrictions, corresponding to $W (x) \equiv 1$ , and second on the case of selective restrictions. This last situation describes an agent system in which restrictions are stronger for individuals exhibiting a larger number of contacts, and it is based on the choice $W (x) = x / {\bar{x}}_{i}$ . By direct computation, we have

C (f_{i}^{\infty}) = {\begin{cases} \frac{1}{ν} (\frac{λ}{λ - 1} \frac{x_{T, i}}{{\bar{x}}_{i}} - 1) & W (x) = 1, \\ \frac{1}{ν} (\frac{1 + λ}{λ} {\bar{x}}_{i} - x_{T, i}) & W (x) = \frac{x}{{\bar{x}}_{i}} . \end{cases}

2.26

Hence, for small penalizations $ν > 0$ of the control, the mean number of connections stabilizes towards the values

{\bar{x}}_{i}^{\infty} = {\begin{cases} x_{T, i} \frac{λ}{λ - 1} > x_{T, i} & W (x) = 1, \\ \frac{λ}{1 + λ} x_{T, i} < x_{T, i} & W (x) = \frac{x}{{\bar{x}}_{i}} . \end{cases}

It is worth noting that, in both situations, the role of the penalization parameter

ν > 0

is limited to a scaling of time, without changing the final goal of the control. Thus, a smaller penalization parameter is reflected in a quicker convergence of the solution towards the target value as expected. We also note that the selective control, i.e. the one especially operating on individuals with larger social nets, allows us to reach for the global population a value below the desired target

x_{T}

. In particular, the mean contact rate is lower than

x_{T}

of an amount directly proportional to the degree of heterogeneity of the population. This can be explained by the fact that individuals with larger social networks are mainly responsible for the creation of the contact dynamics. This situation suggests that it is important to act on this class of people in order to limit the spread of a disease instead of taking measures with uniform impact over the whole population.

From an operative point of view, it is clear that the application of a uniform control, which is based on restrictive measures uniformly valid for the global population, are simpler to apply with respect to measures which provide for heavier restrictions for the categories that show a greater number of social contacts. In this direction, the previous analysis seems to suggest that a mixture of the two types of controls will present the best compromise between containment results and effective applicability of containment measures.

3. Numerical tests

In this section, using (2.23)–(2.24), we compare possible interventions for the containment of an epidemic based on the limitation of the mean number of contacts. These non-pharmaceutical measures are designed to drive the population towards a prescribed target dependent on the relevant epidemiological state of agents, namely $x_{T, i}$ , $i \in C$ in the considered framework.

To this end, we first estimate the relevant epidemiological parameters of the second order macroscopic SEIAR model taking into account the available dataset provided by John Hopkins University.¹ After the calibration of the model with the data at our disposal, we assess the impact of different control strategies leading to the observed epidemic dynamics in different countries where lockdown policies have been implemented.

We will mainly focus on the so-called first wave of infection lasting the first half of 2020 since in this period similar contact restrictions have been proposed in all Western countries. In the following, we concentrate on the cases of Italy and the UK.

(a) Calibration of the second order model

We consider two countries that adopted the stay-at-home policies during early stages of the infection. In table 1, we report the days of implementation of such measures together with a sketch of the timeline of the epidemic. In the second column, we highlight the begin of restrictions at the global level for each different country. Anyway, it is worth mentioning that often local restriction measures, i.e. isolation of a small portion of a territory, have been implemented as well.

Table 1. Timeline of COVID-19 pandemic.
Collapse
country	first detected case	begin restrictions	$β_{I}$	$λ$
Italy	30 January	9 March	0.0042	4.9593
UK	31 January	23 March	0.0043	5.8502

Before proceeding with the estimation of the model parameters some remarks are deserved. First, we stress that the calibration of epidemiological models is generally difficult and requires special care at the numerical level. Furthermore, available data typically represent a lower bound with respect to the true numbers of an infection, especially in its initial phases, due to a limited testing capacity. Here, we do not consider this aspect and we refer to recent attempts to overcome this issue based on tools of uncertainty quantification [4,19] that can be successfully employed in this context. For related approaches, we mention also [31,–34].

For the calibration of the model (2.23)–(2.24), we have adopted a bilevel approach. In particular, in the phase preceding the lockdown, we estimated the epidemic parameters in the unconstrained regime assuming that no restriction of social contacts was activated, i.e. $C \equiv 0$ . We also fixed the known clinical parameters according to values reported in the literature: see [21,35]. We then set $γ_{E} = 1 / 3.32$ , $γ_{I} = 1 / 10$ and $γ_{A} = 2 γ_{I}$ , $β_{A} = 2 β_{I}$ and we also assumed that $ξ = 0.2$ , such that $80 %$ of infections come from asymptomatic or paucisymptomatic individuals, in agreement with the recent serological campaigns promoted in Italy.² The proposed choice for $β_{A}, β_{I} > 0$ is coherent with existing literature, see e.g. [36].

We then solved a least-square problem based on the minimization of the relative $L^{2}$ norm of the difference between reported number of infected ${\hat{N}}_{I}$ and recovered ${\hat{N}}_{R}$ , and the theoretical evolution of the model $N_{I} (t)$ and $N_{R} (t)$ for $t \in [t_{0}, t_{ℓ}]$ with $N_{E} (t_{0}) = N_{I} (t_{0}) = N_{A} (t_{0}) = N_{R} (t_{0}) = 1$ , i.e. we assume that in $t_{0}$ the pandemic was nearly started. We also considered ${\bar{x}}_{S} (t_{0}) = {\bar{x}}_{E} (t_{0}) = {\bar{x}}_{A} (t_{0}) = {\bar{x}}_{R} (t_{0}) = 10$ , in agreement with the experimentally observed mean number of contacts in a Western country before the pandemic [6], and ${\bar{x}}_{I} (t) \equiv 3$ for recognized infected, corresponding to the average number of family contacts. Thus once the known parameters have been fixed, we considered the following minimization problem

min_{β_{I}, λ} [(1 - θ) | | N_{I} (t) - {\hat{N}}_{I} (t) | |_{L^{2} ([t_{0}, t_{ℓ}])} + θ | | N_{R} (t) - {\hat{N}}_{R} (t) | |_{L^{2} ([t_{0}, t_{ℓ}])}],

where

θ \in [0, 1]

| | \cdot | |_{L^{2} ([t_{0}, t_{ℓ}])}

the relative norm over the time horizon

[t_{0}, t_{ℓ}]

and under the constraints

β_{I} \in [0, 1 / 100]

λ \in [0, 50]

. In the third and fourth columns of table 1, we report the estimated epidemiological parameters obtained for

θ = 10^{- 3}

to get a better resolution of the infected dynamics. It is worth noticing that, according to these values, Italy manifests higher heterogeneity linked to the contact distribution compared to the UK.

(b) Assessing the impact of non-pharmaceutical interventions

Once the relevant epidemiological parameters are estimated, we focus on the subsequent lockdown phase where we look for the optimal value $x_{T}$ in the control term which permits us to fit the data. We assume this value $x_{T}$ to be equal for each compartment and we study the two cases of uniform and selective restrictions, respectively obtained by fixing $W (x) = 1$ and $W (x) = x / {\bar{x}}_{i}$ and $C (f^{\infty})$ of the form (2.26) in the dynamics. Hence, we solve an optimization problem in the lockdown time interval $[t_{ℓ} + 1, t_{f}]$ , for a sequence of time steps $t^{n}$ over a time window of one week

min_{x_{T} (t^{n}) \in R^{+}} [(1 - θ) | | N_{I} (t) - {\hat{N}}_{I} (t) | |_{L^{2} ([t^{n} - κ_{ℓ}, t^{n} + κ_{r}])} + θ | | N_{R} (t) - {\hat{N}}_{R} (t) | |_{L^{2} ([t^{n} - κ_{ℓ}, t^{n} + κ_{r}])}],

with

κ_{ℓ} = 3

κ_{r} = 4

In figure 1, we show the obtained trends for the registered cases for Italy and the UK. We may observe how for both controls the macroscopic trends are coherent with the observed ones. Therefore, the adopted fitting procedure allows us to compare the estimated target $x_{T, i}$ for mean contacts in the considered control strategies. In particular, we observe how the uniform control implies a stronger contact reduction with respect to the selective one, while a similar cost, computed as the sum $J_{S} + J_{E} + J_{A} + J_{R}$ , is produced by the two containment strategies. In other words, the adoption of selective-type containment strategies based on the average number of contacts allows us to mitigate the severity of lockdown policies.

In the previous part, we focused on getting the target mean number of contacts for both the control strategies. We remark that the NPIs implemented by both countries have been mainly uniform in the first wave of the epidemic, meaning that, in the present model framework, the average number of contacts were approximately the ones defined by the red line with circle markers in the second row of figure 1. Hence, it is interesting to consider a retrospective analysis where the estimated $x_{T, i}$ , $i \in C$ , is instead implemented in the dynamics with selective control. In figure 2, we present the results of such a study, where the evolution of the disease in the presence of selective control with target given by the uniform control is reported. It is evident how the amplitude of the epidemic peak is strongly reduced and the fraction of asymptomatic cases is under control in both cases using this alternative strategy. Furthermore, it is worth observing that a selective-type control is more effective in reducing the peak of infection in the presence of high heterogeneity like in the case of Italy.

Figure 2. Retrospective analysis of the epidemic dynamics in the presence of selective control strategy with target $x_{T, i}$ , $i \in C$ , estimated in the uniform case. (Online version in colour.)

4. Conclusion

The recent spreading of COVID-19 epidemic has been thwarted in most countries through non-pharmaceutical interventions, intended to diminish transmission rates and, to this end, reduce person-to-person contact via social distancing. Even if the confinement measures appear to be a correct tool to limit the infection spreading, a precise quantification of the effects played by them has still to be carefully understood.

Aiming to investigate the impact of the number of daily contacts in the spread of infectious diseases, in [18], we introduced a mathematical description of the epidemic spreading by integrating the epidemiological dynamics of the classical SIR-type compartmental model with a kinetic term quantifying the population-based contacts.

The kinetic description of the formation of social contacts considered in [18], based on elementary interactions, turns out to be very flexible with respect to external inputs, like a control forcing the number of daily social contacts towards a fixed target. From the modelling point of view, this strategy leads to integrating the compartmental models with social contacts with a new term quantifying in a precise way the effect of the control strategy depending on the population heterogeneity.

The main conclusion of our theoretical analysis, confirmed by numerical experiments, is that a selective control, when compared with a uniform control, and at the same social cost, leads to a marked reduction of the epidemic spreading. This reduction is further linked to the parameter of social heterogeneity, so that the selective control gives a better performance in a population with higher heterogeneity.

Data accessibility

All data used in this work are publicly available at https://github.com/CSSEGISandData/COVID-19.

Authors' contributions

G.D., G.T. and M.Z. conceived and wrote the paper; G.D. and M.Z. designed and performed the experiments. All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Competing interests

We declare we have no competing interests.

Funding

The research was partially supported by the Italian Ministry of Education, University and Research(MIUR: Project ‘Optimal mass transportation, geometrical and functional inequalities with applications’, and Dipartimenti di Eccellenza Program(2018–2022)—Dept. of Mathematics ‘F. Casorati’, University of Pavia. G.D. would like to thank the Italian Ministry of Instruction, University and Research (MIUR) for supporting this research with funds coming from PRIN Project 2017(no. 2017KKJP4X entitled ‘Innovative numerical methods for evolutionary partial differential equations and applications’).