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.
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 the framework, a population of hosts is divided into susceptibles (), infectives (), and recovered (), 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 transits back to the susceptible compartment; hence when waning of immunity is included, the model is called .
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 Section 2, 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
Let and denote the total population of susceptibles and infectives, respectively, at time . The total population shall be assumed to be in balance (hence normalized) , with birth rate equal to the natural death rate and no disease-induced death. We assume that newborns are all susceptible.
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 life-long protection. Immune hosts who experience immunity loss become susceptible again.
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. The value corresponds to the maximal immunity, whereas corresponds to the lowest level of immunity for hosts in the compartment. We assume that individuals who recover at time enter the immune compartment with maximal level of immunity . The level of immunity tends to decay in time and when it reaches the lower threshold , 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 Figure 1).
Given a host with initial immune status , let us denote by the updated immune status, which is achieved after new contact with the pathogen. The updated immune status is modeled as a random variable taking values in with probability density function . The value represents the relative likelihood that a boost from immune level to level 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 (). 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 .
The immunity level decays in time at rate , with positive, smooth, and bounded, which is the same for all immune individuals with immunity level . In the absence of immune system boosting, an infected host who recovered at time becomes again susceptible at time , where
With the above assumptions, we obtain for the following system of equations with initial values , , coupled with a partial differential equation (PDE) for the immune population, for , with boundary condition and initial distribution . 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].
(a) No boosts
(b) Boost to maximal level
(c) Boost to any higher level
2.2.
Before presenting further results, we introduce the basic reproduction number of models (3), (4), and (5), 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 as time and space, respectively. As and do not depend on , but only on time, the initial conditions are
Since for all the boundary condition for is imposed in (5) at the inflow boundary, that is, at . 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 mesh The initial function on this mesh is given by
For an explicit scheme for a hyperbolic system, a necessary condition for stability is the Courant-Friedrichs-Lewy condition, which for our system means
We define the boosting matrix as the numerical discretization of the probability density function on the equidistant mesh . Each entry of the boosting matrix represents the probability that is the updated immune level, given initial immune level , that is, . It follows that , with if and
3. Results and Discussion
3.1. Scenarios for the Immune Response
Let us consider an immune host who has immune status 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 (Figure 2(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
Recall the definition of in (2). The transport equation (10) with boundary condition (5) is solved along characteristics and we obtain 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 and have a classical SIRS model with constant delay (cf. [12]) where is the total immune population at time 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 (Figure 2(b)). This corresponds to the limit case in which the boosting probability function is a Dirac measure with support on the maximal immune level (). As all boosted individuals are transferred to the boundary, and the integral term in (4) moves to the boundary condition, yielding
The BVP (13) and (14) can be solved along characteristics and one obtains for a system of delay equations with constant and distributed delay: with as defined in (2) and with given initial functions , , and , such that , for all . 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, with , for all . The boosting matrix for the numerical implementation is shown in Figure 3(a).
(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 () 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 , the atomic measure of the updated immunity is concentrated on , where 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 () and after () reexposure is governed by the function
De Graaf et al. [14] suggest that a small jump from to 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 and small modifications which allow to have boundedness. In detail, we assume that for very low initial immunity () and for large initial immunity (), the immune status is always boosted to the maximal level , whereas for , the boosting probability is atomic along the graph of the boosting function 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 rate , and the initial conditions for , , and . This means that we assume that of the total population is infectious at the beginning of our observations, while is immune to the pathogen and the level of immunity is equally distributed.
Concerning the disease dynamics, we set and , for all , corresponding to when the model can be reduced to DDE.
We allow the basic reproduction number 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 Section 3.1. We compare in Figure 4 the numerical solution corresponding to different boosting mechanisms. In Figure 4(a), we show the component 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.
(a)
(b)
In Figure 4(b), for each boosting mechanism, we visualize the stationary distribution corresponding to the endemic equilibrium in Figure 4(a). This indicates how immunity is distributed among the -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, , for all . For certain choices of the probability density function, one can calculate the stationary distribution explicitly. For example, in the absence of immune boosts (NOboost), we have where 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 with the basic reproduction number , as in (6). With the parameter values indicated in Section 3.2, and 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 ) 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 , there is no explicit analytical formulation for the endemic equilibrium , which can only be determined numerically (see also [7]). We fix to the value computed on the finest mesh ( 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).
3.4. Sustained Oscillations, Repeated Outbreaks
When 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 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 ( 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 with respect to time and immune status for the four selected boosting mechanisms.
(a) t = 1
(b) t = 5
(c) t = 10
(d) t = 20
(a) NOboost
(b) MAXboost
(c) ANYboost
(d) ODboost
In Figure 9, we visualize the contour plots of the component of the solution for different immune response mechanisms. Entries of the solution matrix are interpreted as heights with respect to the plane; isolines are calculated and displayed using colors corresponding to the colormap on the right.
(a) NOboost
(b) MAXboost
(c) ANYboost
(d) ODboost
In Figure 10, we compare the immune distribution corresponding to different levels of an outbreak, for different boosting mechanisms. Figure 10(a) shows four points on a typical infective curve () 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 () which slowly decays, reaching the level at which most of the population has a low to intermediate immune level () (Figure 10(d)). Just before a new outbreak, most of the hosts have very low () or very high () immunity (Figure 10(e)).
(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 (). 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 ().
The in-host dynamics described in this paper is rather coarse. The immune status of recovered individuals is represented by a single scalar quantity () 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 threshold , 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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
M. V. Barbarossa is supported by the European Social Fund and by the Ministry of Science, Research and the Arts Baden-Württemberg. The work of M. Polner was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences, the EU-funded Hungarian Grants EFOP-3.6.2-16-2017-00015 and NKFI KH 125628. G. Röst was supported by Marie Sklodowska-Curie Grant nos. 748193 and NKFI FK 124016.