A SIR Epidemic Model Allowing Recovery

: 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.


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), (1) 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 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 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 Sections 3-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 Sections 3-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 , 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 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: The first of these is the reciprocal of the more commonly occurring relative removal rate.

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 As is well known and clear from these expressions, 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 ↑ ∞, Differentiation will show (for all possible values of δ) that I ′′ (t) has the same sign as the product then this second factor has a unique zero t i ∈ (0, ∞).We thus conclude: is concave-increasing, and if (7) holds then there is a time of inflection, 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

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 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 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 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 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.
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 This functional equation has the form attributed (incorrectly) to J. Lambert, and we write its solution as 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 Another consequence of ( 12) is that

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: Proof.Simply observe from (2) that and recall implications from Section 3.
Lemma 3. In all cases τ > 0, equivalently, 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 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 In what follows, we need only, and will, consider Cases 1 and 3.It follows from (10) and its limiting form that where Hence, S(t) and R(t) approach their limits at the same rate.In addition, Hence, the state variables approach their limits at the same rate and hence it suffices to determine how quickly I(t) → 0.
Proof.Let ℓ(t) = e τt I(t) and observe from Lemma 2 that where 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 It follows from Lemma 2 and (2) that Next, it follows from (3) that

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 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.

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.
Next, twice differentiating (10) yields 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 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 The following theorem embraces all three state functions.
(i) I(t) is decreasing and it is ultimately convex-decreasing.
(ii) If then I(t) is concave-convex with a unique inflection point, t i , which solves (iv) R(t) is concave-increasing and S(t) is convex-decreasing.
We now consider I(t).Differentiating (1) and using (2) to eliminate first-order derivatives on the right-hand side yields the identity 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 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 Hence, S ′ (t)(R 0 S(t) − S 0 ) > 0, i.e., the left-hand side of (19) increases from 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 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.

R(t)
In addition, the maximum number of infectives is and S(t m ) = S 0 /R 0 .
Remark 1.Note the simplification of ( 21) for the SIR model where a = 0.
The susceptible numbers function behaves in a more complicated manner.
It follows from Theorem 5 that the epidemic curve has its maximum at t = 0 if (22) holds and the maximum is at t si < t m if (23) holds.The values of c * := max t≥0 (−S ′ (t)) are stated in our next result.
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 −1 0 then the epidemic begins more slowly and, in terms of the loss of susceptibles, it gathers pace until the time t si , after which it slows down.Observe too that the condition B(t si ) = 0 is equivalent to It then follows from (10) that and hence the second evaluation in Lemma 5 is made explicit because 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 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 ii ∈ (0, t m ), such that Proof.The proof relies on determining how σ(t) behaves along the trajectory of S(t) in (0, t m ).We know from (8) that dR dS = − γ βS − ρ and that the right-hand side < 0 in (0, ∞).Hence, and the right-hand side is positive for t ∈ (0, t m ).
Remark 4. Condition (23) implies that σ(t si ) < 0, i.e., I(t) is concave-increasing in a neighbourhood of t si and hence t si > t ii .

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).
Proof.Observe first that (28) is equivalent to Dividing ( 11) by γ renders it as Writing y = αR ∞ and expanding the exponential term in powers of y yields, after cancellation, the identity 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 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, ρ.
It also follows from (33), or (38) with c = 0, that , we obtain the approximation which generalises the known result for the SIR model.In particular, we still observe the incremental change 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 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: Theorem 9.If the assumptions of Theorem 8 hold, then where where 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 the same form as (4.4) in [5] (with differing notation) derived for the SIR model.

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.

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.,