Waning and boosting: on the dynamics of immune status

The aim is to describe the distribution of immune status (as captured by antibody level) on the basis of a within-host submodel for continuous waning and occasional boosting. Inspired by Feller’s fundamental work and the more recent delay equation formulation of models for the dynamics of physiologically structured populations, we derive, for given force of infection, a linear renewal equation. The solution is obtained by generation expansion, with the generation number corresponding to the number of times the individual became infected. Our main result provides a precise characterization of the stable distribution of immune status.


Introduction
The immune system defends an individual host against pathogens. After infection with a specific species of pathogen has been cleared, some memory remains, providing protection against future attacks of that same pathogen. Over time the memory wanes, until it is boosted by a new encounter. Thus the immune status is shaped by a combination of the exogenous process of exposure and endogenous processes (fighting the invader to achieve clearance of the infection and subsequent waning).
In practice the immune status of an individual is often quantified by measuring the concentration of specific antibodies in serum. Distributions of such serological measurements are used to assess the immune status of a population, for example in the context of vaccination programs (Wilson et al. 2012). The immune status of a population impacts the risk of outbreaks of an infection and can provide information on incidence of infection, including asymptomatic infection (Metcalf et al. 2016). Longitudinal changes of antibody titers of individuals may show the effects of boosting and waning over time, e.g., for pertussis (Versteegh et al. 2005). In mathematical models for vaccine preventable diseases, immunity is often represented by a dichotomous variable -individuals are either susceptible or immune-even though this distinction is not straightforward in reality. Therefore, it is useful to have a mathematical modeling framework that is capable of describing immunity as a continuous variable subject to waning and boosting over time. Our aim here is to provide a first step towards such a framework. We neglect all subtleties of specific infectious diseases and focus on the processes of waning and boosting in their simplest form.
So we ignore much of the subtlety and complexity of the immune system by postulating that the immune status is fully described by a positive quantity y (antibody titer against pertussis toxin is what we have in mind as a concrete example). Waning is described by the ordinary differential equation dy/dt = g(y) for the decline of y between encounters with the pathogen. Such encounters occur at rate . So is the constant force of infection and is considered a parameter (in Sect. 5 we shall briefly indicate how to formulate a feedback consistency condition for ; this condition involves assumptions about infectiousness and thus increases the number of parameters). We assume that, on the time scale set by g and , the time it takes the immune system to clear infection is negligible. This assumption allows us to introduce as a third model ingredient the instantaneous boosting map f that sends the immune status y just before the infection, to the immune status f (y) just after clearance of infection. In De Graaf et al. (2014) an explicit expression for f was derived from a submodel for the struggle between the pathogen and the immune system; see (Teunis et al. 2016) for a follow-up.
The three ingredients g, and f define a Piecewise Deterministic Markov Process (Davis 1993;Rudnicki and Tyran-Kamińska 2015). Indeed, waning and boosting are both deterministic, the only randomness is in the hitting times of the Poisson process with rate . Let the random variable Y (a), with Y (0) = y b , correspond to the immune status at age a of an immortal individual. Here y b , the immune status at birth, is another parameter (heterogeneity can of course be captured by an assumed distribution of y b ). Provided the pathogen under consideration does not contribute to mortality, we can later on introduce an independent age-specific survival probability.
Having established this mathematical framework, we now want to answer the following questions: (a) Can we compute the distribution of an individual's immune status at age a, given that it starts life with an immune level y b ? (b) If we let the process run for a long time, will the distribution of immune status converge to a stable distribution? Since feedback through on the transmission process is ignored, the stable distribution will describe the distribution of immunity in a population in steady state, if everybody is born with immune status y b .
Let y be an element of (0, ∞) and let be a measurable subset of (0, ∞). Development over time of the distribution of immune status is described by the kernel which, by assumption, does not depend on a. In Sect. 2 we formulate a Renewal Equation (RE) for Q and solve it by generation expansion (using the techniques of Sect. 4 of Diekmann et al. (1998) one can show that Q does have, as it should, the Chapman-Kolmogorov property; in the "Appendix A1" we formulate the more traditional Kolmogorov backward and forward PDE that are associated with Q). In Sect. 3 we introduce the corresponding next-generation operator and construct its fixed point that describes the stable distribution at the generation level. A general result from (Diekmann et al. 1998) relates the stable distribution at the generation level to the stable distribution of the process itself. This result is formulated in Sect. 4.

