key: cord-0530653-95w3qul1
authors: Alutto, Martina; Como, Giacomo; Fagnani, Fabio
title: On SIR epidemic models with feedback-controlled interactions and network effects
date: 2021-05-04
journal: nan
DOI: nan
sha: 40ac3272ab355b9204e5fcac812ffec33b828f58
doc_id: 530653
cord_uid: 95w3qul1

We study extensions of the classical SIR model of epidemic spread. First, we consider a single population modified SIR epidemics model in which the contact rate is allowed to be an arbitrary function of the fraction of susceptible and infected individuals. This allows one to model either the reaction of individuals to the information about the spread of the disease or the result of government restriction measures, imposed to limit social interactions and contain contagion. We study the effect of both smooth dependancies of the contact rate for which we prove the existence of a threshold phenomenon that generalizes the well-known dichotomy associated to the reproduction rate parameter in the classical SIR model, and discontinuous feedback terms, which can be studied using tools from sliding mode control. Finally, we consider network SIR models involving different subpopulations that interact on a contact graph and present some preliminary simulations of modified versions of the classic SIR network.

As a result of the COVID-19 pandemic, there has been a renewed interest in the mathematical modeling, analysis, and control of epidemic spreadings. See, e.g., [1] - [6] . Of the two best known mathematical models of epidemics -the socalled susceptible-infective-succeptive (SIS) and susceptibleinfective-recovered (SIR) models-it is the latter that better approximates the spread of diseases like COVID-19 over a time horizon during which individuals tend not to get infected more than once, either because they have deceased or since they have recovered achieving some degree of immunity. The deterministic SIR epidemic model, as first presented in the pioneering work [7] , is a compartmental model consisting of a nonlinear system of three coupled differential equations describing the evolution of the fractions of susceptible, infected, and removed individuals in a fully mixed closed population.

The main feature of the deterministic SIR model is the existence of a phase transition described in terms of a scalar parameter, known as the reproduction number, whose value can determine two fundamentally different behaviors of the epidemics. Specifically, if the reproduction number does not exceed a unitary value, then the fraction of the infected individuals is bound to remain monotonically decreasing in time, and in fact asymptotically vanishing as time grows large, thus preventing an epidemic outbreak. In contrast, if the reproduction number exceeds the unitary threshold, then the fraction of infected individuals is initially increasing until reaching a peak, after which it starts to decrease monotonically as in the previous case and vanishes asymptotically in the large time limit. Thus, reproduction number values above one are equivalent to the occurrence of an epidemic outbreak. This phase transition is a crucial aspect of the epidemics and motivates the recent focus on a correct estimation of the reproduction number, as well the attempts to design control policies capable to stir the reproduction number below the unitary threshold level.

The original deterministic SIR model relies on the assumption that the rate at which individuals get infected is proportional to the product between the fraction of the susceptible individuals and the fraction of infected individuals. This is an appropriate model if we envision a fully mixed population with constant contact and transmission rates. There are many reasons for considering different interaction terms, as the result of either endogenous or exogenous interactions, e.g.:

• self-isolation policies or precautionary measures such as mask wearing put into place by the individuals who have become aware of the danger of the epidemic spreading; • feedback policies enforcing contact rate reduction by a central controller; • heterogeneities in the population and network effects. In this paper, we consider epidemic models that address the points listed above. Specifically, we study two extensions of the classical SIR model. First, we introduce a modified deterministic SIR model with a general interaction term describing the contact frequency rate as an arbitrary function of the fractions of susceptible and infected individuals in the population. Such function is typically decreasing in the fraction of infected individuals and can be interpreted as an endogenous reaction to the spread of epidemics (people tend to self-isolate), or rather an exogenous feedback control term modeling the action of a central planner actuating some partial social distancing measure. Then, we consider deterministic network SIR models where nodes represent subpopulations, derived by a split by either biological attributes (age, gender, risk) or geographical ones.

