key: cord-0102222-hije3nbh
authors: Lazzizzera, Ignazio
title: An analytic approximate solution of the SIR model
date: 2020-11-15
journal: nan
DOI: nan
sha: f5bb831da7a70f29b6fe618feb2002968a7a96ed
doc_id: 102222
cord_uid: hije3nbh

The SIR(D) epidemiological model is defined through transcendental equations not solvable by elementary functions. In the present paper, those equations are successfully replaced by approximate ones, whose solutions are given explicitly in terms of elementary functions, namely, piece-wisely, generalized logistic functions: they unveil a useful feature, that in fact is also owned by the (numerical) solutions of the exact equations.

The SIR model [1] [2] [3] [4] [5] [6] is a simple compartmental model of infectious diseases developed by Kermack and McKendrick [1] in 1927. It considers three compartments: S, the set of susceptible individuals; I, the set of the infectious (or currently positive) individuals, who have been infected and are capable of infecting susceptible individuals; R, the set of the removed individuals, namely people who recovered (healed, H subset) from the disease or deceased (D subset), the former assumed to remain immune afterwards. The SIR model does not consider at all the sub-compartments H and D; instead the SIRD model simply assumes them to constitute a partition of R, fractionally fixed over time, so that, actually compared to the SIR model, nothing substantially changes in the dynamics of the epidemic progression. It is assumed that births and non-epidemic-related deaths can be neglected in the epidemic timescale and that the incubation period is negligible too. Indicating with letters not in bold the cardinality of each of the compartments, it is taken

where t 0 is an initial time; then the time evolution of the SIR(D) model is defined by the following system of non-linear first order differential equations:

Clearly

so that S(t) + I(t) + R(t) = S(t 0 ) + I(t 0 ) + R(t 0 ) = N .

Usually R(t 0 ) = 0. Both the parameters β and γ have dimension of a frequency. γ in eq.2c is the fractional removal rate (1/I)(dR/dt) of individuals from the infectious compartment. SI in eq.2a is understood as the number of possible contacts among the infectious and the susceptible individuals, so that β/N is the fractional decrease rate −(SI) −1 (dS/dt) of the number of individuals in the susceptible compartment; correspondingly β/N is the fractional increment rate of the number of infected individuals, determining the increment rate in the infectious compartment I , after subtraction of the rate of people entering the removed compartment R: this is in fact what eq.2b states. It is obvious that for the epidemic to spread, the increment rate of the newly infectious individuals must be higher then the increment rate of the newly removed individuals. Thus, introducing the so called basic reproduction ratio α (quite often denoted as R • ), namely

and dividing eq.2a by eq.2c , it must be

As a matter of fact, re-writing eq.2b as

condition 6 is seen to be equivalent to require the time derivative of I to be positive, or I itself increasing. By the way, since all the functions in the model are defined positive, from eq.2a it is seen that S is monotonic decreasing, tending to zero; then from eq.7 it follows that necessarily the number of infectious individuals must decrease after reaching a maximum and tend to zero even before the susceptible compartment may get empty. Taking eq.6 at the very beginning of the epidemic, when R = 0 and I S, thus S ≈ N , one understands that the meaning of α is the number of newly infectious individuals while just one infectious individual gets removed. Of course the average number of new infections from an infectious individual strongly influences the basic reproduction ratio, but it is not exactly the basic reproduction ratio itself. It is convenient to introduce the following non-dimensional variable and functions:

Then the basic equations of the SIR(D) model are written as

and

As is well known, the solutions of the equations of the SIR(D) model depend completely on the basic reproduction number (and the initial conditions), while of course β (indeed not γ !) gives the time scale. From eq.9a and then eq.9c one gets

Using this in 9b one easily finds then

Using this in 9c again, one gets

This is a transcendental equation, whose solutions one cannot give explicitly in closed analytic form by elementary functions. So in the sequel it is replaced by an approximate one, whose solutions are given explicitely. Of course, once r(x) is given, i(x) comes from eq.13 and s(x) from eq.12 .

The functions s(x), i(x) and r(x) are defined positive, so s(x) must be monotonic decreasing according to eq.9a and r(x) monotonic increasing according to eq.9c.

that is the condition for an epidemic to trigger, according to the short discussion above. Due to eq.9b , the function i(x) starts growing to a maximum which is reached at a time