The RE and its solution by generation expansion
Concerning the three model ingredients, g, and f we assume the following H : > 0, to ensure that exposure actually occurs.
In addition, lim y↓0 f (y) = ∞, f (y) > y on (0, ∞) and for some δ ≥ 0 Fig. 1. See "Appendix A2" for an alternative. This form of f (y) is motivated by De Graaf et al. (2014) and implies that for an exposure occurring in an individual with immune status y < y c a large jump of immune level occurs during infection, while the increase in immune level is small if the immune status is higher than y c at exposure. We interpret the threshold y c as the immune level that distinguishes symptomatic and asymptomatic infection. In other words, an immune level y > y c provides protection against symptoms but nevertheless the encounter with the pathogen leads to a slight increase in immune level, whilst y < y c does not provide much protection and leads to a large boost of the immune level. H g : The function g describes the rate of waning of immunity between exposures and should therefore ensure that y(a) is a monotone decreasing (but positive) function. So let g : (0, ∞) → (−∞, 0) be such that the initial value problem dy da = g(y), y(a 0 ) = y 0 > 0 has a unique solution and lim a→∞ y(a) = 0, lim a→−∞ y(a) = ∞.
An alternative formulation of this assumption is that 1 g has a primitive T , say such that T (y) → ∞ for y ↓ 0 and T (y) → −∞ for y → ∞. The relationship between π and T is given by For later use we observe that An example is provided by with w > 0 a parameter. This choice for g(y) models an exponential decline in immune level between boosting events. In that case π(a, y 0 ) = e −wa y 0 (2.6) and In the generation expansion, we distinguish according to the number of hits of the Poisson process. So we start with the possibility of no hit at all.
So Q 0 t (y, ) describes the probability that an individual who has immune level y at age a and survives up to age a + t has had no exposures in the time interval [a, a + t] and has now an immune level in the set (e.g. within a given range [y low , y high ]). Of course, in this case we know exactly what the individual's immune level is, but as soon as an exposure does occur, the immune level at a + t has a range of possible values depending on when exactly the exposure occurred. Therefore, we now formulate an equation for the probability Q t (y, ), which is the probability that the individual has an immune state y in the set at age a + t given that it had immune status y at age a and without any restriction on the number of exposures since then.
If an infection does occur in [a, a + t], there has to be a first infection in this time window. The probability per unit of time that it occurs after exactly time σ equals e − σ . In that case, the immune status jumps to f (π(σ, y)) and there is time t − σ left before the clock reaches t. Accordingly Q should satisfy the RE In order to rewrite (2.9) in a more condensed symbolic form we introduce another kernel. The kernel B 1 t (y, ) describes the "position" on the y-axis immediately after the first jump (infection). Since it may happen that no jump occurs in the time interval of length t under consideration, this is not described by a probability distribution, but by a measure of total size less than one (this is sometimes called a "defective probability distribution"). The precise definition reads after the first infection belongs to |Y (a) = y) (2.10) In addition we define the product of two kernels as follows These definitions allow us to write (2.9) as Successive approximation amounts to substituting the right hand side for Q at the right hand side, yielding and then repeat this procedure again and again. By induction we define kernels B k for k ≥ 2: Explicitly we have y)))) ( )dσ 2 dσ 1 (2.15) and in general The interpretation reads The formulas (2.15) and (2.16) clearly show that all randomness derives from the hitting times of the Poisson process.
Returning to (2.13) and its "successors" obtained by repeating the approximation procedure, we are led to introduce, for k ≥ 1 and to observe that The upshot is that we define Cautionary note on notation: Q differs from Q 1 . Thus we constructed Q a (y b , ·), the distribution of the random variable Y (a), given its value y b at birth, on the basis of the three model ingredients π , and f . If we consider the limit t → ∞ for the kernel B t , the probability that no infection occurs vanishes. So describes the probability distribution immediately after the next infection, given that the current immune status equals y, when we do not restrict the length of the time interval under consideration. Note that indeed for all y we have B 1 ∞ (y, (0, ∞)) = 1, reflecting what was stated above: in an infinite time interval infection occurs with probability 1. With this kernel we associate the map K defined by The measure b describes the probability distribution of immune status of an individual directly after an infection event. The support of b is contained in the range of f , i.e. in [ f (y c ), ∞), since infection results in a new immune status obtained by applying f to the old immune status. Note that b is a probability measure, i.e., is positive and has total measure one. The map K preserves these properties. We call K the next-generation operator. Our aim in this section is to derive conditions on f and g that guarantee that K has a unique fixed point.
In order to derive a more handsome representation of K , we need some definitions. Since for every point in the range of f there are two points in the domain that are mapped to it, the inverse of f is double valued and we need notation to distinguish these two values from each other. For y 1 > y 2 , letτ (y 1 , y 2 ) be the time involved in decreasing, by waning, from immune level y 1 to y 2 . As it is convenient to haveτ also defined when y 1 > y 2 does not hold, we extend by zero.
We also introduce the notation The motivation for our special interest in sets of this form is the following. Measures in the range of K have their support in [ f (y c ), ∞). As our interest is in iterating K , we might as well restrict the domain of K to probability measures concentrated on The values that such measures take on sets of the form (3.6), for arbitrary y > 0, provide full information about the measure. Next note that Since (recall Fig. 1 and Definitions 3.1 and 3.2) we can perform the integration with respect to σ and rewrite (3.7) in the form The first term at the right hand side limits starting points for jumps to the right of y c , the second limits starting points to the left of y c . As the distinction is helpful, we emphasize it by defining and The range of K − is one-dimensional and spanned by an absolutely continuous measure. This reflects that, in order to have immune status less than y c when the next infection hits, the immune status must first wane to y c without any hit (the probability that this happens is e T (η) and determines the contribution to c). But once y c is "safely" reached, any information about η is irrelevant. After arrival in y c , the waiting time till being hit is exponentially distributed with parameter and the waiting time translates to an arrival position after the infection, as detailed in (3.12). Now, let us focus attention on For Combining (3.12) and (3.14) we see that a measure in the range of K has a density on [ f (y c ), f (2) (y c )]. The reason is similar as given above concerning K − : in order to lead to an arrival position in [ f (y c ), f (2) (y c )], the immune status has to wane to below f (y c ) and information about η becomes irrelevant upon passing f (y c ). We can subsequently use the explicit representation (3.15) to conclude that a measure in the range of K 2 has a density on [ f (y c ), f (3) (y c )]. Etcetera. It follows that a fixed point of K necessarily has a density.
Our aim is now to construct a fixed point of K . Apart from a multiplicative constant c, the fixed point is known when we restrict attention to the interval [ f (y c ), f (2) (y c )], cf. (3.12) and (3.14). The idea is to use (3.15) to extend the interval on which we know the fixed point and in the very end use the condition that the construction should yield a probability measure to determine the free constant c.
Assuming that b has a density φ we can, using (3.12) and (3.15), write the fixed The factor ∞ f −1 + (y) e T (η) φ(η)dη at the right hand side of (3.16) involves values of φ beyond y and this makes the use of (3.16) for extending the fixed point from the explicit (apart from an as yet unknown constant c) expressions provided by the right hand sides of (3.12) and (3.14) to larger values of y problematic. The solution is to differentiate (3.16) with respect to y and next combine the two identities to eliminate the factor. This is just a formula manipulation trick that has no biological interpretation whatsoever. We realize that unpleasant technical details are involved, but are not able to avoid these. To facilitate the formulation, we introduce some notation.
Differentation of (3.16) yields 1 which in combination with (3.16) itself leads to Equation (3.20) should hold for y > f (2) (y c ) and is supplemented by (3.17), which we rewrite as Now observe that for any integrable φ and for y Combining these two identities we deduce iii. α(y) Identity iii. allows us to rewrite the rewritten (3.20) as d dy Integrating the initial condition, see (3.17), we find that for y = f (2) (y c ) the quantity in the square brackets equals zero. By uniqueness, it equals zero for all y ≥ f (2) (y c ).
Using identity ii. above we see that this amounts to (3.16).
Standard contraction arguments (or, alternatively, monotone iteration arguments) guarantee that φ defined by (3.21) on [ f (y c ), f (2) (y c )] can be extended via (3.20) to an interval having f (2) (y c ) as its left end point. By continuation we obtain a maximal existence interval. Provided α is non-singular, this interval is [ f (y c ), ∞), since the equation is linear. A key point for us, however, is that φ should in fact be integrable over [ f (y c ), ∞). In other words, we want an a priori bound on the L 1 -norm.
To derive some intuition, let us briefly look at the constant coefficient homogeneous version This equation has a solution of the form φ(y) = e λy provided λ is a root of the characteristic equation (where the right hand side should be interpreted as αδ for λ = 0). Whenever αδ < 1, all roots of (3.23) have negative real part (see e.g. Chapter XI of Diekmann et al. (1995)). If αδ > 1, a positive real root exists. Clearly the corresponding solution is not integrable over the positive real axis. This example illustrates that we need some kind of condition on the behaviour of α(y) and f −1 + (y) for y → ∞. The probabilistic formulation is that we need tightness, meaning that mass is prevented from moving ever higher up the y-axis when we iterate K . The following consideration serves to build up intuition.
On average, the time in between successive jumps equals 1 . The condition f (π(1/ , y)) < y for large y (3.24) states that high-up the y-axis we lose immunity if we jump after the expected time.
As applying a nonlinear map and taking an expectation do not commute, this provides an idea, not a workable argument. Moreover, as we can also end up high on the y-axis by first waning to a very low y-level, we might need an additional and rather different condition to control this route to high immune levels. We now return to (3.20) and present sufficient conditions. Lemma 3.5 Assume that ρ and z exist, with 0 < ρ < 1 and z ≥ f (2) (y c ), such that Then the unique solution of (3.20), (3.21) is integrable over [ f (y c ), ∞).
Proof Integrating (3.20) from z to y > z we obtain Interchanging the two integrals in the first term at the right hand side leads to Since α is positive and, moreover, a derivative (cf. (3.18)) We claim that (3.25) implies that the right hand side is bounded by ρ. To verify this claim, we first write (3.25) as π(ρ/ , y) < f −1 + (y) and next use (2.3) to reformulate this inequality as Since T and T −1 are decreasing, it follows that and, finally Returning to the identity at the start of the proof, we conclude that Now note that, on account of (3.19), since T is positive on (0, y c ), cf. (2.2). Assumption (3.26) therefore guarantees that a constant C exists such that the right hand side, and hence the left hand side, is bounded by C for y ∈ [z, ∞). It follows that y z φ(ξ )dξ converges to a finite number for y → ∞.
For condition (3.26), the behaviour of f for y ↓ 0 matters too. We therefore focus on It follows that and consequently If ( is the object of interest. The above analysis leads to the following conclusion:

The stationary distribution
Throughout this section we assume that the conditions (3.25) and (3.26) appearing in Lemma 3.5 are satisfied. As demonstrated in Theorem 6.1 of Diekmann et al. (1998) (and likewise in Lemma's 3.7 and 3.8 of Diekmann et al. (2003)), there is a one-to-one correspondence between fixed points of the next-generation operator and stationary distributions of the process itself. Here the relationship is rather simple, since the jump rate does not depend on the position in the state space (or, in another jargon, the risk of encountering the pathogen does not depend on the immune level). In fact we have is a fixed point of the next-generation operator and, conversely, if b is a fixed point of the next-generation operator then m defined by is a stationary distribution.
The intuitive argument underlying Theorem 4.1 is that a steady distribution yields a steady stream of jumps with uniform departure rate and that f −1 relates the point of arrival to the point of departure. So the probability per unit of time of "landing" in after a jump equals m( f −1 ( )) and, with b the fixed point of K , this equals cb( ).
Since both b and m are probability measures, c must be equal to . For a formal proof see (Diekmann et al. 1998(Diekmann et al. , 2003 and take the limit t ↓ 0 in (6.10) of Diekmann et al. (1998) or (3.27) of Diekmann et al. (2003).

Corollary 4.2 There exists a unique normalized stationary distribution. This distribution has a density and the density takes strictly positive values on (0, ∞).
Our next aim is to establish the asymptotic stability of the stationary distribution by invoking results of Pichór and Rudnicki (2000). To do so, we need to make some preparations.
In the following we use the symbol m again to denote an arbitrary (so not necessarily a stationary) measure on (0, ∞). We define If m is a positive measure, so is Q t × m. In order to show that Q t × m is a probability measure if m is a probability measure, we need the following observation. Proof We claim that Clearly (4.4) is true for k = 0. Suppose (4.4) has been verified for k = n. Then π(s, z)), (0, ∞))ds and the claim is verified. Finally, If m has a density, does Q t × m have a density? Or, in other words, does the semigroup of operators associated with the kernel Q t leave the subspace of absolutely continuous measures invariant, so that we can associate with this kernel a semigroup of operators on L 1 (0, ∞)? The answer is probably affirmative without severe conditions on g and f . Below we shall derive a stronger result under a condition that guarantees the monotonicity of s → π(t − s, f (π(s, z))) so of the position at time t when starting in z, given that precisely one jump occurs, as a function of the time s at which the jump occurs. But first we pay attention to the situation with no jump at all.

Lemma 4.4
If m has a density φ, then Q 0 t × m has a density T 0 (t)φ defined by where in the last step we used the identity (2.4).
Proof π(s, z)))) and T −1 is decreasing, so it suffices to prove that is decreasing as well. The derivative with respect to s is f (π(s, z))g(π(s, z)) and (4.6) guarantees that this is negative.
The interpretation is that, by postponing the jump, you end up higher on the y-axis. And (4.6) is indeed an infinitesimal version of exactly this condition: Our motivation for this assumption derives from Proof A straightforward calculation reveals that Lemma 4.7 Assume that (4.6) holds. Denote the inverse function of η = η(s) = π(t − s, f (π(s, y))) by S = S(η) (where we have suppressed the dependence of S on t and y in the notation). Then Note that q 1 (t, y, η) ≥ 0 and in fact is strictly positive for π(t, f (y)) < η < f (π(t, y)) which in the limit t → ∞ amounts to η ∈ (0, ∞).
In the terminology of Pichór and Rudnicki (2000) {T (t)} is a Markov semigroup. From Corollary 4.2 we know that {T (t)} has a unique invariant density ψ and that ψ is positive. The fact that T (t) ≥ T 1 (t) and ||T 1 (t)|| = te − t > 0 for t > 0 guarantees that {T (t)} is partially integral as defined in Pichór and Rudnicki (2000). Thus we verified all assumptions of Theorem 2 in Pichór and Rudnicki (2000) (and in fact also of Theorem 2 in the preprint (Gerlach and Glück 2017)) and we obtain Theorem 4.11 Assume (4.6) holds and assume m has a density φ. Let T (t)φ be the density of Q t × m. Then where ψ is the density corresponding to the normalized stationary distribution.
But what happens if m does not have a density? The fact that the singular part of Q t × m converges to zero guarantees that in that case Q t × m converges tom in the total variation norm, wherem is the unique normalized stationary distribution. Indeed, we can write

