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 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 , correspond to the immune status at age a of an immortal individual. Here , the immune status at birth, is another parameter (heterogeneity can of course be captured by an assumed distribution of ). 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:Let y be an element of and let be a measurable subset of . Development over time of the distribution of immune status is described by the kernelwhich, 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.

Can we compute the distribution of an individual’s immune status at age a, given that it starts life with an immune level ?

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 .

The RE and its solution by generation expansion

Concerning the three model ingredients, g, and f we assume the following

, to ensure that exposure actually occurs.

is continuously differentiable and there exists such that for , , for . In addition, , on and for some Moreover, for . See 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 a large jump of immune level occurs during infection, while the increase in immune level is small if the immune status is higher than at exposure. We interpret the threshold as the immune level that distinguishes symptomatic and asymptomatic infection. In other words, an immune level provides protection against symptoms but nevertheless the encounter with the pathogen leads to a slight increase in immune level, whilst does not provide much protection and leads to a large boost of the immune level.

Fig. 1

The graph of f

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 be such that the initial value problem has a unique solution and , .

An alternative formulation of this assumption is that has a primitive T, saysuch that for and for . The relationship between and T is given byFor later use we observe thatAn example is provided bywith a parameter. This choice for g(y) models an exponential decline in immune level between boosting events. In that caseandIn 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.

Definition 2.1

So describes the probability that an individual who has immune level y at age a and survives up to age has had no exposures in the time interval and has now an immune level in the set (e.g. within a given range ). 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 has a range of possible values depending on when exactly the exposure occurred. Therefore, we now formulate an equation for the probability , which is the probability that the individual has an immune state y in the set at age 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 , there has to be a first infection in this time window. The probability per unit of time that it occurs after exactly time equals . In that case, the immune status jumps to and there is time left before the clock reaches t. Accordingly Q should satisfy the REIn order to rewrite (2.9) in a more condensed symbolic form we introduce another kernel. The kernel 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

Definition 2.2

In addition we define the product of two kernels as follows

Definition 2.3

These definitions allow us to write (2.9) asSuccessive approximation amounts to substituting the right hand side for Q at the right hand side, yieldingand then repeat this procedure again and again. By induction we define kernels for :Explicitly we haveand in generalwithThe interpretation readsThe 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 and to observe thatThe upshot is that we defineCautionary note on notation: Q differs from .

Thus we constructed , the distribution of the random variable Y(a), given its value at birth, on the basis of the three model ingredients , and f.

The next-generation operator and its fixed point

In this section we define an operator, which maps the distribution of immune status right after an infection event to the distribution right after the next such event. We call this the next-generation operator and we will show that, under certain conditions, this operator has a unique fixed point (and that the fixed point has a density). In the next section we shall use these results to derive the existence and uniqueness of a stationary distribution of immunity to which, under the influence of repeated episodes of waning interrupted by boosting events, a general distribution converges for t tending to infinity.

If we consider the limit for the kernel , the probability that no infection occurs vanishes. Sodescribes 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 , reflecting what was stated above: in an infinite time interval infection occurs with probability 1. With this kernel we associate the map K defined byThe 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 , 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.

Definition 3.1

Both of these functions are defined on . Both take the value in .

For , let be the time involved in decreasing, by waning, from immune level to . As it is convenient to have also defined when does not hold, we extend by zero.

Definition 3.2

Hence [cf. (2.2) and (2.3)]We also introduce the notationThe motivation for our special interest in sets of this form is the following. Measures in the range of K have their support in . 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 , provide full information about the measure.

Next note thatwithSince (recall Fig. 1 and Definitions 3.1 and 3.2)we can perform the integration with respect to and rewrite (3.7) in the formThe first term at the right hand side limits starting points for jumps to the right of , the second limits starting points to the left of . As the distinction is helpful, we emphasize it by definingwithandwithNote that for , cf. (2.2), so c is well-defined. The range of is one-dimensional and spanned by an absolutely continuous measure. This reflects that, in order to have immune status less than when the next infection hits, the immune status must first wane to without any hit (the probability that this happens is and determines the contribution to c). But once is “safely” reached, any information about is irrelevant. After arrival in , 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 we have and hence for all . ConsequentlyFor , on the other hand, we have for and consequentlyCombining (3.12) and (3.14) we see that a measure in the range of K has a density on . The reason is similar as given above concerning : in order to lead to an arrival position in , the immune status has to wane to below and information about becomes irrelevant upon passing . We can subsequently use the explicit representation (3.15) to conclude that a measure in the range of has a density on . 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 , 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 point problem in the formand with for defined byThe factor 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.

Definition 3.3

Differentation of (3.16) yields1 which in combination with (3.16) itself leads toEquation (3.20) should hold for and is supplemented by (3.17), which we rewrite asBefore analysing (3.20), let us check that no information was lost when deriving (3.20) from (3.16).

Lemma 3.4

Assume that is integrable over , satisfies (3.20) for and (3.21) for . Then (3.16) holds.

Proof

We write (3.20) in the formNow observe that for any integrable and for Combining these two identities we deduceIdentity iii. allows us to rewrite the rewritten (3.20) asIntegrating the initial condition, see (3.17), we find that for the quantity in the square brackets equals zero. By uniqueness, it equals zero for all . 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 can be extended via (3.20) to an interval having as its left end point. By continuation we obtain a maximal existence interval. Provided is non-singular, this interval is , since the equation is linear. A key point for us, however, is that should in fact be integrable over . In other words, we want an a priori bound on the -norm.

