# Temporal Evolution of Immunity Distributions in a Population with Waning and Boosting

Authors: M. V. Barbarossa; M. Polner; G. Röst
Journal: Complexity
 (2018)

---

## Abstract

We investigate the temporal evolution of the distribution of immunities in a population, which is determined by various epidemiological, immunological, and demographical phenomena: after a disease outbreak, recovered individuals constitute a large immune population; however, their immunity is waning in the long term and they may become susceptible again. Meanwhile, their immunity can be boosted by repeated exposure to the pathogen, which is linked to the density of infected individuals present in the population. This prolongs the length of their immunity. We consider a mathematical model formulated as a coupled system of ordinary and partial differential equations that connects all these processes and systematically compare a number of boosting assumptions proposed in the literature, showing that different boosting mechanisms lead to very different stationary distributions of the immunity at the endemic steady state. In the situation of periodic disease outbreaks, the waveforms of immunity distributions are studied and visualized. Our results show that there is a possibility to infer the boosting mechanism from the population level immune dynamics.


---

## Body

## 1. Introduction

The outcome of an infection within an individual host depends on the specific pathogen and the status of the immune system of the host. At a larger scale, the outcome of an epidemic in a population is influenced by the ensemble of individual immunities. There are a number of processes in play that determine how these immunities change in time. Upon recovery from infection, some immune memory remains, which may persist for long time after pathogen clearance. Eventually, memory cells slowly decay, and in the long run, recovered hosts could lose pathogen-specific immunity [1]. Waning immunity is possibly one of the contributing factors which cause, in particular in highly developed regions, recurrent outbreaks of infectious diseases such as chickenpox and pertussis. Immune memory can be boosted due to repeated exposure to the pathogen thus prolonging the time during which immune hosts are protected. Our goal in this paper is to monitor the distributions of immune memories in a population and track their temporal evolution, which provides very important insights about the interplay of individual and population level disease dynamics.In theSIR framework, a population of hosts is divided into susceptibles (S), infectives (I), and recovered (R), and interactions among individuals from the different compartments are considered. Susceptibles are those hosts who either have not contracted the disease in the past or have lost immunity against the disease-causing pathogen. When a susceptible host gets in contact with an infective, the pathogen can be transmitted and the susceptible host may become infective. After the loss of immunity, an individual from compartment R transits back to the susceptible compartment; hence when waning of immunity is included, the model is called SIRS.To account for immune system boosting, we structure the immunes according to their level of immunity. The high complexity of the immune status of an individual is simplified into a single parameter that reflects the strength of immunity in the sense that it indirectly indicates the duration of immunity until waning. The challenge due to secondary or multiple exposures to the pathogen initiates an immune response that results in a higher level of immunity.The details are yet unclear how exactly the immune response and in particular a boost of the immune system work [1]; most likely, there is a range of mechanisms underlying these processes that are specific to each host and pathogen [2]. Laboratory analysis on vaccines tested on animals or humans suggests 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 [2, 3]. In previous mathematical models, it was assumed that a boost restores the maximal immune status, the same as after natural infection [4]. Few authors have assumed that at each new contact with a known pathogen, the immune system is boosted a little, with a small increase in memory cells [5]. Others have considered a combination of both possibilities [6].To investigate the temporal evolution of how the distribution of these immunity levels changes in the population, we propose a mathematical modeling framework along the lines of our previous works [7, 8]. The core of the model is a hybrid system of equations of SIRS type, in which the immune population is structured by the level of immunity, whereas the susceptible and the infective populations are nonstructured. In [8], we investigated the well-posedness of the general model and its basic qualitative properties, whereas in [7], we considered a special case of the hybrid system in form of delay differential equations (DDEs) with constant and distributed delay. Here, we focus on the immune response identifying several possible scenarios for immune system boosts.We systematically compare a number of assumptions on the boosting mechanism that have been used in the literature. The goal is to observe the effects of different immune responses not only at an individual but also at population level, such as the stationary distribution of immunities in the case of an endemic disease and the periodic change of these distributions in case of repeated disease outbreaks. To the best of our knowledge, there are no further studies about the temporal evolution of the distribution of immune statuses in a population, and the only paper that has predicted a specific immunity distribution from a mathematical model is [6].In Section2, we provide the details of the mathematical model and present the numerical methods that are used for its solution. Four scenarios for immune boosting mechanisms are described and compared in Section 3, together with the parametrizations that lead to an either stable endemic state or oscillatory disease outbreaks. A careful and systematic numerical analysis is performed to understand the effect of various boosting mechanisms, and the results are summarized and discussed in Section 4.


## 2. Materials and Methods

### 2.1. The Mathematical Model

LetSt and It denote the total population of susceptibles and infectives, respectively, at time t. The total population shall be assumed to be in balance (hence normalized) Nt≡1, with birth rate equal to the natural death rate d≥0 and no disease-induced death. We assume that newborns are all susceptible.Contact with infectives (at rateβI) induces susceptible hosts to become infective themselves. Infected hosts recover at rate γ>0, that is, 1/γ is the average infection duration. Once recovered from the infection, individuals become immune; however, there is no guarantee for life-long protection. Immune hosts who experience immunity loss become susceptible again.Letrt,z denote the density of immune individuals at time t with immunity level z∈zmin,zmax,0≤zmin<zmax<∞. The total population of immune hosts is given by
(1)Rt=∫zminzmaxrt,zdz.The parameterz describes the immune status and can be related to the number of specific immune cells of the host. The value zmax corresponds to the maximal immunity, whereas zmin corresponds to the lowest level of immunity for hosts in the R compartment. We assume that individuals who recover at time t enter the immune compartment with maximal level of immunity zmax. The level of immunity tends to decay in time and when it reaches the lower threshold zmin, the host becomes susceptible again. However, contact with infectives, or equivalently, exposure to the pathogen, can boost the immune system from z∈zmin,zmax to any higher immune status (see Figure 1).Figure 1
Sketch of the mathematical models (3), (4), and (5). Susceptible hosts become infective after pathogen transmission. Infected hosts who recover enter the immune compartment R which is structured by the level of immunity. Natural infection induces the maximal level of immunity zmax. Immunity decays in time and when the immune status reaches the minimal value zmin; the recovered host becomes susceptible again. Meanwhile, exposure to the pathogen can boost the immune system and prolong protection.Given a host withinitial immune status v∈zmin,zmax, let us denote by Zv∈zmin,zmax the updated immune status, which is achieved after new contact with the pathogen. The updated immune status Zv is modeled as a random variable taking values in Sv≔v,zmax⊆zmin,zmax⊆ℝ with probability density function p·,v. The value pz,v represents the relative likelihood that a boost from immune level v to level z occurs. Secondary exposures to the pathogen might have no effects on the host’s immune system or might restore the immunity level induced by the disease (zmax). In order to capture these particular aspects, in Section 3.1, we shall also consider limit cases in which the probability density function is a Dirac measure centered either on the current immune status or on the maximal level zmax.The immunity level decays in time at rategz, with g positive, smooth, and bounded, which is the same for all immune individuals with immunity level z. In the absence of immune system boosting, an infected host who recovered at time t0 becomes again susceptible at time t0+τ, where
(2)τ=∫zminzmax1gxdx.With the above assumptions, we obtain fort>0 the following system of equations
(3)S′t=d1−St−βStIt+gzminrt,zmin,I′t=βStIt−γ+dIt,with initial values S0=S0>0, I0=I0≥0, coupled with a partial differential equation (PDE) for the immune population,
(4)∂∂trt,z−∂∂zgzrt,z=−drt,z+βIt∫zminzpz,vrt,vdv−rt,z,for z∈zmin,zmax, with boundary condition
(5)gzmaxrt,zmax=γIt,and initial distribution r0,z=ψz,z∈zmin,zmax. A sketch of the full model is given in Figure 1, whereas some possible boosting mechanisms are depicted in Figure 2, and they will be discussed later in detail. The formal derivation of a slightly more general version of models (3), (4), and (5) with variable total population size and disease-induced death is given in [8].Figure 2
Exposure to the pathogen has a boosting effect on the immune system, whereby several possible scenarios are possible: (a) No boosting events (NOboost): the host who recovered at timet becomes susceptible at time t+τ, with τ>0 given in (2). (b) Boosts restore disease-induced immunity, that is, the immune system is always boosted to the maximal level of immunity (MAXboost). (c) Variable boost, to any higher immune level.
(a)
No boosts
(b)
Boost to maximal level
(c)
Boost to any higher level


