Next Article in Journal
Robustness Analysis for Sundry Disturbed Open Loop Dynamics Using Robust Right Coprime Factorization
Previous Article in Journal
Hermite–Hadamard–Mercer Inequalities Associated with Twice-Differentiable Functions with Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A SIR Epidemic Model Allowing Recovery

Department of Mathematics & Statistics, University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia
Axioms 2024, 13(2), 115; https://doi.org/10.3390/axioms13020115
Submission received: 20 December 2023 / Revised: 29 January 2024 / Accepted: 2 February 2024 / Published: 8 February 2024
(This article belongs to the Topic Mathematical Modeling)

Abstract

:
The deterministic SIR model for disease spread in a closed population is extended to allow infected individuals to recover to the susceptible state. This extension preserves the second constant of motion, i.e., a functional relationship of susceptible and removed numbers, S ( t ) and R ( t ) , respectively. This feature allows a substantially complete elucidation of qualitative properties. The model exhibits three modes of behaviour classified in terms of the sign of S ( 0 ) , the initial value of the epidemic curve. Model behaviour is similar to that of the SIS model if S ( 0 ) > 0 and to the SIR model if S ( 0 ) < 0 . The separating case is completely soluble and S ( t ) is constant-valued. Long-term outcomes are determined for all cases, together with determination of the rate of convergence. Determining the shape of the epidemic curve motivates an investigation of curvature properties of all three state functions and quite complete results are obtained that are new, even for the SIR model. Finally, the second threshold theorem for the SIR model is extended in refined and generalised forms.

1. Introduction

