Immunoepidemiology of a population
structured by immune status: a mathematical study of
waning immunity
and immune system boosting
Abstract
When the body gets infected by a pathogen the immune system develops pathogenspecific immunity. Induced immunity decays in time and years after recovery the host might become susceptible again. Exposure to the pathogen in the environment boosts the immune system thus prolonging the time in which a recovered individual is immune. Such an interplay of within host processes and population dynamics poses significant challenges in rigorous mathematical modeling of immunoepidemiology.
We propose a new framework to model SIRS dynamics, monitoring the immune status of individuals and including both waning immunity and immune system boosting. Our model is formulated as a system of two ODEs coupled with a PDE. After showing existence and uniqueness of a classical solution, we investigate the local and the global asymptotic stability of the unique diseasefree stationary solution. Under particular assumptions on the general model, we can recover known examples such as large systems of ODEs for SIRWS dynamics, as well as SIRS with constant delay. Moreover, a new class of SIS models with delay can be obtained in this framework.
KEYWORDS: Immunoepidemiology; Waning immunity; Immune status; Boosting; Physiological structure; Reinfection; Delay equations; Global stability; Abstract Cauchy problem
AMS Classification: 92D30; 35Q91; 47J35; 34K17
1 Introduction
Models of SIRS type are a traditional topic in mathematical epidemiology. Classical approaches present a population divided into susceptibles (S), infectives (I) and recovered (R), and consider interactions and transitions among these compartments. Susceptibles are those hosts who either have not contracted the disease in the past or have lost immunity against the diseasecausing pathogen. When a susceptible host gets in contact with an infective one, the pathogen can be transmitted from the infective to the susceptible and with a certain probability, the susceptible host becomes infective himself. After pathogen clearance the infective host recovers and becomes immune for some time, afterwards he possibly becomes susceptible again.
From the inhost point of view, immunity to a pathogen is the result of either active or passive immunization. The latter one is a transient protection and is due to the transmission of antibodies from the mother to the fetus through the placenta, thanks to which the newborn is immune for several months after birth McLean and Anderson (1988a). Active immunization is either induced by natural infection or can be achieved by vaccine administration Siegrist (2008); Goldsby et al. (2003).
Let us first consider the case of natural infection. A susceptible host, also called naive host, has a very low level of specific immune cells for a pathogen (mostly a virus or a bacterium, but possibly also a fungus). The first response to a pathogen is nonspecific, as the innate immune system cannot recognize the physical structure of the pathogen. The innate immune response slows down the initial growth of the pathogen, while the adaptive (pathogenspecific) immune response is activated. Clonal expansion of specific immune cells (mostly antibodies or CTL cells) and pathogen clearance follow. The population of pathogenspecific immune cells is maintained for long time at a level that is much higher than in a naive host. These are the socalled memory cells and are activated in case of secondary infection (see Fig. 1). Memory cells rapidly activate the immune response and the host mostly shows mild or no symptoms Wodarz (2007); Antia et al. (2005).
Each exposure to the pathogen might have a boosting effect on the population of specific memory cells. Indeed, the immune system reacts to a new exposure as it did during primary infection, thus yielding an increased level of memory cells Antia et al. (2005). Though persisting for long time after pathogen clearance, the memory cell population slowly decays, and in the long run the host might lose his pathogenspecific immunity. Wodarz (2007) writes that it “is unclear for how long hosts are protected against reinfection, and this may vary from case to case.”
Vaccineinduced immunity works similarly to immunity induced by the natural infection. Agents contained in vaccines resemble, in a weaker form, the diseasecausing pathogen and force a specific immune reaction without leading to the disease. If the vaccine is successful, the host is immunized for some time. Vaccinees experience immune system boosting and waning immunity, just as hosts recovered from natural infection do. In general, however, diseaseinduced immunity induces a much longer lasting protection than vaccineinduced immunity does Siegrist (2008). Waning immunity might be one of the factors which cause, also in highly developed regions, recurrent outbreaks of infectious diseases such as measles, chickenpox and pertussis. On the other side, immune system boosting due to contact with infectives prolongs the protection duration.
In order to understand the role played by waning immunity and immune system boosting in epidemic outbreaks, in the recent past several mathematical models were proposed. Few of these models describe only inhost processes which occur during and after the infection Wodarz (2007); Heffernan and Keeling (2008); De Graaf et al. (2014). Many more models, formulated in terms of ordinary differential equations (ODEs), consider the problem only at population level, defining compartments for individuals with different levels of immunity and introducing transitions between these compartments Dafilis et al. (2012); Heffernan and Keeling (2009). Possibly, also vaccinated hosts or newborns with passive immunity are included, and waning of vaccineinduced or passive immunity are observed Rouderfer et al. (1994); Mossong et al. (1999); Moghadas et al. (2008). In some cases, e. g. in the model by Heffernan and Keeling (2009), a discretization of the immune status leads to a very large system of ODEs. Simplified versions were suggested by Glass and Grenfell (2003); Arinaminpathy et al. (2012) who add the class W of hosts with waning immunity. Whosts can receive immune system boosting due to contact with infectives or move back to the susceptible compartment due to immunity loss. Lavine et al. (2011) extend the SIRWS model in Glass and Grenfell (2003), dividing each population into ageclasses of length of 6 months each and further classifying immune hosts (R) and hosts with waning immunity (W) by the level of immunity.
To describe the sole waning immunity process, authors have chosen delay differential equations (DDEs) models with constant or distributed delay Kyrychko and Blyuss (2005); Taylor and Carr (2009); Blyuss and Kyrychko (2010); Bhattacharya and Adler (2012); Yuan and Bélair (2013). The delay represents the average duration of the diseaseinduced immunity. However, neither a constant nor a distributed delay allows for the description of immune system boosting.
Martcheva and Pilyugin (2006) suggest an SIRS model in which infective and recovered hosts are structured by their immune status. In infective hosts the immune status increases over the course of infection, while in recovered hosts the immune status decays at some nonconstant rate. When the immune status has reached a critical level, recovered hosts transit from the immune to the susceptible compartment.
The goal of the present paper is to suggest a general framework for modeling waning immunity and immune system boosting in the context of population dynamics. We want to understand how the population dynamics, and in particular the number of infectives, affects waning immunity and immune system boosting in a recovered host and how, in turn, these two inhost processes influence the population dynamics. In contrast to the models proposed in Heffernan and Keeling (2009); Lavine et al. (2011), we shall maintain the number of equations as low as possible. For the sake of simplicity we do not include vaccination in this model.
We suggest a model in which the immune population is structured by the level of immunity, . Individuals who recover at time enter the immune compartment (R) with maximal level of immunity . The level of immunity tends to decay in time and when it reaches the minimal value , the host becomes susceptible again.
Immune system boosting is included by the mean of the probability , that an individual with immunity level moves to immunity level , when exposed to the pathogen. At the same time we choose the susceptible and the infective populations to be nonstructured. In this way we combine a PDE for the immune population with two ODEs for susceptible and infective hosts. The physiological structure in the immune population is the key innovation of the modeling approach, as it allows to describe at once both waning immunity (as a natural process which occurs as time elapses) and immune system boosting. There is hence no need to include any compartment for hosts with waning immunity as it was done, e. g. in Glass and Grenfell (2003); Arinaminpathy et al. (2012).
The paper is organized as follows. In Sect. 2 we carefully derive the model equations for the SIRS dynamics. This system will be referred to as model (M1). ODEs for susceptibles and infectives are easily introduced, whereas the PDE for the structured immune population is derived from balance laws. Though at a first glance the resulting immune hosts equation might resemble a sizestructured model, there is a crucial difference. In sizestructured models, indeed, transitions occur only in one direction, i. e., individuals grow in size but never shrink. On the other hand, our physiologically structured population is governed by a differential equation which includes transport (decay of immune status) and jumps (boosting to any higher immune status), describing movements to opposite directions. This makes the model analysis particularly challenging.
In Sect. 3 we consider fundamental properties of solutions for the model (M1). Writing the system as an abstract Cauchy problem in the statespace
,
we can guarantee existence of a unique classical solution. The proof requires a quite long computation which is postponed to the appendix. Further, using an appropriate comparison system, we show nonnegativity of solutions and determine an invariant subset in .
In Sect. 4 we consider stationary solutions. The proof of existence of stationary solution requires connection with the theory of Volterra integral equations.
We show that there exists a unique disease free equilibrium (DFE) and investigate its local and global stability. Investigation of endemic equilibria is considered in a following work.
Finally in Sect. 5 and Sect. 6 we show how to obtain from the general model (M1) various systems of ODEs, such as those in Glass and Grenfell (2003); Arinaminpathy et al. (2012); Heffernan and Keeling (2009); Lavine et al. (2011), and systems of equations with constant delay, such as those in Taylor and Carr (2009); Kyrychko and Blyuss (2005). Moreover we obtain as a special case a new class of SIS systems with delay.
2 A new class of models
In this section we derive the model equations for the SIRS dynamics. We start with the ODEs for susceptibles and infectives, which can be easily introduced, and continue with a PDE for the structured immune population, which we obtain from a discrete approach. The result is the general model (M1). We conclude the session with remarks on the total immune population.
2.1 Susceptibles and infectives
Let and denote the total population of susceptibles, respectively infectives, at time . In contrast to previous models which include shortterm passive immunity Rouderfer et al. (1994); McLean and Anderson (1988b); Moghadas et al. (2008), we assume that newborns are all susceptible to the disease. To make the model closer to real world phenomena we suppose that newborns enter the susceptible population at rate , dependent on the total population size . In general one could choose the natural death rate to be a function of the total population, too. However, for simplicity of computation, we consider the case of constant death rate . The death rate is assumed to be the same for all individuals, but infectives might die due to the infection, too.
Assumption 1.
Let with , be a nonnegative function, with .
To guarantee that there exists a nontrivial equilibrium , such that (see Fig. 2), we require the following.
Assumption 2.
There is an such that for we have that , whereas for we have . Further it holds that
When we include diseaseinduced death at rate , an equilibrium satisfies
(1) 
Contact with infectives (at rate ) induces susceptible hosts to become infective themselves. Infected hosts recover at rate , that is, is the average infection duration. Once recovered from the infection, individuals become immune, however there is no guarantee for lifelong protection. Immune hosts who experience immunity loss become susceptible again.
The dynamics of and is described by the following equations:
The term , which represents transitions from the immune compartment to the susceptible one, will be specified below together with the dynamics of the immune population.
2.2 Immune hosts
Let denote the density of immune individuals at time with immunity level . The total population of immune hosts is given by
The parameter describes the immune status and can be related to the number of specific immune cells of the host (cf. Sect. 1). The value corresponds to maximal immunity, whereas corresponds to low level of immunity.
We assume that individuals who recover at time enter the immune compartment (R) with maximal level of immunity . The level of immunity tends to decay in time and when it reaches the minimal value , the host becomes susceptible again. However, contact with infectives, or equivalently, exposure to the pathogen, can boost the immune system from to any higher immune status, see Fig. 3. It is not straightforward to determine how this kind of immune system boosting works, as no experimental data are available. Nevertheless, laboratory analysis on vaccines tested on animals or humans suggest that the boosting efficacy might depend on several factors, among which the current immune status of the recovered host and the amount of pathogen he receives Amanna et al. (2007); Luo et al. (2012). Possibly, exposure to the pathogen can restore the maximal level of immunity, just as natural infection does (we shall consider this special case in Sect. 3.1).
Let denote the probability that an individual with immunity level moves to immunity level , when exposed to the pathogen. Due to the definition of we have
and , for all .
As we effectively consider only immunity levels in the interval , we set
Then we have
Exposure to the pathogen might restore exactly the immunity level induced by the disease (). In order to capture this particular aspect of immune system boosting, we write the probability as the combination of a continuous () and atomic measures (Dirac delta):
where
 ,