### 2.2.R0

Before presenting further results, we introduce the basic reproduction numberR0 of models (3), (4), and (5),
(6)R0=βγ+d,which indicates the average number of secondary infections generated in a fully susceptible population by one infected host over the course of his infection. The basic reproduction number is a reference parameter in mathematical epidemiology used to understand if, and in which proportion, the disease will spread among the population.


### 2.3. Numerical Solution of the Hybrid System

In the following, we outline the numerical method used to solve the hybrid system (3), (4), and (5).Consider (3), (4), and (5) as a one-dimensional, first-order nonlinear PDE system. In this section, we refer to the independent variables t,z as time and space, respectively. As S and I do not depend on z, but only on time, the initial conditions are
(7)S0,z≡S0≡S0>0,I0,z≡I0≡I0≥0,r0,z=ψz,z∈zmin,zmax.Sincegz>0 for all z∈zmin,zmax, the boundary condition for r is imposed in (5) at the inflow boundary, that is, at z=zmax. All together we have an initial-boundary value problem (IBVP). For the numerical integration of the IBVP, we employ the MATLAB code hpde [9], developed to solve IBVPs for first-order systems of hyperbolic PDEs in one space variable and time. The hpde routine implements Richtmyer’s two-step variant of the Lax-Wendroff method [10, 11], which is well established to solve hyperbolic PDEs. This scheme is explicit and second-order accurate in both space and time. To compute the numerical solution of our IBVP, we have modified the code hpde so that the integral term in (4) is efficiently implemented, preserving the second-order accuracy of the Lax-Wendroff scheme. In Section 3.3, we have verified the spatial order of accuracy for a number of computations.In the computation of the numerical solutions of the IBVP, we used an equidistant meshzmin=z^0<z^1<⋯<z^M=zmax,Δz^=z^i−z^i−1,i=1,…,M. The initial function on this mesh is given by S0,I0,r0,z^i=S0,I0,ψz^i.For an explicit scheme for a hyperbolic system, a necessary condition for stability is the Courant-Friedrichs-Lewy condition, which for our system means(8)ΔtΔz^maxz^∈z^0,…,z^Mgz^<1.We define theboosting matrix Π as the numerical discretization of the probability density function p on the equidistant mesh Δz^. Each entry πij=Πi,j of the boosting matrix represents the probability that z^i is the updated immune level, given initial immune level z^j, that is, πij=pz^i,z^j. It follows that πij∈0,1, with πij=0 if j>i and
(9)∑i=1Mπij=∑i=jMπij=1,forallj=1,…,M.


## 2.1. The Mathematical Model

LetSt and It denote the total population of susceptibles and infectives, respectively, at time t. The total population shall be assumed to be in balance (hence normalized) Nt≡1, with birth rate equal to the natural death rate d≥0 and no disease-induced death. We assume that newborns are all susceptible.Contact with infectives (at rateβI) induces susceptible hosts to become infective themselves. Infected hosts recover at rate γ>0, that is, 1/γ is the average infection duration. Once recovered from the infection, individuals become immune; however, there is no guarantee for life-long protection. Immune hosts who experience immunity loss become susceptible again.Letrt,z denote the density of immune individuals at time t with immunity level z∈zmin,zmax,0≤zmin<zmax<∞. The total population of immune hosts is given by
(1)Rt=∫zminzmaxrt,zdz.The parameterz describes the immune status and can be related to the number of specific immune cells of the host. The value zmax corresponds to the maximal immunity, whereas zmin corresponds to the lowest level of immunity for hosts in the R compartment. We assume that individuals who recover at time t enter the immune compartment with maximal level of immunity zmax. The level of immunity tends to decay in time and when it reaches the lower threshold zmin, the host becomes susceptible again. However, contact with infectives, or equivalently, exposure to the pathogen, can boost the immune system from z∈zmin,zmax to any higher immune status (see Figure 1).Figure 1
Sketch of the mathematical models (3), (4), and (5). Susceptible hosts become infective after pathogen transmission. Infected hosts who recover enter the immune compartment R which is structured by the level of immunity. Natural infection induces the maximal level of immunity zmax. Immunity decays in time and when the immune status reaches the minimal value zmin; the recovered host becomes susceptible again. Meanwhile, exposure to the pathogen can boost the immune system and prolong protection.Given a host withinitial immune status v∈zmin,zmax, let us denote by Zv∈zmin,zmax the updated immune status, which is achieved after new contact with the pathogen. The updated immune status Zv is modeled as a random variable taking values in Sv≔v,zmax⊆zmin,zmax⊆ℝ with probability density function p·,v. The value pz,v represents the relative likelihood that a boost from immune level v to level z occurs. Secondary exposures to the pathogen might have no effects on the host’s immune system or might restore the immunity level induced by the disease (zmax). In order to capture these particular aspects, in Section 3.1, we shall also consider limit cases in which the probability density function is a Dirac measure centered either on the current immune status or on the maximal level zmax.The immunity level decays in time at rategz, with g positive, smooth, and bounded, which is the same for all immune individuals with immunity level z. In the absence of immune system boosting, an infected host who recovered at time t0 becomes again susceptible at time t0+τ, where
(2)τ=∫zminzmax1gxdx.With the above assumptions, we obtain fort>0 the following system of equations
(3)S′t=d1−St−βStIt+gzminrt,zmin,I′t=βStIt−γ+dIt,with initial values S0=S0>0, I0=I0≥0, coupled with a partial differential equation (PDE) for the immune population,
(4)∂∂trt,z−∂∂zgzrt,z=−drt,z+βIt∫zminzpz,vrt,vdv−rt,z,for z∈zmin,zmax, with boundary condition
(5)gzmaxrt,zmax=γIt,and initial distribution r0,z=ψz,z∈zmin,zmax. A sketch of the full model is given in Figure 1, whereas some possible boosting mechanisms are depicted in Figure 2, and they will be discussed later in detail. The formal derivation of a slightly more general version of models (3), (4), and (5) with variable total population size and disease-induced death is given in [8].Figure 2
Exposure to the pathogen has a boosting effect on the immune system, whereby several possible scenarios are possible: (a) No boosting events (NOboost): the host who recovered at timet becomes susceptible at time t+τ, with τ>0 given in (2). (b) Boosts restore disease-induced immunity, that is, the immune system is always boosted to the maximal level of immunity (MAXboost). (c) Variable boost, to any higher immune level.
(a)
No boosts
(b)
Boost to maximal level
(c)
Boost to any higher level


