An Analytic Approximate Solution of the SIR Model

I. Lazzizzera^1,2

¹Department of Physics, University of Trento, via Sommarive 14 Povo (Trento), Italy.
²Trento Institute for Fundamental Physics and Applications (INFN), via Sommarive 14 Povo (Trento), Italy.

DOI: 10.4236/am.2021.121005                

Abstract

The SIR(D) epidemiological model is defined through a system of 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, originating, piece-wisely, from generalized logistic functions: they ensure exact (in the numerical sense) asymptotic values, besides to be quite practical to use, for example with fit to data algorithms; moreover they unveil a useful feature, that in fact, at least with very strict approximation, is also owned by the (numerical) solutions of the exact equations. The novelties in the work are: the way the approximate equations are obtained, using simple, analytic geometry considerations; the easy and practical formulation of the final approximate solutions; the mentioned useful feature, never disclosed before. The work’s method and result prove to be robust over a range of values of the well known non-dimensional parameter called basic reproduction ratio, that covers at least all the known epidemic cases, from influenza to measles: this is a point which doesn’t appear much discussed in analogous works.

Keywords

SIR Epidemic Model, Kermack-McKendrick Model, Epidemic Dynamics, Approximate Analytic Solution

1. Introduction

The SIR model [1] - [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 due to the disease (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

$S\left(t_0\right)+I\left(t_0\right)+R\left(t_0\right)=N\mathrm{,}$ (1)

where $t_0$ is an initial time, usually with $R\left(t_0\right)=0$.

The model introduces two parameters, β and γ, having dimension of a frequency. Saying t the time variable, γ is defined as the fractional removal rate $\left(1/I\right)\left(\text{d}R/\text{d}t\right)$ of individuals from the infectious compartment. Since SI is understood as the number of possible contacts among the infectious and the susceptible individuals, β/N is defined as the fractional decrease rate $-{\left(SI\right)}^{-1}\left(\text{d}S/\text{d}t\right)$ of the number of individuals in the susceptible compartment: it expresses therefore the fractional increment rate of the number of infectious individuals, that is the increment rate of the infectious compartment I, after subtraction of the rate of people entering the removed compartment R.

Usually one introduces the following non-dimensional variable and new functions:

$\alpha :=\frac{\beta }{\gamma },x:=\gamma t,s\left(x\right):=\frac{S\left(t\right)}{N},i\left(x\right):=\frac{I\left(t\right)}{N},r\left(x\right):=\frac{R\left(t\right)}{N},$ (2)

$\alpha$ called basic reproduction ratio. Then the basic equations given by Kermack and McKendrick [1] are written as

$\frac{\text{d}s}{\text{d}x}\left(x\right)=-\alpha i\left(x\right)s\left(x\right)$ (3a)

$\frac{\text{d}i}{\text{d}x}\left(x\right)=i\left(x\right)\left(\alpha s\left(x\right)-1\right)$ (3b)

$\frac{\text{d}r}{\text{d}x}\left(x\right)=i\left(x\right)$ (3c)

with

$s\left(x\right)+i\left(x\right)+r\left(x\right)=s\left(x_0\right)+i\left(x_0\right)+r\left(x_0\right)=\mathrm{1,}$ (4)

and

$s_0:=s\left(x_0\right),i_0:=i\left(x_0\right),r_0:=r\left(x_0\right)\equiv 0.$ (5)

Using Equation (3c) in Equation (3a) and formally integrating, one gets $s\left(x\right)=s_0{\text{e}}^{-\alpha r\left(x\right)}$ ; using this and Equation (3c) again, from Equation (3b) one easily finds $i\left(x\right)=1-s_0{\text{e}}^{-\alpha r\left(x\right)}-r\left(x\right)$ ; then from Equation (3c) she/he will obtain

$\frac{\text{d}r}{\text{d}x}\left(x\right)=1-s_0{\text{e}}^{-\alpha r\left(x\right)}-r\left(x\right)\mathrm{.}$ (6)

This is a transcendental equation, whose solutions one cannot give explicitly in closed analytic form by elementary functions. In their original paper Kermack and McKendrick themselves [1] gave approximate solutions, however without any exhaustive discussion of applicability for various values of the basic reproduction ratio. Quite recently various authors have approached the problem in different ways, but with the same incompleteness [7] [8] [9] [10] [11]. In the sequel, on the basis of simple, analytic geometry considerations, a novel method is introduced, producing approximate but accurate solutions, given explicitly, piece-wisely, from generalized logistic function (see [12] for a description of the origin of the logistic function and its adoption in bio-assay); due attention is paid for the method to be robust over the whole range of possible known values of $\alpha$, from just above 1 as for influenza, to 1.4 - 3.9 as for Covid-19, to 3 - 5 as for SARS, to 5 - 7 as for polio, to 10 - 12 as for varicella, to 12 - 18 as for measles (see for instance [13] and references therein).

2. Getting the Key Differential Equation

For the epidemic to spread, the increment rate of the newly infectious individuals must be higher than the increment rate of the newly removed individuals. Dividing Equation (3a) by Equation (3c), one finds that it must be

$1<-\frac{\text{d}s}{\text{d}r}\left(t\right)=\alpha s\left(t\right).$ (7)

As a matter of fact this condition implies that $i\left(t\right)$ increases over time due to Equation (3b). The functions $s\left(x\right)$, $i\left(x\right)$ and $r\left(x\right)$ are all defined positive and less or equal to 1; consequently it must be $\alpha >1$ for the epidemic to spread and $s\left(x\right)$ turns to be monotonic decreasing according to Equation (3a), while $r\left(x\right)$ monotonic increasing according to Equation (3c). It follows that the function $i\left(x\right)$ starts growing due to (7), reaching necessarily a maximum at a time $t_M=x_M/\gamma$ such that

$\alpha s\left(x_M\right)=\mathrm{1,}$ (8)

then asymptotically decreasing to zero. This implies that the bounded monotonically increasing function $r\left(x\right)$ must exhibit a point of inflection at $t_M$, after which it bends, increasing slower and slower, finally flattening to some limiting value

$r_{\infty }\equiv r\left(+\infty \right)\le 1.$ (9)

So one must have

$0=\underset{x\to +\infty }{\lim }\frac{\text{d}r}{\text{d}x}\left(x\right)=1-s_0{\text{e}}^{-\alpha r_{\infty }}-r_{\infty },$ (10)

thus getting a transcendental equation for $r_{\text{}\infty }$.

Conveniently for the following developments, a new function is introduced, namely

$w\left(x\right)=1-s_0{\text{e}}^{-\alpha r\left(x\right)}\mathrm{,}$ (11)

in terms of which Equation (6) is re-written as

$\frac{\text{d}w}{\text{d}x}=F\left[w\right]\mathrm{,}$ (12a)

$F\left[w\right]:=\left(1-w\right)\left[ϵ+\alpha w+\ln \left(1-w\right)\right],$ (12b)

$ϵ=-\ln \left(s_0\right)=-\ln \left(1-i_0\right)\mathrm{.}$ (12c)

Clearly

$\hat{w}:=\underset{x\to +\infty }{\lim }w\left(x\right)=1-s_0{\text{e}}^{-\alpha r_{\infty }}=r_{\infty }$ (13)

must be solution of the equation

$F\left[\hat{w}\right]=0,$ (14)

for Equation (10) and the fact that

$\frac{\text{d}w}{\text{d}x}=s_0\alpha {\text{e}}^{-\alpha r\left(x\right)}\frac{\text{d}r}{\text{d}x},$

so that

$\frac{\text{d}w}{\text{d}x}=0\iff \frac{\text{d}r}{\text{d}x}=0.$

The functional $F\left[w\right]$ is null in $w=1$, but $\hat{w}$ cannot be 1 because $0\le r\left(x\right)\le 1$ and $s_0$ is not null (see Equation (11)); thus $\hat{w}$ must be solution of the equation

$ϵ+\alpha \hat{w}+\ln \left(1-\hat{w}\right)=0,$ (15)

which is nothing but Equation (10), as can be easily verified. Equation (15) is transcendental and is to be solved numerically; the interval $\left[\mathrm{0,}\hat{w}\right]$ is the range of $w\left(x\right)$ as x runs from $x_0$ to $+\infty$.

The second derivative of F, namely

$\frac{{\text{d}}^{2}F}{\text{d}{w}^2}\left[w\right]=-2\alpha +\frac{1}{1-w}$ (16)

starts and remains negative from $w=0$, until it reaches the point of inflection $w_{f\text{}l\text{}x}$, given by

$w_{f\text{}l\text{}x}=1-\frac{1}{2\alpha };$ (17)

then it becomes positive: thus $F\left[w\right]$ starts and remains concave until $w=w_{f\text{}l\text{}x}$ ; then it becomes convex. Of course, in an interval around its inflection point, $F\left[w\right]$ is nearly straight. Figure 1 shows how $\hat{w}$ and $w_{f\text{}l\text{}x}$ vary as a function of $\alpha$ : for $\alpha <{\alpha }_{\text{cr}}\simeq 1.75$ one has $\hat{w}<w_{f\text{}l\text{}x}$ and consequently $F\left[w\right]$ is always concave in the domain $\left[\mathrm{0,}\hat{w}\right]$ ; otherwise it changes from concave to convex

Figure 1. Point of inection and $\hat{w}$ as a function of $\alpha$.

after $w=w_{f\text{}l\text{}x}$. It is worth noting that as $\alpha$ increases, $\hat{w}$ (together with $w_{f\text{}l\text{}x}$ ) approaches more and more the limiting value 1, namely a region where the log term in $F\left[w\right]$ 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\left[w\right]$ 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\ll 1$, it is taken

$\left(1-w\right)\ln \left(1-w\right)\approx -w\left(1-\frac{1}{2}w\right)\mathrm{,}$ (18)

so that

$\frac{\text{d}w}{\text{d}x}\approx ϵ+\left(\alpha -1-ϵ\right)w-\left(\alpha -\frac{1}{2}\right)w^2=F^{\left(1\right)}\left[w\right]\mathrm{.}$ (19)

Figure 2 shows on the left, in red, this $F^{\left(1\right)}\left[w\right]$ segment against $F\left[w\right]$ (black curve) for $\alpha =2.74$ and (consequently) $\hat{w}\simeq 0.92$, extending to its maximum point, which is rather close to the maximum of $F\left[w\right]$. Clearly $F^{\left(1\right)}\left[w\right]$ is a parabola with axis along the ordinate line, so that the maximum is its vertex.

Denoting by $w_1^{\left(1\right)}$ and $w_2^{\left(1\right)}$ the roots of $F^{\left(1\right)}\left[w\right]$, one can write

$F^{\left(1\right)}\left[w\right]=-A\left(w-w_1^{\left(1\right)}\right)\left(w-w_2^{\left(1\right)}\right)\mathrm{,}$ (20a)

$A:=\alpha -\frac{1}{2}\mathrm{,}$ (20b)

with

$w_{1/2}^{\left(1\right)}=\frac{\alpha -1-ϵ\pm \sqrt{{\left(\alpha -1-ϵ\right)}^2+2\left(2\alpha -1\right)ϵ}}{2\alpha -1}\mathrm{.}$ (21)

The vertex is located in

$w_M=\frac{w_1^{\left(1\right)}+w_2^{\left(1\right)}}{2}\mathrm{.}$ (22)

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

$F^{\left(2\right)}\left[w\right]=-Z^{\star }\left(w-w_1^{\left(2\right)}\right)\left(w-w_2^{\left(2\right)}\right),$ (23a)

$\frac{w_1^{\left(1\right)}+w_2^{\left(1\right)}}{2}=w_M=\frac{w_1^{\left(2\right)}+w_2^{\left(2\right)}}{2},$ (23b)

$-A\left(w_M-w_1^{\left(1\right)}\right)\left(w_M-w_2^{\left(1\right)}\right)=-Z^{\star }\left(w_M-w_1^{\left(2\right)}\right)\left(w_M-w_2^{\left(2\right)}\right),$ (23c)

$F^{\left(2\right)}\left[w\right]=F\left[w\right],$ (23d)

$\frac{\delta F^{\left(2\right)}}{\delta w}\left[w\left(x\right)\right]=\frac{\delta F}{\delta w}\left[w\left(x\right)\right].$ (23e)

Equations (23b) and (23c) 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^{\left(2\right)}\left[w\right]$ to be tangent to $F\left[w\right]$. It is convenient expressing $Z^{\star }$, appearing in Equation (23c), in terms of the unknown tangency point $w^{\star }$ using Equation (23e), so that consequently one solves Equation (23d) for $w^{\star }$.

Namely, introducing

$\delta w^{\left(1\right)}:=\frac{w_1^{\left(1\right)}-w_2^{\left(1\right)}}{2}\mathrm{,}$ (24a)

$\delta w^{\left(2\right)}:=\frac{w_1^{\left(2\right)}-w_2^{\left(2\right)}}{2}\mathrm{,}$ (24b)

due to Equation (23c) one can write

$Z^{\star }{\left(\delta w^{\left(2\right)}\right)}^2=A{\left(\delta w^{\left(1\right)}\right)}^2\mathrm{,}$ (25)

while from Equation (23e) and Equation (23d) one has

$\left(1-w^{\star }\right)\left[ϵ+\alpha w^{\star }+\ln \left(1-w\right)\right]=-Z^{\star }{\left(w-w_M\right)}^2+A{\left(\delta w^{\left(1\right)}\right)}^2\mathrm{,}$ (26a)

$Z^{\star }=\frac{1+ϵ+2\alpha w^{\star }+\ln \left(1-w^{\star }\right)-\alpha }{2\left(w^{\star }-w_M\right)}\mathrm{.}$ (26b)

Using this expression for $Z^{\star }$ in Equation (26a), one obtains a transcendental ordinary equation for $w^{\star }$, to be solved numerically:

$\begin{array}{l}2ϵ+\left(\alpha -ϵ-1\right)w_M-2A{\left(\delta w^{\left(1\right)}\right)}^2+\left(\alpha -ϵ-2\alpha w_M+1\right)w^{\star }\\ +\left(2-w^{\star }-w_M\right)\ln \left(1-w^{\star }\right)=0.\end{array}$ (27)

Using $w^{\star }$ so obtained, one gets $Z^{\star }$ from Equation (26b) and finally $w_1^{\left(2\right)}$ and $w_2^{\left(2\right)}$ via Equation (25) and Equation (23b). In Figure 2, on the left, the second segment for $\alpha =2.6$ is shown in blue, extending from $w_M$ to the point of tangency of the successive approximation segment still to be chosen.

Figure 2. Examples of two cases, with three approximation stretches on the left (red, blue, green) and four approximation stretches on the right (red, blue, green, magenta).

With reference to the discussion before the end of Section 2, it should be noted that $F\left[w\right]$ remains concave up to $w=\hat{w}$ when $\alpha \le {\alpha }_{\text{cr}}$, while it happens that the root $w_1^{\left(2\right)}$ of $F^{\left(2\right)}\left[w\right]$ (see Figure 3) remains very close to $\hat{w}$ : this suggests in that range of $\alpha$ values replacing the above $F^{\left(2\right)}\left[w\right]$ by a different arc of parabola $f^{\left(2\right)}\left[w\right]$, keeping its vertex in common with $F^{\left(1\right)}\left[w\right]$ as $F^{\left(2\right)}\left[w\right]$ does, but just ending in $\hat{w}$, thus imposing the constraint $w_1^{\left(2\right)}=\hat{w}$ instead of the tangency to $F\left[w\right]$.

Then for $\alpha \le {\alpha }_{\text{cr}}$

$f^{\left(2\right)}\left[w\right]=-Z\left(w-w_1^{\left(f\right)}\right)\left(w-w_2^{\left(f\right)}\right),$ (28a)

$\frac{w_1^{\left(1\right)}+w_2^{\left(1\right)}}{2}=w_M=\frac{w_1^{\left(f\right)}+w_2^{\left(f\right)}}{2},$ (28b)

$-A\left(w_M-w_1^{\left(1\right)}\right)\left(w_M-w_2^{\left(1\right)}\right)=-Z\left(w_M-w_1^{\left(f\right)}\right)\left(w_M-w_2^{\left(f\right)}\right),$ (28c)

$Z=\left(\alpha -\frac{1}{2}\right)\frac{{\left(\delta w^{\left(1\right)}\right)}^2}{{\left(\hat{w}-w_M\right)}^2},$ (28d)

$w_1^{\left(f\right)}=\hat{w},w_2^{\left(f\right)}=2w_M-\hat{w},\delta w^{\left(f\right)}=\hat{w}-w_M.$ (28e)

For ${\alpha }_{\text{cr}}<\alpha \le 6$ $F\left[w\right]$ is almost always concave, ending roughly as a straight line when approaching $\hat{w}$. In this range of α’s one keeps $F^{\left(2\right)}\left[w\right]$, but completes the approximation through a new parabola, requiring it to be tangent to $F^{\left(2\right)}\left[w\right]$ and to reach $\hat{w}$ along the tangent to $F\left[w\right]$ in $\hat{w}$ ; an alternative is the ray originating in $\hat{w}$, tangent to $F^{\left(2\right)}\left[w\right]$. The latter is settled by

$L\left[w\right]:=-2u{Z}^{\star }\left(w-\hat{w}\right)\mathrm{,}$ (29a)

$∍\text{'}\left\{L\left[w\right]-F^{\left(2\right)}\left[w\right]=0\wedge \Delta \left(L\left[w\right]-F^{\left(2\right)}\left[w\right]\right)=0\right\}\mathrm{,}$ (29b)

Figure 3. $w_1^{\left(2\right)}$ and $\hat{w}$ as functions of $\alpha$.

where $\Delta \left(L\left[w\right]-F^{\left(2\right)}\left[w\right]\right)$ is the discriminant of the second order algebraic equation $L\left[w\right]-F^{\left(2\right)}\left[w\right]=0$, set to zero to assure $L\left[w\right]$ to be tangent to $F^{\left(2\right)}\left[w\right]$. The appropriate solution for u is

$u_{-}=\hat{w}-w_M-\sqrt{{\left(\hat{w}-w_M\right)}^2-{\left(\delta w^{\left(2\right)}\right)}^2}\mathrm{.}$ (30)

The problem with this approximation is that, looking for instance at the function $r\left(x\right)$ obtained from $w\left(x\right)$, it gets unacceptably overestimated in the region where it bends to reach the asymptotic value as $x\to +\infty$ : this is because $L\left[w\right]$ necessarily remains below $F\left[w\right]$ due to the concavity of the latter.

The quadratic alternative is defined by

$F^{\left(3\right)}\left[w\right]:=-2\lambda \left(w-\hat{w}\right)+\sigma {\left(w-\hat{w}\right)}^2\mathrm{,}$ (31a)

$\lambda ={{\left(F^{\left(2\right)}\left[w\right]\right)}^{\prime }|}_{w=\hat{w}}=\frac{1-\alpha \left(1-\hat{w}\right)}{2},$ (31b)

$∍\text{'}\left\{F^{\left(3\right)}\left[w\right]-F^{\left(2\right)}\left[w\right]=0\wedge \Delta \left(F^{\left(3\right)}\left[w\right]-F^{\left(2\right)}\left[w\right]\right)=0\right\}\mathrm{,}$ (31c)

where “prime” stands for derivative and $\Delta \left(F^{\left(3\right)}\left[w\right]-F^{\left(2\right)}\left[w\right]\right)$ is the discriminant of the second order algebraic equation $F^{\left(3\right)}\left[w\right]-F^{\left(2\right)}\left[w\right]=0$, set to zero so to assure $F^{\left(3\right)}\left[w\right]$ to be tangent to $F^{\left(2\right)}\left[w\right]$. In this case, however, with respect to using $L\left[w\right]$, one has the opposite effect on $r\left(x\right)$, because the given choice for $\lambda$ forces $F^{\left(3\right)}\left[w\right]$ to stay somewhat above $F\left[w\right]$.

The solution is to keep the quadratic alternative, but replacing the previous value of $\lambda$ by a compromise one, defined through

${\lambda }^{\circ }:=\tan \left(\arctan \left(-2\lambda \right)\right)+\tan \left(\frac{\arctan \left(-2\lambda \right)-\arctan \left(-2u_{-}{Z}^{\star }\right)}{2}\right).$ (32)

Then the parameter $\sigma$ in (31a) is set by means of the condition (31c):

$\sigma =\frac{Z^{\star }h-g^2}{2\hat{w}g-h-Z^{\star }{\hat{w}}^2}\mathrm{,}$ (33a)

$g=Z^{\star }{w}_M+{\lambda }^{\circ },h=Z^{\star }{w}_1^{\left(2\right)}{w}_2^{\left(2\right)}+2{\lambda }^{\circ }\hat{w},$ (33b)

$w^{\circ }=\frac{\sigma \hat{w}+g}{\sigma +Z^{\star }}\mathrm{,}$ (33c)

where $w^{\circ }$ is the tangency point of $F^{\left(3\right)}\left[w\right]$ to $F^{\left(2\right)}\left[w\right]$.

So, for ${\alpha }_{\text{cr}}<\alpha \le 6$ the third and last approximation segment is given by (31a), with $\lambda$ replaced by ${\lambda }^{\circ }$, extending from $w^{\circ }$ to $\hat{w}$.

For $w>6$ the convexity trait of $F\left[w\right]$, following the almost straight stretch around $w_{flx}$, gets more and more included in the domain $\left[\mathrm{0,}\hat{w}\right]$, because $\hat{w}$ increases with $\alpha$. Then, the solution adopted is to introduce a linear segment $T\left[w\right]$ parallel to the tangent in $w_{flx}$ to $F\left[w\right]$ and tangent to $F^{\left(2\right)}\left[w\right]$ in a point that will be denoted $\tilde{w}$ ; this linear segment will be continued by a new parabola $F^{\left(4\right)}\left[w\right]$, which is similar to $F^{\left(3\right)}\left[w\right]$, thus ending in $\hat{w}$, but tangent to $T\left[w\right]$. Namely

$T\left[w\right]:=-2\tilde{f}w+\tilde{I}\mathrm{,}$ (34a)

$-2\tilde{f}:={F^{\prime }\left[w\right]|}_{w=w_{flx}}=\ln \left(2\alpha \right)-\alpha -ϵ\mathrm{,}$ (34b)

$∍\text{'}\left\{T\left[w\right]-F^{\left(2\right)}\left[w\right]=0\wedge \Delta \left(T\left[w\right]-F^{\left(2\right)}\left[w\right]\right)=0\right\}\mathrm{,}$ (34c)

giving

$\tilde{I}=Z^{\star }\left[{\tilde{w}}^2-w_M^2+{\left(\delta w^{\left(2\right)}\right)}^2\right]$ (35a)

$\tilde{w}=w_M+\frac{\tilde{f}}{Z^{\star }}\mathrm{.}$ (35b)

Then the $F^{\left(4\right)}\left[w\right]$ approximation stretch, constrained to end in $\hat{w}$ and to be tangent to $T\left[w\right]$ in a point $w^{\text{u}}$ chosen by trial and error optimization, is given by:

$F^{\left(4\right)}\left[w\right]:=-2{\lambda }^{\text{u}}\left(w-\hat{w}\right)+{\sigma }^{\text{u}}{\left(w-\hat{w}\right)}^2,$ (36a)

$w^{\text{u}}:=\left(1-z\right)w_{\text{}f\text{}l\text{}x}+z\hat{w},z=0.575,$ (36b)

${∍}^{\prime }\left\{F^{\left(4\right)}\left[w\right]-T\left[w\right]=0\wedge \Delta \left(F^{\left(4\right)}\left[w\right]-T\left[w\right]\right)=0\right\},$ (36c)

giving

${\lambda }^{\text{u}}=\tilde{f}+\frac{2\hat{w}\tilde{f}-\tilde{I}}{w^{\text{u}}-\hat{w}},$ (37a)

${\sigma }^{\text{u}}=\frac{2\hat{w}\tilde{f}-\tilde{I}}{{\left(w^{\text{u}}-\hat{w}\right)}^2}.$ (37b)

4. The Approximate Analytic Solution

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

$\frac{\text{d}w}{\text{d}x}\left(x\right)=\mathcal{F}\left[w\left(x\right)\right]\mathrm{,}$ (38)

where $\mathcal{F}\left[w\right]$ is one of $F^{\left(i\right)}\left[w\right]$ $\left(i=1,2,3,4\right)$ or $f^{\left(2\right)}\left[w\right]$ or $T\left[w\right]$, with given $\alpha$ and β parameters (or βand γ) and initial conditions. For $\mathcal{F}\left[w\right]=F^{\left(1\right)}\left[w\right]$, from the definition in Equation (11), the initial condition is $w\left(x_0\right)=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 $\mathcal{F}\left[w\right]$ is at most a second order polynomial, Equation (38) is indeed quite trivially solved, giving a generalized logistic function.

For $\mathcal{F}\left[w\right]=F^{\left(1\right)}\left[w\right]$ :

$w^{\left(1\right)}\left(x\right)=\frac{w_1^{\left(1\right)}+w_2^{\left(1\right)}k{\text{e}}^{-x/\gamma {\tau }_1}}{1+k{\text{e}}^{-x/\gamma {\tau }_1}}\mathrm{,}$ (39a)

$k=\frac{w_1^{\left(1\right)}-i_0}{i_0-w_2^{\left(1\right)}},{\tau }_1=\frac{1}{\gamma \cdot \left(\alpha -1/2\right)\left(w_1^{\left(1\right)}-w_2^{\left(1\right)}\right)}.$ (39b)

For $\mathcal{F}\left[w\right]=f^{\left(2\right)}\left[w\right]$, thus $\alpha \le {\alpha }_{\text{cr}}$ :

$w^{\left(f\right)}\left(x\right)=\frac{\hat{w}+\left(2w_M-\hat{w}\right){\text{e}}^{-\left(x-x_M\right)/\gamma {\tau }_f}}{1+{\text{e}}^{-\left(x-x_M\right)/\gamma {\tau }_f}}\mathrm{,}$ (40a)

$x_M=\gamma {\tau }_1\ln \left(k\right)∍\text{'}{w}^{\left(1\right)}\left(x_M\right)=w_M\mathrm{,}{\tau }_f=\frac{\delta w^{\left(f\right)}}{\delta w^{\left(1\right)}}{\tau }_1>{\tau }_1\mathrm{.}$ (40b)

For $\mathcal{F}\left[w\right]=F^{\left(2\right)}\left[w\right]$, thus $\alpha >{\alpha }_{\text{cr}}$ :

$w^{\left(2\right)}\left(x\right)=\frac{w_1^{\left(2\right)}+w_2^{\left(2\right)}{\text{e}}^{-\left(x-x_M\right)/\gamma {\tau }_2}}{1+{\text{e}}^{-\left(x-x_M\right)/\gamma {\tau }_2}}\mathrm{,}$ (41a)

$x_M=\gamma {\tau }_1\ln \left(k\right)∍\text{'}{w}^{\left(1\right)}\left(x_M\right)=w_M\mathrm{,}{\tau }_{\text{2}}=\frac{\delta w^{\left(2\right)}}{\delta w^{\left(1\right)}}{\tau }_{\text{1}}>{\tau }_{\text{1}}\mathrm{.}$ (41b)

For $\mathcal{F}\left[w\right]=F^{\left(3\right)}\left[w\right]$, thus ${\alpha }_{\text{cr}}<\alpha \le 6$ :

$w^{\left(3\right)}\left(x\right)=\frac{\hat{w}-\left(\hat{w}+2{\lambda }^{\circ }/\sigma \right){\varphi }^{\circ }{\text{e}}^{-\left(x-x^{\circ }\right)/\gamma {\tau }_{\text{3}}}}{1-{\varphi }^{\circ }{\text{e}}^{-\left(x-x^{\circ }\right)/\gamma {\tau }_{\text{3}}}}\mathrm{,}$ (42a)

$x^{\circ }\mathrm{<mo>\; =\; </mo>}\gamma x_M+\gamma {\tau }_{\text{2}}\ln \left(\frac{w^{\circ }-w_2^{\left(2\right)}}{w_1^{\left(2\right)}-w^{\circ }}\right)∍\text{'}{w}^{\left(2\right)}\left(x^{\circ }\right)=w^{\circ }\mathrm{,}$ (42b)

${\varphi }^{\circ }=\frac{\hat{w}-w^{\circ }}{\hat{w}-w^{\circ }+\frac{2{\lambda }^{\circ }}{\sigma }}\mathrm{,}{\tau }_{\text{3}}=\frac{1}{2{\lambda }^{\circ }\gamma }\mathrm{.}$ (42c)

For $\mathcal{F}\left[w\right]=T\left[w\right]$, thus $\alpha >6$ (see (34) and (35)):

$w^{\left(T\right)}\left(x\right)=\frac{1}{2\tilde{f}}\left[\tilde{I}-\left(\tilde{I}-2\tilde{w}\tilde{f}\right){\text{e}}^{-\left(x-\tilde{x}\right)/\gamma \tilde{\tau }}\right]\mathrm{,}$ (43a)

$\tilde{\tau }=\frac{1}{2\tilde{f}\gamma }\mathrm{,}\tilde{x}=\gamma x_M+\gamma {\tau }_{\text{2}}\ln \left(\frac{\tilde{w}-w_2^{\left(2\right)}}{w_1^{\left(2\right)}-\tilde{w}}\right)∍\text{'}{w}^{\left(2\right)}\left(\tilde{x}\right)=\tilde{w}\mathrm{.}$ (43b)

Finally for $\mathcal{F}\left[w\right]=F^{\left(4\right)}\left[w\right]$ thus $\alpha >6$ :

$w^{\left(4\right)}\left(x\right)=\frac{\hat{w}-\left(\hat{w}+2{\lambda }^{\text{u}}/{\sigma }^{\text{u}}\right){\varphi }^{\text{u}}{\text{e}}^{-\left(x-x^{\text{u}}\right)/\gamma {\tau }_4}}{1-{\varphi }^{\text{u}}{\text{e}}^{-\left(x-x^{\text{u}}\right)/\gamma {\tau }_4}},$ (44a)

$x^{\text{u}}=\gamma \tilde{x}+\gamma \tilde{\tau }\ln \left(\frac{\tilde{I}-2\tilde{f}\tilde{w}}{\tilde{I}-2\tilde{f}{w}^{\text{u}}}\right){∍}^{\prime }{w}^{\left(T\right)}\left(x^{\text{u}}\right)=w^{\text{u}}$ (44b)

${\varphi }^{\text{u}}=\frac{\hat{w}-w^{\text{u}}}{\hat{w}-w^u+\frac{2{\lambda }^{\text{u}}}{{\sigma }^{\text{u}}}},{\tau }_4=\frac{1}{2{\lambda }^{\text{u}}\gamma }.$ (44c)

It is convenient to introduce

$\stackrel{\vee }{r}\left(t\right):=r\left(\gamma t\right),\stackrel{\vee }{i}\left(t\right):=i\left(\gamma t\right),\stackrel{\vee }{s}\left(t\right):=s\left(\gamma t\right),\stackrel{\vee }{w}\left(t\right):=w\left(\gamma t\right),\text{etc}.,.$ (45)

Then, from Equation (11) one has

$\stackrel{\vee }{r}\left(t\right)=\frac{1}{\alpha }\ln \frac{1-i_{\text{0}}}{1-w\left(\gamma t\right)}\mathrm{,}$ (46)

so that

$\stackrel{\vee }{i}\left(t\right)=\frac{\text{d}\stackrel{\vee }{r}}{\text{d}t}\left(t\right)=\frac{1}{\alpha }{\left[\frac{1}{1-w\left(x\right)}\frac{\text{d}w}{\text{d}x}\left(x\right)\right]}_{x=\gamma t}\mathrm{.}$

On the other hand Equation (12) implies

$\frac{1}{1-w}\frac{\text{d}w}{\text{d}x}=\alpha w-\ln \frac{1-i_{\text{0}}}{1-w}$

and consequently (see Equation (46))

$\stackrel{\vee }{i}\left(t\right)={\left[w\left(x\right)-\frac{1}{\alpha }\ln \frac{1-i_0}{1-w\left(x\right)}\right]}_{x=\gamma t}=\stackrel{\vee }{w}\left(t\right)-\stackrel{\vee }{r}\left(t\right).$ (47)

Finally, of course, due to (4):

$\stackrel{\vee }{s}\left(t\right)=1-\stackrel{\vee }{i}\left(t\right)-\stackrel{\vee }{r}\left(t\right)=1-\stackrel{\vee }{w}\left(t\right).$ (48)

In the case of the SIRD model one defines

$\stackrel{\vee }{r}=\stackrel{\vee }{h}+\stackrel{\vee }{d}\mathrm{,}$ (49a)

$\gamma \to \gamma +\mu \text{so}\text{that}\stackrel{\vee }{h}=\frac{\gamma }{\gamma +\mu }\stackrel{\vee }{r}\text{and}\stackrel{\vee }{d}=\frac{\mu }{\gamma +\mu }\stackrel{\vee }{r}\mathrm{.}$ (49b)

Figure 4 shows a comparison between the numerical “exact” solutions of the SIRD model and the approximate solutions of this work with $\beta =0.25$ and $\alpha =1.6,2.5,4.5,8.3$.

^{1 $\left(a\le b\right)$} ? then $c=f$ : otherwise $c=g$.

Imitating a formal expression typical of computing languages¹, the result for w can be summarized as follows:

Figure 4. Comparison of “exact” numerical solutions and approximate solutions for the SIRD model.

for $\alpha \le {\alpha }_{\text{cr}}$ (50a)

$\stackrel{\vee }{w}\left(t\right)=\left(t\le t_M\right)\mathrm{<mo>\; ?\; </mo>}\stackrel{\vee }{w}{}^{\left(1\right)}\left(t\right)\mathrm{<mo>\; :\; </mo>}\stackrel{\vee }{w}{}^{\left(f\right)}\left(t\right)$ (50b)

for ${\alpha }_{\text{cr}}<\alpha \le 6$ :

$\stackrel{\vee }{w}\left(t\right)=\left(t\le t_M\right)\mathrm{<mo>\; ?\; </mo>}\stackrel{\vee }{w}{}^{\left(1\right)}\left(t\right)\mathrm{<mo>\; :\; </mo>}\left(\left(t\le t^{\circ }\right)\mathrm{<mo>\; ?\; </mo>}\stackrel{\vee }{w}{}^{\left(2\right)}\left(t\right)\mathrm{<mo>\; :\; </mo>}\stackrel{\vee }{w}{}^{\left(3\right)}\left(t\right)\right)$ (50c)

for $\alpha >6$ :

$\begin{array}{c}\stackrel{\vee }{w}\left(t\right)=\left(t\le t_M\right)\mathrm{<mo>\; ?\; </mo>}\stackrel{\vee }{w}{}^{\left(1\right)}\left(t\right)\mathrm{<mo>\; :\; </mo>}\\ \left(\left(t\le \tilde{t}\right)\mathrm{<mo>\; ?\; </mo>}\stackrel{\vee }{w}{}^{\left(2\right)}\left(t\right)\mathrm{<mo>\; :\; </mo>}\left(\left(t\le t^u\right)\mathrm{<mo>\; ?\; </mo>}\stackrel{\vee }{w}{}^{\left(T\right)}\left(t\right)\mathrm{<mo>\; :\; </mo>}\stackrel{\vee }{w}{}^{\left(4\right)}\left(t\right)\right)\right)\mathrm{.}\end{array}$ (50d)

Similarly for $\stackrel{\vee }{s}\left(t\right)$, $\stackrel{\vee }{i}\left(t\right)$, $\stackrel{\vee }{h}\left(t\right)$ and $\stackrel{\vee }{d}\left(t\right)$.

In practice one does:

· solve numerically the transcendental ordinary Equation (15) to get $\hat{w}$ ;

· use Equation (21) and Equation (39) to get $w^{\left(1\right)}\left(x\right)$ as in Equation (39);

· for $\alpha \le {\alpha }_{\text{cr}}$ use Equation (22), Equations (28d) and (28e) to get $w^{\left(f\right)}\left(x\right)$ as in Equation (40);

· for $\alpha >{\alpha }_{\text{cr}}$ use Equation (27), Equation (26b), Equation (22), Equation (24b), Equation (25) and Equation (41) to get $w^{\left(2\right)}\left(x\right)$ as in Equation (41);

· for ${\alpha }_{\text{cr}}<\alpha \le 6$ use Equation (32), Equation (33a) and Equation (33b), Equation (33c) and finally Equation (42) to get $w^{\left(3\right)}\left(x\right)$ as in Equation (42);

· for $\alpha >6$ use Equation (34b), Equation (35) and Equation (43) to get $w^{\left(T\right)}\left(x\right)$ as in Equation (43);

· for $\alpha >6$ use Equation (17), Equation (36b), Equation (37) and Equation (44) to get $w^{\left(4\right)}\left(x\right)$ as in Equation (44);

· eventually use Equation (46), Equation (47), Equation (48), Equation (49).

The four plots in Figure 4 are produced by a C++ code implementing the above steps, then sending the produced analytic function to the graphing utility “gnuplot”: the C++ code could be re-used easily to fit-study data.

5. A Useful Feature

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

$\frac{1}{w^{\left(1\right)}{}^{{}^2}}\frac{\text{d}\text{}{w}^{\left(1\right)}}{\text{d}x}=-\frac{A}{w^{\left(1\right)}{}^{{}^2}}\left(w^{\left(1\right)}-w_2^{\left(1\right)}-\delta w^{\left(1\right)}\right)\left(w^{\left(1\right)}-w_2^{\left(1\right)}\right)\mathrm{.}$ (51)

Using the explicit solution Equation (39), one has

$w^{\left(1\right)}-w_2^{\left(1\right)}=\frac{2\delta w^{\left(1\right)}}{{\left[1+k{\text{e}}^{-\left(x-x_{\text{0}}\right)/\gamma {\tau }_1}\right]}^2}\mathrm{.}$ (52)

and consequently

$\frac{1}{w^{\left(1\right)}{}^{{}^2}}\frac{\text{d}\text{}{w}^{\left(1\right)}}{\text{d}x}=4Ak\frac{{\left(\delta w^{\left(1\right)}\right)}^2}{w_1^{\left(1\right)}{}^{{}^2}}\frac{{\text{e}}^{-\left(x-x_{\text{0}}\right)/\gamma {\tau }_1}}{{\left[1+\frac{w_1^{\left(1\right)}}{w_2^{\left(1\right)}}k{\text{e}}^{-\left(x-x_{\text{0}}\right)/\gamma {\tau }_1}\right]}^2}.$ (53)

Typically

$\left|\frac{w_1^{\left(1\right)}}{w_2^{\left(1\right)}}\right|\ll 1\text{and}\left|\frac{w_1^{\left(1\right)}}{w_2^{\left(1\right)}}k\right|\lesssim \mathrm{1,}$ (54)

but anyway with $t-t_0$ greater than some ${\tau }_1$ ’s, in the end one can write

$\ln \left(\frac{1}{{\left(\stackrel{\vee }{w}{}^{\left(1\right)}\right)}^2}\frac{\text{d}\stackrel{\vee }{w}{}^{\left(1\right)}}{\text{d}t}\right)\left(t\right)\simeq \ln \left(4A\gamma k\right)-\frac{t-t_{\text{0}}}{{\tau }_1}\mathrm{.}$ (55)

Analogous results hold for all the approximation stretches in the different $\alpha$ intervals as summarized in Equation (50); for instance, with $t-t^{\circ }$ greater enough then ${\tau }_3$, one has

$\ln \left(\frac{1}{{\left(\stackrel{\vee }{w}{}^{\left(3\right)}\right)}^2}\frac{\text{d}\stackrel{\vee }{w}{}^{\left(3\right)}}{\text{d}t}\right)\left(t\right)\simeq \ln \left[4\sigma \gamma {\varphi }^{\circ }{\left(\frac{2{\lambda }^{\circ }}{\sigma \hat{w}}\right)}^2\right]-\frac{t-t^{\circ }}{{\tau }_{\text{3}}}\mathrm{.}$ (56)

These piecewise linear behaviors can be seen in Figure 5 for $\alpha =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 Equation (47)) that $w\left(x\right)=r\left(x\right)+i\left(x\right)$, so that w is directly related to the data. The plot on the right shows that the function $\stackrel{\vee }{r}\left(t\right)$ of the removed individuals exhibits an analogous behavior: since in the SIRD model the $\stackrel{\vee }{d}\left(t\right)$ function is a fraction of $\stackrel{\vee }{r}\left(t\right)$, then one has the analogous behavior for the function of the deceased individuals. Figure 6 refers to the data of the deceased individuals during the winter-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 $\tau$, are completely determined by the parameters $\alpha$ and $\beta$ (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.

Figure 5. A couple of checks of the particular piece-wise linearity disclosed by the approximation approach of the present work against the numerical solution of the exact equations.

Figure 6. The special piece-wise linearity disclosed in the present work as exhibited by the real data for Covid-19 originated deaths during the winter-spring 2020 wave in Italy.

6. Conclusions

In this paper the equations of the SIR(D) epidemiological model are replaced by approximate ones, whose solutions are totally defined uniquely by the basic reproduction ratio α and the fractional removal rate γ (alternatively by $\beta =\gamma /\alpha$ ). These solutions are continuous (with the first derivative) chains of two or three or four generalized logistic related functions, the number depending on the value of $\alpha$ only; they are summarized in Equation (50) and easily implementable and usable, for instance, to fit-study data.

The analytic geometry based approximation method used here is novel and set stably at least over the range of the measured values of the basic reproduction ratio for several known pandemic diseases. A useful feature of the SIR(D) model, never disclosed before, is also given.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Kermack, W.O. and McKendrick, A.G. (1927) Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Socoety A, 115, 700-721. https://doi.org/10.1098/rspa.1927.0118
[2]	Murray, J. (1993) Mathematical Biology. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-662-08542-4
[3]	Daley, D. and Gani, J. (1999) Epidemic Modelling. Cambridge University Press, Cambridge.
[4]	Brauer, F. (2017) Mathematical Epidemiology: Past, Present, and Future. Infectious Disease Modelling, 2, 113-127. https://doi.org/10.1016/j.idm.2017.02.001
[5]	Martcheva, M. (2015) An Introduction to Mathematical Epidemiology. Springer, Berlin. https://doi.org/10.1007/978-1-4899-7612-3
[6]	Brauer, F., Castillo-Chavez, C. and Feng, Z. (2019) Mathematical Models in Epidemiology. Springer, Berlin. https://doi.org/10.1007/978-1-4939-9828-9
[7]	Ozyapici, A. and Bilgehan, B. (2018) Generalized System of Trial Equation Methods and Their Applications to Biological Systems. Applied Mathematics and Computation, 338, 722-732. https://doi.org/10.1016/j.amc.2018.06.020
[8]	Barlow, N.S. and Weinstein, S.J. (2020) Accurate Closed-Form Solution of the SIR Epidemic Model. Physica D: Nonlinear Phenomena, 408, Article ID 132540. https://doi.org/10.1016/j.physd.2020.132540
[9]	Kroeger, M. and Schlickeiser, R. (2020) Analytical Solution of the SIR-Model for the Temporal Evolution of Epidemics. Part A: Time-Independent Reproduction Factor. Journal of Physics A: Mathematical and Theoretical, 53, Article ID 505601. https://doi.org/10.1088/1751-8121/abc65d
[10]	Pakes, A.G. (2015) Lambert’s W Meets Kermack-McKendrick Epidemics. IMA Journal of Applied Mathematics, 80, 1368-1386. https://doi.org/10.1093/imamat/hxu057
[11]	Fowler, A.C. and Hollingsworth, T.D. (2015) Simple Approximations for Epidemics with Exponential and Fixed Infectious Periods. Bulletin of Mathematical Biology, 77, 1539-1555. https://doi.org/10.1007/s11538-015-0095-3
[12]	Cramer, J.S. (2002) The Origins of Logistic Regression, TI 2002-119/4. Tinbergen Institute Discussion Paper. https://doi.org/10.2139/ssrn.360300
[13]	Heesterbeek, J.A.P. (2002) A Brief History of R0 and a Recipe for Its Calculation. Acta Biotheoretica, 50, 189-204. https://doi.org/10.1023/A:1016599411804