Turning points of analytical solutions to the linear,
constant-population SIR model for campylobacteriosis and cryptosporidiosis
Graham McBride and Sharleen Harper, NIWA, Hamilton, NZ, March 2014
;
Under assumptions stated in McBride et al. ("Projected changes in reported campylobacteriosis and cryptosporidiosis rates as a function of climate change: a New Zealand study", SERRA, submitted), the governing differential equations for a constant-size population (N) for age-dependent but time-invariant (static) SIR model compartments (Susceptibles, Ill,Recovered) are:
dS/da = b – αS – cK1K2S – cK1(1 – K2)S + δR,
dI/da = cK1K2S – (α+ γ)I,
dR/da = γI + cK1(1 – K2)S – (α + δ)R,
where:
a denotes age, Sa=0 = N, Ia=0 = 0, Ra=0 = 0 and N = S + I + R.Parameters are:
b = specific immigration rate of susceptibles (# per annum);
c = specific rate of contact with pathogen (per annum)
K1 = probability of infection given contact;
K2 = probability of illness given infection;
γ= reciprocal of the shedding period (per annum);
α = specific death rate in population (per annum);
δ = specific immunity loss rate (per annum).
Note that summing these equations gives us b = αN (recruitment rate = death rate), which is used to remove b from the first equation.
The analytical solutions presented by McBride & French (2006) refer.[†]They are cited in square brackets. A summary table of findings is given near the end of this document.
We consider three cases, based on the sign of the model's discriminant (): (i) = 0; (ii) > 0; (iii) < 0, where ν2 –y (see below).
Case (i): = 0
The first derivative of the illness rateequation [16] is
(1)
where
(2)
with
(3)
and
(4)
We will also need
(5)
and
(6)
These relationships will also be useful
(7)
(8)
(9)
(10)
(11)
(12)
Therefore, noting that = 0 y = , we have
(13)
Now we need a useful expression for the steady states in terms of the variables defined above, i.e.,
(14)
(15)
and to get r it seems simplest to use the fact that r = 1 – (s + i),[‡] giving
(16)
Therefore
(17)
a pretty result. So the solution is
(18)
and so
(19)
This has extrema(where di/da = 0 in the a > 0 domain) at
(20)
This result suggests that will play a significant role in the form of the solutions to be obtained.[§]
The extrema predicted by (20) will be a maximum iff d2i/da2 < 0 at a = aext, and it will be a flat line if d2i/da2 = 0 (and all higher derivatives also vanish). So let's see. The second derivative is
(21)
and so
(22)
which is strictly non-positive, QED.
Susceptible curve
We have (from [15])
(23)
So, noting that
(24)
we have
(25)
and so
(26)
This equation is always negative because, from (7), we have = + cK1.
Note that at the extrema
(27)
Recovery curve
Noting that dr/da = –(ds/da + di/da), we have
(28)
and so
(29)
As a check, we can evaluate the solution directly, i.e., from [17]
(30)
Now – = + cK1 – – g = cK1(1 – K2) 0. Therefore (30) is always positive.
Also at the extrema we have
(31)
Check:
(32)
Case (ii): > 0
The first derivative of the illness rate equation [13] is
(33)
where
(34)
and
(35)
The steady-states are now
(36)
(37)
and to get r we again use the fact that r = 1 – (s + i), giving
(38)
We can now calculate
(39)
another pretty result! So by defining
,(40)
the derivative is
(41)
and so
(42)
Note that this solution collapses to the special case [when = 0, i.e., (19)].[**]
The infected proportion of the population, i(a), has maxima/minima when Equation (42) is equal to zero. This occurs as a, and also when/if the terms in the bracket vanish.[††] The latter occurs only if a value of a can be found such that ( + ) – ( – )ea = 0, that is
(43)
which isdefined only if . Let’s examine this term.
The numerator is always positive, so we need only look at the denominator. Noting that is strictly non-negative—but is not—we simply need.
(44)[‡‡]
and so the turning points ofi(a) are:
(45)
Now if < 0 and ||, aext will be negative and that is not in the solution domain of interest (which is a > 0). We therefore need to ensure that ( + )/( – ) > 1. Therefore (45) is modified to
(46)
Now consider the case (>0). The illness rate at aext will be a maximumiff d2i/da2 < 0 at that age—so let's see. The second derivative is
(47)
and so
(48)
thus
(49)
which, given that (> 0), is strictly negative, QED.
Let's now consider the gradients of s and r with age.
Susceptible curve
We have (from [12])
(50)
where
(51)
Now consider the expression in braces in (50); call it J. Therefore
(52)
The term in brackets can be multiplied out to give
(53)
so the solution is
(54)
This equation is notstrictly non-positive (a), because .
Check: as 0 the term in brackets in (54) becomes
(55)
Substituting (55) back into (54) we obtain
(56)
Which is equation (26), QED.
Finally, note that
(57)
and so
(58)
Recovery curve
Noting that dr/da = –(ds/da + di/da) we have
and so
(59)
This equation is not strictly non-positive (a), because .
Check: as 0 the term in brackets becomes
(60)
Substitution into (59) recovers (29), QED.
Finally, note that
(61)
and so
(62)
Case (iii): < 0
In this case, . Therefore
(63)
so that
(64)
(65)
and therefore, using r = 1 – (s + r),
(66)
Now the first derivative of the illness rate equation [19] is
(67)
Inserting (65)–(66) for the steady states, we have the cosine term [J1 =( + )i + gr]
(68)
and for the sine term [J2 =( – )i – gr]
(69)
and so we have
(70)
The extrema of i(a) occur when Equation (70) is equal to zero. This occurs at as a, and also when cos(a) = sin(a), for example, when
(71)
Note that, because of the periodic nature of the tangent function, there are many possible values for aext. In fact,
[§§](72)
The derivative of (70) is
(73)
Now in order to examine the sign of (73) at the extrema given by (72), we will need these identities[***]
for x > 0(74)
for x < 0, and(75)
for all real x,(76)
and also these identities
(77)
We must consider four possibilities for the sign of (73) at the extrema given by (72):
1. < 0 and n odd
and so
(78)
which is strictly negative, therefore are maxima.
2. < 0 and n even
and so
(79)
which is strictly positive, therefore are minima.
3.0 and n odd
and so
(80)
which is strictly positive, therefore are minima.
4. > 0 and n even
and so
(81)
which is strictly negative, therefore are maxima.
Recovered curve
We have (from [20])
(82)
The coefficient of the cosine term [J1 = i +g( – )r] can be written as
(83)
The coefficient of the sine term [J2 = i +g( + )r] can be written as
(84)
Substituting (83) (84) into (82) we obtain
and so
(85)
Also, using the relation that dr/da = –(ds/da + di/da), we have
(86)
Note that as 0, then (85) & (86) collapse to the zero discriminant cases (27) & (31).
Furthermore, using (74–(77), we have:
1.For < 0 and n odd
and so
.(87)
Also,
and so
.(88)
2.For < 0 and n even
and so
.(89)
Also,
and so
.(90)
3.For 0 and n odd
(91)
and so
.(92)
Also,
(93)
and so
.(94)
4. For > 0 and n even
and so
.(95)
Also,
and so
.(96)
Speculation:
The maximum value of i predicted by [19] occurs using (72) with n = 0 (when > 0) and with n = 1 (when < 0). That this is so is now shown.
The maxima/minima of i(a) are located at
For < 0, the maxima are given by aext,n with n odd, and for > 0 the maxima are given by aext,n with n even.
1.For < 0 and n odd,
and so
(97)
Since is strictly positive, decreases with increasing n. As a result, the smallest aext,n which is positive will correspond to the absolute maximum of in the solution domain. The smallest positive aext,n is given by the lowest whole odd number n satisfying
.(98)
Now, therefore . As a result, , and so from (98). Consequently, the lowest whole odd number satisfying (98) is n = 1, and therefore aext,1 corresponds to the absolute maximum of i(a) in the solution domain.
2.For > 0 and n even,
and so
(99)
As above, since is strictly positive, decreases with increasing n. As a result, the smallest aext,n which is positive will correspond to the absolute maximum of in the solution domain. The smallest positive aext,n is given by the lowest whole even number n satisfying (98).
In this case, therefore . As a result, , and so from (98). Consequently, the lowest whole even number satisfying (98) is n = 0, and therefore aext,0 corresponds to the absolute maximum of i(a) in the solution domain.
Note that i = g( – )/(2 + 2), so (97) & (99) show that i(aext) > i.
Special case: = 0
(i)If = 0, then from (19), (20), (22), (26) and (29),
aext, di/da = d2i/da2 = ds/da = dr/da = 0.
(ii)If 0, then from (42), (46), (49), (54) and (59),
aext, di/da = d2i/da2 = ds/da = dr/da = 0.
(iii)If < 0, then using (70) in the a > 0 domain, di/da = 0, when aext = (2n + 1)/n = 0, 1, 2,…,
, n = 0, 1, 2,…(100)[†††]
and so from (73)
(101)
which is strictly non-positive, QED.
Finally, using (85) & (86) we have
(102)
(103)
The value of i at finite aext
(i)For aext = 1/. Noting that i = g(–)/2 and i – r = –/, we have (from [16]),
(104)
giving
(105)
Note that i = g( – )/2, so (105) shows that i(aext) > i.
(ii)For aext = ln[(+)/(–)]/(). Noting that i = g(m – )/(2– 2) and i – r = –(–)/(–, by rewriting [13] at aext as
(106)
and noting that exp(aext) – 1 = /( – ), we obtain
(107)
and so
(108)
Note that i = g( – )/(2 + 2), so (108) shows that i(aext) > i.
(iii)For we merely appeal to (97) and (99).
1
O:\ESR12204\Working\Papers (CCID)\Resubmit drafts\Resubmitted\SERR_13_0025 McBride_Tait_Slaney_SIR_CC_zoonosis_paper_SI_resubmitted (24 March 2014) No Track Changes.doc
SUMMARY OF SOLUTIONS (in the positive a domain)
Variable / Case (i): = 0 / Case (ii): > 0 / Case (iii): < 0aext / / /
i(aext) / / /
The parameter groups are: = + ( + cK1 + )/2; = ( – cK1 – )/2; = ( + cK1 – )/2; g = cK1K2; = + g; = ||, where = 2 – y and y = g; = 2. The age-at-peak-illness-rate is aext, the value of a at which di/da = 0(which can be demonstrated from the above). Note that aext:<0,=0 = /.
1
APPENDIX: Special cases
If K2 = 0, we have (from McBride & French 2006)
(A.1)
(A.2)
(A.3)
(A.4)
(A.5)
So defining s = S/N and r = R/N, (A.1) becomes
(A.6)
where
(A.7)
and so
(A.8)
(A.9)
(A.10)
If = 0 and = 0, we have
(A.11)
and so
(A.12)
(A.13)
(A.14)
where the Methuselah states are now
(A.15)
(A.16)
(A.17)
Analytical solutions to the linear SIR model
1
[†] McBride, G.B.; French, N.P. (2006). Accounting for age-dependent susceptibility and occupation-dependent immune status: a new linear analytical model. WSEAS Transactions on Mathematics 11(5): 1241–1246. They are repeated on the last page of these notes.
[‡] It is simplest because the numerator of rin [21] is pcK1 – g, and just how this can be translated into a useful combination of the variables defined on the first page isn't immediately obvious.
[§] The result for 0 in (20) follows by noting that in (19)g(1 – a)e–a = g(1 – a)/ea = g(1 – a)/[1 + (a)/1! +(a)2/2! +(a)3/3! + …] = (g/a – )/[1/a + +a/2 +a2/6 + …] 0 as a.
[**] As 0 we have as the approximate right-hand-side of (42): g[ – a]e–a/(noting that e1 + as 0). This collapses to the right-hand-side of(19), QED [i.e., to g(1 – a)e–a].
[††] Note that (42) can be written as di/da = (g/)[( + )e–(+)a – ( – )e–(–)a]. For di/da 0 as a we therefore require that and this is always true [noting that = (p + q)/2, = (p – q)/2 and = |2 – y| where y = cK1K2 > 0].
[‡‡]Equation (44) can be shown to be equivalent to requiring that (1 – K2). To see that,we note from the preceding equation that we want 22. Using (6) we therefore require that 2 +2cK1 + (cK1)22 – cK1K2 2 + cK12. Then, using (5), – cK1 – + cK1K2, thus (1 – K2).
[§§] Only those extrema which are positive lie in the solution domain, andaext,n > 0 iffn > –[tan–1(/)]/. Note too that if = 0 and di/da = 0, (70) requires that cos(aext) = 0 [so sin(aext) = 1]. To keep the second derivative negative (at aext) we must have sin(aext) > 0 in (73), giving aext= (2n+ 1)/(n = 0, 2, 4,…) . But we must also maximize [19] and so the smallest age must be take, i.e., aext = /.
[***]Given at
[†††] Note that tan–1(x) = /2 – 1/x + 1/(3x3) – 1/(5x5) + … (x 1), so as → 0 then / (≡ x) → and so Eq. (71) → /(2).