## 2.2.R0

Before presenting further results, we introduce the basic reproduction numberR0 of models (3), (4), and (5),
(6)R0=βγ+d,which indicates the average number of secondary infections generated in a fully susceptible population by one infected host over the course of his infection. The basic reproduction number is a reference parameter in mathematical epidemiology used to understand if, and in which proportion, the disease will spread among the population.


## 2.3. Numerical Solution of the Hybrid System

In the following, we outline the numerical method used to solve the hybrid system (3), (4), and (5).Consider (3), (4), and (5) as a one-dimensional, first-order nonlinear PDE system. In this section, we refer to the independent variables t,z as time and space, respectively. As S and I do not depend on z, but only on time, the initial conditions are
(7)S0,z≡S0≡S0>0,I0,z≡I0≡I0≥0,r0,z=ψz,z∈zmin,zmax.Sincegz>0 for all z∈zmin,zmax, the boundary condition for r is imposed in (5) at the inflow boundary, that is, at z=zmax. All together we have an initial-boundary value problem (IBVP). For the numerical integration of the IBVP, we employ the MATLAB code hpde [9], developed to solve IBVPs for first-order systems of hyperbolic PDEs in one space variable and time. The hpde routine implements Richtmyer’s two-step variant of the Lax-Wendroff method [10, 11], which is well established to solve hyperbolic PDEs. This scheme is explicit and second-order accurate in both space and time. To compute the numerical solution of our IBVP, we have modified the code hpde so that the integral term in (4) is efficiently implemented, preserving the second-order accuracy of the Lax-Wendroff scheme. In Section 3.3, we have verified the spatial order of accuracy for a number of computations.In the computation of the numerical solutions of the IBVP, we used an equidistant meshzmin=z^0<z^1<⋯<z^M=zmax,Δz^=z^i−z^i−1,i=1,…,M. The initial function on this mesh is given by S0,I0,r0,z^i=S0,I0,ψz^i.For an explicit scheme for a hyperbolic system, a necessary condition for stability is the Courant-Friedrichs-Lewy condition, which for our system means(8)ΔtΔz^maxz^∈z^0,…,z^Mgz^<1.We define theboosting matrix Π as the numerical discretization of the probability density function p on the equidistant mesh Δz^. Each entry πij=Πi,j of the boosting matrix represents the probability that z^i is the updated immune level, given initial immune level z^j, that is, πij=pz^i,z^j. It follows that πij∈0,1, with πij=0 if j>i and
(9)∑i=1Mπij=∑i=jMπij=1,forallj=1,…,M.


## 3. Results and Discussion

### 3.1. Scenarios for the Immune Response

Let us consider an immune host who has immune statusv at the moment of reexposure to the pathogen. We shall investigate the following possible scenarios for immune response in case of secondary exposure to the pathogen.


#### 3.1.1. NOboost

Assume that the immune system does not respond to reinfection, that is, we observe only waning of immunity after recovery (Figure2(a)). This corresponds to the limit case in which the boosting probability function is simply a Dirac measure with support on the initial immune level. It follows that (4) reduces to
(10)∂∂trt,z−∂∂zgzrt,z=−drt,z.Recall the definition ofτ>0 in (2). The transport equation (10) with boundary condition (5) is solved along characteristics and we obtain
(11)gzminrt,zmin=γIt−τe−dτ,meaning that τ time after recovery immune hosts who did not die become susceptible again. In turn, we find a delay term in the equation for S and have a classical SIRS model with constant delay (cf. [12])
(12)S′t=d−βStIt−dSt+γIt−τe−dτ,I′t=βStIt−γ+dIt,R′t=γIt−γIt−τe−dτ−dRt,where Rt is the total immune population at time t as in (1).


#### 3.1.2. MAXboost

Assume that at any new encounter with the pathogen the immune system of a recovered host is boosted in such a way that the disease-induced (maximal) immunity is restored (Figure2(b)). This corresponds to the limit case in which the boosting probability function is a Dirac measure with support on the maximal immune level (zmax). As all boosted individuals are transferred to the boundary, and the integral term in (4) moves to the boundary condition, yielding
(13)∂∂trt,z−∂∂zgzrt,z=−drt,z−βItrt,z,(14)gzmaxrt,zmax=γIt+βItRt.The BVP (13) and (14) can be solved along characteristics and one obtains for t≥0 a system of delay equations with constant and distributed delay:
(15)S′t=d1−St−βItSt+It−τγ+βRt−τe−dτ−β∫−τ0It+udu,I′t=βItSt−γ+dIt,R′t=−dRt+γIt−It−τγ+βRt−τe−dτ−β∫−τ0It+udu,with τ>0 as defined in (2) and with given initial functions ϕSt≥0, ϕIt≥0, and ϕRt≥0, such that ϕSt+ϕIt+ϕRt≡1, for all t∈−τ,0. The formal derivation of system (15) is given in [7, 8].


#### 3.1.3. ANYboost