Our contribution is three-fold. First, for the extended scalar SIR model described above, we prove in Section II that as long as the contact rate is a smooth function of the state that is nondecreasing in the fraction of susceptible individuals, we retrieve the same threshold behavior of the SIR model: either the infection dies out monotonically, or it first increases, reaches a peak, and then decreases monotonically until vanishing asymptotically. Notice that in both cases the fraction of infected individuals remains a unimodal function of time, i.e., it adimts an unique local maximum that occurs at t = 0 when the the reproduction number does not exceed the unitary value and at t > 0 if the reproduction number value is larger than 1. Second, in Section III, we consider the generalized scalar SIR model where the contact rate is a discontinuous feedback term. For the especially important case of piecewise constant feedback controls, we show conditions under which a new sliding motion phenomenon can arise. In fact, we prove that in such models the fraction of infected individuals can remain constant at its maximum level, for a trivial interval of time before starting its monotone convergence to 0. Finally, in Section IV we study the deterministic network SIR model and show the effect of introducing two different types of control within the classical model. The former is a limitation of interactions between different nodes, while the latter is a general reduction of any kind of interaction.

We consider the following system of ODE's

where f : R 2 → R + is a state-dependent contact rate and γ > 0 is a constant recovery rate. Notice how, in the special case f (x, y) = β > 0, i.e., when the contact rate is a positive constant, (1) reduces to the classical deterministic SIR model

first introduced and studied in [7] .

For now we assume f to be a C 1 function. As it happens for the SIR model, the three equations are dependent, yielding that x(t) + y(t) + z(t) is constant. Throughout the paper, we assume this constant to be one, so to interpret x(t), y(t), and z(t) as fractions. The same considerations than in the SIR model, moreover, lead to the fact that solutions of (1) are globally defined (and unique) and that they keep invariant the simplex

. From now on we assume to always pick an initial condition that lays in S. Considering that x(t) and z(t) are monotone functions of time, the former decreasing and the latter increasing, we obtain that they always converge to a limit as t grows large. Invariance of the simplex S implies that the all solution vectors converges asymptotically and that in the limit y(t) → 0. Now, extending the contribution in [8] , we prove that, under mild assumptions on the function f , the threshold behavior of the SIR is retrieved in this model.

To this aim, define, for every solution (x(t), y(t), (t)) of (1), the function

We now show that R(t) plays the same role than the usual reproduction number for the SIR model (to which it reduces when f is constant). The following invariant result holds. Proposition 1: Assume that f is of class C 1 and is such that ∂f /∂x ≥ 0 in every point of S. For every solution of the ODE (1), if R(0) < 1, then R(t) < 1 for all t ≥ 0.

Proof: We first notice that if x(0) = 0, then x(t) = 0 at all time and consequently R(t) = 0 at all time. We now consider the case when x(0) > 0 that yields (by uniqueness of the solution) x(t) > 0 at all time. In this case, we prove the result by contradiction. If not, by continuity, there exists t * > 0 such that R(t * ) = 1 and R(t) < 1 for all t < t * . We now compute the time derivative of R(t):

We study the sign ofṘ(t * ). Because of the sign ofẋ, the standing fact that x > 0 at all time and the assumption on the function f , we have that the first two addends in the righthand side of (4) are always negative. Finally, the third term is 0 at t * sinceẏ(t * ) = 0. We conclude thatṘ(t * ) < 0. By continuity, we can state

and stays strictly below 1 for all t < t * for the assumption made, it follows that also R(t * ) < 1 contrarily to what we had assumed. This yields the result.

We can now state and prove the following result that shows how R(t) plays the same exact role than the reproduction number for the SIR model.

Theorem 1: Assume that f is of class C 1 and is such that ∂f /∂x ≥ 0 in every point of S. Given an initial condition (x 0 , y 0 , z 0 ) ∈ S and called the corresponding solution of (1) as (x(t), y(t), z(t)), the following facts hold (i) If R(0) < 1, then y(t) converges to 0 monotonically. (ii) If R(0) > 1 and y(0) > 0, then there exists t * > 0 such that

