A singularly perturbed age structured SIRS model with fast recovery

Age structure of a population often plays a significant role in the spreading of a disease among its members. For instance, childhood diseases mostly affect the juvenile part of the population, while sexually transmitted diseases spread mostly among the adults. Thus, it is important to build epidemiological models which incorporate the demography of the affected populations. Doing this we must be careful as many diseases act on a time scale different from that of the vital processes. For many diseases, e.g. measles, influenza, the typical time unit is one day or one week, whereas the proper time unit for the vital processes is the average lifespan in the population; that is, 10-100 years. In such a case, the epidemiological model with vital dynamics becomes a multiple time scale model and thus it often can be significantly simplified by various asymptotic methods. The presented paper is concerned with an SIRS type disease spreading in a population with a continuous age structure modelled by the McKendrick-von Foerster equation. We consider a disease with a quick recovery rate in a large population. Though it is not too surprising that in such a model the introduced disease quickly vanishes, the result is mathematically interesting as the error estimates are uniform on the whole infinite time interval, in contrast to the typical results based on the Tikhonov theorem and classical asymptotic expansions.

Abstract. Age structure of a population often plays a significant role in the spreading of a disease among its members. For instance, childhood diseases mostly affect the juvenile part of the population, while sexually transmitted diseases spread mostly among the adults. Thus, it is important to build epidemiological models which incorporate the demography of the affected populations. Doing this we must be careful as many diseases act on a time scale different from that of the vital processes. For many diseases, e.g. measles, influenza, the typical time unit is one day or one week, whereas the proper time unit for the vital processes is the average lifespan in the population; that is, 10-100 years. In such a case, the epidemiological model with vital dynamics becomes a multiple time scale model and thus it often can be significantly simplified by various asymptotic methods. The presented paper is concerned with an SIRS type disease spreading in a population with a continuous age structure modelled by the McKendrick-von Foerster equation. We consider a disease with a quick recovery rate in a large population. Though it is not too surprising that in such a model the introduced disease quickly vanishes, the result is mathematically interesting as the error estimates are uniform on the whole infinite time interval, in contrast to the typical results based on the Tikhonov theorem and classical asymptotic expansions.
1. Introduction. In many diseases the transmission rates vary significantly with age. In fact, for the exanthematic diseases (measles, scarlet fever) the transmission mainly involves early ages, while for sexually transmitted diseases the principal mechanism of infection only involves mature individuals. Thus, it can be expected that the interplay of the demographic processes with the infection mechanism will produce a nontrivial dynamics.
In this paper we consider the SIRS model introduced in [9]. In this model, the total population can be modelled by the linear McKendrick-von Foerster model describing the evolution in time of the density of the population with respect to age a ∈ [0, ω], ω < ∞, denoted by n(a, t). The demography of the population is driven by death and birth processes with vital rates µ(a) and β(a), respectively. Due to the epidemics, we split the population into susceptibles, infectives and recovered n(a, t) = s(a, t) + i(a, t) + r(a, t), so that the scalar McKendrick equation for n splits into the system ∂ t s(a, t) + ∂ a s(a, t) + µ(a)s(a, t) = −Λ(a, i(·, t))s(a, t) + δ(a)i(a, t), ∂ t i(a, t) + ∂ a i(a, t) + µ(a)i(a, t) = Λ(a, i(·, t))s(a, t) − (δ(a) + γ(a))i(a, t), ∂ t r(a, t) + ∂ a r(a, t) + µ(a)r(a, t) = γ(a)i(a, t), where γ(a) and δ(a) are age specific removal and recovery rates, respectively (here the removal refers to the recovery with permanent immunity). The function Λ is the infection rate (or the force of infection). In this paper we consider the so-called intercohort model where K is a nonnegative bounded function which accounts for the age dependence of the infections. For instance, for typical childhood diseases K should be large for small a and a and close to zero for large a or a . System (1) is supplemented by the boundary conditions where p, q ∈ [0, 1] are the vertical transmission parameters of infectiveness and immunity, respectively. Finally, we prescribe the initial conditions However, in building (1) some care should be taken not to mix different time scales.
In fact, if we consider a human population, then the death and birth rates are measured in units 1/year or even the time unit often is taken to be the average lifespan in the population, [14]. Then, in a human population with average life-span of, say, 70 years, averaged µ would be 1 and averaged β also would be 1 (2 children per couple in the lifetime).
On the other hand, if we model the disease such as flu or measles, then the recovery rate is measured in units 1/day (average duration of the disease, which is the inverse of the recovery rate, is 2-7 days). Thus, if we are to use 70 years as the unit of time in (1), then the numerical values of δ and γ, used in the literature for the SIRS compartmental model, should be multiplied by 70 × 365. Certainly, for other diseases, such as HIV/AIDS, the duration of the disease is at the same time scale as the vital dynamics and the above argument cannot be applied.
As far as the scaling of the infection force Λ is concerned, the situation can vary. For a simple compartmental SIR model Λ = λi, where λ is constant. Then, [8], in the population with vital dynamics (that is, with the constant size of the population) we have λ = γR 0 /N , where R 0 is the so-called basic reproduction rate and N is the size of the population. The constant R 0 for, say, measles, is 18, so the relative magnitude of λ depends on the size of the population: if the population is large, then this term can be included in the 'small' terms, while if it is small, it is a 'large' term. In this paper we shall be concerned with the former case and, for the sake of (mathematical) completeness, we allow the population to vary without affecting the value of λ. Thus we consider and is a small parameter reflecting the ratio of the typical time scales of the vital and epidemiological processes. Further, the bounded operator B : The main aim is to analyse the behaviour of the solution (s , i , r ) as → 0. The well-posedness of this problem is discussed in Section 2.
To formulate the main result, we have to introduce some notation.
being the positive cone of X. The norm · will refer to the relevant norm in L 1 . If necessary, we indicate the norm in a specific space X by writing · X . For any measurable function φ on [0, ω] we introducē Then we assume A1: µ ≥ 0 is a measurable function on [0, ω] and µ > 0; with δ > 0, γ > 0. We note that for the model to be realistic, the death rate µ must satisfy additional condition (15). Since this condition does not have any bearing on the formulation of main results, it will be discussed later. Next, we note that if (s , i , r ) satisfies (5) then, by adding the equations and the side conditions, we find that the total population n(t, a) = s (t, a) + i (t, a) + r (t, a) satisfies ∂ t n(a, t) = −∂ a n(a, t) − µ(a)n(a, t), n(0, t) = ω 0 β(a)n(a, t)da, n(a, 0) =: and is independent of . It is known, [9], that there is a dominant eigenvalue λ µ ≤β − µ, defined as the unique real solution to where, see [1,9,14], µ(s)ds (9) is the probability of survival of an individual till age a, and a constant m such that where λ µ is negative, zero or positive if and only if the net reproduction rate R µ = ω 0 β(a)Π µ (a)da is, respectively, smaller, equal or greater than 1. As we mentioned earlier, the considered scaling is realistic for models in which there are no significant changes of the total population. Hence, throughout the paper we will assume though some results will be formulated for a general R µ . Further, let γ # (a) = γ(a) γ(a) + δ(a) , δ # (a) = δ(a) γ(a) + δ(a) (12) and letw be the solution to the McKendrick problem ∂ tw (a, t) = −∂ aw (a, t) − µ(a)w(a, t), Then we have • r) be such that (s , i , r ) is a classical solution to (5). Then, for any ϑ < δ + γ there are constants C 1 , C 2 , C 3 , depending only on the coefficients of the problem and the D(S) ∩ D(M) norm of the initial conditions, such that for all sufficiently small > 0 uniformly for t ∈ R + . Furthermore, if R µ < 1, then C i can be replaced by C i e λµt , i = 1, 3, so that the error decays to zero as t → ∞.
We emphasize that this result shows that the solution of the nonlinear problem (5) can be approximated by the solution of two scalar linear problems and explicitly given initial layer corrector, uniformly for any time.
2. Notation, assumptions and well-posedness results. Problems like (1)- (4) have been relatively well-researched, even in a much more complex setting involving nonlinear dependence of µ and β on the total population, see [15,Section 5.3] or [7,13] and references therein. Since, however, the results are scattered and refer to two distinct cases (with ω < +∞ and ω = ∞), we summarize basic facts in the form relevant to the problem at hand, see [11]. Though the model most resembles that discussed in [15], the main difference is that in op.cit. the maximum age is infinite but the death rate is bounded, which simplifies some arguments. The boundedness of ω affects the analysis of the linear part of (1)-(4). We note that the biologically realistic assumption is that ω < ∞; that is, that no individual can live indefinitely. This requires the survival probability, given by (9), to satisfy Π µ (ω) = 0 which, in turn, yields Hence, µ cannot be bounded as a → ω − . This is in contrast with the case ω = ∞, where commonly it is assumed that µ is a bounded function on R + (see however [12]) and introduces another unbounded operator, the multiplication by µ, in (1). In some papers, e.g. [10], this was circumvented by assuming that there is a maximum reproductive age a r < ω, so that the birth rate satisfies β(a) = 0 for a > a r , and hence ignoring the post-reproductive population by performing the analysis for a ∈ [0, a r ]. Note that in this way we lose conservativeness of the model. The analysis of the model without any simplifying assumption in the scalar linear case was done in [9] by reducing it to an integral equation along characteristics. It is known that the solutions obtained in this way generate a strongly continuous semigroup on L 1 ([0, ω]), see [2,6,15] (though in [15] it is assumed that ω = ∞.) The technical details needed to handle the case with unbounded µ on [0, ω] with ω < ∞ can be found in [11]. In particular, it follows that the realization of the operator A := S + M on the domain generates a positive C 0 -semigroup, denoted by e tA t≥0 . Since, for a fixed , (5) is a quadratic perturbation of the linear system, a standard argument, see [5,11], shows (5). This solution becomes a classical solution if This estimate can be improved. Indeed, as noted earlier, if (s (t), i (t), r (t)) is a classical solution to (5), then n(t, a) = s (t, a) + i (t, a) + r (t, a) is a classical solution to (7). We denote by (A, D(A)) the generator of the solution semigroup e tA t≥0 for (7), with the domain defined analogously to (16). Hence, since , r (t)) ≥ 0, we see that each component of u is controlled by the -independent solution n of (7): for each t ≥ 0 and almost every a ∈ [0, ω] and thus satisfy inequality (10).
3. Formal asymptotic expansion. Here we justify the terms of the approximation which appear in Theorem 1.1. Following the general approach of the asymptotic analysis, see e.g. [4], we are looking for the so-called hydrodynamic space V of the singularly perturbed equation (5) which, in this case is given by the null-space of C.
Since here a is a parameter, we can carry out the calculations for a fixed a. Then it is easy to see that V is a two dimensional subspace where e 1 , e 2 are an appropriate basis for V . The complementary spectral space W, called the kinetic space, corresponding to the eigenvalue λ(a) = −(δ(a) + γ(a)), is spanned by e 3 (a) = (−δ # (a), 1, −γ # (a)), see (12). To find the decomposition u = c 1 e 1 + c 2 e 2 + c 3 e 3 into the hydrodynamic and kinetic space, we have to find the coefficients c i , i = 1, 2, 3 which will give us the aggregated variables. Using standard results from linear algebra, see e.g. [2], c i = f i · u/f i · e i where f i are left eigenvectors of C corresponding to the same eigenvalue as e i , i = 1, 2, 3. Here, the situation is slightly more complicated as V is two dimensional and we have to select its basis in a convenient way. Since f 1 = (1, 1, 1) is a left eigenvector of C corresponding to the zero eigenvalue and f 1 · u = u 1 + u 2 + u 3 , it should play an important role in the analysis due to (7). Then, for convenience of calculations, we take e 2 = (−1, 0, 1) ∈ V, which is orthogonal to f 1 . We have some freedom in selecting f 2 : taking f 2 = (0, γ # (a), 1) yields e 1 = (1, 0, 0). Finally f 3 = (0, 1, 0) so that (19) Next, we use this decomposition to change variables in (5). Accordingly, we define n = s + i + r , w = γ # i + r and leave i unchanged. This yields the system ∂ t n = −∂ a n − µn, Thus, the pair (n, w ) belongs to the hydrodynamic space and i to the kinetic space. Since the total population n decouples from the system, there is no need to approximate it and, in the last two equations, it can be treated as a known function. Thus we shall focus on the pair (w , i ) and, following the Chapman-Enskog procedure, we only expand the kinetic part: Inserting it into the second equation of (20) and denoting the bulk approximation of w by w, we have Comparing coefficients at like powers of and using Λ(0) = 0 we find i k = 0 for all k ≥ 0. A similar behaviour of the bulk expansion recently has been observed in an alignment model, see [3]. Hence, denoting by i the bulk approximation of i , we arrive at the (formal) bulk approximation (n, w, i) = (n,w, 0), where n and w satisfy (7) and (13), respectively. Note that (13) is obtained by substituting the approximation i ≈ i = 0 in the equations for w in (20). To check how good the obtained approximation is, we have to find the error of the approximation. Denoting the error and inserting it to the equations, we find We see that there is at least one reason why the error cannot be of order -the initial condition is of order 1. Thus we have to introduce the initial layer correction by blowing up time according to τ = t/ and looking for the approximation were we anticipate that there is no need to introduce the initial layer for w as w satisfies the exact initial condition. Rescaling time in (20) we see that due to ∂ t = −1 ∂ τ , the only equation for the terms at the −1 level is which, subject to the initial conditionĩ 0 (0) = • i, gives the initial layer in the form The new error is given by The error equations for E can be obtained from (23) by expressing e w and e i in terms of e w and e i , according to (26). Since e w (a, 0) = 0, e i (a, 0) = 0.
We see that on the boundary we still have terms which are of an order lower than . Fortunately, both inhomogeneities on the boundary have an exponential decay in t → ∞ and → 0. Thus, we anticipate that we do not need the boundary layer, compare [6], but only the corner layer which is obtained by simultaneous rescaling of time and age according to τ = t/ , α = a/ . Unfortunately, the standard approach to the corner layer, [6], will not suffice here as the corner layer equation will not incorporate the unbounded operator M. Thus the classical corner layer will not belong to D(M) and we would not be able to substitute the error terms into the equation, as in (27). To alleviate the problem, we define the corner corrector to be the solution to w(a, 0) = φ a ,ȋ(a, 0) = 0, where φ is a function to be determined. The reason for introducing this function is that for the classical solvability of the problem with inhomogeneous boundary conditions one needs the equality of the boundary condition and the initial condition at (a, t) = (0, 0). This is easy to see here, first forȋ, as the equation forȋ is decoupled and can be solved separately. Then, solving it by reducing it to a Volterra equation, as in [9], shows that if the inhomogeneity is zero at t = 0, then the total birth rate is 0 at t = 0 which, together withȋ(a, 0) = 0, ensures the continuity of the solution across the diagonal t = a, see [11]. In our case we have where c is given by Since ω 0 β(a)e − a da → 0 as → 0, c is bounded as → 0.
The necessary estimates of the corner layer solution will be provided later. Here, assuming that it is sufficiently regular, we consider the new approximation w ≈ w +w and i ≈ĩ 0 +ȋ. The corresponding error, satisfies the system then there are constants K and B such that Proof. The proof is very similar to that in [9] for proving the AEG property of the McKendrick problem (with ψ = 0), hence we only provide a brief sketch of it.
Using the linearity, we can write the total birth rate B(t) = φ(0, t) as B(t) = B• φ (t) + B ψ (t) and focus on the case where K M (a) = β(a)Π M (a). To use the uniform notation for the Laplace transform, all functions defined on [0, ω] are considered to be extended by 0 outside this interval.
Taking the Laplace transform of (36), we obtain The difference with the calculations of [9], which will become relevant later, is that ψ is not an entire function. It is, however, analytic in λ > −c, hence the only singularities ofB ψ are due to the zeroes of 1−K M or are contained in the half-plane { λ ≤ −c}. SinceK M is an entire function, its zeroes are isolated of finite order (thus giving rise to poles ofB ψ which do not have finite accumulation point). As mentioned earlier, λ M is the only real root of multiplicity 1 of (8) and in each strip a < λ < b, there is at most a finite number of other roots. By (33), λ M is in the domain of analyticity ofψ. Let us consider the second term in the last formula of (37),Ĥ :=ψK M /(1−K M ). As in [9], on any line {σ + iy; y ∈ R} which does not meet any root of (8), we have InvertingĤ(λ), we have To estimate H(t) we note that, by properties ofĤ, we can shift the line of integration to {σ 1 + iy; y ∈ R}, where ζ < σ 1 < λ M , ζ = max{ λ 1 , −c} and λ 1 is the first eigenvalue such that λ 1 < λ M . Then, using the Cauchy theorem, we can write and H 2 satisfies the estimate where B 1 satisfies B 1 ≤ C k σ1 L2(R) φ σ1 L2(R) for some constant C. Thus, we arrived at the representation or, taking into account that σ 1 < λ M , and denoting B 2 = |B 0 | + |B 1 |, we can write Using the standard representation of the solution in terms of B and the estimates of [9, p. 22] for B• φ , we arrive at (35).
We want to specify (35) to the context of (28). First, let us denote α(a) = γ(a) + δ(a) and consider an arbitrary constant α * < α, see (6). We do not assume (11) here. Proposition 1. Letȋ andw be the solution to (28). Then there are constants C i and C w , depending on the coefficients of the problem and on the W 1 1 (R + ) norm of the initial condition • i, such that for any t ∈ R + , Proof. Let us begin with some general comments and notation. First, we observe that if M ≤M then, since λ M and λM are the unique real solutions to (8), from monotonicity it follows λM ≤ λ M .
(44) Let λ M θ be the (dominant) eigenvalue corresponding to M θ = µ + θ/ for any function of a, or a constant, θ. Note that λ M0 = λ µ , the dominant eigenvalue corresponding to M 0 = µ; that is, of the problem (7). Hereafter, we shall use that later notation. Then, by the definition of α * , we have Let λ be the solution to (8) with the largest real part λ = σ that is smaller than λ µ . Then λ − α * / is the solution to (8) specified to M α * ; that is, of the equation First we consider the equation forȋ. The problem with the direct application of (35) is that the function ψ in this case satisfies for some constant C 1 . Thus,ψ is analytic for λ > −α/ and, by (45), in general λ Mα will not belong to the domain of analyticity ofψ if λ µ ≤ 0. Thus, in this case, we have to settle on a slightly weaker estimate. Namely, by (47), we write and, if λ µ < 0, we assume that ≤ 0 with 0 < −(α − α * )/λ µ . Then λ M α * is in the domain of analyticity ofφ and we carry out the estimates for φ α * . We have Similarly, Hence, we can see that in this case for some constant C. Next, we need to estimate B 1 . To simplify calculations for σ < 0, we further assume that ≤ 0 where 0 < −(α − α * )/σ , otherwise can be arbitrary. Then both λ M α * and λ − α * / belong to the domain of analyticity ofψ. First we observe that inf y∈R |1 −K M α * (σ + iy)| ≥ m > 0 independently of , where σ is between λ M α * and the real part of the next solution to (46). In fact, by the considerations above, we can take σ =σ − α * / , where σ <σ < λ µ is independent of . Then and the claim follows. To complete the estimate of B 1 we calculate where C 3 is bounded independently of . Let us estimateȋ from the relevant equation in (28). By (35), we have Hence, taking into account thatȋ(0) = 0, (48) and (49) give ȋ (t) ≤ C 4 e λµt− α * t + C 5 e − α * t which yields (42). Next we use this estimate forw. Here In this case λ M = λ µ and, using the above, Hence, |B 2 | ≤ C 14 . Therefore, by (29), (35) and (50) with α * = 0 and < α * /λ µ if λ µ > 0, we arrive at (43).
Indeed, for 0 ≤ t ≤ 2ω the left hand side is smaller or equal than 1 and the right hand side is larger or equal than 1. Then, for t ≥ 2ω we have