We assume that the immune system is boosted to any better immunity level with uniform probability. That is,pz,v=p¯>0 with p¯=1/zmax−v, for all z∈v,zmax. The boosting matrix for the numerical implementation is shown in Figure 3(a).Figure 3
Nonmonotonic immune response mechanisms. (a, c, d) Examples of different boosting matricesΠ for the numerical implementation of pz,v, the probability that individuals with immune level v at exposure are boosted to level z≥v after exposure. The immune status is represented in zmin,zmax=0,1 and discretized on an equidistant mesh with 30 points. (b) The dotted curve represents the boosting mechanism inspired by [6], the solid curve the one inspired by [14].
(a)
ANYboost
(b)
Nonmonotonic boosting functions
(c)
ODboost
(d)
HKboost


#### 3.1.4. HKboost and ODboost

Previous works have proposed nonmonotonic functions for modeling the immune response to secondary exposure. In [6, 13], a mathematical model for in-host dynamics during measles infection was proposed. The model explicitly considers the relation between the immune level at time of exposure and the level of memory cells after exposure. Heffernan and Keeling [6, 13] suggest that this relation is nonmonotonic. Indeed, although the level of immunity after exposure is always greater than the level at the time of exposure, the boosted level starts high, decreases, and then increases again. We shall denote by HKboost the boosting function from [6, 13]. For the numerical implementation of HKboost, we first extract values from Figure 2(d) in [6], which shows the relation between the initial and the updated immune status. Then, we normalize the immune status interval (zmin,zmax=0,1) and extend the boosting function to [0.75, 1], assuming that for high initial immune status, the immune system is always boosted to the maximal level (see Figure 3(b)). In terms of our model coefficients, for any initial immunity v, the atomic measure of the updated immunity is concentrated on fv, where f is the boosting function represented by the dotted curve in Figure 3(b). The corresponding boosting matrix is shown in Figure 3(d).Immune boosting in pertussis was considered in [14], where a nonmonotonic boosting function similar to the one in [6, 13] was proposed. The relation between the level of memory cells before (v) and after (z) reexposure is governed by the function
(16)z=fv=v1+15v5,withv∈0,5.De Graaf et al. [14] suggest that a small jump from v to z corresponds to a mild infection that causes no harm but only boosts the antibody level of the individual, whereas a large jump corresponds to a severe infection that may also cause disease. We shall use the boosting function in (16) for simulations, previous normalization to zmin,zmax=0,1, and small modifications which allow to have boundedness. In detail, we assume that for very low initial immunity (0≤v<0.075) and for large initial immunity (0.8<v≤1), the immune status is always boosted to the maximal level z=zmax=1, whereas for v∈0.075,0.8, the boosting probability is atomic along the graph of the boosting function f in (16), as above. The resulting boosting function is given by the solid curve in Figure 3(b), and the corresponding boosting matrix is shown in Figure 3(c).The two mechanisms for immune response HKboost and ODboost are governed by rather similar boosting mechanisms, and we shall see that the solutions behave accordingly.


### 3.2. Parameter Values

For the numerical computations below, we set the birth rate, natural death rated=0.02, and the initial conditions ψ0z=0.05 for z∈zmin,zmax=0,1, I0=0.01, and S0=0.94. This means that we assume that 1% of the total population is infectious at the beginning of our observations, while 5% is immune to the pathogen and the level of immunity is equally distributed.Concerning the disease dynamics, we setγ=3 and gz=0.5, for all z∈zmin,zmax, corresponding to τ=2 when the model can be reduced to DDE.We allow the basic reproduction numberR0 as given in (6) to vary, in order to show both solutions that converge to an endemic equilibrium and solutions that produce periodic oscillations.


### 3.3. Stable Endemic Equilibrium

When the basic reproduction number is sufficiently small, the solution converges to an endemic equilibrium for all the boosting mechanisms presented in Section3.1. We compare in Figure 4 the numerical solution corresponding to different boosting mechanisms. In Figure 4(a), we show the component I of the solution, which indicates how the number of infective evolves in time and rapidly converges to an equilibrium. Changing the boosting mechanism has small effects on the infective population: numerical solutions show the same qualitative behavior, with endemic equilibria deviating less than 1% from each other. In particular, we see that the nonmonotonic boosting mechanisms, ODboost, HKboost, and ANYboost, yield quantitatively equivalent solutions for the infective population. The latter solution curves are bounded from below by the solution of MAXboost and from above by the solution of the NOboost problem.Figure 4
Comparison of different boosting mechanisms in case of a stable endemic equilibrium. (a) The infective population in case of NOboost (black curve), MAXboost (green), ANYboost (red dotted), ODboost (blue), and HKboost (magenta). (b) The stationary distribution in the immune populationr¯z corresponding to (a).
(a)
(b)In Figure4(b), for each boosting mechanism, we visualize the stationary distribution r¯z corresponding to the endemic equilibrium in Figure 4(a). This indicates how immunity is distributed among the R-population at the endemic equilibrium. It is visible that the choice of the boosting mechanism importantly affects the stationary distribution. To understand why this happens, let us consider the case of a constant immune decay rate, gz=g^, for all z∈zmin,zmax. For certain choices of the probability density function, p one can calculate the stationary distribution r¯z explicitly. For example, in the absence of immune boosts (NOboost), we have
(17)r¯z=γI∗g^e−d/g^zmax−z,where I∗ is an endemic equilibrium of the systems (3) and (10), with boundary condition (5).In case of boosts to maximal immune level (MAXboost), from (13) and (14), we have
(18)r¯z=I∗g^γ+β1−1R0−I∗e−d+βI∗/g^zmax−z,with the basic reproduction number R0=β/γ+d, as in (6). With the parameter values indicated in Section 3.2, and R0=1.5, the curves in (18) and (17) are almost straight lines, as it can be observed in Figure 4(b).We use the analytic stationary solution (18) of systems (3), (4), (13), and (14) to verify the order of “spatial” accuracy (i.e., with respect to the variable z) of the numerical discretization. On different meshes, we compute numerically the stationary solution of the MAXboost systems (3), (4), (13), and (14). Although we know the explicit solution for the stationary distribution r¯, there is no explicit analytical formulation for the endemic equilibrium I∗, which can only be determined numerically (see also [7]). We fix I∗ to the value computed on the finest mesh (M=14000 equidistant mesh points), and we insert this value into the (18) of the stationary distribution. Then, we compute the maximum error between the analytical and the numerical solutions on different mesh refinements and observe that our method preserves the second-order accuracy of the Lax-Wendroff scheme (see Figure 5).Figure 5
Spatial accuracy of the numerical discretization, with error represented in a log-log plot. The method is second-order accurate.


### 3.4. Sustained Oscillations, Repeated Outbreaks