is monotonically decreasing in [t * , +∞[ and converges to 0 Proof: Concerning (i), it follows from Proposition 1 that R(t) < 1 for all t. This condition is equivalent to saying thatẏ(t) < 0 at all times t and proves (i).

Concerning (ii), notice that there must exist t > 0 such that R(t) < 1. Indeed, if not,ẏ(t) > 0 ∀t and y(t) would not possibly tend to 0. If we now define

we have that by construction the properties expressed in (ii) hold true. The proof is now complete.

Remark 1: In the case when f (x, y) = β is constant, we are back in the classical SIR model and Theorem 1 retrieves, in this case, the well known result on the behavior of curve of infected in this model. More details on the solutions can be obtained in this case. We briefly recall them below, as they will be needed in the next section. It will be convenient to set up the notation ρ = γ/β. From the first and third equation in (1), we notice that the function Γ(x, y, z) = ρ ln x + z = ρ ln x − x + 1 − y is motion invariant. In particular, given the initial condition x 0 = 1 − , y 0 = , and z 0 = 0, we obtain that the solution will lay in the manifold

When R(0) = (1 − )/ρ > 1, the maximum value reached by the infection happens in correspondence of x = ρ, as for this termẏ = 0, and is thus given by

In most cases of interest, the feedback term f is only function of y. Indeed, it is natural to imagine that a reaction both endogenous or exogenous be correlated to the extent of the current level of infection in the population. We notice that our result do not put any constraint on the way f may depend on y. Natural feedback terms will however be decreasing in y. In Figure 1 we compare the evolution of the classical SIR model in an unstable case with two modified versions having f (x, y) = h 2 (y) with respectively h(y) = 1 − y and h(y) = (1 − y) 2 .

III. SIR MODEL WITH THRESHOLD TERMS When we are modeling a control action of a central planner, it is of interest to study the case when f (y) has discontinuities. Indeed, it is not feasible to imagine a policy that varies with continuity, rather it is more natural a policy that changes when the infection reaches certain thresholds. In this section, we study in detail the case when f (y) is piecewise constant. In this case, the analysis carried on in previous section in general fails because of the discontinuities of the right hand side of the ODE (1). Classical solutions may not exist and, in this context, we will use the concept of solution according to Filippov [9] .

We consider the ODE (1) with an interaction term f (y) as defined below We interpret β as a sort of intrinsic interaction/contagion term that, in the absence of control measures, describes the rate at which infection propagates. When the infection gets above the threshold k, a (partial) lockdown policy takes place and brings this term to a smaller valueβ < β. In accord to Remark 1 we use the notation ρ = γ/β andρ = γ/β. To analyze this model, it is convenient to focus on the first two equations of (1). The right hand side is a discontinuous vector field that present a so called sliding manifold:

Indeed, in a sufficiently small neighborhood of Ω the vector field of the ODE points in the direction of the manifold. This is illustrated in Figure 2 . Trajectories in that region will eventually hit the manifold Ω and then they will remain on it sliding in the direction of decreasing x till the point (ρ,ȳ). From that point, the trajectory will remain in the region below (uncontrolled) and it will coincide with the trajectory of the uncontrolled SIR with interaction rate term equal to β.

In the following result we gather more detailed information on the nature of the sliding phenomenon and the conditions on the parameters for it to happen. To distinguish the two regimes determined by the value of the interaction term, we refer to, respectively, the β and theβ-SIR model.

Theorem 2: Assume that f is as defined in (7) and that the initial condition is of type (x(0), y(0), z(0)) = (1 − , , 0). Assume that (1 − )/ρ > 1 and put (c) If ≤ k ≤ m( , ρ,ρ), then there exists a time instant t * > 0 such that

• y(t) is monotonically decreasing in [t * , +∞[ and converges to 0. The maximum value in this third regime is given by the expression

where x(k) is the abscissa when the solution crosses the manifold y = k and is explicitly described by the relation

Proof: In the regime described in (a), the solution for the classical β-SIR model remains always in the region y < y. Consequently, it is also a solution of the controlled SIR model with control (7).

To study the other cases, we reduce to the xy plane and we only consider the first two equations in (1). We denote by (x(t), y(t)) the solution of the uncontrolled β-SIR model. We first note that, using the relation (5), whenρ < 1 − the expression m( , ρ,ρ) coincides with the value y of the solution corresponding to x =ρ. In other terms, there exists t 1 ≥ 0 such that x(t 1 ) =ρ and y(t 1 ) = m( , ρ,ρ).

If we are in the regime described by (b), we notice that for sure y(0) = < k so initially the solution leaves in the (uncontrolled) region y < k. Indicate by t * > 0 the first time when the solution (x(t), y(t)) hits the threshold level k. Such an instant must exist since the maximum value reached by y(t) is above k. Whenρ < 1 − notice that necessarily t 1 < t * as y(t 1 ) < k and y(t) is increasing till it reaches its maximum value. This implies that

When insteadρ > 1− , we have that x(t) <ρ at all times t so that (11) remains true. This says that, in any case, the solution hits, at time t * the sliding manifold Ω. Considering that the derivative of x is always negative, according to the definition of Filippov solution, the solution of the controlled SIR model from instant t * 1 on will be sliding on Ω till it reaches the point (ρ, k). This is reached at some time t * * > t * . From time t * * on the solution coincides again with the solution of the classical β-SIR model and will be decreasing in the component y and will converge to 0.

Consider now the regime (c). The only interesting case is whenρ < 1 − . Consider again (x(t), y(t)) the solution of the uncontrolled β-SIR model. By the considerations above, we have that at time t 1 > 0 when x(t 1 ) =ρ we have that y(t 1 ) > k. This implies that the solution (x(t), y(t)) has hit the line y = k at some previous time t * for which x(t * ) >ρ.

As (x(t * ), k) is out of the sliding manifold Ω, the solution of the controlled SIR-model will continue with a just a jump in the first derivative and since then it will coincide with the solution of an unstable classicalβ-SIR model. The component y will reach a peak for x =ρ and will then decrease and hit again the line y = k at some further time t 2 > t * . Depending on whether x(t 2 ) < ρ or x(t 2 ) > ρ the solution, in the first case, will undergo another jump in the first derivative and converge as a solution of a β-SIR model, while in the second case, will first slide along Ω to the point (ρ, k) and then will converge again as a solution of a β-SIR model.

Finally, the values for the maximum reached by the fraction of infected are simply obtained through the formula (6) that computes the maximum value of the infected in classical SIR models.

In the regime (c) of Theorem 2, it can happen that the solution, after reaching the peak, exhibits a sliding motion during the decreasing phase. We have not explicitly indicated this in the statement, as this phenomenon does not modify the maximum value reached globally by the component y(t) of the solution. It was however noticed in the proof. In Figure 3 we show all possible behaviors of the trajectory in the plane xy.

Notice, moreover, that the expression (9) for the maximum value reached by the fraction of infected is composed of two term. The first one is the value it would reach under the assumption that k = , namely that the controlled regime is active since the initial time. The true value y max is obtained from this adding a positive extra term that depends on k and accounts for the fact that for a while the epidemics has growth with no control. The estimation of this term can be relevant in the decision process of policy to adopt.

In this section, we analyse versions of the network SIR model.

Let G = (V, E, A) be a weighted di-graph with finite set of nodes V = {1, 2, . . . , n}, set of directed links E ⊆ V × V, and adjacency/weight matrix A in R n×n + . The different nodes i in V represent different subpopulations whereas the positive entries A ij > 0 of the weight matrix are in one-toone correspondence with the links (i, j) in E and measure the contact frequency of members of subpopulation i with members of subpopulation j. Throughout, we shall assume that the diagonal of A is strictly positive, i.e., that A ii > 0 for all i = 1, . . . , n.

For given infection rate β > 0 and recovery rate γ > 0, the network SIR epidemic model on a graph G = (V, E, A) is the dynamical system

for i = 1, . . . , n, where x i , y i , and z i represent respectively the fractions of susceptible, infected, and recovered individuals in population i. Notice that (12) may be more compactly rewritten aṡ

(13) This model has been studied in [10] and [11] . In particular, it is known that all solutions converge to an the equilibrium point of the form (x * , 0, z * ) in R 3n + such that x * + z * = 1 and that the locally asymptotically stable equilibrium points are those such that

where λ max (M ) stands for the dominant eigenvalue of a nonnegative matrix, which coincides with its spectral radius thanks to the Perron-Frobenius Theorem. In fact, under the assumption that the graph G is strongly connected, [10, Theorem 7] shows that the quantity

is decreasing along solutions and it plays a role similar to the one played by the reproduction number in the scalar SIR model. Specifically, if R(0) ≤ 1 then the weighted average v(0) y(t) will monotonically decrease to 0 as t grows large; on the other hand if R(0) > 1, then the weighted average v(0) y(t) will be initially increasing (epidemic outbreak) and there exists some τ > 0 such that R(τ ) ≤ 1 and the weighted average v(τ ) y(t) will be decreasing to 0 for t in the interval [τ, +∞).

Here v(t) stands for the nonnegative leading eigenvector of the matrix diag (x(t))A. Notice that the results summarized above concern the average behavior of the infection curve, with no implication on its behavior at individual nodes.

We now study two versions of the SIR network model, in which we introduce a social interaction mitigation function and, as in the previous section, we assume it is dependent only on the fraction of infected and decreasing with respect to it. This control function can be introduced as a modification of interactions between different nodes or any type of interaction, both inter-nodal and within the same subpopulation. If we consider nodes as geographically distinct subpopulations, the first type of contact limitation will lead to a kind of distancing and isolation per area with movements restriction. The interpretation of nodes in this model as a subdivision of the population into age groups will instead result in a limitation of interactions between people of different ages. This may be justified by an attempt to avoid contact between stronger people and people in age groups more vulnerable to disease, for example. We will then consider the following model

where the term f ij (t) concerns the limitation of contacts between the individuals of population i and those of population j. Obviously this term f ij (t) could depend on the infection level of both populations. For this reason, we can start by assuming that this term is equal to the product between the individual lockdown terms within the populations, that is f ij = f i f j ∀i, j = 1, ..., n where the individual lockdown measures considered are the following functions f i (t) = 1 − y i (t) ∀i = 1, ..., n A second alternative is to consider a control over all types of interaction within the network. In this case the studied model will be instead

where individual lockdown measures f i are assumed as in the previous case.

In Figure 4 we show simulations of the network SIR model with n = 2 subpopulations in the case of epidemic outbreak in both nodes and the two modified versions with the introduction of a internodal and a total control. Regarding the first node, these modifications of the model lead to an attenuation of the infection peak, while for the second node the introduction of the control causes the disappearance of an increasing trait for the curve of the infected. 

We have studied extensions of the classical Kermack and McKendrick's SIR epidemic model [7] that account for feedback-dependant contact rates and network effects. In particular, we have shown that discontinuous piecewise constant feedback rates may give rise to sliding motions, while for network models, simulations of possible modified versions are shown. Future research includes extension of these results, in particular for the network SIR model.

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

Controlling epidemic spread: reducing economic losses with targeted closures

A Simple Planning Problem for COVID-19 Lock-down, Testing, and Tracing

Optimal Targeted Lockdowns in a Multi-Group SIR model

Optimal epidemic suppression under an ICU constraint

Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models

A contribution to the mathematical theory of epidemics

A generalization of the Kermack-McKendrick deterministic epidemic model

Differential Equations with Discontinuous Righthand Sides. Mathematics and its Applications

On the dynamics of deterministic epidemic propagation over networks

Analysis and control of epidemics: A survey of spreading processes on complex networks