To derive some intuition, let us briefly look at the constant coefficient homogeneous versionThis equation has a solution of the form provided is a root of the characteristic equation(where the right hand side should be interpreted as for ). Whenever , all roots of (3.23) have negative real part (see e.g. Chapter XI of Diekmann et al. (1995)). If , 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 and for . 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 . The conditionstates 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 and , such thatandThen the unique solution of (3.20), (3.21) is integrable over .

Integrating (3.20) from z to we obtainInterchanging the two integrals in the first term at the right hand side leads toSince 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 and next use (2.3) to reformulate this inequality asSince T and are decreasing, it follows thatand, finallyReturning to the identity at the start of the proof, we conclude thatand henceNow note that, on account of (3.19),since T is positive on , 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 . It follows that converges to a finite number for .

The conditions (3.25) and (3.26) are somewhat implicit, so let us formulate more explicit conditions that imply them.

Lemma 3.6

Assume that , v and z exist, with , and , such thatandwith as introduced in , then (3.25) holds.

Let be such thatIncreasing z if necessary, we may assume thatFrom (3.27) it follows thatHence

Corollary 3.7

If we can satisfy (3.25) for any in and any choice of .

For condition (3.26), the behaviour of f for matters too. We therefore focus onwhich was derived in De Graaf et al. (2014).

Lemma 3.8

Let f be given by (3.29).

If , condition (3.26) is fulfilled.

If, for some , for and condition (3.26) is fulfilled.

For f given by (3.29) we haveIt follows thatFor small enough y we haveand consequentlyIf (3.30) holds, the exponent is less than and boundedness of suffices to fulfil (3.26). If for all y, for and the fact that the exponent is negative guarantees that (3.26) holds.

Before summarizing the results of this section, let us elucidate the role of the constant c in (3.20), (3.21). To compute the density of the fixed point of K, first solve (3.20) and (3.21) with . When (3.25) and (3.26) hold, the integralis well-defined. The renormalized functionis the object of interest. The above analysis leads to the following conclusion:

Theorem 3.9

Provided the conditions (3.25) and (3.26) are satisfied, the next-generation operator K has a unique fixed point. This fixed point has a density (3.31), where is the solution of (3.20), (3.21) corresponding to the choice , or any other choice. Hence for .

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

Theorem 4.1

If m is a stationary distribution then b defined byis a fixed point of the next-generation operator and, conversely, if b is a fixed point of the next-generation operator then m defined byis 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 relates the point of arrival to the point of departure. So the probability per unit of time of “landing” in after a jump equals and, with b the fixed point of K, this equals . Since both b and m are probability measures, c must be equal to . For a formal proof see (Diekmann et al. 1998, 2003) and take the limit 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 .

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 . We defineIf m is a positive measure, so is . In order to show that is a probability measure if m is a probability measure, we need the following observation.

Lemma 4.3

, .

We claim thatClearly (4.4) is true for . Suppose (4.4) has been verified for . Thenand the claim is verified. Finally,

If m has a density, does have a density? Or, in other words, does the semigroup of operators associated with the kernel leave the subspace of absolutely continuous measures invariant, so that we can associate with this kernel a semigroup of operators on ? 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 ofso 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 has a density defined by

where in the last step we used the identity (2.4).

The condition alluded to above reads(We provide the biological interpretation of (4.6) after the next lemma.)

Lemma 4.5

If (4.6) holds, the functionis, for every given and , an increasing function on [0, t].

and is decreasing, so it suffices to prove thatis decreasing as well. The derivative with respect to s isand (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

Lemma 4.6

Withandcondition (4.6) holds.

A straightforward calculation reveals that

Lemma 4.7

Assume that (4.6) holds. Denote the inverse function ofby (where we have suppressed the dependence of S on t and y in the notation). Thenwith

Note that and in fact is strictly positive forwhich in the limit amounts to .

Corollary 4.8

Assume (4.6) holds. has densityIn particular, if m has a density , then has density defined by(we say that is a kernel operator).

Lemma 4.9

If then for where is inductively defined by

So if (4.10) holds, then

By comparing (4.10) and (4.4) we deduce that

Corollary 4.10

Assume (4.6) holds. If m has a density , then has a density, say . Moreover, for each the map from to is continuous and (recall (4.5)) is a kernel operator.

Since translation is continuous in , we deduce from (4.5) that as . Combining (4.10) and (4.11) we find that as (in fact the stronger result as holds). Since is a semigroup, continuity at implies continuity.

In the terminology of Pichór and Rudnicki (2000) is a Markov semigroup. From Corollary 4.2 we know that has a unique invariant density and that is positive. The fact that and for guarantees that 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 be the density of . Thenwhere 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 converges to zero guarantees that in that case converges to in the total variation norm, where is the unique normalized stationary distribution. Indeed, we can writewhere is the absolutely continuous part of and the singular part. Then and hence . Moreover, for . Henceand by choosing s sufficiently large the last two terms can be made arbitrarily small, while the first term can be made arbitrarily small by choosing t sufficiently large.

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. DefineThen S should satisfy the equationThe point is that one needs S in order to compute the expected time between two passages of (indeed, the time till the first jump after passing is exponentially distributed and the time of the jump determines from which the waning back towards starts). If one manages to prove that the expected time between two passages of 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: In a stationary population the stable age distribution has density given byThe parameter should satisfyand, since we are interested in ,where the dependence of the right hand side on is hidden in the notation, but derives from the fact that Q depends on .

the probability distribution of immune status at birth

the probability to survive till at least age a

the contribution to the force of infection when becoming infected while having immune status y

The modeling approach presented here provides a basis for refining sero-epidemiological 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.