WhenR0=5 and all other parameter values are as indicated in Section 3.2, the systems (3), (4), and (5) show periodic oscillatory behavior, independent on the boosting mechanism. In Figure 6, we compare the I component of the solution of the systems (3), (4), and (5) for five different immune responses. Further, for all kinds of immune response, we compare in Figure 7 the distribution of immune status among the population for fixed times (t=1,5,10,20 years) after the beginning of the observation. Results illustrated in Figures 6 and 7 and Table 1 indicate that also in the case of periodic oscillations, the solutions of problems governed by ODboost and HKboost are qualitatively and quantitatively equivalent (minor differences might be due to computational errors and we consider them negligible). Therefore, in the following, we show only results related to NOboost, ANYboost, MAXboost, and ODboost, which characterize four different immune responses with qualitatively different results. Figure 8 shows the solution rt,z with respect to time and immune status for the four selected boosting mechanisms.Figure 6
Infective populationIt for different boosting mechanisms. Also in case of periodic oscillations, it is evident that the solution of problems governed by the boosting mechanisms ODboost and HKboost are quantitatively equivalent.Figure 7
Comparison of the effects of different boosting mechanisms on the distribution of immunityrt∗,z at fixed times t∗=1,5,10,20 years after the beginning of the observations.
(a)
t = 1
(b)
t = 5
(c)
t = 10
(d)
t = 20Table 1
Quantitative comparison of the periodic oscillatory solutions (I component) for changing immune boost mechanism. Oscillation amplitudes and periods, as well as lowest and highest incidence values, are reported.
Amplitude
Period
Min (I)
Max (I)
NOboost
0.3088
2.766
0.0169
0.3257
MAXboost
0.3588
3.4
0.0024
0.3612
ANYboost
0.3430
3.2
0.0036
0.3466
ODboost
0.4105
3.35
0.0023
0.4129
HKboost
0.3926
3.35
0.0024
0.3950Figure 8
Thert,z solution for different boosting mechanisms. (a) Waning immunity, absence of immune boosts. (b) Boost to maximal level of immunity. (c) Uniform probability of boosting to any higher immunity level. (d) Nonlinear boost according to Figure 3(c).
(a)
NOboost
(b)
MAXboost
(c)
ANYboost
(d)
ODboostIn Figure9, we visualize the contour plots of the r component of the solution for different immune response mechanisms. Entries of the solution matrix r are interpreted as heights with respect to the z,t plane; isolines are calculated and displayed using colors corresponding to the colormap on the right.Figure 9
Contour plot of ther solution for different boosting mechanisms. (a) Waning immunity, no boosts. (b) Boost to maximal level of immunity. (c) Uniform probability of boosting to any higher immunity level. (d) Nonlinear boosting function as in Figure 3(c).
(a)
NOboost
(b)
MAXboost
(c)
ANYboost
(d)
ODboostIn Figure10, we compare the immune distribution rt,z corresponding to different levels of an outbreak, for different boosting mechanisms. Figure 10(a) shows four points on a typical infective curve (I) which we shall consider for comparison. Point A corresponds to the infection peak, that is, the time point at which the infective population is at its maximal value. Point B and point D correspond to intermediate levels in the infective population, just after, respectively, just before an outbreak peak. Point C corresponds to the time point at which the infective population is at its minimal level. We see in Figures 10(b)–10(e) that the different immune responses yield qualitatively similar immune distributions with immunity waves moving from right to left. At the infective peak (Figure 10(b)), following primary infection and reexposure to the pathogen, most of the population has a very high immune status (z∈0.8,1) which slowly decays, reaching the level at which most of the population has a low to intermediate immune level (z∈0,0.6) (Figure 10(d)). Just before a new outbreak, most of the hosts have very low (z∈0,0.2) or very high (z∈0.8,1) immunity (Figure 10(e)).Figure 10
Comparison of immune distribution for different boosting mechanisms at different levels of an outbreak. (a) TheI solution and the four reference points to be considered for comparison: (A) the infection peak, (B) an intermediate level after the outbreak, (C) the minimal infection level, and (D) just before a new outbreak. (b–e) Immune distribution corresponding to the points (A)–(D).
(a)
(b)
Point (A)
(c)
Point (B)
(d)
Point (C)
(e)
Point (D)


## 3.1. Scenarios for the Immune Response

Let us consider an immune host who has immune statusv at the moment of reexposure to the pathogen. We shall investigate the following possible scenarios for immune response in case of secondary exposure to the pathogen.


### 3.1.1. NOboost

Assume that the immune system does not respond to reinfection, that is, we observe only waning of immunity after recovery (Figure2(a)). This corresponds to the limit case in which the boosting probability function is simply a Dirac measure with support on the initial immune level. It follows that (4) reduces to
(10)∂∂trt,z−∂∂zgzrt,z=−drt,z.Recall the definition ofτ>0 in (2). The transport equation (10) with boundary condition (5) is solved along characteristics and we obtain
(11)gzminrt,zmin=γIt−τe−dτ,meaning that τ time after recovery immune hosts who did not die become susceptible again. In turn, we find a delay term in the equation for S and have a classical SIRS model with constant delay (cf. [12])
(12)S′t=d−βStIt−dSt+γIt−τe−dτ,I′t=βStIt−γ+dIt,R′t=γIt−γIt−τe−dτ−dRt,where Rt is the total immune population at time t as in (1).


### 3.1.2. MAXboost

Assume that at any new encounter with the pathogen the immune system of a recovered host is boosted in such a way that the disease-induced (maximal) immunity is restored (Figure2(b)). This corresponds to the limit case in which the boosting probability function is a Dirac measure with support on the maximal immune level (zmax). As all boosted individuals are transferred to the boundary, and the integral term in (4) moves to the boundary condition, yielding
(13)∂∂trt,z−∂∂zgzrt,z=−drt,z−βItrt,z,(14)gzmaxrt,zmax=γIt+βItRt.The BVP (13) and (14) can be solved along characteristics and one obtains for t≥0 a system of delay equations with constant and distributed delay:
(15)S′t=d1−St−βItSt+It−τγ+βRt−τe−dτ−β∫−τ0It+udu,I′t=βItSt−γ+dIt,R′t=−dRt+γIt−It−τγ+βRt−τe−dτ−β∫−τ0It+udu,with τ>0 as defined in (2) and with given initial functions ϕSt≥0, ϕIt≥0, and ϕRt≥0, such that ϕSt+ϕIt+ϕRt≡1, for all t∈−τ,0. The formal derivation of system (15) is given in [7, 8].


### 3.1.3. ANYboost

