key: cord-0111363-kypbb3br
authors: Colangeli, Matteo; Muntean, Adrian
title: Towards a quantitative reduction of the SIR epidemiological model
date: 2021-03-11
journal: nan
DOI: nan
sha: 9f2a8ea1c28784bb4abfd342e4fa74556e029abf
doc_id: 111363
cord_uid: kypbb3br

Motivated by our intention to use SIR-type epidemiological models in the context of dynamic networks as provided by large-scale highly interacting inhomogeneous human crowds, we investigate in this framework possibilities to reduce the classical SIR model to a representative evolution model for a suitably chosen observable. For selected scenarios, we provide practical {em a priori} error bounds between the approximate and the original observables. Finally, we illustrate numerically the behavior of the reduced models compared to the original ones.

The quest of a reduced description from a microscopic dynamics characterized by a large number of degrees of freedom is one of the classical problems of statistical and many-body physics, where model reduction and coarse-graining techniques proved to be a central tool underpinning renormalization group methods [24, 30] . Recently, model reduction techniques also found relevant applications in meteorology [16] and in physical and chemical kinetics [21] . One example is represented by the derivation of the hydrodynamic laws, described in terms of a restricted set of fields (e.g. density, momentum and temperature), from a kinetic description based on an extended set of moments or on the Boltzmann equation [10, 11, 25, 12] . With this occasion, we shall discuss the application of 2 Basic SIR model and its quantitative reduction We focus our attention on the structure of the celebrated SIR model. We refer the reader, for instance, to [7] for a nice description of the modelling ideas behind SIR as well as to [31, 32] for a number of qualitative properties of the solutions to SIR, SIRD, SEIR, and closely related models.

We consider three populations of individuals belonging to a larger population whose total number of individuals is fixed. We shall denote by S, I, and R the fraction of susceptible, infected, and removed individuals, respectively, such that S + I + R = 1.

The model equations entering the SIR model are:

where the initial conditions are prescribed as S(0) = S 0 , I(0) = I 0 , R(0) = R 0 such that S 0 + I 0 + R 0 = 1. Furthermore, b and γ are here strictly positive parameters that refer to an averaged infection rate constant and an averaged recovery rate constant, respectively. As a natural consequence, we see that if S 0 + I 0 + R 0 = 1 holds, then we have that also the mass conservation law S(t) + I(t) + R(t) = 1 holds for any t ∈ (0, T ), where T > 0 is arbitrarily fixed. We refer to the system of ODEs (1)-(3) as the original dynamics.

In this framework, we will offer a couple of other reduced variants of this SIR model. In all cases, we rely on the existence and uniqueness of classical positive solutions to the used models.

To obtain a reduced description from the original dynamics, we make the following ansatz: we assume that a suitable time scale exists, in which the time evolution of I(t) and R(t) is driven by the dynamics of the leading observable S(t).

Under this assumption, we identify the so-called driven observables. For the SIR model we have essentially two options: eitherŜ(t) andR(t) orÎ(t) and R(t), i.e. they correspond tô

The discussion of the reduction method based on Eq. (4) is discussed in Sec. 3, while the analysis of (5) is deferred to Sec. 4.3.

The time evolution of the observablesŜ(t),Î(t) andR(t) is dictated by the original dynamics, Eqs. (1)-(3), complemented by the ansatz (4). In particular, the dynamics ofÎ(t) reads

Equation (7) corresponds to the desired reduced description of the original SIR model, in which an expression for the constitutive law Φ[Î] is yet to be found. The dynamics of the driven observablesŜ(t) andR(t) is given by

The initial value problem for the system (6)-(8) is defined by fixing the valueŝ

We look for the exact expression of the functional Φ[Î], or, at least, a good approximate version thereof. To this aim, using the ansatz (4), we may also write the time derivative of the observableŜ(t) by relying on the chain rule, namely we have

with Φ [Î] := dΦ[Î]/dÎ, whereas the time derivative ofÎ is given by (6) .

The IM reduction method stipulates the equality of the two expressions of the time derivative ofŜ(t) given in Eqs. (7) and (9) . This procedure thus leads to the invariance equation [19, 20] , which reads

While Eq. (10) attains an exact, although not explicit, solution in terms of the Lambert W function [18] , we wish to follow here another route, and look for approximate, possibly explicit, solutions to Eq. (10). We can then integrate the latter by separation of variables, thus obtaining