The very well known general or SIR model of the spread of a disease in a large population lends itself to a fairly complete qualitative understanding owing to the fact that its fundamental system of three differential equations (DEs) possess two constants of motion. See [1] or [2] for accounts of this topic. The first of these is common to very many models: the size of the population is constant in time. The second constant of motion occurs because the time variable can be eliminated between the DEs for the number of susceptible and removed individuals, S ( t ) and R ( t ) , respectively, to give a functional relation between them that holds along each system trajectory.
Recall that the SIR model describes progression of susceptible individuals to possible infection and then ultimate removal with no chance of infecting any in the susceptible subpopulation. The mathematically simpler SIS model ([1,2]) allows for recovery of infected individuals back to full susceptibility. These models are frequently named, respectively, the general and the simple epidemic. In this paper, we consider an extension that interpolates between these models. This extension allows for infected individuals, I ( t ) in number, to either be removed at rate γ or to recover with immediate susceptibility at rate ρ . We assume that γ + ρ > 0 . The immediacy of recovery to a susceptible state can, perhaps, be justified on the assumption that any period of immunity is so brief as to be negligible in comparison with typical periods of infectiousness and susceptibility and that final states are essentially achieved by times that are short relative to demographic changes.
So, assuming infection occurs from homogeneous mixing of infected and susceptible individuals with rate β > 0 , the governing DEs are
S ( t ) = β S ( t ) I ( t ) + ρ I ( t ) ,
I ( t ) = β S ( t ) I ( t ) γ I ( t ) ρ I ( t ) ,
and
R ( t ) = γ I ( t ) .
Initial values are denoted by S 0 > 0 , I 0 > 0 and R 0 = 0 . Thus, if all rate parameters are positive, transitions between compartments occur as for the SIR model with an additional flux of individuals at rate ρ from the infected compartment back to the susceptible compartment.
This model could be denoted by SIR∨S, pronounced ‘serves’, and this is, to the author’s knowledge, due to Mulkern and Nosrati [3] (whose notation for the parameters differs). They pursue the analytic consequences of the system equations to the extent of noting the second constant of motion (with a minor error) and illustrate its behaviour with some numerical calculation and graphical displays. The aim of this paper is to derive as far as is possible the analytical consequences of the above system.
Before describing the structure of the paper, we recall the so-called first threshold theorem for the SIR model. It is known for this model that lim t 0 I ( t ) = 0 and I ( t ) either decreases or initially increases and then decreases. These modes are taken to represent, respectively, a minor or major epidemic. Their occurrence is determined by the sign of the initial rate of infection, I ( 0 ) , specifically I ( 0 ) > 0 for a major epidemic. These initial conditions are characterised by the first threshold theorem asserting the obvious conclusion from (2) that I ( 0 ) > 0 iff S 0 > γ / β . The lower bound is called the relative removal rate and it is the threshold that separates a minor/major outbreak. The threshold criterion is usually expressed in terms of the basic reproduction ratio R 0 = β S 0 / γ , i.e., a major outbreak occurs iff R 0 > 1 .
The basic reproduction ratio is similarly defined for other epidemic models in terms of the sign of I ( 0 ) , and for the SIR∨S model it is clear from (2) that
R 0 = β S 0 γ + ρ .
Restricted to the SIS or SIR models, R 0 can be interpreted as the average number of secondary infections caused by introducing a single infective into a large susceptible population. This extends to the SIR∨S model, as follows. The number of infections, I t , by time, t, caused by a single infective arriving at t = 0 equals N ( β S 0 t ) , where N ( · ) is a unit-rate Poisson process. The duration, T, of infectivity is modelled as a competing risks situation, i.e., T = min ( ϵ γ , ϵ ρ ) , where ϵ γ and ϵ ρ are independent random variables having exponential distributions with parameters γ and ρ and they are independent of N ( · ) . Hence, T has an exponential distribution with parameter γ + ρ . It thus follows that
E ( I T ) = E N ( β S 0 T ) = β S 0 E ( T ) = β S 0 γ + ρ .
Many of the formulae derived for the SIR∨S model have γ > 0 appearing as a denominator. Hence, the results for the SIS model cannot simply be read off from the general results established in Section 3, Section 4, Section 5 and Section 6. In addition, if γ > 0 , then the fundamental property of the SIR model that I ( t ) 0 subsists for the SIR∨S model (Lemma 1), whereas the SIS model admits an endemic level of infection under some circumstances. Consequently, we begin in Section 2 by recalling (mostly) known facts about the SIS model. In addition to behaviour determined by the sign of R 0 1 , where R 0 = β S 0 / ρ , there are curvature behaviours of I ( t ) that qualitatively differ according to the sign of the difference δ : = β N ρ .
In Section 3, Section 4, Section 5 and Section 6, we assume that γ , ρ > 0 , the full SIR∨S model, and derive (in Section 3) the second constant of motion, a functional relation between S ( t ) and R ( t ) , and its consequence for the limiting values S and R of susceptible and removed numbers. This reveals three cases (denoted 1, 2 and 3, respectively) according as β S 0 > ρ , = ρ or < ρ , where S 0 = S ( 0 ) . It is clear from (1) that these cases correspond, respectively, to S ( 0 ) < 0 , S ( 0 ) = 0 and S ( 0 ) > 0 . All three of these can occur for the SIS model and only the first is possible for the SIR model. Case 2 is a degenerate case having a complete elementary solution. We also establish the existence of the long-term limits.
We show in Section 4 that the three state functions converge exponentially quickly to their limiting values, with a common time constant, τ , that depends on the three parameters and the final number S of susceptibles.
In Section 5, we examine the curvature properties of the state functions. To appreciate why, recall that the epidemic curve is defined as the production rate of new infectives, S ( t ) . Thus, the maximum c * of the epidemic curve will occur in ( 0 , ) iff S ( t ) has an infection point therein (Theorem 5 and Lemma 5). This invites a more thorough investigation of second-order properties. Observe too that, e.g., S ( t ) measures the instantaneous acceleration of susceptible numbers. In particular, the speeding up, or slowing down, of the state functions corresponds to regions of their convexity, or concavity. Inflection points delineate, e.g., transition times from acceleration to deceleration. Much of this is new, even for the SIR model. Case 3 exhibits the simplest behaviour—there are no inflection points and, quite unlike the SIR model, S ( t ) is increasing. Case 1 shares behaviour with the SIR model: S ( t ) decreases and both it and R ( t ) have at most one inflection and I ( t ) has one or two inflections.
The second threshold theorem for the SIR model is an assertion of the final size, R (i.e., total number of removed infected individuals), in the case that S 0 is just a little larger than the threshold, γ / β . More precisely, if S 0 = γ / β + ϵ , where 0 < ϵ 1 , and if I 0 ϵ 2 , then R 2 ϵ . Equivalently, S γ / β ϵ , i.e., final susceptible numbers end as far below the threshold as they began above it. See [1,2]. This fundamental result, going back to the pioneering work of Kermack and McKendrick, is usually derived by using a Maclaurin expansion of the exponential function occurring in the functional equation relating S and R . See, e.g., p. 84 in [1] (where the constraint on I 0 is stated a little obscurely—see the equation just after (6.9)) or p. 28 in [2] (an imprecise assertion in that I 0 is specified as being ‘small relative to ϵ ’) and the proof on pp. 30, 31 where the assumption is I 0 ϵ 2 .
In Theorem 8, we follow this mode of proof for the SIR∨S model with the precise specification that I 0 = c ϵ 2 , where 0 c < . The consequences of the general outcome show that the general form of the second threshold theorem subsists for the SIR∨S model, even if c > 0 . If c = 0 , then the above outcome for the SIR model is unchanged, i.e., the parameter ρ is absent from the first-order of approximation. Lemma 6 asserts a consequence of the proof of Theorem 9 in which we allow a larger initial number of infectives: I 0 C ϵ 1 + ω , where ω [ 0 , 1 ) . Finally, by using a branch-point expansion of the Lambert function we extend, in Theorem 9, the approximation (34) by estimating the quadratic and cubic components of the O ( ϵ 2 ) term.
A concluding section ends the paper. The following derived constants frequently occur:
α = β / γ and a = ρ / γ .
The first of these is the reciprocal of the more commonly occurring relative removal rate.

2. The SIS Model