We assume that the immune system is boosted to any better immunity level with uniform probability. That is,pz,v=p¯>0 with p¯=1/zmax−v, for all z∈v,zmax. The boosting matrix for the numerical implementation is shown in Figure 3(a).Figure 3
Nonmonotonic immune response mechanisms. (a, c, d) Examples of different boosting matricesΠ for the numerical implementation of pz,v, the probability that individuals with immune level v at exposure are boosted to level z≥v after exposure. The immune status is represented in zmin,zmax=0,1 and discretized on an equidistant mesh with 30 points. (b) The dotted curve represents the boosting mechanism inspired by [6], the solid curve the one inspired by [14].
(a)
ANYboost
(b)
Nonmonotonic boosting functions
(c)
ODboost
(d)
HKboost


### 3.1.4. HKboost and ODboost

Previous works have proposed nonmonotonic functions for modeling the immune response to secondary exposure. In [6, 13], a mathematical model for in-host dynamics during measles infection was proposed. The model explicitly considers the relation between the immune level at time of exposure and the level of memory cells after exposure. Heffernan and Keeling [6, 13] suggest that this relation is nonmonotonic. Indeed, although the level of immunity after exposure is always greater than the level at the time of exposure, the boosted level starts high, decreases, and then increases again. We shall denote by HKboost the boosting function from [6, 13]. For the numerical implementation of HKboost, we first extract values from Figure 2(d) in [6], which shows the relation between the initial and the updated immune status. Then, we normalize the immune status interval (zmin,zmax=0,1) and extend the boosting function to [0.75, 1], assuming that for high initial immune status, the immune system is always boosted to the maximal level (see Figure 3(b)). In terms of our model coefficients, for any initial immunity v, the atomic measure of the updated immunity is concentrated on fv, where f is the boosting function represented by the dotted curve in Figure 3(b). The corresponding boosting matrix is shown in Figure 3(d).Immune boosting in pertussis was considered in [14], where a nonmonotonic boosting function similar to the one in [6, 13] was proposed. The relation between the level of memory cells before (v) and after (z) reexposure is governed by the function
(16)z=fv=v1+15v5,withv∈0,5.De Graaf et al. [14] suggest that a small jump from v to z corresponds to a mild infection that causes no harm but only boosts the antibody level of the individual, whereas a large jump corresponds to a severe infection that may also cause disease. We shall use the boosting function in (16) for simulations, previous normalization to zmin,zmax=0,1, and small modifications which allow to have boundedness. In detail, we assume that for very low initial immunity (0≤v<0.075) and for large initial immunity (0.8<v≤1), the immune status is always boosted to the maximal level z=zmax=1, whereas for v∈0.075,0.8, the boosting probability is atomic along the graph of the boosting function f in (16), as above. The resulting boosting function is given by the solid curve in Figure 3(b), and the corresponding boosting matrix is shown in Figure 3(c).The two mechanisms for immune response HKboost and ODboost are governed by rather similar boosting mechanisms, and we shall see that the solutions behave accordingly.


## 3.1.1. NOboost

Assume that the immune system does not respond to reinfection, that is, we observe only waning of immunity after recovery (Figure2(a)). This corresponds to the limit case in which the boosting probability function is simply a Dirac measure with support on the initial immune level. It follows that (4) reduces to
(10)∂∂trt,z−∂∂zgzrt,z=−drt,z.Recall the definition ofτ>0 in (2). The transport equation (10) with boundary condition (5) is solved along characteristics and we obtain
(11)gzminrt,zmin=γIt−τe−dτ,meaning that τ time after recovery immune hosts who did not die become susceptible again. In turn, we find a delay term in the equation for S and have a classical SIRS model with constant delay (cf. [12])
(12)S′t=d−βStIt−dSt+γIt−τe−dτ,I′t=βStIt−γ+dIt,R′t=γIt−γIt−τe−dτ−dRt,where Rt is the total immune population at time t as in (1).


## 3.1.2. MAXboost

Assume that at any new encounter with the pathogen the immune system of a recovered host is boosted in such a way that the disease-induced (maximal) immunity is restored (Figure2(b)). This corresponds to the limit case in which the boosting probability function is a Dirac measure with support on the maximal immune level (zmax). As all boosted individuals are transferred to the boundary, and the integral term in (4) moves to the boundary condition, yielding
(13)∂∂trt,z−∂∂zgzrt,z=−drt,z−βItrt,z,(14)gzmaxrt,zmax=γIt+βItRt.The BVP (13) and (14) can be solved along characteristics and one obtains for t≥0 a system of delay equations with constant and distributed delay:
(15)S′t=d1−St−βItSt+It−τγ+βRt−τe−dτ−β∫−τ0It+udu,I′t=βItSt−γ+dIt,R′t=−dRt+γIt−It−τγ+βRt−τe−dτ−β∫−τ0It+udu,with τ>0 as defined in (2) and with given initial functions ϕSt≥0, ϕIt≥0, and ϕRt≥0, such that ϕSt+ϕIt+ϕRt≡1, for all t∈−τ,0. The formal derivation of system (15) is given in [7, 8].


## 3.1.3. ANYboost

We assume that the immune system is boosted to any better immunity level with uniform probability. That is,pz,v=p¯>0 with p¯=1/zmax−v, for all z∈v,zmax. The boosting matrix for the numerical implementation is shown in Figure 3(a).Figure 3
Nonmonotonic immune response mechanisms. (a, c, d) Examples of different boosting matricesΠ for the numerical implementation of pz,v, the probability that individuals with immune level v at exposure are boosted to level z≥v after exposure. The immune status is represented in zmin,zmax=0,1 and discretized on an equidistant mesh with 30 points. (b) The dotted curve represents the boosting mechanism inspired by [6], the solid curve the one inspired by [14].
(a)
ANYboost
(b)
Nonmonotonic boosting functions
(c)
ODboost
(d)
HKboost


## 3.1.4. HKboost and ODboost