We note in passing that, by proceeding in the same manner with the observ-ableR, we obtain the corresponding invariance equation

Using now (10), we can rewrite (12) in the form

which yields the expression

Finally, letting Φ 0 =Ŝ 0 and Ξ 0 =R 0 (15) and using (11) and (14), we obtain the consistency relation

which shows that the conservation of the total number of individuals in the population is inherited by the reduced description.

Next, in order to determine an explicit approximate expression of the functional Φ[Î], we seek approximate solutions of Eq. (11) obtained via an iteration method. As a possible choice of iteration method, we propose

For instance, by setting the initial condition for the recurrence equation (17) equal to Φ (0) [Î(t)] =Î(t), we find, at the the iteration level i = 1, the solution

A natural question thus arises: What do we learn from approximate solutions? It will turn out that the better we can approximate the exact solution of Eq. (11), given in terms of the Lambert W function, the better we can reduce the SIR model; see Claim 3 later on.

We now aim at estimating the nearness of the solution S(t) to Eq. (1) with initial datum S(0) = S 0 , and the constitutive law

It will turn out that our bound is meaningful only for a small time interval of observation and for a convenient parameter regime.

For an arbitrarily fixed value δ > 0 with t * ∈ (0, δ), we can derive for any t ∈ (0, δ) the next upper bound:

where c * > 0 is a constant depending on the parameters of the model, as well as on a priori uniform bounds on S, I and on the smoothness of Φ.

The last term comes form the Taylor expansion of Φ[I(t)], while the term involving the point evaluation in t * is the result of the application of the meanvalue theorem. Based on the rough estimate (19) , we observe that the quantity

can be made small if δ > 0 is sufficiently small, Φ is at least twice differentiable and S(t) and I(t) are bounded positive continuous functions. Mind, however, that the smallness of e(t) strongly depends on the choice of the parameters b and γ. For instance, in the limit of large values of γ, most likely the quantity e(t) will grow. On the other hand, if γ takes moderate values, then we expect e(t) ∼ O(δ) for sufficiently small t and e(t) ∼ O(δT ) for t ∈ (0, T ). (ii) We expect that instead of estimating from above the quantity e(t), it is more practical to bound the quantity

The dynamics of I(t) is governed by

where I(0) = I 0 . The starting point of this discussion is the fact that besides

we may also consider