is a continuously differentiable function and describes the probability that, due to contact with infectives, a host with immunity level boosts to the maximal level of immunity .
 ,

is a continuously differentiable function and describes the probability that, due to contact with infectives, a host with immunity level boosts to any other level , according to the continuous probability .
 ,

describes the probability that getting in contact with infectives, the host with immunity level does not experience immune system boosting.
The immunity level decays in time at some rate which is the same for all immune individuals with immunity level , that is,
Assumption 3.
Let be continuously differentiable.
The positivity of is given by the fact that, if for some value , there would be no change of the immunity level at , contradicting the hypothesis of natural decay of immune status.
In absence of immune system boosting, an infected host who recovered at time becomes again susceptible at time , where
The above equality becomes evident with an appropriate change of variables. Setting , we find . For the integration boundaries we have and . It follows that
To obtain a correct physical formulation of the equation for , we start from a discrete approach, as it could be done for agestructured or sizestructured models Webb (2008); Ellner (2009). The number of immune individuals in the immunity range is .
We track how many individuals enter and exit a small immunity interval in a short time .
Starting at a time with individuals, we want to describe the number of individuals at time .
Recall that the immunity level tends to decay, hence individuals enter the interval from and exit from (i. e., transition occurs the other way around than in agestructured or sizestructured models). As contact with infective individuals might boost the immune system, an immune host with immunity level in can move from this interval to any higher immunity level. For the same reason, immune individuals with any lower immunity level can “jump” to the interval . We assume that contacts between infectives and immune hosts occur at the same rate, , as between infectives and susceptibles.
Given any , denote by the partition of the interval into small intervals of length . Observe for all .
The total number of immune individuals with immunity level in at time is given by
Perform division by ,
and compute the limit (observe that also , as and are elements of the same partition) to get
Finally, we divide by and let . From the discrete approach derivation it becomes clear that the quantity , initially introduced in the equation at p. 2.1 to represent the hosts who experience immunity loss, is given by the number of immune hosts who reach the minimal level of immunity.
Altogether for we have the system of equations
(2)  
with initial values , coupled with a partial differential equation for the immune population,
(3) 
where , with boundary condition
(4) 
and initial distribution .
At a first glance equation (3) might resemble a sizestructured model, but there is an important difference. In sizestructured models transitions occur only in one direction, i. e., individuals grow in size but never shrink. The PDE (3) for the immune population is governed by a transport process (decay of immune status) and jumps (boosting to any higher immune status). The model analysis results much more challenging than the one of classical sizestructured models.
In the following we refer to the complete system (2) – (4) as to model (M1). An overview on the model notation is provided in Table 1.
Symbol  Description 