Previous works have proposed nonmonotonic functions for modeling the immune response to secondary exposure. In [6, 13], a mathematical model for in-host dynamics during measles infection was proposed. The model explicitly considers the relation between the immune level at time of exposure and the level of memory cells after exposure. Heffernan and Keeling [6, 13] suggest that this relation is nonmonotonic. Indeed, although the level of immunity after exposure is always greater than the level at the time of exposure, the boosted level starts high, decreases, and then increases again. We shall denote by HKboost the boosting function from [6, 13]. For the numerical implementation of HKboost, we first extract values from Figure 2(d) in [6], which shows the relation between the initial and the updated immune status. Then, we normalize the immune status interval (zmin,zmax=0,1) and extend the boosting function to [0.75, 1], assuming that for high initial immune status, the immune system is always boosted to the maximal level (see Figure 3(b)). In terms of our model coefficients, for any initial immunity v, the atomic measure of the updated immunity is concentrated on fv, where f is the boosting function represented by the dotted curve in Figure 3(b). The corresponding boosting matrix is shown in Figure 3(d).Immune boosting in pertussis was considered in [14], where a nonmonotonic boosting function similar to the one in [6, 13] was proposed. The relation between the level of memory cells before (v) and after (z) reexposure is governed by the function
(16)z=fv=v1+15v5,withv∈0,5.De Graaf et al. [14] suggest that a small jump from v to z corresponds to a mild infection that causes no harm but only boosts the antibody level of the individual, whereas a large jump corresponds to a severe infection that may also cause disease. We shall use the boosting function in (16) for simulations, previous normalization to zmin,zmax=0,1, and small modifications which allow to have boundedness. In detail, we assume that for very low initial immunity (0≤v<0.075) and for large initial immunity (0.8<v≤1), the immune status is always boosted to the maximal level z=zmax=1, whereas for v∈0.075,0.8, the boosting probability is atomic along the graph of the boosting function f in (16), as above. The resulting boosting function is given by the solid curve in Figure 3(b), and the corresponding boosting matrix is shown in Figure 3(c).The two mechanisms for immune response HKboost and ODboost are governed by rather similar boosting mechanisms, and we shall see that the solutions behave accordingly.


## 3.2. Parameter Values

For the numerical computations below, we set the birth rate, natural death rated=0.02, and the initial conditions ψ0z=0.05 for z∈zmin,zmax=0,1, I0=0.01, and S0=0.94. This means that we assume that 1% of the total population is infectious at the beginning of our observations, while 5% is immune to the pathogen and the level of immunity is equally distributed.Concerning the disease dynamics, we setγ=3 and gz=0.5, for all z∈zmin,zmax, corresponding to τ=2 when the model can be reduced to DDE.We allow the basic reproduction numberR0 as given in (6) to vary, in order to show both solutions that converge to an endemic equilibrium and solutions that produce periodic oscillations.


## 3.3. Stable Endemic Equilibrium

When the basic reproduction number is sufficiently small, the solution converges to an endemic equilibrium for all the boosting mechanisms presented in Section3.1. We compare in Figure 4 the numerical solution corresponding to different boosting mechanisms. In Figure 4(a), we show the component I of the solution, which indicates how the number of infective evolves in time and rapidly converges to an equilibrium. Changing the boosting mechanism has small effects on the infective population: numerical solutions show the same qualitative behavior, with endemic equilibria deviating less than 1% from each other. In particular, we see that the nonmonotonic boosting mechanisms, ODboost, HKboost, and ANYboost, yield quantitatively equivalent solutions for the infective population. The latter solution curves are bounded from below by the solution of MAXboost and from above by the solution of the NOboost problem.Figure 4
Comparison of different boosting mechanisms in case of a stable endemic equilibrium. (a) The infective population in case of NOboost (black curve), MAXboost (green), ANYboost (red dotted), ODboost (blue), and HKboost (magenta). (b) The stationary distribution in the immune populationr¯z corresponding to (a).
(a)
(b)In Figure4(b), for each boosting mechanism, we visualize the stationary distribution r¯z corresponding to the endemic equilibrium in Figure 4(a). This indicates how immunity is distributed among the R-population at the endemic equilibrium. It is visible that the choice of the boosting mechanism importantly affects the stationary distribution. To understand why this happens, let us consider the case of a constant immune decay rate, gz=g^, for all z∈zmin,zmax. For certain choices of the probability density function, p one can calculate the stationary distribution r¯z explicitly. For example, in the absence of immune boosts (NOboost), we have
(17)r¯z=γI∗g^e−d/g^zmax−z,where I∗ is an endemic equilibrium of the systems (3) and (10), with boundary condition (5).In case of boosts to maximal immune level (MAXboost), from (13) and (14), we have
(18)r¯z=I∗g^γ+β1−1R0−I∗e−d+βI∗/g^zmax−z,with the basic reproduction number R0=β/γ+d, as in (6). With the parameter values indicated in Section 3.2, and R0=1.5, the curves in (18) and (17) are almost straight lines, as it can be observed in Figure 4(b).We use the analytic stationary solution (18) of systems (3), (4), (13), and (14) to verify the order of “spatial” accuracy (i.e., with respect to the variable z) of the numerical discretization. On different meshes, we compute numerically the stationary solution of the MAXboost systems (3), (4), (13), and (14). Although we know the explicit solution for the stationary distribution r¯, there is no explicit analytical formulation for the endemic equilibrium I∗, which can only be determined numerically (see also [7]). We fix I∗ to the value computed on the finest mesh (M=14000 equidistant mesh points), and we insert this value into the (18) of the stationary distribution. Then, we compute the maximum error between the analytical and the numerical solutions on different mesh refinements and observe that our method preserves the second-order accuracy of the Lax-Wendroff scheme (see Figure 5).Figure 5
Spatial accuracy of the numerical discretization, with error represented in a log-log plot. The method is second-order accurate.


## 3.4. Sustained Oscillations, Repeated Outbreaks