where Φ (n) is the solution to our iterative method at the step n ∈ N. We provide also the information on the initial data I(0),Î(0), andĨ(0). Claim 1: The iteration method works such that for any r ∈ [0,

with lim n→∞ n = 0. Here ||·|| ∞ denotes the standard uniform norm on C[0, T ].

Proof of Claim 1. Previous work done in the existing literature on the rigorous numerical approximation of the Lambert function gives trust in this Claim; see, for instance, [18] and references cited therein. Claim 2: Under the assumptions for which Claim 1 holds, there exist constantsĉ 1 > 0 andĉ 2 > 0 such that

holds for any t ∈ (0, T ). Hereĉ 1 ,ĉ 2 are independent of t, T .

Proof of Claim 2. We observe firstly that

Subtracting (22) from (21) gives:

Noting that

Now, using the Grönwall's inequality (cf. e.g. Appendix B in [17] ) gives

Choosing nowĉ 1 := γ + ||Φ (n) || ∞ |(1 + ||I|| ∞ ) andĉ 2 := ||Φ || ∞ leads to the desired estimate proposed by Claim 2.

In this section, we aim to bound from above the error produced by the reduction method proposed within this framework. Claim 3: Assume the hypothesis of Claim 2 to be true. Let τ > 0 be arbitrarily fixed. Then there exist strictly positive constantsĉ 1 ,ĉ 2 , andĉ 3 such that the following a priori estimate holds:

where I andÎ satisfy (2) , and respectively (51).

Proof of Claim 3. We have

By a direct manipulation of the structure of the equations (2) and (51), we obtain:

As the first term on the right-hand side of the last inequality can be bounded above by (γ + ||S|| ∞ b) 1 β τ 0 |I(s)−Î(s)|ds (with 0 < β ≤ I(t)), and respectively, the last term by 2 b |I(τ ) −Î(τ )| + |I(0) −Î(0)| , the claim is now proven by choosing correspondingly the constantsĉ 1 ,ĉ 3 , andĉ 3 . 

This is obtained by adding and subtracting anĨ in each term on the righthand side of the estimate provided by Claim 3 and employing conveniently the statement of Claim 2. Note that the estimate (27) is quite practical. It basically tells that if I(0) =Î(0), then the quality of the reduction method depends mostly on the quality of the numerical approximation of the constitutive law Φ n .

The analysis of the SIR model, developed in Sec. 2, has shed light on the conditions under which one may hope to quantitatively capture the features of the original dynamics by using a reduced description based on the IM method.

We shall now turn our attention to another model, amenable to an analytical solution, which will also clarify the strengths and limitations of the IM method.

The Mori-Zwanzig method, in its essence, performs a partition of the dynamical variables into two subsets corresponding to the "relevant" and the "irrelevant" variables. Suitable projection operators are then employed to project the original dynamics onto the subspace of the relevant variables [29] . The simplest case in which this method can be discussed is a linear system of coupled first order ODEs:ẋ = L 11 x + L 12 y (28)

with x(0) = x 0 and y(0) = y 0 , and where L ij , i, j = 1, 2, are real parameters. Here x(t) and y(t) are regarded, respectively, as the "relevant" and the "irrelevant" variables. Using the set-up of Sec. 2, we may also regard y(t) as the dynamical variable whose time evolution is driven by x(t).

In this case, the original dynamics, given by Eqs. (28)- (29) , is amenable to an analytical solution. Namely, we can first solve (29) for y(t):

which represents the exact constitutive law linking y(t) to x(t). We can plug, next, (30) into (28) to obtain a closed ODE for x(t), which reads

Equation (31) represents the exact reduced description in terms of the relevant variable x(t). The price we paid to get a reduced description, in this example, amounts to the presence of a memory term in the evolution equation for x(t), which echoes the dynamics of the irrelevant variables.

We point out that, in the modelling of multiscale phenomena, the choice of the relevant dynamical variables is not always supported by guiding thermodynamic principles [22] . In fact, an improper choice of the relevant variables may not lead, eventually, to a successful reduced description [27] . On the other hand, a meaningful selection of the relevant variables proved to be extremely important, in statistical mechanics, to establish general results such as the Fluctuation-Dissipation Relations and the Fluctuation Relations in nonequilibrium systems [13, 14] .

We shall now discuss the application of the IM method to the system (28)- (29) . To get started, we focus on a perturbative method known as the Chapman-Enskog expansion, which stems from the geometrical theory of singular perturbations [23] . Namely, we introduce a singular perturbation into the equation of y(t) as follows:ẏ

where ε > 0 is a small parameter. Proceeding as in Sec. 2, we introduce the new variableŷ such thatŷ

wherex(t) is defined via its ODĖ

The dynamics ofŷ(t) can be described equivalently using the following two equations

We then impose the equality of two expressions of the time derivative ofŷ(t), (35) and (36), thus obtaining the invariance equation

In the Chapman-Enskog method, the solution of Eq. (37) is sought by expanding the variableŷ in a form of a series in powers of the number ε, i.e.

Thus, the Eq. (37) takes the form

which must be solved order by order by equating terms on both sides of the equation. At the lowest orders of ε, one finds the following sequence of constitutive laws order ε −1 : [20] . Therefore, we seek for a constitutive law (33) endowed with the linear structure

where A(L 11 , L 12 , L 21 , L 22 ) is an unknown function of the parameters L ij , i, j = 1, 2, yet to be determined. Hence (37) takes the form

Equation (45) becomes the quadratic equation

whose roots are

Thus, the IM method leads to the following constitutive laŵ which should be compared with the exact constitutive law (30) .

Using (48) we can now rewrite (34) as followṡ

where we have set A * = L 11 + L 12 A * . Figure 1 shows the behavior of x(t), obtained by integration of Eq. (31), and of the real part (x(t)) of the solution of the reduced system (49), with initial datum x 0 =x 0 = 1, for two different values of the coupling parameter L 21 . In particular, Fig. 1 evidences that fixing a larger value of L 21 , while keeping the other parameters fixed, leads to a better performance of the IM reduction method. The reason is that a larger value of L 21 , in (29) , makes the time derivative of y(t) more strongly affected by the behavior of x(t). This, hence, fits nicely with the ansatz (33), which requires the dynamics ofŷ(t) to be driven byx(t).

In this concluding section, we return to the SIR model of Sec. 2 and discuss the application of the IM method by using the ansatz (5) . It will turn out that, in this case, the IM method may yield an exact reduced description.

The dynamics of the fieldŜ(t) now reads 

Next, as in Sec. 3, we also write the time derivative ofÎ(t) by using the chain rule, i.e.

We thus obtain the invariance equation:

We can then integrate Eq. (53) by separation of variables, thus obtaining In a similar fashion, we find

Finally, by setting

and by summing up (54) and (55), we obtain

which yields, again, the conservation of the total number of individuals in the reduced dynamics. In Fig. 2 (50), equipped with the constitutive law (54). Figure 2 shows that the behavior ofŜ(t) recovers with striking accuracy that of S(t). Moreover, the two panels of Fig. 3 show the parametric plots of the I(t) vs. S(t) (left panel) and R(t) vs. S(t) (right panel) for the original SIR model and for the reduced description. We would like, however, to point out that the nice agreement between original and reduced dynamics outlined above might easily be lost when considering time dependent parameters b = b(t) and γ = γ(t), for t > 0. A careful rewriting of the invariance equation is demanded to handle such case. We will discuss this scenario elsewhere.

We have succeeded to identify constitutive laws to reduce either the presence of the population of susceptibles or of the infectives in the standard SIR model. The reduced descriptions, obtained using the IM method, agree via numerical simulations and practical a priori error bounds with what is expected from the original SIR dynamics.

Our work opens the possibility to use the reduced SIR dynamics for readingoff data available, for instance, on demonstrated Covid-19 infections and deaths and, based on a parameter identification approach done at this level, produce a new forecast on the effects of the pandemic evolution.

From a long-term research perspective, the method discussed in these notes indicates new routes to be exploited to obtain reduced descriptions in yet un-charted, or only partially explored, territories, such as the mathematical modelling of crowd dynamics [4, 5, 6, 15] and uphill diffusions [3, 8, 9] , in the framework of interacting particle systems.

Final size and convergence rate for an epidemic in heterogeneous populations

A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2

Stationary uphill currents in locally perturbed zero-range processes

Effects of communication efficiency and exit capacity on fundamental diagrams for pedestrian motion in an obscure tunnel -a particle system approach

A lattice model for active -Passive pedestrian dynamics: A quest for drafting effects

When diffusion faces drift: Consequences of exclusion processes for bi-directional pedestrian flows

Mathematical modeling of infectious diseases dynamics

Microscopic models for uphill diffusion

Nonequilibrium two-dimensional Ising model with stationary uphill diffusion

From hyperbolic regularization to exact hydrodynamics for linearized Grad's equations

Hyperbolicity of exact hydrodynamics for three-dimensional linearized Grad's equations

Boltzmann equation and hydrodynamic fluctuations

Fluctuation-dissipation relation for chaotic non-Hamiltonian systems

Equilibrium, fluctuation relations and transport for irreversible deterministic dynamics

Modelling interactions between active and passive agents moving through heterogeneous environments, Modeling and Simulation in Science, Engineering and Technology

Strategies for model reduction: Comparing different optimal bases

Partial Differential Equations

An algorithm for the matrix Lambert W function

Invariant Manifolds for Physical and Chemical Kinetics

Hilbert's 6th problem: Exact and approximate manifolds for kinetic equations

Model Reduction and Coarse-Graining Approaches for Multiscale Phenomena

Role of thermodynamics in multiscale physics

Geometric Singular Perturbation Theory, Dynamical systems

Scaling laws for Ising models near T c

Exact linear hydrodynamics from the Boltzmann equation

COVID-19 dynamics in an Ohio prison

Beyond Equilibrium Thermodynamics

A spatial Markov chain cellular automata model for the spread of viruses

Nonequilibrium Statistical Mechanics

Renormalization group and critical phenomena I. Renormalization group and the Kadanoff scaling picture

On build-up of epidemiologic models-Development of a SEI 3 RSD model for the spread of SARS-CoV-2, Zeitschrift für angewandte

Mathematische Bemerkungen zu einem epidemiologischen Modell im Zusammenhang mit der Corona-Pandemie, preprint