number of susceptibles at time  
number of infective hosts at time  
number of immune individuals at time  
density of immune individuals with immunity level at time  
total population () at time  
recruitment rate  
natural death rate  
diseaseinduced death rate  
transmission rate  
recovery rate  
natural decay of immunity  
probability that boosting raises immunity level to level  
maximal level of immune status in immune hosts  
minimal level of immune status in immune hosts  
initial distribution for 
2.3 Total number of immune individuals
From the PDE (3) – (4) we obtain the total number of immune individuals at time . Integrating the lefthand side of (3) in , we find
Similarly, integrating the righthand side of (3) we have
Altogether:
(5)  
When an immune host comes in contact with infectives, his immune system gets boosted so that either the maximal level of immunity or any other higher (or equal) level of immunity is restored. This means that the terms in the square bracket in (5) cancel out and the total immune population at time satisfies
(6) 
In other words, inflow at time into the immune class occurs by recovery of infected hosts, while outflow is either due to death or to immunity loss. Observe that the result (6) holds also in the special case of no boosting ( and , and in the case in which boosting always restores the maximal level of immunity ( and ).
3 Existence, uniqueness and positivity
Consider model (M1), with given initial data and for all . Let with its usual norm, and let be defined by
For given values of and , we define the map by
(7) 
with . In this notation the initial condition for is given by .
The state space for model (M1) is with the norm defined by
Observe that for wellposedness of the model (M1) the zero is not an element of . For all , we define . Then the model (M1) can be formulated as a perturbation of a linear abstract Cauchy problem,
(8)  
where and are respectively a linear and a nonlinear operator on , defined by
and
with .
The next result guarantees the existence and uniqueness of a classical solution. We remark that it is called mild solution of (8) the continuous solution of the integral equation
where is the semigroup on generated by . On the other side, a function is a classical solution of (8) on if is continuous on , continuously differentiable on , is in the domain of for all and (8) is satisfied on (see Pazy (1983)).
Theorem 3.1.
Proof. To have existence of a unique (classical) solution on we need to show:
(i) that is the generator of a semigroup on and (ii) that is continuously differentiable in (see (Pazy, 1983, Chap. 6)).
(i) The first hypothesis is easily verified as corresponds to the linear homogeneous part of the system and its domain is
Similar linear operators arising from population dynamics were considered in Webb (2008); Calsina and Farkas (2012).
(ii) Continuous differentiability of can be shown in two steps. First, for all , we determine the existence of the operator defined by
Second, we show that the operator is continuous in , that is
where is the operator norm. These two steps require the long computation included in the Appendix. ∎
Note that to have the existence of a classical solution one has to show continuous differentiability of . To have existence and uniqueness of a mild solution, as well as continuous dependence on the initial data, it is sufficient to prove Lipschitz continuity of (see Pazy, 1983).
From now on, we shall assume that all hypotheses of Theorem 3.1 hold.
Next we show that, given nonnegative initial data, model (M1) preserves nonnegativity. We proceed in two steps. First we introduce model (M2) in which boosting restores the maximal immune status. The PDE in model (M2) has the same characteristic curves as the PDE in model (M1), allowing us to use it as a comparison system in the proof of nonnegativity. There are other reasons to consider model (M2) separately. On the one side, the assumption of boosting restoring maximal immunity has been used in previous models, e.g. in Arinaminpathy et al. (2012). On the other side, in Sect. 6.2 we shall use model (M2) to obtain a new SIRS system with delay.
3.1 Model (M2): Boosting restores the maximal level of immunity
Let us assume that whenever an immune host comes in contact with the pathogen, his immune system is boosted in such a way that the maximal immunity level is restored. This means that
or equivalently,
This assumption modifies the equation (3) and the boundary condition (4) in model (M1) as follows
(9)  
(10) 
The equations for and in model (M1) remain unchanged. We shall refer to the system
(2), (9) – (10) as to model (M2).
Just for a moment, assume that for some the values and are known, and recall . For all let us define
and
Definition 1 (cf. Calsina and Saldaña (1995)).
We introduce the characteristic curve which identifies the cohort of individuals who recovered (hence have maximal level of immunity) at time . In this way we distinguish between those individuals who recovered before time and are already immune at the initial time, and those who recovered later than . Then the problem (9)–(10) can be solved along the characteristics (see, e.g Calsina and Saldaña, 1995; Webb, 2008) and we have the solution
(11) 
where the time is implicitly given by
As the death rate is bounded, we can extend the solution to all positive times (see also Calsina and Saldaña, 1995, Sec. 3.3). It is obvious that the solution is nonnegative for all .
3.2 Nonnegative solutions of model (M1)
Define the set
Theorem 3.2.
The cone is positively invariant for the model (M1).
Proof. We start with the infective population. Let be given. If for some , we have , hence the component is always nonnegative. Further we have that
It follows that, given positive initial values, the total population is larger than zero for all .
The equation for includes the term , which is given by the solution of the PDE (3)–(4). Let be given. Assume for all , hence (recall by Assumption 3). We see that implies
Hence, if , also does not leave the nonnegative cone.
To conclude the proof, we have to show that for all , .
First we show that strictly positive initial data lead to a strictly positive solution for all , . Assuming the contrary, there is a time such that
(i) for all ;
(ii) for all , whereas , for some .
In other words, is the first time at which for some value the solution of (3)–(4) is zero. Note that the PDE (3) in model (M1) and the PDE (9) have the same characteristic curves. By assumptions (i)(ii), along the characteristics we have the estimate
We can use the equations (9)–(10) as comparison system for (3)–(4) along characteristics. From (11) it is evident that the solution of (9)–(10) is going to be zero if and only if the characteristic curve associated to the point has a starting value (either or , ) equal to zero, contradicting to (i)(ii). We conclude that, given strictly positive initial data , the solution of (3)–(4) is strictly positive.
To complete the proof, we extend the result to nonnegative initial data. Let . We introduce a value and repeat the same argument as above for initial data . Finally we let . From the continuous dependence on initial data (cf. Theorem 3.1 and Pazy, 1983, Chap. 6), it follows that for . ∎
4 The disease free equilibrium
For investigation of stationary solutions of model (M1) we set the time derivative equal to zero and consider the problem
(12)  
(13)  