WhenR0=5 and all other parameter values are as indicated in Section 3.2, the systems (3), (4), and (5) show periodic oscillatory behavior, independent on the boosting mechanism. In Figure 6, we compare the I component of the solution of the systems (3), (4), and (5) for five different immune responses. Further, for all kinds of immune response, we compare in Figure 7 the distribution of immune status among the population for fixed times (t=1,5,10,20 years) after the beginning of the observation. Results illustrated in Figures 6 and 7 and Table 1 indicate that also in the case of periodic oscillations, the solutions of problems governed by ODboost and HKboost are qualitatively and quantitatively equivalent (minor differences might be due to computational errors and we consider them negligible). Therefore, in the following, we show only results related to NOboost, ANYboost, MAXboost, and ODboost, which characterize four different immune responses with qualitatively different results. Figure 8 shows the solution rt,z with respect to time and immune status for the four selected boosting mechanisms.Figure 6
Infective populationIt for different boosting mechanisms. Also in case of periodic oscillations, it is evident that the solution of problems governed by the boosting mechanisms ODboost and HKboost are quantitatively equivalent.Figure 7
Comparison of the effects of different boosting mechanisms on the distribution of immunityrt∗,z at fixed times t∗=1,5,10,20 years after the beginning of the observations.
(a)
t = 1
(b)
t = 5
(c)
t = 10
(d)
t = 20Table 1
Quantitative comparison of the periodic oscillatory solutions (I component) for changing immune boost mechanism. Oscillation amplitudes and periods, as well as lowest and highest incidence values, are reported.
Amplitude
Period
Min (I)
Max (I)
NOboost
0.3088
2.766
0.0169
0.3257
MAXboost
0.3588
3.4
0.0024
0.3612
ANYboost
0.3430
3.2
0.0036
0.3466
ODboost
0.4105
3.35
0.0023
0.4129
HKboost
0.3926
3.35
0.0024
0.3950Figure 8
Thert,z solution for different boosting mechanisms. (a) Waning immunity, absence of immune boosts. (b) Boost to maximal level of immunity. (c) Uniform probability of boosting to any higher immunity level. (d) Nonlinear boost according to Figure 3(c).
(a)
NOboost
(b)
MAXboost
(c)
ANYboost
(d)
ODboostIn Figure9, we visualize the contour plots of the r component of the solution for different immune response mechanisms. Entries of the solution matrix r are interpreted as heights with respect to the z,t plane; isolines are calculated and displayed using colors corresponding to the colormap on the right.Figure 9
Contour plot of ther solution for different boosting mechanisms. (a) Waning immunity, no boosts. (b) Boost to maximal level of immunity. (c) Uniform probability of boosting to any higher immunity level. (d) Nonlinear boosting function as in Figure 3(c).
(a)
NOboost
(b)
MAXboost
(c)
ANYboost
(d)
ODboostIn Figure10, we compare the immune distribution rt,z corresponding to different levels of an outbreak, for different boosting mechanisms. Figure 10(a) shows four points on a typical infective curve (I) which we shall consider for comparison. Point A corresponds to the infection peak, that is, the time point at which the infective population is at its maximal value. Point B and point D correspond to intermediate levels in the infective population, just after, respectively, just before an outbreak peak. Point C corresponds to the time point at which the infective population is at its minimal level. We see in Figures 10(b)–10(e) that the different immune responses yield qualitatively similar immune distributions with immunity waves moving from right to left. At the infective peak (Figure 10(b)), following primary infection and reexposure to the pathogen, most of the population has a very high immune status (z∈0.8,1) which slowly decays, reaching the level at which most of the population has a low to intermediate immune level (z∈0,0.6) (Figure 10(d)). Just before a new outbreak, most of the hosts have very low (z∈0,0.2) or very high (z∈0.8,1) immunity (Figure 10(e)).Figure 10
Comparison of immune distribution for different boosting mechanisms at different levels of an outbreak. (a) TheI solution and the four reference points to be considered for comparison: (A) the infection peak, (B) an intermediate level after the outbreak, (C) the minimal infection level, and (D) just before a new outbreak. (b–e) Immune distribution corresponding to the points (A)–(D).
(a)
(b)
Point (A)
(c)
Point (B)
(d)
Point (C)
(e)
Point (D)


## 4. Conclusion

Understanding the role of immune system boosting due to the interplay of in-host and between hosts dynamics is a central point in immunoepidemiology of infectious diseases. Employing the mathematical models (3), (4), and (5), in this work, we have systematically investigated the effects of different (in-host) immune responses at population level. The boosting mechanisms studied here are in part taken or inspired by previous literature, for example, on measles [6, 13], pertussis [14], and influenza [15].Our results indicate that the in-host immune response, that is, the individual boosting mechanism, has important effects on the qualitative dynamics of the epidemics. The numerical results show that observing the distribution of immune level among recovered individuals, one can reconstruct the underlying boosting mechanism. When the trajectories of the systems (3), (4), and (5) approach an endemic equilibrium (Figure 4(b)), the distribution of immunity uniquely identifies the boosting mechanism governing the immune response in case of secondary infections. Also for diseases with repeated outbreaks, it is possible to infer the immune boosting mechanism from the temporal evolution of immunity, though in this case, it is convenient to look at the dynamic of the immune population over time (Figures 8 and 9), rather than at the immune distribution (Figure 10).This paper presents a novel mathematical study of temporal evolution of immune status in a population. Despite the large number of publications in mathematical epidemiology or mathematical immunology, to the best of our knowledge, there have been only few authors who have effectively combined in-host dynamics and processes at population level [6, 13, 16, 17], though only in [6, 13] the temporal evolution of immunity has been considered. If we compare the stationary distribution of immunity of HKboost in Figure 4(b) with Figure 3(c) in [6], we observe that the results are rather different, the curve in [6] being smoother. We conjecture that the differences are on the one hand due to the choice of the parameters in the system, on the other hand to the different nature of the mathematical models. In [6], a large system of ordinary differential equation is presented, including exposed hosts and immune structure in all compartments, not only in the recovered ones as it is in our work. Here, we focused on a simple model for susceptible-infective-recovered dynamics in order to understand the role of immune boosting and waning immunity. We plan to extend the models (3), (4), and (5), by means of physiologically structured susceptible and exposed populations.Parameter values used for all numerical simulations presented in this work are not inspired by a specific infectious disease, as we purely focused on the effects of waning and boosting. Nevertheless, the model is suitable to describe several infectious diseases, prior proper choice of the parameters. Certain studies, for example [5, 18], introduced a further parameter κ to express a relative force of infection for reexposure, so they represent the rate of boosting by the term (κβI). This is reasonable if we think that during contacts, boosting occurs with a different chance compared to primary infection. In this paper, we have assumed that the force of infection in secondary exposure is the same as in primary exposure (κ=1).The in-host dynamics described in this paper is rather coarse. The immune status of recovered individuals is represented by a single scalar quantity (z) which does not specifically represent any kind of cells of the immune system. We are aware of the fact that the immune system is indeed very complex, and careful mathematical modeling should also take into consideration nonlinear in-host processes and interactions among the different players of the innate and adaptive immune system. Nevertheless, such a fine description would quickly lead to a complicated mathematical model, making it hard to achieve any reliable analytical or numerical result at population level.A further limitation of the model is the assumption on the sharp thresholdzmin, which defines the criterion for transition from immune to susceptible compartment. It is plausible that this transition does not occur in such an on-off manner, but it is rather a continuous process; that is, recovered individuals with a certain critically low level of immunity are less likely to get immune boost and more likely to experience a new infection, as suggested in [14]. Moreover, it is known that the immune system is differently reactive at different life stages, in particular immune boosts are weaker in children and elderly than in adults (see, e.g., the case of pertussis considered in [19]). We plan to extend the model including age heterogeneity.Immune boosts occur not only after natural infection but also in case of vaccine-induced immunity. Current vaccination schedules include “boosters” whose goal is the prolongation of protection against a certain pathogen [20]. The mathematical models (3), (4), and (5) can be extended to include vaccination and natural boosts induced in vaccinated hosts by contact with infectives (see [21]). Similar investigation as proposed in this paper can be repeated in the context of vaccination.Our results indicate that immune system boosts from natural secondary infections importantly affect the immune dynamics at population level. We hope that our findings will stimulate future research in controlling infectious diseases via vaccination and vaccine design.



---

*Source: 9264743-2018-08-30.xml*