Discussion
On the basis of a given force of infection , a given waning rate g and a given boosting map f , we derived a constructive description of how the distribution of immune status of an individual changes in time. Under the not-so-very-restrictive conditions (3.25), (3.26) and (4.6) we established that the distribution converges to a unique stationary distribution that has an (almost) everywhere positive density. Our approach builds on a long tradition as exposed in the seminal work (Lasota and Mackey 1994) of Lasota and Mackey. See (Rudnicki and Tyran-Kamińska 2015;Mackey et al. 2013;Lasota et al. 1992;Banasiak et al. 2012;Hille et al. 2016;Pichór and Rudnicki 2016). Also see the recent (Gerlach and Glück 2017), which builds on (and simplifies) work of Greiner (1982) and the somewhat older work of Heijmans (1986a, b) who uses spectral methods.
It might be interesting to explore yet another approach. Define  π(σ, y))))dσ (5. 2) The point is that one needs S in order to compute the expected time between two passages of y c (indeed, the time till the first jump after passing y c is exponentially distributed and the time of the jump determines from which y > y c the waning back towards y c starts). If one manages to prove that the expected time between two passages of y c is finite, a general result, viz. Theorem V.1.2, page 126, of Asmussen (1987) guarantees existence, uniqueness and asymptotic stability (with respect to the weak * topology) of a stationary distribution. So far we considered the force of infection as a parameter. We now formulate a consistency condition that should satisfy. It involves three new modeling ingredients: μ: the probability distribution of immune status at birth F(a): the probability to survive till at least age a ξ(y): the contribution to the force of infection when becoming infected while having immune status y In a stationary population the stable age distribution has densityF given bỹ The parameter should satisfy and, since we are interested in > 0, where the dependence of the right hand side on is hidden in the notation, but derives from the fact that Q depends on . The modeling approach presented here provides a basis for refining seroepidemiological methods. These methods exploit serum antibodies as biomarkers for infection and yield powerful tools for inferring the frequency of asymptomatic infections, those that cannot be observed directly (De Melker et al. 2006;Kretzschmar et al. 2010). Our work was motivated by questions arising in the epidemiology of pertussis, where serological data suggest that exposure to pertussis occurs much more frequently than documented in notification data of national surveillance systems. Whether this is due to underascertainment and underreporting, or to the fact that many infections are asymptomatic or mild and therefore not diagnosed, is unclear. To design vaccination strategies that protect young infants, it would be very useful to understand better how the boosting and waning of immunity in a population is reflected in serological surveys of the population and what the distribution of serological markers can tell us about the risk of exposure to the pathogen.
Current estimates of incidence of seroconversions only account for backward recurrence time of the most recent infection, ignoring any previous history of infections (Teunis et al. 2012). Carry-over from prior infections, leading to elevated baseline antibody levels, influence the current seroconversion, through the function f . The function f also determines whether infection causes a "small" or a "large" jump in antibody level, putatively corresponding to asymptomatic or symptomatic infection, respectively (De Graaf et al. 2014).
Thus, the relation between an infection-history dependent degree of protection against symptomatic infection and the incidence with which infection events occur may be exploited to better characterize transmission in the population, not only for pertussis but for any infectious pathogen with a measurable serum antibody response.
The statistical methods that have been developed to estimate the incidence of infection from serological data might be adapted and extended in order to incorporate a continuum immunity status as in the current paper. While thinking about the possibilities, one quickly realizes that our assumption of an age-independent force of infection is doubtful. And although incorporation of an age-specific force of infection in the formalism does not lead to great difficulties, such an extension of course leads to major complications when it comes to parametrization and identification. Here u(t, ·) is a continuous function of the form u(t, ·) = (0,∞) Q t (·, dz)ψ(z), while m(t, ·) is a density such that m(t, y)dy = (0,∞) φ(z)Q t (z, )dz for all ⊂ (0, ∞). One derives (A1.2) from (A1.1) by taking adjoints and restricting to absolutely continuous measures. This entails the assumption that one starts at t = 0 with a density. Indeed, as the Q 0 t (y, ) contribution to Q t (y, ) shows, the solution has a (diminishing) Dirac component when one starts with a Dirac initial condition.

A2 An alternative type of function f
The motivation for a graph of f as depicted in Fig. 1 derives from the within-host submodel introducted in De Graaf et al. (2014). The behaviour for y ↓ 0 and for y ↑ ∞ is questionable. We conclude with a graphical representation, Fig. 2, of an alternative. Actually the proofs given above remain valid when f (0) < ∞ but the behaviour of f (y) for y → ∞ is as described in H f .