If γ = 0 , then R 0 = β S 0 / ρ and the DE (2) reduces to a logistic DE because in that case S ( t ) + I ( t ) N . Defining δ = β N ρ , its solution is
I ( t ) = I 0 e δ t 1 + ( I 0 β / δ ) e δ t 1 > 0 if   δ 0 , I 0 1 + β t if   δ = 0 .
As is well known and clear from these expressions,
lim t I ( t ) = I = N ρ / β if   δ > 0 , 0 if   δ 0 .
However, the finer aspects of I ( t ) (and S ( t ) ) also depend on cases analogous to 1–3 specified above.
We observe first that if δ 0 (implying R 0 < 1 ), then it is easy to check that I ( t ) is convex-decreasing to zero, and hence S ( t ) is concave-increasing to N.
The case δ > 0 exhibits richer behaviour, as follows (but omitting algebraic detail). Clearly, R 0 = β S 0 / ρ > 1 implies that δ > 0 and that this latter can hold if R 0 1 . The above limit result for infective numbers takes the monotone form: As t ,
I ( t ) I if   I 0 > I , i . e . , R 0 < 1 ; I if   I 0 = I , i . e . , R 0 = 1 ; I if   I 0 < I , i . e . , R 0 > 1 .
Differentiation will show (for all possible values of δ ) that I ( t ) has the same sign as the product
I I 0 × I I 0 I 0 e δ t .
Hence, if δ > 0 and R 0 < 1 , i.e., I 0 > I , then I ( t ) is convex-decreasing to I .
If R 0 > 1 , then the second factor in (6) decreases from I 2 I 0 to . Hence, if 2 I 0 < I , equivalently,
R 0 > 1 + β I 0 / ρ ,
then this second factor has a unique zero t i ( 0 , ) . We thus conclude: If 1 < R 0 1 + β I 0 / ρ then I ( t ) is concave-increasing, and if (7) holds then there is a time of inflection,
t i = δ 1 log R 0 1 β I 0 / ρ > 0 ,
such that
I ( t ) i s c o n v e x i n c r e a s i n g i n ( 0 , t i ) , c o n c a v e i n c r e a s i n g i n ( t i , ) .
These results give a precise expression to the idea that the infection initially accelerates if R 0 is sufficiently larger than unity before necessarily slowing as it approaches its endemic level, I . Recalling the interpretation of R 0 as the initial relative per-capita rate of increase of infectives, the dual quantity R ^ 0 = β I 0 / ρ is the initial per-capita relative rate of increase of recoveries. Hence, we have a threshold result for the SIS model asserting that I ( t ) is initially convex-increasing iff R 0 R ^ 0 > 1 .

3. The Second Constant of Motion