then asymptotically it decreases to zero. Consequently, for eq.9c the bounded monotonically increasing function r(x) must exhibit a point of inflection at t M , after which it bends, increasing slower and slower, finally flattening to some limiting value

So one must have

thus getting a transcendental equation for r ∞ . Conveniently for the following developments, a new function is introduced, namely

in terms of which eq.14 is re-written as

must be solution of the equation

for eq.18 and the fact that

The functional F [w] is null in w = 1, but ∧ w cannot be 1 because 0 ≤ r(x) ≤ 1 and s 0 is not null (see eq.19 ); thus ∧ w must be solution of the equation

which is nothing but eq.18 , as can be easily verified. The ordinary eq.23 is transcendental and is to be solved numerically; the interval [0,

starts and remains negative from w = 0, until it reaches the point of inflection w flx , given by It is worth noting that as α increases, ∧ w (together with w flx ) approaches more and more the limiting value 1, namely a region where the log term in F [w] becomes important: this fact is relevant here because such log term, with its argument approaching zero, rises complications in searching for an effective approximation.

3. Approximating the key differential equation The idea is to approximate F [w] by few stretches of up to second order polynomials, joining continuously each other with the first derivative. Then in each stretch the obtained approximate differential equation becomes analytically and explicitly solvable by a generalized logistic function. For w 1 , it is taken [w] is a parabola with axis along the ordinate line, so that the maximum is its vertex. Denoting by w (1) 1 and w (1) 2 the roots of F (1) [w], one can write

The vertex is located in

A new parabola is chosen as the second approximation stretch, tangent to F [w] on its descending side, with axis along the ordinates and the vertex coincident with that of the first segment F (1) [w]:

Equations 31b and 31c impose that the two stretches have in common their vertexes, located in w = w M ; the system of the last two equations states the conditions for F (2) [w] to be tangent to F [w]. It is convenient expressing Z in terms of the unknown tangency point w using eq.31e , so that then one solves eq.31d for w . Introducing

due to eq.31c one can write

while from eq.31e and eq.31d one has

Using this expression for Z in eq.34a , one obtains a transcendental ordinary equation for w , to be solved numerically:

Using w so obtained, one gets Z from eq.34b and finally w (2) 1 and w (2) 2 via eq.33 and eq.31b . In fig.2 , on the left, the second segment for α = 2.6 is shown in blue, extending from w M to the point of tangency of the successive approximation segment still to be chosen. With reference to the discussion at the end of Section 2 , it should be noted that F [w] remains concave up to w = ∧ w when α ≤ α cr , while it happens that the root w (2) 1 of F (2) [w] (see fig.3 ) remains very close to ∧ w : this suggests in that range of α values replacing the above F (2) [w] by a different arc of parabola f (2) [w] , keeping its vertex in common with F (1) [w] as F (2) [w] does, but just ending in 

For α cr < α ≤ 6 F [w] is almost always concave, ending roughly as a straight line when approaching ∧ w. In this range of α's one keeps F (2) [w] and completes the approximation through a new parabola, requiring it to be tangent to F (2) [w] and to reach ∧ w along the tangent to F [w] in ∧ w; an alternative is the ray tangent to F (2) [w], extending from the point of tangency to ∧ w . The latter is settled by

[w] is the discriminant of the second order algebraic equa-

[w] = 0 , set to zero to assure L[w] to be tangent to F (2) [w] . The appropriate solution for u is

The problem with this approximation is that, looking for instance at the function r(x) obtained from w(x), it gets unacceptably overestimated in the region where it bends to reach the symptotic value as x → +∞: this is because L[w] necessarily remains below F [w] due to the concavity of the latter. The quadratic alternative is defined by

[w] − F (2) [w] = 0 ∧ ∆ F

[w] − F (2) [w] = 0 ,

where "prime" stands for derivative and ∆ F

[w] − F (2) [w] is the discriminant of the second order algebraic equation F

[w] − F (2) [w] = 0 , set to zero so to assure F (3) [w] to be tangent to F (2) [w] . In this case, however, with respect to using L[w] , one has the opposite effect on r(x), because the given choice for λ forces F (3) [w] to stay somewhat above F [w]. The solution is to keep the quadratic alternative, but replacing the previous value of λ by a compromise one, defined through

Then the parameter σ in 39a is set by means of the the condition 39c :

where w • is the tangency point of F (3) [w] to F (2) [w]. So, for α cr < α <= 6 the third and last approximation segment is given by 39a , with λ repalced by λ • , extending from w • to ∧ w. For w > 6 the convexity trait of F [w], following the almost straight stretch around w flx , gets more and more included in the domain [0, ∧ w] , because ∧ w increases with α. Then, the solution adopted is to introduce a linear segment T [w] parallel to the tangent in w flx to F [w] and tangent to F (2) [w] in a point that will be denotedw; this linear segment will be continued by a new parabola F (4) [w], which is similar to F (3) [w], thus ending in ∧ w, but tangent to T [w]. Namely

Then the already mentioned F (4) [w] approximation streatch is constrained to end in ∧ w and to be tangent to T [w] in a point w u chosen by trial and error optimization:

giving

(45b)

For each of the above approximation segments, a differential equations is defined of the type dw dx

where F[w] is one of F (i)

[w] (i = 1, 2, 3, 4) or f (2) [w] or T [w], with given α and β parameters (or β and γ) and initial conditions. For F[w] = F (1) [w] , from the definition in eq.19 the initial condition is w(x 0 ) = 1 − s 0 = i 0 (x 0 = 0 without loss of generality), while for each of the remaining approximation segments it is given by the value of the respective preceding segment at the junction point. Since F[w] is at most a second order polynomial, eq.46 is indeed quite trivially solved, giving a generalized logistic function.

For F[w] = F (1) [w]:

[w] , thus α ≤ α cr :

[w] , thus α > α cr :

[w] , thus α cr < α ≤ 6 :

, thus α > 6 (see 42 and 43):

[w] thus α > 6 :

Then, from eq.19 one has ∨ r(t) = 1 α ln

so that

On the other hand eq.20 implies

Finally, of course, due to 10,: In the case of the SIRD model one defines

(57b) can be summarized as follows:

The equation of the first approximation segment can be re-written as 

Using the explicit solution eq.47 , one has

and consequently

Typically

but anyway with t − t 0 greater then some τ 1 's, in the end one can write

Analogous results hold for all the approximation stretches in the different α intervals as summarized in eq.s58; for instance, with t − t • greater enough then τ 3 , one has

These piecewise linear behaviors can be seen in fig.5 for α = 2.6 . The plot on the left shows the numerical solution of the exact equation, compared with the corresponding approximate analytic solution: it is worth recalling (see eq.55 ) that w(x) = r(x) + i(x), so that w is directly related to the data. The plot on the right shows that the function ∨ r(t) of the removed individuals exhibits an analogous behavior: since in the SIRD model the spring 2020 first wave of Covid-19 in Italy: it remarkably confirms this model feature. One important point here is that the slopes of the straight segments, that are inversely proportional to the related time constants τ , are completely determined by the parameters α and β (besides the initial conditions) and so is the angle between such straight segments: consequently one can compare that angle with the theoretically predicted one and argue about the effects of social measures to reduce the pandemic, of course within the trustworthiness of the model.

In this paper the equations of the SIR(D) epidemiological model are replaced by approximate ones, whose solutions are totally defined by the basic reproduction ratio α and the fractional removal rate γ of individuals from the infectious compartment (alternatively by β = γ/α). These solutions are chains of two or three or four generalized logistic functions, depending on the value of α only; they are summarized in eq.s 58 .

In practice, to get them one does:

• solve numerically the transcendental ordinary eq.23 to get ∧ w;

• use eq.29 and eq.s47 to get w (1) (x) as in eq.47;

• for α ≤ α cr use eq.30 , eq.36d and 36e to get w (f ) (x) as in eq.48 ;

• for α > α cr use eq.35 , eq.34b , eq.30 , eq.32b , eq.33 and eq.s49 to get w (2) (x) as in eq.49 ;

• for α cr < α <= 6 use eq.40 , eq.41a and eq.41b , eq.41c and finally eq.s50 to get w (3) (x) as in eq.50 ;

• for α > 6 use eq.42b , eq.s43 and eq.s51 to get w (T ) (x) as in eq.51 ;

• for α > 6 use eq.25 , eq.s44b , eq.s45 and eq.s52 to get w (4) (x) as in eq.52 ;

• eventually use eq.54 , eq.55 , eq.56 , eq.57 .

Having such explicit solution would help, for instance, to study the data through easy fits.

Contribution to the mathematical theory of epidemics

Epidemic Modelling

Mathematical epidemiology: Past, present, and future