We assume in the sequel that γ > 0 and ρ 0 . Observe that, exactly as for the SIR model, R ( t ) > 0 for all t and hence R ( t ) R , say, and R N . Next, dividing (3) into (1) yields
d S d R = β S ρ γ ,
whose solution for t > t > 0 is
log β S ( t ) ρ β S ( t ) ρ = β γ R ( t ) R ( t ) .
This is valid for all parameter combinations because (3) implies that R ( t ) is strictly increasing and hence the argument of the logarithm function must be less than unity. Assuming that R ( 0 + ) = 0 , then letting t 0 and exponentiating the result we obtain for t 0 that
β S ( t ) ρ = β S 0 ρ e ( β / γ ) ( R ( t ) .
It follows that the sign of β S ( t ) ρ is the same as the sign of β S 0 ρ , i.e., the sign of S ( 0 ) . This invariance distinguishes three cases.
Case 1:
β S 0 > ρ 0 . Here, susceptibles are initially being infected at a faster rate than they are recovering. Recalling (5), the solution (9) becomes
β S ( t ) = ρ + ( β S 0 ρ ) e α R ( t ) .
It is clear that susceptible numbers preserve the classical SIR behaviour: S ( t ) S , say, as t . Setting ρ = 0 yields the long-known relation for the SIR model.
Case 2:
β S 0 = ρ . In this case, we have the degenerate outcome that β S ( t ) ρ because, independently of the value of γ , the rate of infection is exactly balanced by the rate of recovery to the susceptible state. It then follows that (2) reduces to I = ρ I , whence
I ( t ) = I 0 e γ t
and (3) yields
R ( t ) = I 0 1 e γ t .
Hence, infected numbers decrease convexly to zero and the final size of the epidemic is R = I 0 .
Case 3:
β S 0 < ρ , implying that ρ > 0 . Here, (7) subsists with the result that S ( t ) S as t .
We observe in passing that (25) in [3] is a rearrangement of (10), with a minor error in that the first occurrence of γ 2 on their right-hand side should be γ 1 . Observe too that it is not obvious what should result from (10) if γ 0 because then α and, of course, R ( t ) 0 in the actual limit value γ = 0 described in Section 2.
We wish to let t in (10). As previously remarked, R = lim t R ( t ) exists and hence so does I = lim t I ( t ) . The following fundamental result generalises the well-known, long-term outcome for the SIR model.
Lemma 1.
If γ > 0 , then the limit I = 0 .
Proof. 
If I > 0 , it follows from (3) that there exists t > 0 such that R ( t ) > 1 2 I if t > t . This implies that R ( t ) , a clear contradiction. □
This result shows that even the slightest degree of removal prevents an endemic level of infection. In particular, S + R = N , so taking the limit in (10) yields the functional equation
β S ρ = ( β S 0 ρ ) e α ( N S ) .
Defining x = S ρ / β and x 0 = S 0 ρ / β , and recalling (5), identity (11) can be recast as
α x e α x = α x 0 e a α N .
This functional equation has the form attributed (incorrectly) to J. Lambert, and we write its solution as
β S ρ = γ W ( α S 0 a ) e a α N ,
where W ( · ) denotes the principal branch of the Lambert W-function [4]. This is the appropriate choice because, if not, then the left-hand side would be an unbounded function of the argument of W ( · ) . See [5] for the version with ρ = 0 and the references there.
Observe that we can let γ 0 in (11) because α , which yields S ρ / β , agreeing with the analytical outcome for the SIS model. The explicit solution (12) is also consistent in this regard. To see this, observe for the SIS case that ρ < β N , implying that the argument of W ( · ) tends to zero as γ 0 . In fact, we have the approximation
S = ρ / β + ( β S 0 ρ ) e ( ρ β N ) / γ ( 1 + o ( 1 ) ) ,
as γ 0 .
Another consequence of (12) is that
S > ρ / β in Case 1 , < ρ / β in Case 3 .

4. The Rate of Convergence

We begin with a lemma that introduces the convergence exponent, which we denote by τ .
Lemma 2.
The function I ( t ) is log-concave in Case 1 and it is log-convex (hence convex) in Case 3. In addition, the following limit exists:
τ : = lim t I ( t ) / I ( t ) = ( β S γ ρ ) .
Proof. 
Simply observe from (2) that
d 2 d t 2 log I ( t ) = β S ( t )
and recall implications from Section 3. □
Next, we establish that τ > 0 .
Lemma 3.
In all cases τ > 0 , equivalently,
β S γ + ρ < 1 .
Proof. 
For Cases 2 and 3, it follows from (11) that β S ρ 0 and hence that τ = γ ( β S ρ ) > 0 . Assume now that Case 1 holds and recast (12) as
β S ρ + γ W ( α S 0 a ) e a α N = 0 .
The argument of the Lambert function is negative, and hence W ( · ) 1 , implying that 0 β S ρ γ , i.e., τ 0 . This inequality is strict provided that the positive-valued term T = ( α S 0 a ) e a α N < e 1 , but S 0 < N and hence T < ( α N a ) e ( α N a ) e 1 . □
In what follows, we need only, and will, consider Cases 1 and 3. It follows from (10) and its limiting form that
β ( S t S ) = ( β S 0 ρ ) e α R ( t ) e α R = ( β S 0 ρ ) e α R e α ( R ( t ) R ) 1 = ( β S ρ ) ( α ( R ( t ) R ) ( 1 + o ( 1 ) ) ,
i.e.,
S ( t ) S A ( R R ( t ) ,
where
0 A = β S ρ γ = 1 τ / γ < 1 .
Hence, S ( t ) and R ( t ) approach their limits at the same rate. In addition,
I ( t ) = N S ( t ) R ( t ) ( 1 A ) ( R R ( t ) ) .
Hence, the state variables approach their limits at the same rate and hence it suffices to determine how quickly I ( t ) 0 .
Theorem 1.
If γ > 0 , then
0 < lim t e τ t I ( t ) < .
Proof. 
Let ( t ) = e τ t I ( t ) and observe from Lemma 2 that
( t ) = τ ( t ) + I ( t ) e τ t = ε ( t ) ( t ) ,
where
ε ( t ) = τ + I ( t ) / I ( t ) 0
as t .
We infer from the monotone nature of S ( t ) and (2) that I ( t ) has at most one critical point in ( 0 , t ) , and if such exists then it must be at a global maximum. Hence, there exists t > Argmax ( I ( t ) ) such that I ( t ) < 0 if t > t . The first equality in (14) is a differential equation whose solution is
I ( t ) = ( t ) exp τ t + t t ε ( v ) d v .
It follows from Lemma 2 and (2) that
ε ( t ) = β ( S ( t ) S ) c o n s t . I ( t ) , ( t ) .
Next, it follows from (3) that 0 I ( t ) d t = R / γ < and hence that t | ε ( t ) | d t < . The assertion now follows from (15). □
Approaching the convergence rate through I ( t ) seems to be the simplest route. To appreciate this, we could begin by defining δ ( t ) = R R ( t ) and observe that (3) yields
δ ( t ) = γ ( N S ( t ) R ( t ) ) = γ ( S S ( t ) + R R ( t ) ) τ δ ( t ) ,
where we have appealed to (13) and Lemma 3. The problem thus reduces to estimating τ + δ ( t ) / δ ( t ) . This looks to be less straightforward than the above proof.

5. Shape Properties

We now investigate the convexity/concavity properties of the state functions, beginning with the simplest case.
Theorem 2.
Assume Case 3, β S 0 < ρ . Then the functions S ( t ) and R ( t ) are concave-increasing and I ( t ) is convex decreasing.
Proof. 
Knowing that S ( t ) S in Case 3, it follows from the definition of τ and Lemma 3 that β S ( t ) γ ρ < 0 and hence that I ( t ) < 0 for all t > 0 . We conclude from (3) that R ( t ) = γ I ( t ) < 0 .
Next, twice differentiating (10) yields
β S ( t ) = α R ( t ) + R ( t ) ρ β S ( t ) .
Each factor on the right-hand side is positive, and hence S ( t ) < 0 . The convexity assertion for I ( t ) is now obvious. □
We now turn to Case 1. Recall that this case can occur if R 0 1 and if R 0 > 1 . The next result concerns the first of these possibilities. We define the constant
R = a 1 + a = ρ γ + ρ < 1
Recalling the stochastic model of infection from Section 1, we observe that R = P ( ϵ ρ < ϵ γ ) , the probability of recovery occurring before removal. In relation to parts (ii) and (iii) of the next result, note that
R 0 R = β S 0 ρ γ + ρ > 0 .
The following theorem embraces all three state functions.
Theorem 3.
Assume Case 1, β S 0 > ρ , and that R 0 1 . The following hold.
(i) I ( t ) is decreasing and it is ultimately convex-decreasing.
(ii) If
I 0 S 0 ( 1 R 0 ) 2 R 0 ( R 0 R ) ,
then I ( t ) is convex-decreasing in ( 0 , ) .
(iii) If
I 0 > S 0 ( 1 R 0 ) 2 R 0 ( R 0 R ) ,
then I ( t ) is concave–convex with a unique inflection point, t i , which solves
R 0 S ( t ) S 0 2 = R 0 I ( t ) R 0 S ( t ) a S 0 .
(iv) R ( t ) is concave-increasing and S ( t ) is convex-decreasing.
Proof. 
Observe first that, since S ( t ) is decreasing and R 0 1 ,
β S ( t ) γ ρ < β S 0 γ ρ 0 ,
i.e., I ( t ) < 0 for all t > 0 . The concavity of R ( t ) follows because, from (3), R ( t ) = γ I ( t ) < 0 . A double differentiation of (10) yields the identity
S ( t ) = ( β S 0 ρ ) e α R ( t ) α R ( t ) + ( α R ( t ) ) 2 > 0 ,
and (iv) follows.
We now consider I ( t ) . Differentiating (1) and using (2) to eliminate first-order derivatives on the right-hand side yields the identity
I ( t ) = I ( t ) β S ( t ) γ ρ 2 β I ( t ) ( β S ( t ) ρ ) .
The right-hand side is asymptotically equal to β S γ ρ 2 I ( t ) > 0 , and Assertion (i) follows.
It follows from (19) that I ( t ) is convex or concave within intervals determined by the sign of the right-hand side of (19). In addition, such intervals are separated by points of inflection that are solutions of (18), and this equation is obtained by dividing (19) by γ + ρ and equating the result to zero. We now show that (18) has at most one solution. As a function of t, the right-hand side of (18) is strictly decreasing to zero. So, its maximum value occurs at t = 0 .
The derivative of the left-hand side is
R 0 S ( t ) R 0 S ( t ) S 0 = R 0 S 0 S ( t ) γ + ρ β S ( t ) γ ρ .
However, S ( t ) is strictly decreasing for Case 1, hence S ( t ) < 0 . In addition, the second factor on the right-hand side decreases from
β S 0 γ ρ = ( γ + ρ ) ( R 0 1 ) 0 .
Hence, S ( t ) ( R 0 S ( t ) S 0 ) > 0 , i.e., the left-hand side of (19) increases from S 0 2 ( R 0 1 ) 2 up to ( R 0 S S 0 ) 2 > 0 . It follows that (19) either has no solution in ( 0 , ) , i.e., I ( t ) > 0 , or it has exactly one solution in ( 0 , ) . The latter case holds if (17) holds. All assertions now are evident. □
It follows from Theorem 3 (iii) that, although infected numbers decrease, if I 0 is sufficiently large then this decrease is initially quite slow—a concave decrease—before gathering speed and then slowing again. Observe that for the SIR model, ρ = 0 , Case 1 is the only possibility, and then, e.g., (17) takes the form
I 0 > S 0 R 0 1 1 2 .
Hence, the concave–convex behaviour is achieved if R 0 is sufficiently close to unity.
We now consider the case R 0 > 1 (implying Case 1), i.e., I ( t ) is increasing for small t. The next result extends a fundamental result for the SIR model.
Theorem 4.
If R 0 > 1 , then the function I ( t ) has a unique maximum, I max , at the point t m > 0 , a solution of β S ( t ) = γ + ρ , i.e.,
R ( t ) : = β S ( t ) γ + ρ = 1 .
In addition, the maximum number of infectives is
I max = N S 0 R 0 1 + log ( 1 + a ) R 0 a 1 + a
and
S ( t m ) = S 0 / R 0 .
Proof. 
It follows from (2) that I ( t ) = 0 for any t solving (20). In addition, R ( 0 ) = R 0 > 1 by hypothesis, and since S ( t ) S , it follows from Lemma 3 that R ( t ) β S / ( γ + ρ ) < 1 . Hence, (20) has a unique positive solution, t m . Evaluation (21) follows from I max = N S ( t m ) R ( t m ) , (20), (10) and the definition of R 0 . The final assertion follows from (2). □
Remark 1.
Note the simplification of (21) for the SIR model where a = 0 .
The recovered numbers Equation (3) implies that R ( t ) = γ I ( t ) , resulting in the following consequence.
Lemma 4.
If R 0 > 1 , then the removed numbers function,
R ( t ) i s c o n v e x i n c r e a s i n g i n ( 0 , t m ) c o n c a v e i n c r e a s i n g i n ( t m , ) .
The susceptible numbers function behaves in a more complicated manner.
Theorem 5.
(a) If R 0 > 1 , then the function S ( t ) is:
(a.i) log-concave and decreasing in ( 0 , t m ) , and
(a.ii) convex-decreasing in ( t m , ) .
(b) In addition:
(b.i) If
I 0 1 R 0 1 S 0 ,
then S ( t ) is convex-decreasing in ( 0 , ) ; and
(b.ii) If
I 0 < 1 R 0 1 S 0 ,
then there is an inflection point t s i ( 0 , t m ) , such that
S ( t ) i s c o n c a v e d e c r e a s i n g i n ( 0 , t s i ) , c o n v e x d e c r e a s i n g i n ( t s i , ) .
Proof. 
Twice differentiating the logarithm of each side of (10) gives the identity
d 2 d t 2 log S ( t ) = ϕ 1 ( t ) R ( t ) ϕ 2 ( t ) R ( t ) 2 ,
where the ϕ i ( t ) are positive-valued functions in ( 0 , ) . Assertion (a.i) follows from Lemma 4.
Twice differentiating (10) yields the identity
β S ( t ) = α ( β S 0 ρ ) e α R ( t ) B ( t ) , where B ( t ) = R ( t ) α R ( t ) 2 .
Assertion (a.ii) follows again from Corollary 5.1.
Next, it follows from (2) and (3) that
B ( t ) = γ I ( t ) β S ( t ) γ ρ β I ( t ) .
The right-hand side is decreasing in ( 0 , t m ) . Moreover,
B ( t m ) = β γ I 2 ( t ) < 0 ,
implying by continuity and (24) that S ( t ) is convex-decreasing in an interval ( t , ) where 0 t < t m .
Now,
S ¯ 0 γ ρ β I 0 = ( γ + ρ ) R 0 1 R 0 I 0 / S 0 .
The right-hand side 0 if (22) holds, and hence B ( t ) < 0 in ( 0 , t m ) and Assertion (b.i) follows.
If (23) holds, then B ( 0 ) < 0 and, since B ( t m ) < 0 and B ( t ) is decreasing, it follows that there exists exactly one number t s i ( 0 , t m ) such that S ( t s i ) = 0 , and (b.ii) follows. □
It follows from Theorem 5 that the epidemic curve has its maximum at t = 0 if (22) holds and the maximum is at t s i < t m if (23) holds. The values of c * : = max t 0 ( S ( t ) ) are stated in our next result.
Lemma 5.
If (22) holds, then
c * = S ( 0 ) = I 0 ( β S 0 ρ ) ,
and if (23) holds, then
c * = I ( t s i ) ( β S ( t s i ) ρ ) = I ( t s i ) ( γ β I ( t s i ) ) .
These outcomes make intuitive sense in as much as they imply that if I 0 / S 0 is large (which is unlikely) then susceptible numbers initially fall at maximum speed, whereas if I 0 / S 0 is less than the critical number 1 R 0 1 then the epidemic begins more slowly and, in terms of the loss of susceptibles, it gathers pace until the time t s i , after which it slows down. Observe too that the condition B ( t s i ) = 0 is equivalent to
R ( t s i ) = N S 0 / R 0 .
It then follows from (10) that
β S ( t s i ) = ρ + ( β S 0 ρ ) e α ( N S 0 / R 0 )
and hence the second evaluation in Lemma 5 is made explicit because
I ( t s i ) = N S ( t s i ) R ( t s i ) = S 0 / R 0 S ( t s i ) .
Next, still assuming R 0 > 1 , we look at the shape of the infective curve in the interval ( 0 , t m ) . Recall that the sign of I ( t ) is the same as the sign of
σ ( t ) : = ( β S ( t ) γ ρ ) 2 β I ( t ) ( β S ( t ) ρ ) .
Theorem 6.
Assume that R 0 > 1 and 0 t t m .
(i) If σ ( 0 ) 0 , then I ( t ) is concave-increasing in ( 0 , t m ) .
(ii) If σ ( 0 ) > 0 , then there is a single point of inflection, t i i ( 0 , t m ) , such that
I ( t ) i s c o n v e x i n c r e a s i n g i n ( 0 , t i i ) , c o n c a v e i n c r e a s i n g i n ( t i i , t m ) .
Proof. 
The proof relies on determining how σ ( t ) behaves along the trajectory of S ( t ) in ( 0 , t m ) . We know from (8) that
d R d S = γ β S ρ
and that the right-hand side < 0 in ( 0 , ) . Hence,
d I d S = 1 + γ β S ρ = β S γ ρ β S ρ
and the right-hand side is positive for t ( 0 , t m ) .
It follows from (25) that
d σ d S = β 3 ( β S γ ρ ) β I .
Differentiating again yields
d 2 σ d S 2 = β 2 3 + β S γ ρ S ρ ,
which is positive for t [ 0 , t m ] , i.e., σ ( t ) is convex in S ( t ) as t traverses [ 0 , t m ] , but σ ( t m ) = β 2 I max < 0 . Assertions (i) and (ii) follow. In particular, if σ ( 0 ) > 0 , then t i i is the unique number in ( 0 , t m ) at which σ ( t ) vanishes. □
Remark 2.
We note that ( d σ / d S ) | t = t m = β 2 I max < 0 .
Remark 3.
The condition (22) for S ( t ) to be convex-decreasing is equivalent to β I 0 β S 0 γ ρ . Hence,
( β S 0 ρ ) β I 0 ( β S 0 ρ ) ( β S 0 γ ρ ) > ( β S 0 γ ρ ) 2 ,
implying that I ( t ) is concave-increasing in ( 0 , t m ) .
Remark 4.
Condition (23) implies that σ ( t s i ) < 0 , i.e., I ( t ) is concave-increasing in a neighbourhood of t s i and hence t s i > t i i .
Finally, consider the shape of I ( t ) in ( t m , ) . We have shown it is concave around t m and ultimately convex-decreasing. Hence, there is at least one inflection point in ( t m , ) . Recalling that β S ( t ) < γ + ρ if t > t m , it follows from (26) that d σ / d S < 0 , i.e., that σ ( t ) is increasing in ( t m , ) . In addition, it follows from (25) that σ ( t m ) < 0 and σ ( ) = τ 2 > 0 . Hence, σ ( t ) vanishes at exactly one point t i i > t m . This yields our final shape result.
Theorem 7.
If R 0 > 1 , then there exists a unique critical point t i i ( t m , ) , such that
I ( t ) i s c o n c a v e d e c r e a s i n g i n ( t m , t i i ) c o n v e x d e c r e a s i n g i n ( t i i , ) .

6. Threshold Theorems

The following result extends the well-known ‘second’ threshold theorem (p. 84 in [1] or p. 29 in [2]), which dates back to the pioneering work of Kermack and McKendrick in the 1920s. Our version is expressed in a slightly more precise manner than in [2]. Let r = ( γ + ρ ) / β , the relative rate of removal or recovery and recall notation (4).
Theorem 8.
Let ϵ ( 0 , 1 ) ,
S 0 = r + ϵ
(implying R 0 = 1 + ϵ / r ) and let I 0 c ϵ 2 ( ϵ 0 ), where c [ 0 , ) is a constant. Then the final size R satisfies
R A ( ϵ ) : = ϵ 1 + α ϵ 1 + 1 + ( 2 c / α ) ( 1 + α ϵ ) .
Proof. 
Observe first that (28) is equivalent to
α S 0 a = 1 + α ϵ .
Dividing (11) by γ renders it as
α S a = α ( S 0 + I 0 R ) a = ( 1 + α ϵ ) e α R .
Writing y = α R and expanding the exponential term in powers of y yields, after cancellation, the identity
( 1 + α ϵ ) y 2 2 α ϵ y 2 α c ϵ 2 + O ( y 3 ) = 0 .
This implies that y = O ( ϵ ) , which, neglecting the O ( y 3 ) term and solving the resulting quadratic equation, yields the exact result R = A ( ϵ ) + O ( ϵ 3 ) . □
Taking account of the last step of the proof, it follows from (29) that
R = 1 + 1 + 2 c γ / β ϵ + O ( ϵ 2 ) .
Not surprisingly, if c > 0 then the number removed exceeds the initial excess, ϵ , of susceptibles. What is surprising is that this dominant contribution is independent of the recovery rate, ρ .
If c = 0 , then (29) simplifies to the approximation
R 2 ϵ 1 + α ϵ = 2 ϵ α S 0 a .
Using (28) and defining γ ¯ = γ / β and ρ ¯ = ρ / β , respectively, the relative removal and recovery rates, then this approximation takes the form
R 2 γ ¯ 1 γ ¯ S 0 ρ ¯
which reduces to a classical form when ρ = 0 ; e.g., (2.3.14) in [2].
It also follows from (33), or (38) with c = 0 , that
R = 2 ϵ + O ( ϵ 2 )
and, since S = N R = S 0 R + O ( ϵ 2 ) , we obtain the approximation
S = r ϵ + O ( ϵ 2 )
which generalises the known result for the SIR model. In particular, we still observe the incremental change S 0 S = R 2 ϵ .
The following result admitting a relatively larger value of I 0 than in Theorem 8 is a corollary of its proof. It asserts a substantially stronger progression of the epidemic than seen in Theorem 8.
Lemma 6.
Let ϵ ( 0 , 1 ) and (28) hold. If I 0 C ϵ 1 + ω , where C > 0 and ω [ 0 , 1 ) , then
R 2 C α · ϵ 1 2 ( 1 + ω ) .
Refinements to (35) can in principle be obtained by using more terms in the expansion of e y , giving rise to polynomial equations of an order exceeding two. A more direct approach is using (12) and observing that outcome (35) corresponds to expanding the Lambert function around its branch point at e 1 . We specify the following constants:
c 1 = 1 + 2 c γ / β , d 1 = 2 α 3 c 1 2 and d 2 = α 2 2 c 1 2 .
Theorem 9.
If the assumptions of Theorem 8 hold, then
R = A 1 ϵ + A 2 ϵ 2 + A 3 ϵ 3 + O ( ϵ 4 ) ,
where
A 1 = 1 + c 1 , A 2 = 1 3 c α α c 1 a n d A 3 = α 2 36 7 c 1 + c 1 3 + 8 .
Proof. 
The required Lambert function expansion is
W e 1 z 2 / 2 = 1 σ ( z ) ,
where
σ ( z ) = z z 2 3 + z 3 36 + z 4 270 + .
See #4.13.6 in [4], and [5] for a survey of Lambert function expansions. Write (12) in the form
α R = α N a + W ( α S 0 a ) e a α N .
Recalling (30), we have
a α N = a α S 0 α I 0 = 1 α ϵ α I 0 ,
and hence in (38) we choose (as expected)
z 2 = 2 α ϵ + α I 0 log ( 1 + α ϵ ) = O ( ϵ 2 ) .
Recalling the assumptions and notation of Theorem 8, expanding the log-term yields
z 2 = ( α c 1 ϵ ) 2 1 d 1 ϵ + d 2 ϵ 2 + O ( ϵ 3 ) .
It follows that
z = α c 1 ϵ 1 1 2 d 1 ϵ 1 2 ( d 2 d 1 2 / 4 ) ϵ 2 + O ( ϵ 3 ) ,
whence
z 3 = α c 1 ϵ 3 + O ( ϵ 4 ) .
Collecting powers of ϵ and substituting into (39) using (37) and (38) leads to Assertion (36). □
Observe that the leading term corresponds to (38). In the classical case, c = 0 , (36) simplifies to
R = 2 ϵ 2 α 3 ϵ 2 + 2 α 3 2 ϵ 3 + O ( ϵ 4 ) ,
the same form as (4.4) in [5] (with differing notation) derived for the SIR model.

7. Concluding Remarks

We have examined the SIR model modified to allow recovery from the infected to the susceptible state at per-capita rate ρ . If γ = 0 , then this is the SIS model that exhibits an endemic level of infection if the derived parameter δ > 0 . On the other hand, if γ > 0 , then I ( t ) 0 —the epidemic ultimately fades. The extended model behaves similarly to the subcritical SIS model ( δ 0 ) if ρ is sufficiently large (i.e., ρ > β S 0 , implying R 0 < 1 ), in that S ( t ) and R ( t ) both increase to positive-valued limits (Theorem 2).
Behaviour similar to the SIR model is manifested if ρ is small, i.e., 0 ρ < β S 0 . If also R 0 1 , then I ( t ) decreases, possibly with an inflection. The other state functions are monotone with no inflections (Theorem 3). If R 0 > 1 , then curvature properties are more complicated, as seen in Theorems 4–7. In particular, I ( t ) has a single positive mode at t m with exactly one inflection in ( t m , ) , and it may, or may not, have an inflection in ( 0 , t m ) .
Elucidating these shape properties depends on the existence of two constants of motion—the total population size is constant the the relation (10) holds. The model can be generalised at the expense of one or both of these invariants by allowing for birth and death. In this case, it is usual to assume that newborns are susceptible. The first invariant is preserved under the assumption of balanced growth—the birth and death rates are equal. If they are not, then the population size either grows or diminishes exponentially fast. In either case, a successful analysis probably will be limited to identifying equilibria and determining their stability using results from the qualitative theory of differential equations. The case of balanced growth for the SIS and SIR models is surveyed in [6], and some more general models are treated there. The review paper [7] lists papers extending classical models to allow varying total population sizes.

Funding

This research received no external funding.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Bailey, N. The Mathematical Theory of Infectious Diseases, 2nd ed.; Charles Griffin: London, UK, 1975. [Google Scholar]
  2. Daley, D.; Gani, J. Epidemic Modelling: An Introduction; C.U.P.: Cambridge, UK, 1999. [Google Scholar]
  3. Mulkern, R.V.; Nosrati, R. A tutorial on common differential equations and solutions useful for modeling epidemics like COVID-19: Linear and non-linear compartment models. J. Appl. Math. Phys. 2022, 10, 3053–3071. [Google Scholar] [CrossRef]
  4. Olver, F.; Lozier, D.; Boisvert, R.; Clark, C. NIST Handbook of Mathematical Functions; C.U.P.: Cambridge, UK, 2010. [Google Scholar]
  5. Pakes, A.G. Lambert’s W meets the Kermack-McKendrick epidemics. IMA J. Appl. Math. 2015, 80, 1368–1386. [Google Scholar] [CrossRef]
  6. Hethcote, H.W. Qualitative analysis of communicable disease models. Math. Biosci. 1976, 28, 335–356. [Google Scholar] [CrossRef]
  7. Hethcote, H.W. The mathematics of infectious diseases. SIAM Rev. 2000, 42, 599–653. [Google Scholar] [CrossRef]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Pakes, A.G. A SIR Epidemic Model Allowing Recovery. Axioms 2024, 13, 115. https://doi.org/10.3390/axioms13020115

AMA Style

Pakes AG. A SIR Epidemic Model Allowing Recovery. Axioms. 2024; 13(2):115. https://doi.org/10.3390/axioms13020115

Chicago/Turabian Style

Pakes, Anthony G. 2024. "A SIR Epidemic Model Allowing Recovery" Axioms 13, no. 2: 115. https://doi.org/10.3390/axioms13020115

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop