Next Article in Journal
To Solve Forward and Backward Nonlocal Wave Problems with Pascal Bases Automatically Satisfying the Specified Conditions
Next Article in Special Issue
Efficient Numerical Solutions to a SIR Epidemic Model
Previous Article in Journal
Online Support Vector Machine with a Single Pass for Streaming Data
Previous Article in Special Issue
Exploring Radial Kernel on the Novel Forced SEYNHRV-S Model to Capture the Second Wave of COVID-19 Spread and the Variable Transmission Rate
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Multi-Scale Model for Cholera Outbreaks

1
Center for Mathematics, Technische Universität München, 85748 Garching, Germany
2
Department of Mathematics, Physics and Computing, Moi University, Eldoret P.O. Box 3900, Kenya
3
Institute for Computational Biology, Helmholtz Center Munich, 85764 Neuherberg, Germany
4
Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(17), 3114; https://doi.org/10.3390/math10173114
Submission received: 28 July 2022 / Revised: 22 August 2022 / Accepted: 24 August 2022 / Published: 30 August 2022
(This article belongs to the Special Issue Mathematical Methods and Models in Epidemiology)

Abstract

:
Cholera, caused by the pathogenic Vibrio cholerae bacteria, remains a severe public health threat. Although a lot of emphasis has been placed on the population-level spread of the disease, the infection itself starts within the body. As such, we formulated a multi-scale model that explicitly connects the within-host and between-host dynamics of the disease. To model the within-host dynamics, we assigned each susceptible individual with a pathogen load that increases through the uptake of contaminated food and water (booster event). We introduced minimal and maximal times when the booster events happen and defined a time since the last booster event. We then scaled the within-host dynamics to the population where we structured the susceptible population using the two variables (pathogen load and time since the last booster event). We analyzed the pathogen load’s invariant distribution and utilized the results and time scale assumptions to reduce the dimension of the multi-scale model. The resulting model is an SIR model whose incidence function has terms derived from the multi-scale model. We finally conducted numerical simulations to investigate the long-term behavior of the SIR model. The simulations revealed parameter regions where either no cholera cases happen, where cholera is present at a low prevalence, and where a full-blown cholera epidemic takes off.

1. Introduction

Cholera, caused by the pathogenic Vibrio cholerae bacteria, is associated with the consumption of contaminated food and water. Three different cases are observed when the infection is introduced in a region: either no cholera cases are observed, even though the pathogen is present in the environment, or there are persistently few cases or a full epidemic occurs [1]. Most models have difficulties capturing these three outcomes. In the present work, we take up recent ideas and return to first principles to derive a refined model of the within and between-host dynamics. It turns out that this model qualitatively explains the observations in a natural way.
Since the seminal papers of Kermack and McKendrick [2,3], epidemic models have described the time course of infections primarily at the population level. Although that has proved to be an effective approximation, the transmission dynamics of an infection span multiple scales: within-host and population-level among others [4,5]. Unlike the population scale, where the interest is in the spread of the infection between individuals, the focal point of the within-host scale is the evolution of the infection inside a single individual [6]. It is therefore necessary to bridge the multiple scales for a comprehensive description of the entire infection process. An appropriate formulation of the incidence term can potentially be gained from such models; so far, the incidences used are mainly mass action [7], standard incidence [8] or various forms of nonlinear incidences [9].
Models that link the two scales of infection (multi-scale models) have been advanced in recent years [10,11,12,13]. An important aspect in the formulation of such models is the linking of the two scales. Many models keep the within-host and the between-host dynamics fairly separated, and extract information from the within-host model to feed into the population-level model [13]. Others [12,14] go a step further and structure the population according to the state of individuals, in line with what is proposed in [15,16]. This creates a better connection between the two levels. However, most papers mentioned above keep the strict distinction between susceptible and infected persons, which appears natural in the tradition of epidemic modeling; however, as we discuss next, on a second glance and at least in certain cases, it is somewhat arbitrary.
In the present work, we use cholera as a base to propose a further step in connecting the within-host dynamics to the disease outbreak at the population level; that is, the model bridges two different scales. The disease, as a severe public health threat, has been targeted by the modeling community for a long time. Population-level models for the infection have been advanced in [1,17,18,19], some of which structure the population according to the age [20,21]. A multi-scale model for the disease proposed in [22] makes use of the interaction between the environmental vibrios and human vibrios (inside the body) to link the within- and between-host dynamics. The ideas are extended in [23] where the within-host model includes the interaction with the immune system and the connection of the two scales is through the pathogen evolution in the environment. In [24], the epidemic model is structured using a function defined as the immune status, which is derived from the within-host model.
Similarly, our model seeks to advance the multi-scale modeling of the disease. On one hand, we aimed to introduce a novel model structure that is closer to the immunological and environmental processes. On the other hand, we aimed to provide a rigorous way to analyze this novel model structure. Unlike the previously mentioned models, the epidemic model is structured using the pathogen load such that the distinction between a susceptible and infected individual is relaxed. We formulate the model by explicitly addressing the pathogen level in susceptible individuals. This sounds contradictory at first glance since susceptible individuals are usually thought to be pathogen free. However, the bacteria causing cholera, Vibrio cholerae, exist in the environment [25]. Contaminated food and water, enhanced by improper sanitation and hygiene practices, is the main cause of infection. In some instances, the consumption of contaminated food and water does not immediately trigger infection. This is because a cholera infection only occurs when the pathogen load of an individual exceeds a critical threshold [26,27] (which varies from individual to individual, e.g., caused by differences in the gut microbiome [28,29]). That is, individuals take up V. cholerae, but the innate immune system fights off and eliminates the bacteria as long as the critical pathogen threshold is not exceeded. Furthermore, laboratory experiments on mice show that a moderate increase in the infectious dose leads to an increase in the pathogen burden over time (with a time scale of approximately 12 h) [30], and not to infection. We are thus led to a model, where the susceptible are structured by their pathogen load, and the transition into the infected class happens at a rate dependent on the pathogen load. This step allows for better connection of the within- and between-host dynamics, as the somewhat arbitrary distinction between a susceptible person and an infected person is softened.
Mathematically, the model structure resembles that of fragmentation–aggregation equations: we have a time-continuous process of pathogen clearance (we assume that pathogens are degraded at a constant rate within a host), and at randomly distributed times, by food uptake, we have booster events that instantaneously increase the pathogen load. Accordingly, the analysis of the invariant distribution for the pathogen load follows the theory developed for aggregation–fragmentation equations, particularly for cell division, given in the 1980s by Heijmans, Gyllenberg and others [31,32], which was later extended, e.g., by Doumic and others [33,34,35,36]. The understanding of the pathogen load’s invariant distribution and the assumption of a time scale separation then allows the reduction of the rather complex multi-scale model to a population-level only model, where the incidence term assumes a form derived from the multi-scale model. We then analyze the resultant model numerically.
The paper is structured as follows: In Section 2, we introduce the model. In Section 3, we focus on the within-host pathogen dynamics and show that these dynamics yield, in the long run, an invariant pathogen distribution. In Section 3, we perform a rigorous mathematical analysis of the model. The results are intuitive while the content is rather technical; a reader who is interested in the application side rather than the mathematical theory should consider skipping this section. In Section 4, we derive a kind of SIR-model from our structured model and discuss the long term behavior. The results are put into context in the discussion in Section 5.

2. Model

The human population is subdivided into three compartments that represent the susceptible (S), infected (I) and recovered (R) individuals. Additionally, we have the bacterial concentration in the environment B. The pathogens (Vibrio cholerae) can survive and reproduce in the environment without interaction with human hosts. We simply assume an environmental pathogen production rate a and a pathogen death rate σ . As a detailed model of the pathogen’s life cycle is out of the scope of the present article, we aim to focus on the human population.
In particular, we aimed at a model connecting the within-host level with the population level. Thereto, we assign a pathogen load P to each susceptible person. Although transmission can take place directly through person-to-person contact, we only considered the case of indirect transmission since direct transmission is a rare occurrence [37]. That is, contact with contaminated food boosters the pathogen load of an individual by ψ B . As food uptake does not take place randomly, there is a minimal and a maximal time between these events, we also take the time since the last booster event τ into account. That is, we have S = S ( t , τ , P ) , the susceptible class is structured according to their pathogen load, and the time since the last uptake event. It turns out, that this structure is not only more realistic but also mathematically convenient. The timing of booster events is modeled by ρ ( τ ) , which we will discuss below in more detail.
It is well known, that a critical pathogen load is required for an infection to occur [26,38]. For a sub-critical level, the innate immune system is able to control the pathogens. We take that into account by defining a clearance rate γ for the within-host pathogens. That is, within a single susceptible individual, we have an interplay between the booster event and the degradation of pathogens, which establishes a stochastic process (see Figure 1). Experiments on mice with a moderate pathogen level show the exact kind of decline in the pathogen load addressed by this sub-model for the susceptible [30]. Furthermore, the infection rate β is a function depending on the pathogen load, β = β ( P ) . Here, β ( P ) = 0 if the pathogen load is sub-critical. β only becomes positive if the pathogen load P exceeds some value and persons with such a high pathogen load are transferred into the infected class. An infected individual sheds pathogens at rate ξ into the environment, and recovers at a rate of α . We assume that the recovered individuals become immune and stay immune: the time scale of our model covers months (time scale of one cholera epidemic), and not years, where recovered individuals become susceptible again.
We add one more assumption: we expect the dynamics on the population level to be slow, while the within-individual pathogen dynamics are fast. A typical cholera epidemic has a time scale of weeks or months, while the time between two booster events is rather hours, and accordingly, the time scale of the native immune event to handle the ingested pathogens is the same as the time scale for the uptake events (otherwise the pathogens accumulate and eventually an infection would be inevitable). The slow time scale is expressed by the rate ϵ .
All in all, the model reads
t S ( t , τ , P ) + τ S ( t , τ , P ) + P ( γ P S ( t , τ , P ) ) = ρ ( τ ) S ( t , τ , P ) ϵ β ( P ) S ( t , τ , P ) S ( t , 0 , P ) = τ ˇ τ ^ ρ ( τ ) S ( t , τ , P ψ B ) d τ S ( t , τ , P ) = 0 for P < 0 I ˙ = 0 τ ˇ τ ^ ϵ β ( P ) S ( t , τ , P ) d τ d P ϵ α I R ˙ = ϵ α I B ˙ = ϵ ( a + ξ I σ B ) S ( 0 , τ , P ) = S 0 ( τ , P ) , I ( 0 ) = I 0 , R ( 0 ) = R 0 , B ( 0 ) = B 0 .
The boundary condition S ( t , 0 , P ) defines the pathogen density in the susceptible right after a booster event. The time since the last booster event τ is zero in this case, while the delay term P ψ B implies that the boostered individuals had a lower pathogen load before the jump. After defining the model equations, we turn to discuss a central aspect, the timing of booster events coded by ρ ( τ ) . The idea is that, beneath our deterministic model, there is a stochastic process describing the uptake of pathogens. We will later use this idea to formulate a numerical method to obtain a stationary solution of our model. The time between two booster events are i.i.d. as a random variable T with probability density φ ( t ) C 1 ( R + ) . That is, P ( T < τ ) = 0 τ φ ( τ ) d t . We assume that there is a minimal time τ ˇ > 0 between two booster events, P ( T < τ ˇ ) = 0 , as well as a maximal time τ ^ > τ ˇ , P ( T < τ ^ ) = 1 . Therefore,
φ ( τ ) = 0 f o r τ [ 0 , τ ˇ ) ( τ ^ , ) , τ ˇ τ ^ φ ( τ ) d τ = 1 .
These assumptions are sensible given the application, and, as we will find out, are also convenient in the analysis of the model. The parameter ρ ( τ ) in the model is now the hazard rate defined by φ ( τ ) ,
ρ ( τ ) = φ ( τ ) 1 0 τ φ ( s ) d s .
We superimpose several rather technical assumptions on the hazard rate.
Assumption 1.
The hazard rate ρ ( τ ) = φ ( τ ) 1 0 τ φ ( s ) d s stems from a distribution φ ( τ ) C 1 ( R + ) , where φ satisfies (2). We also assume that
( a ) lim τ τ ^ 0 τ ρ ( τ ) d τ = , ( b ) sup τ [ 0 , τ ^ ] ρ ( τ ) e 0 τ ρ ( s ) d s < , lim τ τ ^ ρ ( τ ) e 0 τ ρ ( s ) d s = 0 .
( c ) τ ˇ τ ^ ρ 2 ( τ ) e 0 τ ρ ( s ) d s d τ < , ( d ) τ ˇ τ ^ | ρ ( τ ) | e 0 τ ρ ( s ) d s d τ < .
In Appendix B.1, we show that these assumptions are satisfied by a wide range of random variables T (resp. their distribution φ ).

3. Pathogen Distribution for Constant Environmental Pathogen Load

It follows that the analysis of the pathogen distribution S ( t , τ , P ) under the condition that the environmental pathogen load B is constant and no infections are occurring, we take ϵ to zero in (1), is central in the reduction in the model. That is, we consider
t S ( t , τ , P ) + τ S ( t , τ , P ) + P ( γ P S ( t , τ , P ) ) = ρ ( τ ) S ( t , τ , P ) S ( t , 0 , P ) = g ( t , P ) = τ ˇ τ ^ ρ ( τ ) S ( t , τ , P ψ B ) d τ S ( 0 , τ , P ) = S 0 ( τ , P ) .
for B > 0 given, fixed. The aim of this section is to show that there is a stationary solution of (3) (the invariant distribution of the underlying stochastic model), and that any non-negative initial condition eventually tends to this solution. Thereto, we first show that (3) defines a strongly continuous and eventually compact semigroup, and then inspect the spectrum of the infinitesimal generator of this semigroup. As we will find out, the generator has an eigenvalue 0 (caused by mass conservation), which is dominating. The non-local term in the boundary condition (jump of B by a booster event) is not completely straightforward to handle.
Mathematically, this semigroup is close to systems describing cell division, or more generally, aggregation–fragmentation equations. While in aggregation–fragmentation equations, entities are growing and are decreased by sudden non-local fragmentation events, in our case, the continuous degradation of the immune system decreases P and the sudden disruptive non-local events (booster events) increase P. This is the difference between aggregation–fragmentation equations and our model. However, for the analysis, we mostly follow the strategy used for those equations developed by Heijmans, Gyllenberg and others [31,32,39], and advanced in recent years [33,34].
We adapt the methods above to address a technical problem that arises in the analysis: the lack of strong positivity. That is, we find a compact interval [ P * , P * ] , with P * and P * being the lower and upper bound of the pathogen, respectively, for the pathogen load right after a booster event happens. However, the probability mass is zero at the boundary of this interval. To solve this, we need to first regularize the problem before using the standard arguments. Additional effort is necessary if we pass to the limit and de-regularize the operators again, to ensure that the desired results remain valid in the limit.
Most part of the present section focuses on the mathematical theory and is rather technical. While it is rather intuitive that the results we aim at are true, a reader who is mainly interested in the epidemic modeling aspect might consider skipping this section and jumping directly to Section 4, where we use the results derived in the present section to obtain a reduced model that will be analyzed under the application point of view. For an introduction to concepts of function spaces, compactness, semigroups and spectral theory, we refer the reader to [40,41].

3.1. Semigroup

Below, using the methods of characteristics, we show that there is a bounded region Ω , such that the semigroup defined by (3) has the space of continuous functions C 0 ( Ω ) as an invariant function space. We thus work with the state space E = C 0 ( Ω ) . For now, we take Ω for granted, and define the operator A : D ( A ) E E as
A ϕ ( τ , P ) = τ ϕ ( τ , P ) P ( γ P ϕ ( τ , P ) ) ρ ( τ ) ϕ ( τ , P ) ϕ D ( A ) D ( A ) = { ϕ ( τ , P ) | ϕ , A ϕ E , ϕ ( 0 , τ ) = τ ˇ τ ^ ρ ( τ ) ϕ ( τ , P ψ B ) d τ } .
We can rewrite (3) as the abstract Cauchy equation (with the understanding that S ( t ) E ).
d d t S ( t ) = A S ( t ) , S ( 0 ) = S 0 .
We next show that the operator A is the infinitesimal generator of a strongly continuous semigroup { T ( t ) | t 0 } on E. As it is usual for this kind of equation [42], the construction of the existence of the semigroup is based on the method of characteristics. The characteristic curves of (3) are given by
d t d s = 1 , d τ d s = 1 , d P d s = γ P , d z d s = ( ρ ( τ ) γ ) z .
For t < τ < τ ˇ , we obtain the solution by a pure transport of the initial conditions along the characteristics,
t = τ , τ = s + τ 0 , P = P 0 e γ s , z = S 0 ( τ 0 , P 0 ) e 0 s ρ ( s + τ 0 ) γ d s . S ( t , τ , P ) = S 0 ( τ t , P e γ τ ) e τ t τ ρ ( τ ) γ d τ .
Therewith, for t < τ ˇ , we obtain the boundary values S ( t , 0 , P ) by an appropriate integral over S 0 ,
g ( t , P ) = τ ˇ τ ^ ρ ( τ ) S ( t , τ , P ψ B ) d τ = τ ˇ τ ^ ρ ( τ ) S 0 ( τ t , ( P ψ B ) e γ τ ) e τ t τ ρ ( τ ) γ d τ d τ
such that we are also able to define S ( t , τ , P ) for τ < t < τ ˇ .
t = τ + t 0 , τ = s , P = P 0 e γ s , z = S ( t 0 , 0 , P 0 ) e 0 s ρ ( τ ) γ d τ S ( t , τ , P ) = S ( t τ , 0 , P e γ τ ) e 0 τ ρ ( τ ) γ d τ = g ( t τ , P e γ τ ) e 0 τ ρ ( τ ) γ d τ = τ ˇ τ ^ ρ ( τ ) S 0 ( τ + τ t , ( P e γ τ ψ B ) e γ τ ) e τ + τ t τ ρ ( τ ) γ d τ d τ e 0 τ ρ ( τ ) γ d τ .
Combining the two observations, we are able to explicitly define the semigroup for 0 t τ ˇ ,
T ( t ) S 0 ( τ , P ) = S 0 ( τ t , P e γ τ ) e τ t τ ρ ( τ ) γ d τ t τ τ ˇ τ ^ ρ ( τ ) S 0 ( τ + τ t , ( P e γ τ ψ B ) e γ τ ) e τ + τ t τ ρ ( τ ) γ d τ d τ e 0 τ ρ ( τ ) γ d τ t < τ .
We have an explicit representation of the semigroup for t [ 0 , τ ˇ ] . The condition for strong continuity immediately follows from this representation. We extend that definition to t 0 . For t = n τ ˇ + δ (where 0 δ < τ ˇ and n N 0 ), we define T t S 0 = T τ ˇ n T δ S 0 . The following theorem is a consequence of the results above.
Theorem 1.
Equation (3) defines a strongly continuous semigroup T t on C 0 ( [ 0 , τ ^ ] × R + ) .
As indicated above, the state space ( τ , P ) [ 0 , τ ^ ] × R + is possible but too large, as for P 1 , the degradation is faster than the boosting, and hence the mass of the solution will eventually collect in a bounded region Ω , which turns out to be invariant under the within-host pathogen dynamics.
We obtain an upper bound P * for the pathogen load right after a booster event ( τ = 0 ): starting at ( τ , P ) = ( 0 , P * ) , the characteristic originating in this point is given by P ( τ ) = P * e γ τ . The minimal time interval to the next booster event might take place is τ ˇ , such that P * = P * e γ τ ˇ + ψ B and
P * = ψ B 1 e γ τ ˇ .
We also introduced a minimal pathogen load P * right after a booster event in the same way: now we consider booster events happening at the maximal time span τ ^ , such that
P * = ψ B 1 e γ τ ^ .
Therewith, we can define the invariant region Ω , bounded above and below by the characteristic starting in ( τ , P ) = ( 0 , P * ) , respectively, in ( τ , P ) = ( 0 , P * ) , and extending to the right up to τ ^ (see Figure 2).
Corollary 1.
Let
Ω = { ( τ , P ) | 0 τ τ ^ , P * e γ τ P P * e γ τ } .
If S 0 C 0 ( [ 0 , τ ^ ] × R + ) , s u p p ( S 0 ) Ω , then for all t 0   s u p p ( T t S 0 ) Ω , such that (in slight abuse of notation) T t ( C 0 ( Ω ) ) C 0 ( Ω ) .
Note that in the representation of T t given by Equation (6) the integral extends over values ( τ , P ) outside of Ω , where the solution S ( t , τ , P ) consequently is zero. From now on, we consider the semigroup to be acting on E = C 0 ( Ω ) .
We prove the result of the following theorem in Appendix B.2.
Theorem 2.
If t > 3 τ ^ , the semigroup T t is compact.

3.2. Stationary Solution and Spectral Gap

After establishing the semigroup, we now turn to the spectrum of the infinitesimal generator. For this, we follow once more the standard approach [31,32,39]: we first convert the eigenvalue problem into a fixed point equation, then analyze the fixed point operator (we derive suitable a priori estimates to show compactness), and use the theory of positive operators (Krein–Rutman) to obtain information about eigenvalues and particularly the dominant eigenvalue. The problem herein is that we will not find a strongly positive fixed point operator, and as such, we will have to first regularize our operator before applying the Krein–Rutman theory and related ideas. Afterwards, we will check that the results hold for the de-regularized operator.
However, the existence of a stationary solution is simply a consequence of mass conservation, which implies that f 0 = 1 is an adjoint eigenfunction for eigenvalue 0. Compactness properties imply that there is also an eigenfunction for eigenvalue zero [40]. The involving part is to establish a spectral gap leading to a spectral decomposition of the underlying Banach space and the stability of the stationary solution [43].

3.2.1. Eigenvalues and Fixed Point Operator

To check for the existence of solutions, we consider the eigenvalue problem associated with (3).
τ S + P ( γ P S ) = ( λ + ρ ( τ ) ) S S ( 0 , P ) = g ( P ) = 0 ρ ( τ ) S ( τ , P ψ B ) d τ .
In the same way as above (Equation (5)), we use the method of characteristics to transform Equation (10) into a linear integral operator. We obtain
S ( 0 , P ) = τ ˇ τ ^ S ( 0 , ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ ( λ + ρ ( s ) ) d s d τ .
Note that we know that Ω is invariant for the semigroup, such that s u p p ( S ( 0 , P ) ) [ P * , P * ] . We are led to the definition of an operator K λ . The kernel of this operator e γ τ ρ ( τ ) e 0 τ ( λ + ρ ( s ) ) d s is integrable due to Assumption 1(b). We also know that Ω is invariant for the semigroup, such that s u p p ( K λ [ S ( 0 , P ) ] ) [ P * , P * ] . The integral itself extends over points outside of Ω , where consequently S ( 0 , P ) = 0 .
Definition 1.
Let K λ : C 0 [ P * , P * ] C 0 [ P * , P * ] be defined by
K λ [ g ] ( P ) = τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ .
The integral bounds of K are τ ˇ and τ ^ , with the understanding that g = 0 if ( P ψ B ) e γ τ is outside of [ P * , P * ] (see also Figure 2). To be more precise, we can define τ ̲ ( P ) and τ ¯ ( P ) (see also Proposition A2 in Appendix B.3) such that
K λ [ g ] ( P ) = τ ̲ ( P ) τ ¯ ( P ) g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ .
Particularly, τ ¯ ( P ) = τ ̲ ( P ) for P { P * , P * } , such that for all g C 0 ( [ P * , P * ] )
K λ [ g ] ( P * ) = K λ [ g ] ( P * ) = 0
which leads to technical difficulties as this equation shows that K λ is not strictly positive.
Corollary 2.
An eigenfunction of (10) for eigenvalue λ is a fixed point of K λ . Particularly, a stationary solution of (3) is a fixed point of K 0 .

3.2.2. A Priori Estimates

The first statement of the next proposition is basically a fact that the semigroup is mass-preserving.
Proposition 1.
(a) For g L ( P * , P * ) , g 0 , we find
K 0 [ g ] L 1 ( P * , P * ) = g L 1 ( P * , P * ) .
(b) If g C 0 [ P * , P * ] , then K λ [ g ] C 0 [ P * , P * ] , and there is c = c ( λ ) > 0 such that
K λ [ g ] C 0 [ P * , P * ] c g C 0 [ P * , P * ] .
In Appendix B.3, we give the proof of Proposition 1 and use the a priori estimates to show the compactness of the operator K λ in Proposition A2.

3.2.3. Regularized Operator

The Perron–Frobenius theory of positive operators will be useful in the proof of the existence of a dominant eigenvalue. For the convenience of the reader, we refer to [32,44] for the definition of basic terms. The Krein–Rutman theorem below (Theorems 3 and 4) is an extension of the theorem to infinite dimensional Banach spaces.
Theorem 3
(Krein–Rutman Theorem, [44]). Let S be a total cone, K : S S be a compact positive linear operator and r ( K ) > 0 . Then, r ( K ) is an eigenvalue of S that corresponds to a positive eigenvector Ψ S + .
Theorem 4
([44]). Let S be a solid cone and K : S S be a compact strongly positive linear operator. Then,
  • r ( K ) > 0 , r ( K ) is a simple eigenvalue with an eigenvector in the non-empty interior S o and no other eigenvalue has a positive eigenvector;
  • λ < r ( K ) ∀ eigenvalues, λ r ( K ) .
Furthermore, the following theorem is well known (e.g., [45], Proposition 1.4).
Theorem 5.
If K i are linear positive operators on the same Banach space and K 1 K 2 , then K 1 K 2 and r ( K 1 ) r ( K 2 ) .

Regularized Operator K λ ε

We regularize the integral operator by replacing K λ with a convex combination of K λ and a strictly positive rank 1 operator, which also preserves the L 1 norm of positive functions.
Definition 2.
Let Q : C 0 [ P * , P * ] C 0 [ P * , P * ] be the rank-1 operator
Q [ g ] = 1 P * P * P * P * g ( τ ) d τ
For ε [ 0 , 1 ] , introduce K λ ε : C 0 [ P * , P * ] C 0 [ P * , P * ] by
K λ ε [ g ] = ( 1 ε ) K λ [ g ] + ε Q [ g ] .
Furthermore, we introduce
Λ ε = { λ C | g C ( [ P * , P * ] , C ) : K λ ε [ g ] = g } .
From the definition and our knowledge about K λ , we immediately obtain the following corollary.
Corollary 3.
(a) For ε [ 0 , 1 ] , λ C the operator K λ ε is compact, and for ε ( 0 , 1 ) , λ R strongly positive;
(b) The map ( λ , ε ) K λ ε is continuous with regard to the operator norm;
(c) Furthermore, for g 0 , we have
K 0 ε [ g ] L 1 ( P * , P * ) = g L 1 ( P * , P * ) .
The adjoint eigenfunction of K 0 ε for eigenvalue 1 is f 0 ε = 1 , independently on ε [ 0 , 1 ] .
Note that Λ 0 coincides with the point spectrum of the infinitesimal generator.

Eigenvalues

Later, we will show that the semigroup has a non-negative stationary solution, that is, that K 0 has a fixed point. To prepare for that result, we show that our regularized operator K 0 ε already has a positive fixed point.
Theorem 6.
For λ = 0 , K λ ϵ has a fixed point for all ε [ 0 , 1 ] , which is positive for ε ( 0 , 1 ] .
Proof. 
Corollary 3 indicates that there is an adjoint eigenfunction of K 0 ϵ for eigenvalue 1, which implies using the compactness of the operator that there also is an eigenfunction for this eigenvalues. For ε > 0 , the Krein–Rutman theorem implies the positivity of the eigenfunction. □

3.2.4. De-Regularization and Spectral Gap

Uniqueness of Fixed Point for λ = 0

Note that the next theorem looks like the Krein–Rutman Theorem. However, as K 0 is not strictly positive, we need to prove that theorem via approximations of K λ by K λ ε . Particularly, for λ = 0 , we already know that there is an eigenfunction. The additional information in this special case is the non-negativity of the eigenfunction.
Theorem 7.
For λ R , the spectral radius of K λ is an eigenvalue with a positive eigenfunction. Particularly, K 0 has a non-negative fixed point and spectral radius 1.
Proof. 
There are positive functions g λ ε C 0 [ P * , P * ] with K λ [ g λ ε ] = g λ ε (Theorem 6), that is,
g λ ε ( P ) = ( 1 ε ) K λ [ g λ ε ] ( P ) + ε Q [ g λ ε ] .
Let λ R be arbitrary, fixed. We normalize the eigenfunction to g λ ε L 1 ( P * , P * ) = 1 . Due to Assumption 1(b), the function ρ ( τ ) e 0 τ ρ ( s ) d s is bounded (supremum finite), such that
K λ [ g λ ε ] ( P ) C τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e ( γ λ ) τ d τ C g ε L 1 ( P * , P * ) = C .
and
g λ ε C 0 [ P * , P * ] = K λ [ g λ ε ] C 0 [ P * , P * ] ( 1 ε ) + ε .
The C 0 norm of the family B = { g λ ε | ε ( 0 , 1 ) } is uniformly bounded. Then, the set K 0 [ B ] is also bounded in C 0 , 1 and relatively compact in C 0 , such that we find a subsequence g λ ε n in C 0 , ε n 0 for n 0 , that converges in C 0 to a function g λ C 0 . As g λ ε 0 , also g λ 0 . We need to exclude that g λ = 0 .
As the topology C 0 is stronger than the L 1 topology, the sequence also converges in L 1 , and hence, (the L 1 norm of g λ ε is 1) also g λ L 1 ( P * , P * ) = 1 , such that g λ ¬ 0 . Additionally, from continuity, K λ [ g λ ] = g λ . We establish the existence of a non-negative fixed point.
The function g λ is also an eigenfunction for the spectral radius of K λ : K λ is the limit of the family K λ ε of compact operators depending (with regard to the operator norm) continuously on ε ; as K λ ε [ g λ ] = r ( K λ ε ) g λ ε , this equation carries over to ε = 0 . For λ = 0 we have r ( K 0 ε ) = 1 and thus also r ( K 0 ) = 1 . □
Then, we ensure that the fixed point of the K 0 solution still is unique, as it is for K 0 ε for ε > 0 . Since we used the limit ε 0 to construct a fixed point, it is not clear if possibly two eigenvalues merged in the eigenvalue 1, such that we have a higher dimensional eigenspace. Ultimately, we use the knowledge that the underlying stochastic process mixes well enough to prevent a non-unique invariant measure.
Proposition 2.
If the space of fixed points of K 0 has at least dimension 2, then there is a fixed point that changes sign.
Proof. 
Assume there are two different fixed points g 1 , g 2 which are independent (no α , β R , | α | + | β | > 0 and α g 1 + β g 2 = 0 ). Both functions are non-zero; without restriction, both are non-negative (otherwise, we are done). Furthermore, again without restriction, g i L 1 = 1 .
Consider g = g 1 g 2 ¬ 0 as g 1 g 2 . g is again a fixed point. If g assumes positive and negative values, we are done. Otherwise, either g 0 or g 0 .
If g 0 , then
g 1 g 2 0 .
As g i 0 , g i C 0 , and the L 1 -norm of g 1 and g 2 are equal, we conclude that g 1 = g 2 , which is a contradiction. The second case g 0 gives us a contradiction by the parallel argument. □
The next proposition is a way to express that the underlying Markov process is well mixing.
Lemma 1.
Let g 0 , p 0 s u p p ( g ) . Then, for n N and λ R
s u p p ( K λ n [ g ] ) [ ( 1 e k γ τ ^ ) P * + e k γ τ ^ p 0 , ( 1 e k γ τ ˇ ) P * + e k γ τ ˇ p 0 ] .
Proof. 
Let p 0 s u p p ( g ) . Then, { P [ P * , P * ] : τ [ τ ˇ , τ ^ ] : ( P ψ B ) e γ τ = p 0 } s u p p ( K λ [ g ] ) . That is,
[ ψ B + e γ τ ^ p 0 , ψ B + e γ τ ˇ p 0 ] [ P * , P * ] s u p p ( K λ [ g ] ) .
With the same argument, we find
[ ψ B + ψ B e γ τ ^ + e 2 γ τ ^ p 0 , ψ B + ψ B e γ τ ˇ + e 2 γ τ ˇ p 0 ] [ P * , P * ] s u p p ( K λ 2 [ g ] )
and, if we iterate k times with operator K,
[ a k , b k ] s u p p ( K λ k [ g ] )
where
a k = = 0 k 1 ψ B e γ τ ^ + e k γ τ ^ p 0 = ψ B 1 e k γ τ ^ 1 e γ τ ^ + e k γ τ ^ p 0 = ( 1 e k γ τ ^ ) P * + e k γ τ ^ p 0 b k = = 0 k 1 ψ B e γ τ ˇ + e k γ τ ˇ p 0 = ψ B 1 e k γ τ ˇ 1 e γ τ ˇ + e k γ τ ˇ p 0 = ( 1 e k γ τ ˇ ) P * + e k γ τ ˇ p 0 .
The boundaries of the interval we obtained is a convex combination between p 0 and P * (resp. P * ). We find that the support of any positive function expands under iteration with K, and becomes [ P * , P * ] after an infinite number of iterations. Unfortunately, the operator is not strictly positive, as K λ [ g ] ( P * ) = K λ [ g ] ( P * ) = 0 , and thus, for point measures μ with s u p p ( μ ) { P * , P * } , the pairing μ , K λ n [ g ] = 0 for all n N . In contrast, the proposition indicates that K λ is semi-supporting in the L 2 -setting.
Theorem 8.
The eigenspace of K 0 for eigenvalue 1 is one-dimensional.
Proof. 
If this is not the case, we have an eigenfunction g C 0 that changes sign (Proposition 2). That is, we find two non-negative functions g ± C 0 , both not identically zero, with
g = g + g , s u p p ( g + ) s u p p ( g ) = , s u p p ( g + ) , s u p p ( g ) .
As K 0 is linear,
K 0 [ g + ] K 0 [ g ] = K [ g ] = g = g + g .
Furthermore, as g ± 0 , we know that P * P * K 0 [ g ± ] ( P ) d P = P * P * g ± ( P ) d P .
Let us focus on g + . We know that the support of g + is strictly smaller than [ P * , P * ] , as s u p p ( g ) . We also know that (Lemma 1) s u p p ( K 0 [ g + ] ) is strictly larger than that of g + . As the integral of g + is preserved by K 0 , it is not possible that g + K 0 [ g + ] . There is x 0 s u p p ( g + ) with
g + ( x 0 ) > K 0 [ g + ] ( x 0 ) .
Note that g ( x 0 ) = g + ( x 0 ) as x 0 s u p p ( g + ) . Therewith,
g + ( x 0 ) > K 0 [ g + ] ( x 0 ) K 0 [ g + ] ( x 0 ) K 0 [ g ] ( x 0 ) = K 0 [ g ] ( x 0 ) = g ( x 0 ) = g + ( x 0 )
which is a contradiction. □

Spectral Gap

In the last part of this section, we establish the spectral gap. We do that in two steps: First, we establish the dominance of the eigenvalue λ d = 0 in Λ 0 , then we exclude that the real parts of a sequence of elements in Λ 0 can approximate λ d = 0 .
For general real λ , we first show that the spectral radius of K λ is always larger zero.
Lemma 2.
For λ R arbitrary fixed, we find that r ( K λ ) > 0 .
Proof. 
We use Theorem 5, and construct a positive operator K ˜ which yields a lower bound for r ( K λ ) . Thereto, we rewrite K λ as
K λ [ g ] ( P ) = τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ = 1 γ ( P ψ B ) max { ( P ψ B ) e γ τ ˇ , P * } min { ( P ψ B ) e γ τ ^ , P * } g ( x ) ρ ( τ ( x ; P ) ) e 0 τ ( x ; P ) λ + ρ ( s ) d s d x
with
x = ( P ψ B ) e γ τ τ = τ ( x ; P ) = 1 γ ln x P ψ B .
Now we check for points in [ P * , P * ] that are in the integration region
[ max { ( P ψ B ) e γ τ ˇ , P * } , min { ( P ψ B ) e γ τ ^ , P * } ]
and where the integral weight is strictly positive, that is, τ ( x ; P ) ( τ ˇ , τ ^ ) . Simple computations show that for P * < P < P * , both conditions are given. Choose P ¯ = ( P * + P * ) / 2 as a reference point. Then, there are δ 1 , δ 2 > 0 such that
K λ [ g ] ( P ) P ¯ δ 1 P ¯ + δ 1 g ( x ) δ 2 d x = : K ˜ [ g ] ( P ) .
Hence, r ( K λ ) r ( K ˜ λ ) . As K ˜ [ g ] is a positive rank one operator (compact) with eigenfunction g ( x ) = χ [ P ¯ δ 1 , P ¯ + δ 1 ] ( P ) (this is the only eigenfunction for an eigenvalue 0 ), the spectral radius (and only positive eigenvalue) is given by r ( K ˜ ) = 2 δ 1 δ 2 > 0 . □
Proposition 3.
For λ R , the operator K λ has a fixed point if and only if λ = 0 . Furthermore, r ( K λ ) is strictly increasing in λ.
Proof. 
We already know that r ( K 0 ) = 1 . Furthermore, the eigenfunction g λ and the adjoint eigenfunction f λ of K λ for eigenvalue r ( K λ ) are non-negative. Due to Lemma 1, g λ > 0 in the open interval ( P * , P * ) . If the support of f λ is a subset of { P * , P * } , then there are a , b R such that
f λ , K λ [ f ] ( P ) = a K λ [ g ] ( P * ) + b K λ [ g ] ( P * ) = 0
according to Equation (12). That is, in this case, f λ is an adjoint eigenfunction for eigenvalue 0, which contradicts the fact that ρ ( K λ ) > 0 (Lemma 2). Thus, the intersection of ( P * , P * ) and the support of f λ is non-void, and hence f λ , g μ > 0 for all λ , μ R .
Now, we use an argument by Heijmans [31]: Let λ , μ R , λ > μ and g non-negative.
K μ [ g ] ( P ) = τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ μ + ρ ( s ) d s d τ e ( λ μ ) τ ˇ τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ = e ( λ μ ) τ ˇ K λ 0 [ g ] ( P )
Taking g = g μ , we obtain r ( K μ 0 ) g μ ( P ) = K μ 0 [ g μ ] ( P ) e ( λ μ ) τ ˇ K λ 0 [ g μ ] ( P ) . If we take the duality pairing with the positive eigenfunctional f λ on both sides, and use that f λ , g μ > 0 , we obtain
r ( K n 0 ) e ( λ n ) τ ˇ r ( K λ 0 )
Hence, λ r ( K λ 0 ) is continuous and strictly decreasing on R . □
Theorem 9.
If λ Λ 0 , λ λ d = 0 , then R e λ < λ d = 0 .
Proof. 
We basically adapt the argument ([31], Theorem 6.13) Suppose that λ Λ 0 and there is a corresponding eigenfunction g λ C 0 ( [ P * , P * ] , C ) such that K λ [ g λ ] = g λ . Then
| g λ | = | K λ [ g λ ] | = | τ ˇ τ ^ g λ ( ( P ψ B ) e γ τ ) ρ ( τ ) e γ τ e 0 τ λ + ρ ( s ) d s d τ | τ ˇ τ ^ | g λ ( ( P ψ B ) e γ τ ) | e γ τ ρ ( τ ) e 0 τ ( λ ) + ρ ( s ) d s d τ = K ( λ ) [ | g λ | ] .
If we iterate with K ( λ ) , we obtain for n N
K ( λ ) n [ | g λ | ] K ( λ ) n + 1 [ | g λ | ] .
Now, we know that there is a non-negative eigenfunctional f ( λ ) corresponding to the eigenvalue r ( K ( λ ) ) of K ( λ ) . We know (by the argument in the proof of Proposition 3) that the support of f ( λ ) has a non-zero intersection with ( P * , P * ) .
f ( λ ) , K ( λ ) n [ | g λ | ] f ( λ ) , K ( λ ) n + 1 | g λ | = r ( K ( λ ) ε ) f ( λ ) , K ( λ ) n [ | g λ | ] .
For a sufficiently large n, due to the expansion property of the support under iteration with K ( λ ) (Lemma 1), we can ensure that f ( λ ) , K ( λ ) n [ | g λ | ] > 0 , and thus, r ( K ( λ ) ε ) 1 . Since r ( K λ ε ) , λ R is a non-increasing function (Lemma 1) and r ( K λ d ε ) = 1 , it implies that ( λ ) λ d .
Suppose λ = λ d + i η Λ 0 , we show that η = 0 . As ( λ ) = λ d , we proceed
| g | = | K λ [ g ] | = | τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ | τ ˇ τ ^ | g ( ( P ψ B ) e γ τ ) | e γ τ ρ ( τ ) e 0 τ ( λ ) + ρ ( s ) d s d τ = K λ d ε [ | g | ] ,
that is, K λ d [ | g | ] | g | . Assume that K λ d [ | g | ] > | g | . As we know that the adjoint eigenvalue of K λ d (recall λ d = 0 ) for eigenvalue 1 is f 0 = 1 (identical 1 on [ P * , P * ] , see Corollary 3), we have
1 , | g | = 1 , K λ d [ | g | ] > 1 , | g |
which is a contradiction. Thus, K λ d [ | g | ] = | g | .
If we let g d to be the eigenfunction corresponding to the eigenvalue r ( K λ d ) = 1 , we can write that | g | = c g d for some constant c which we may assume to be one (the eigenspace of K λ d for eigenvalue 1 is one-dimensional, Theorem 8). Hence g ( P ) = g d ( P ) e i ζ ( p ) for some real-valued function ζ ( P ) . If we substitute this relation into | K λ [ g ] | = | g | = g p = K λ d [ g d ] , we obtain
τ ˇ τ ^ g d ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ d + ρ ( s ) d s d τ = | τ ˇ τ ^ g d ( ( P ψ B ) e γ τ ) e γ τ e i ζ ( p ) ρ ( τ ) e 0 τ λ d + i η + ρ ( s ) d s d τ | .
From ([32], Lemma 6.12), there exists a constant b C , | b | = 1 such that ζ ( P ) η τ = b . Substituting this relation in K λ [ g ] = g , we obtain
e i b τ ˇ τ ^ ρ ( τ ) g d [ ( P ψ B ) e γ τ ] e γ τ e 0 τ λ d + ρ ( s ) d s d τ = g d ( p ) e i ζ ( p ) .
Thus, e i b K λ d [ g d ] ( P ) = g d ( P ) e i ζ ( P ) , which in turn implies b = ζ ( P ) , such that η = 0 . Hence, there is no element with real part λ d apart from λ d itself. □
We then show that the spectral gap follows from an argument by Gyllenberg and Heijmans [39].
Theorem 10.
There exists δ > 0 , such that for all λ Λ 0 , λ λ d = 0 , it is true that
( λ ) < λ d δ = δ .
Proof. 
We know that the desired inequality is true for δ = 0 . If no such δ > 0 exists, then there is a sequence of elements λ n Λ 0 with ( λ n ) λ d = 0 . For each λ n , there is a fixed point g n of K λ n , which in turn defines an eigenfunction S n ( τ , P ) of (4). This eigenfunction is also an eigenfunction of T t for eigenvalue e λ n t . As the real values of λ n converge to λ d = 0 , the spectrum of the operator T t has an accumulation point. Since T t is compact for t > 3 τ ^ (Theorem 2), that is impossible. □
From here, as it is standard for the analysis of aggregation–fragmentation equations [39,46,47], a result by Webb [43,47] immediately implies information about the asymptotic behavior of the semigroup: the underlying Banach space can be decomposed by a spectral projection into the eigenspace of the dominant eigenvalue 0, and a remaining part. If Σ is the stationary solution, then (as 1 is the adjoint eigenfunction), the spectral projector is given by Π [ f ] = 1 , f Σ . The remaining part ( I Π ) [ T t S 0 ] will tend to zero for t . To sum it up, we find:
Theorem 11.
Let Σ C 0 ( Ω ) denote the non-negative stationary solution of (3), normalized to Σ L 1 ( Ω ) = 1 . Consider T t S 0 for a non-negative, non-trivial initial condition S 0 C 0 ( Ω ) . Then, with R ˜ ( t ) = ( I Π ) [ T t S 0 ] ,
T t S 0 = Σ 1 , S 0 + R ˜ ( t )
and R ˜ ( t ) 0 for t exponentially fast in C 0 .

4. Reduced Model

We use the theory from Section 3 above to reduce the dimension of the model, and investigate the behavior of the resulting equations by numerical simulation.

4.1. Fast-Slow Analysis

We intend to use the singular perturbation theory. For an intuitive introduction, see [48], while a formal approach is given in [49]. As above, we rewrite the semigroup as d d t S = A S . Then, there is a stationary solution Σ ( τ , P ) C 0 ( Ω ) (an eigenfunction of A for eigenvalue 0). Please note that the generator A depends on the environmental pathogen load B (which we take to be fixed for the moment), such that the eigenfunction Σ also depends on B. The adjoint eigenfunction is simply Σ * ( τ , P ) = 1 , as the semigroup is mass-preserving,
d d t Ω Σ * ( τ , P ) S ( t , τ , P ) d ( τ , P ) = d d t Ω S ( t , τ , P ) d ( τ , P ) = 0 .
Then, Σ and Σ * are the right and left eigenfunctions of A for eigenvalue 0. Note that Σ * ( τ , P ) = 1 is independent of B.
We define the spectral projector Π f ( τ , P ) = Σ ( τ , P ; B ) Ω f ( τ , P ) d ( τ , P ) = Σ 1 , f . Then, A Π = Π A and
d d t Π S = A Π S = A Σ 1 , S = 0 , d d t ( I Π ) S = A ( I Π ) S .
The definition of the projector Π implies that
Π S ( t , τ , P ) = Σ ( τ , P ; B ) s ( t ) w i t h s ( t ) = Ω S ( t , τ , P ) d ( τ , P ) .
We use these two projectors to define a new coordinate system that disentangles the slow and fast dynamics, S = Π S + ( I Π ) S = Σ s + ( I Π ) S . Until now, the considerations have been made under the assumption that B is fixed. Even if B (slowly) varies with time, we can still use these new coordinates, but now the projectors Π and ( I Π ) do not commute with the time derivative. In the case of I Π , we simply use the chain rule and the fact that B scales with O ( ϵ ) to obtain
( I Π ) d d t S = d d t ( ( I Π ) S ) + O ( ϵ ) .
Hence, multiplying the first model equation d d t S = A S ϵ β ( P ) S by ( I Π ) from the left yields
d d t ( I Π ) S = A ( I Π ) S ϵ ( I Π ) [ β ( P ) ( Σ s + ( I Π ) S ) ] + O ( ϵ ) .
In the case of Π , we look slightly more closely, using the fact that the left eigenfunction Σ * = 1 does not depend on B, such that
Π d d t S ( t , τ , P ) = Σ ( τ , P ) Ω d d t S ( t , τ , P ) d ( τ , P ) = Σ ( τ , P ) d d t s ( t ) .
If we integrate this equation over Ω , we have
Ω Π d d t S ( t , τ , P ) d ( τ , P ) = d d t s ( t ) .
That is, multiplying d d t S = A S ϵ β S by Π from the left, and integrating over gives us an equation for s ( t ) (recall Π A = 0 )
d d t s ( t ) = ϵ s ( t ) Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P ) ϵ 1 , β ( P ) ( I Π ) S .
All in all, our model becomes
d d t s ( t ) = ϵ s ( t ) Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P ) ϵ 1 , β ( P ) ( I Π ) S d d t ( I Π ) S = A ( I Π ) S ϵ ( I Π ) [ β ( P ) ( Π S + ( I Π ) S ) ] + O ( ϵ ) d d t I = ϵ Ω β ( P ) ( Σ ( τ , P ; B ) s ( t ) + ( I Π ) S ) d P d τ ϵ α I d d t R = ϵ α I d d t B = ϵ ( a + ξ I σ B ) .
That is, only ( I Π ) S is the fast variable, while all other variables are slow.
Fast system. If we take ϵ = 0 , we find for, the fast variable,
d d t ( I Π ) S = A ( I Π ) S .
Due to the spectral gap, we know that ( I Π ) S 0 . Therefore, ( I Π ) S = 0 forms the slow manifold.
Slow system. Now we use the slow time T = ϵ t , and the result for the slow manifold, and obtain the reduced model
d d T s = s Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P )
d d T I = s Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P ) α I
d d T R = α I
d d T B = a + ξ I σ B .
The model suggests that the distribution of the susceptible population is always in its (quasi) steady state, and hence the force of infection becomes Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P ) .

Behavior of the Reduced Model: A Simulation Study

We show the dependency of Σ on B.
Lemma 3.
If Σ ( τ , P ) is a stationary solution of (3) for ψ B = 1 with Ω Σ d ( τ , P ) = 1 , then
S ( τ , P ; B ) = 1 ψ B Σ ( τ , P / ( ψ B ) )
is a stationary solution for a given value of ψ B with Ω S d ( τ , P ) = 1 .
Proof. 
We suppress the multiplicative factor 1 ψ B that only ensures that the norm is preserved. Let P ˜ = P / ( ψ B ) .
τ S + P ( γ P S ) = τ Σ ( τ , P / ( ψ B ) ) + P ( γ P Σ ( τ , P / ( ψ B ) ) ) = τ Σ ( τ , P ˜ ) + P ˜ ( γ P ˜ Σ ( τ , P ˜ ) ) = ρ ( τ ) Σ ( τ , P / ( ψ B ) ) = ρ ( τ ) S ( τ , P ; B )
and for the boundary value we obtain
τ ˇ τ ^ ρ ( τ ) S ( τ , ( P ψ B ) ; B ) d τ = τ ˇ τ ^ ρ ( τ ) Σ ( τ , ( P ψ B ) / ψ B ) d τ = Σ ( 0 , P ˜ ) = S ( 0 , P ; B ) .
As we are not interested in S ( τ , P ) , but only in the marginal distribution S ( τ , P ) d τ , it is convenient to use the underlying stochastic process to obtain a numerical approximation of the function (kind of Monte-Carlo integration): we determine the realizations of the random variable T that is distributed according to the time between two booster events. To keep things simple, we use a uniform distribution on [ τ ˇ , τ ^ ] . With this aspect, it is straightforward to compute a realization of the time course of the pathogen load P t (also see Figure 1). After dismissing a burn-in phase, we sample at discrete time points the values of P t ; the histogram of those values is proportional to S ( τ , P ) d τ (see Figure 3).
The transition from susceptible to infected corresponds to a massive replication of pathogens, which (at a short time scale), cannot be controlled any more by the immune system. A branching process, going from the sub- to the supercritical parameter range, can suit as a toy model for that process. Inspired by that idea, we use the probability to take off for a branching process as β ( P ) and define β ( P ) = β 0 max { 0 , 1 π 0 / P } (Figure 3). With these two components and equipped with Lemma 3, we are able to determine the incidence’s dependency on B. As we have an SIR model, we do not have stationary solutions to address. However, if we assume that the class of susceptible is reduced relatively slowly, we can identify parameter combinations where infections take off: given that a and I, the environmental pathogen load B asymptotically tends to
B = ( a + ξ I ) / σ .
If we feed this pathogen level into the incidence, we obtain for the r.h.s. of I
F ( a , I ) = s Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P ) | ( a + ξ I ) / σ α I .
The incidence grows if F ( a , I ) > 0 , and decreases if F ( a , I ) < 0 . We need to emphasize that, here, we use a kind of quasi-steady state for B. In more realistic cases, B and I will change on similar time scales, such that F ( I , a ) mainly yields a heuristic about the behavior to expect, and not a rigorous threshold argument.
Although cholera is the infection that inspired this model, we are rather interested in the model structure than in realistic parameters. Therefore, we do not even try to determine parameters suited for cholera but focus on the discussion of the model for a fairly arbitrary set of parameters (see Appendix A), which then yield Figure 4. On the left side of the figure, we find a structure resembling typical bistable behavior. This is not a coincidence: If a is small, without additional shedding by infected individuals, the pathogen load in the environment is too small to trigger an epidemic. However, for a certain parameter range of a, the positive feedback (infected shed pathogens which in turn additionally infect further individuals) is able to trigger an epidemic: If I is small, the incidence stays low, if I is large, the number of infections increases further. Only if a is large, such that the pathogen load in the environment becomes supercritical, can an epidemic happen without a considerable number of initial infections.
Accordingly, for the three scenarios shown on the right in Figure 4, all parameters are fixed except a, and I ( 0 ) = B ( 0 ) = 0 in all cases. For the chosen parameter values, we observe that there is a threshold value of a between 20.5 and 20.8, which determines whether an epidemic can take place. For example, for smaller a (a = 15), the pathogen load tends to a positive equilibrium, while the incidence is identically zero. The pathogen cannot accumulate sufficiently high in the susceptible host to trigger an infection. If we further increase a ( a = 20.5 ), then we again have a positive equilibrium (approximately) for the bacterial compartment, but now there is also a small but positive incidence. The pathogen concentration is large enough to trigger infections once in a while, but we are still in a region, where this small incidence is not able to cause a larger outbreak. Only if we further increase a ( a = 20.8 ), then we reach the situation where an epidemic takes place: in the initial time interval, the bacteria reach a plateau, where they are able to increase the prevalence to a certain level. Then, the combined pathogen load due to natural sources (described by a) and shedding is able to trigger the positive feedback, and we obtain a steep increase in the prevalence.
This model behavior can be compared with the incidence of cholera: due to the life cycle of V. cholerae, their abundance undergoes seasonal changes, which can be modeled by the different values of a. If they are barely present, the incidence of cholera is zero. However, there are also situations, where we have a low, fluctuating incidence and there are distinct outbreaks of larger cholera epidemics. The reason for these different figures is under discussion, but the main idea is the change of V. cholerae abundance, as our model also seems to hint.
It might also be interesting that at least some cholera epidemics share in the onset the sharp increase in cases, followed by a low decline, e.g., the outbreak in Katsina (Nigeria) in 1982 [50]. An interpretation might be that the pathogen load rapidly increases on a high level, also increasing the prevalence. In contrast to SIR-type models, where the increase is stopped by the depletion of the susceptible class, we rather find here a saturation in the incidence function: also, if B becomes arbitrarily large, the force of infection Ω β ( P ) Σ ( τ , P ; B ) d ( τ , P ) stays bounded ( β ( P ) is bounded). This property of our model’s force of infection is in line with other model approaches [1,18].

5. Discussion

The main contribution of this work is that we have developed a novel framework for linking the within-host and between-host dynamics of cholera. We have also provided a rigorous mathematical analysis of this novel structure that can be applied to other multi-scale disease models. We did this by assigning each susceptible individual with a pathogen load that grows through the uptake of contaminated food and water (booster event) from the environment and decreases between two booster events through elimination by the immune system. The transition from susceptible to infected took place at a certain pathogen-load dependent rate. This rate was only positive if a critical pathogen threshold was surpassed. We further took the population dynamics to occur on a slower time scale in comparison to the within-host dynamics. We analyzed the model on the fast time scale and showed the existence of an invariant solution. This was performed by the construction of a semigroup and analysis of the spectrum of its infinitesimal generator. We then used the results obtained from the spectral analysis to reduce the dimension of the original multi-scale model on the slow time scale to an SIR model and performed numerical simulations on the resulting model.
Our work provided a mathematical framework for utilizing the methods used in aggregation–fragmentation models [32,35] in the analysis of multi-scale disease models. The interplay between the booster events and the immune responses when the pathogen load is sub-critical can be compared to the mechanisms of integrate and fire neural models, which have been used before to address the spread of infections [51,52]. The underlying mechanism of our model is a velocity-jump process, where the jumps are fostered by the booster events. Velocity-jump processes are close to age structured models, as both consist of a transport equation with a nonlocal term. This nonlocal term, however, is focused at age 0 in age-structured models (the boundary condition—which can be newborns or new infections in other cases—addresses age zero, see, e.g., [20,21])), but is distributed over all states in velocity-jump processes. On that, velocity-jump processes are often more complex to analyze [33,36].
The numerical simulations showed three different parameter regions for the infection dynamics: In the first region, the V. cholerae bacteria were present in the environment but there were no cholera cases (incidence is zero); in the next region, few cholera cases would occur once in a while (the incidence was low such that larger outbreaks could not occur), and in the third parameter region, a full outbreak of the disease could be observed. Through the results of these simulations, we can qualitatively explain the different outcomes that occur when the infection is introduced in a region. The abundance of the bacteria in the environment was seen as a driving force in the occurrence of an epidemic which was in line with other studies [1] that highlight the role of environmental reservoirs in the infection process. Nevertheless, we have to emphasize that our set of parameters is not specific to cholera and therefore we only obtain an overview of what could happen. The incidence term derived from our model is close to other saturated incidence terms used for the disease [1,18,19] which have been found to be more realistic. However, our incidence can drop to zero if the pathogen levels are low which, contrasts with the other models that maintain positive values.
Our study had several limitations: To reduce the complexity of the analysis, we did not explicitly account for the role of immune responses which would have made our within-host model more robust. We also focused on structuring the susceptible class even though more information can be gained from a postulation of the same for the infected class. Finally, the aim of this study was to gain insights on the model structure and the mathematical theory and therefore we did not use parameters specific to cholera in the simulations, nor did we perform parameter estimation. Parameters derived from empirical findings could be considered to make the model more sophisticated. We intend to investigate the aforementioned issues in our future work.

Author Contributions

Conceptualization, B.M. and J.M.; methodology, B.M., J.M. and Z.F.; formal analysis, B.M. and J.M.; writing—original draft preparation, B.M. and J.M.; writing—review and editing, B.M., J.M. and Z.F. All authors have read and agreed to the published version of the manuscript.

Funding

B.M. is funded by the German Academic Exchange Service (DAAD); J.M. and B.M. are supported by a grant from the Deutsche Forschungsgemeinschaft (DFG) through the TUM International Graduate School of Science and Engineering (IGSSE), GSC 81, within the project GENOMIE-QADOP.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors thank Daniel Matthes and the TUM GENOMIE_QADOP project members for their helpful discussions of the work.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

We again emphasize that we did not try to obtain realistic parameter values for cholera, as the aim of the present paper was to make a proposal for a modeling approach, and the analysis of this approach. For Figure 4, we used the following, rather arbitrary parameter values.
Table A1. Parameter values.
Table A1. Parameter values.
ParameterValue
τ ˇ 2
τ ^ 10
γ 0.2
α 12
ζ 2000
σ 10
β ( P ) β 0 max { 0 , 1 π 0 / P }
β 0 0.01
π 0 3
s(0)100
I(0), B(0)0

Appendix B

In this Appendix, we provide the proofs of some of the theorems and propositions.

Appendix B.1

In this section, we show that Assumption 1 is satisfied by a wide range of random variables.
Proposition A1.
Let φ C 1 ( [ 0 , τ ^ ) ) , s u p p ( φ ) [ τ ˇ , τ ^ ] , 0 τ ^ φ ( s ) d s = 1 , and assume that there is z C 1 ( R + ) and m > 0 such that φ ( τ ) = z ( τ ) ( τ ^ τ ) m . Then, the hazard rate ρ ( τ ) satisfies Assumption 1.
Proof. 
The conditions (2) are obviously satisfied. We now show (a)–(d) from Assumption 1.
(a)
0 τ ρ ( τ ) d τ = 0 τ φ ( τ ) 1 0 τ φ ( s ) d s d τ = ln 1 0 τ φ ( s ) d s
implies that lim τ τ ^ 0 τ ρ ( τ ) d τ = .
(b) Note that the function ρ ( τ ) e 0 τ ρ ( s ) d s is continuous for τ < τ ^ . If we can show that this function converges to 0 for τ τ ^ , it is already bounded ( sup ( ) < ). We use the formula derived in (a) and obtain
lim τ τ ^ ρ ( τ ) e 0 τ ρ ( s ) d s = lim τ τ ^ φ ( τ ) 1 0 τ φ ( s ) d s e 0 τ d d s ln 1 0 s φ ( τ ) d τ d s = lim τ τ ^ φ ( τ ) 1 0 τ φ ( s ) d s 1 0 τ φ ( s ) d s = lim τ τ ^ φ ( τ ) = 0 .
(c)
τ ˇ τ ^ ρ 2 ( τ ) e 0 τ ρ ( s ) d s d τ = τ ˇ τ ^ φ ( τ ) 1 0 τ φ ( s ) d s 2 1 0 τ φ ( s ) d s d τ = τ ˇ τ ^ φ 2 ( τ ) 1 0 τ φ ( s ) d s d τ .
We discuss the asymptotic of the integrand for τ τ ^ for the case φ ( τ ) = z ( τ ) ( τ ^ τ ) m for some smooth z ( τ ) . For the denominator, we obtain
1 0 τ φ ( s ) d s = τ ˇ τ ^ φ ( s ) d s = τ ˇ τ ^ z ( τ ) ( τ ^ s ) m d s = O ( ( τ ^ s ) m + 1 ) .
Therewith,
φ 2 ( τ ) 1 0 τ φ ( s ) d s = O ( ( τ ^ s ) 2 m ) O ( ( τ ^ s ) m + 1 ) = O ( ( τ ^ s ) m 1 )
which is integrable, as we assume m > 0 .
(d) Note that
τ ˇ τ ^ ρ ( τ ) e 0 τ ρ ( s ) d s d τ = τ ˇ τ ^ | d d τ φ ( τ ) 1 0 τ φ ( s ) d s | 1 0 τ φ ( s ) d s d τ τ ˇ τ ^ | φ ( τ ) | d τ + τ ˇ τ ^ φ 2 ( τ ) 1 0 τ φ ( s ) d s d τ
which is integrable, as seen previously. □

Appendix B.2

In this section, we show the proof of Theorem 2. The Arzel a –Ascoli theorem provides a criterion for checking the compactness of a subset in the space of continuous functions [40]. To prove the compactness of the semigroup, we first iterate the boundary condition and express it as integral over boundary values at earlier times. We then check whether the C 1 estimates of the transformed integral are bounded in C 0 . If this is satisfied, it implies compactness by the Arzel a –Ascoli theorem.
Proof. 
If t > τ ^ , the booster event has already occurred, we have
S ( t , τ , P ) = g ( t τ , P e γ τ ) e 0 τ ρ ( τ ) γ d τ
such that,
S ( t , τ , P ψ B ) = g ( t τ , ( P ψ B ) e γ τ ) e 0 τ ρ ( τ ) γ d τ
and therefore
g ( t , P ) = τ ˇ τ ^ ρ ( τ ) S ( t , τ , P ψ B ) d τ = τ ˇ τ ^ g ( t τ , ( P ψ B ) e γ τ ) ρ ( τ ) e 0 τ ρ ( τ ) γ d τ d τ .
Assuming t > 2 τ ^ , we iterate twice to express the boundary value g ( t , P ) as an integral over the boundary value at earlier times. We already know that g ( t , P ) = τ ˇ τ ^ g ( t τ , ( P ψ B ) e γ τ ) ρ ( τ ) e 0 τ ρ ( s ) γ d s d τ and thus,
g ( t τ , ( P ψ B ) e γ τ ) = τ ˇ τ ^ g ( t τ τ , ( ( P ψ B ) e γ τ ψ B ) e γ τ ) ρ ( τ ) e 0 τ ρ ( s ) γ d s d τ
which yields
g ( t , P ) = τ ˇ τ ^ τ ˇ τ ^ g ( t τ τ , ( ( P ψ B ) e γ τ ψ B ) e γ τ ) ρ ( τ ) e 0 τ ρ ( s ) γ d s d τ ρ ( τ ) e 0 τ ρ ( s ) γ d s d τ = τ ˇ τ ^ τ ˇ τ ^ g ( t τ τ , ( ( P ψ B ) e γ τ ψ B ) e γ τ ) ρ ( τ ) e 0 τ ρ ( s ) γ d s ρ ( τ ) e 0 τ ρ ( s ) γ d s d τ d τ .
To transform the integral, we let
u = t τ τ v = ( ( P ψ B ) e γ τ ψ B ) e γ τ = ( P ψ B ) e r ( τ + τ ) ψ B e r τ
thus
v = ( P ψ B ) e r ( t u ) ψ B e r τ
e r τ = ( P ψ B ) e r ( t u ) v ψ B , τ = τ ( u , v ; t , P ) = 1 γ ln ( P ψ B ) e r ( t u ) v ψ B
and therefore
τ = τ ( u , v ; t , P ) = t u τ ( u , v ; t , P ) .
Then, we find the derivatives of τ ( u , v ; t , P ) and τ ( u , v ; t , P ) with respect to t and P.
t τ ( u , v ; t , P ) = 1 r γ ψ B ( P ψ B ) e r ( t u ) ψ B ( ( P ψ B ) e r ( t u ) v ) = ( P ψ B ) e r ( t u ) ( P ψ B ) e r ( t u ) v .
Since ( P ψ B ) e r ( t u ) v = e r τ , t τ ( u , v ; t , P ) will be bounded for all feasible values of u , v , t , P .
P τ ( u , v ; t , P ) = 1 r ψ B e r ( t u ) ψ B ( ( P ψ B ) e r ( t u ) v ) = 1 γ e r ( t u ) ( P ψ B ) e r ( t u ) v ) .
This derivative is also bounded. Thus, τ ( u , v ; t , P ) is bounded in C 1 with regard to t and P.
Since τ ( u , v ; t , P ) = t u τ ( u , v ; t , P ) , the function τ ( u , v ; t , P ) = t u τ ( u , v ; t , P ) is C 1 with regard to t and P.
Transformation of the integral domain: The integral domain reads ( τ , τ ) [ τ ˇ , τ ^ ] × [ τ ˇ , τ ^ ] . Due to (A1)
t 2 τ ˇ u t 2 τ ^ .
Furthermore, for a given v, we have according to (A2)
v [ ( P ψ B ) e γ ( t u ) ψ B e γ τ ˇ , ( P ψ B ) e γ ( t u ) ψ B e γ τ ^ ]
which implies that the integral boundaries transform to
τ ˇ τ ^ τ ˇ τ ^ . . . . d τ d τ = t 2 τ ˇ t 2 τ ^ ( P ψ B ) e γ ( t u ) ψ B e γ τ ˇ ( P ψ B ) e γ ( t u ) ψ B e γ τ ^ . . . . | det ( τ , τ ) ( u , v ) | d u d v .
Transformation of the infinitesimal: From (A3) and (A4),
τ ( u , v ; t , P ) = t u τ ( u , v ; t , P ) τ ( u , v ; t , P ) = 1 γ ln ( P ψ B ) e r ( t u ) v ψ B .
We find the Jacobian of τ ( u , v ; t , P ) and τ ( u , v ; t , P ) with regard to u and v. Therefore, we first consider their derivatives.
u τ ( u , v ; t , P ) = ( P ψ B ) e γ ( t u ) ( P ψ B ) e γ ( t u ) v v τ ( u , v ; t , P ) = 1 r ( P ψ B ) e γ ( t u ) v u τ ( u , v ; t , P ) = 1 u τ ( u , v ; t , P ) u τ ( u , v ; t , P ) = u τ ( u , v ; t , P ) .
All four derivatives are C 1 with regard to t and P. Thus, the determinant of the Jacobian also
| det J ( u , v ; t , P ) | = | det u τ ( u , v ; t , P ) v τ ( u , v ; t , P ) u τ ( u , v ; t , P ) v τ ( u , v ; t , P ) |
is in C 1 with regard to t and P. Taking all the elements together we obtain
g ( t , P ) = τ ˇ τ ^ τ ˇ τ ^ ρ ( τ ) g ( t τ τ , ( ( P ψ B ) e γ τ ψ B ) e γ τ ) e 0 τ ρ ( s ) γ d s ρ ( τ ) e 0 τ ρ ( s ) γ d s d τ d τ = t 2 τ ˇ t 2 τ ^ ( P ψ B ) e γ ( t u ) ψ B e γ τ ˇ ( P ψ B ) e γ ( t u ) ψ B e γ τ ^ ρ ( τ ( u , v ; t , P ) ) g ( u , v ) e 0 τ ( u , v ; t , P ) ρ ( s ) γ d s ρ ( τ ( u , v ; t , P ) ) e 0 τ ( u , v ; t , P ) ρ ( s ) γ d s | det J ( u , v ; t , P ) | d u d v .
As the derivative of t and P no longer acts on g ( . ) and all the expressions are C 1 with regard to t and P, for t ¯ > 2 τ ^ , we can obtain an estimate of the form
g ( t , P ) C 1 ( [ 2 τ ^ , t ¯ ] × [ 0 , P * ] ) C ( t ¯ ) g ( t , P ) C 0 ( [ 0 , t ¯ ] × [ 0 , P * ] ) C ( t ¯ ) T t S 0 C 0 ( [ 0 , t ¯ ] , C 0 ( Ω ) ) C ( t ¯ ) S 0 C 0 ( Ω ) .
That is, the boundary conditions are smooth for t > 2 τ ^ , and—as the transport of the boundary conditions into the region Ω is also smooth—this inequality implies for t > 3 τ ^ that
T t S 0 C 0 ( [ 3 τ ^ , t ¯ ] , C 1 ( Ω ) ) C ( t ¯ ) S 0 C 0 ( Ω ) .
Note that for t > 3 τ ^ , the integral extends to S 0 . Therefore, the semigroup is eventually compact by Arzel a –Ascoli. □

Appendix B.3

In this section, we provide the proof for the Proposition 1 on a priori estimates. We further use these estimates to show the compactness of operator K λ in Proposition A2.
Proof. 
(Proposition 1). (a) Note that the integrand is non-negative for g 0 . Thus,
P * P * K 0 [ g ] ( P ) d P = P * P * τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ ρ ( s ) d s d τ d P = τ ˇ τ ^ P * P * g ( ( P ψ B ) e γ τ ) d P e γ τ ρ ( τ ) e 0 τ ρ ( s ) d s d τ
where we exchanged the integrals with the understanding (as above) that g becomes zero outside of [ P * , P * ] . We now change the integration variable, x = ( P ψ B ) e γ τ . For any τ ( τ ˇ , τ ^ ) , we have
( P * ψ B ) e γ τ = P * e γ ( τ ^ τ ) < P * a n d ( P * ψ B ) e γ τ = P * e γ ( τ ˇ τ ) > P *
such that P * P * g ( ( P B ) e γ τ ) e γ τ d P = P * P * g ( x ) d x and we proceed
P * P * K 0 [ g ] ( P ) d P = τ ˇ τ ^ P * P * g ( x ) d x ρ ( τ ) e 0 τ ρ ( s ) d s d τ = g L 1 ( P * , P * ) τ ˇ τ ^ ρ ( τ ) e 0 τ ρ ( s ) d s d τ = g L 1 ( P * , P * ) .
(b) If g C 0 [ B , P * ] , then K ^ λ [ g ] is continuous for any P [ P * , P * ] , as the integral kernel is smooth. The a priori estimate follows as before,
| K λ [ g ] ( P ) | τ ˇ τ ^ | g ( ( P ψ B ) e γ τ ) | e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ C ( λ ) τ ˇ τ ^ e γ τ ρ ( τ ) e 0 τ ρ ( s ) d s d τ g C 0 [ P * , P * ] c ( λ ) g C 0 [ P * , P * ] .
Based on this observation, we obtain a proper C 0 , 1 estimate from which we deduce compactness by the Arzel a –Ascoli theorem.
Proposition A2.
There is a c > 0 such that K λ [ g ] C 0 , 1 [ P * , P * ] c g C 0 [ P * , P * ] . The operator K λ : C 0 [ P * , P * ] C 0 [ P * , P * ] is compact.
Proof. 
We already know that K λ [ g ] C 0 [ P * , P * ] c ( λ ) g C 0 [ P * , P * ] . We now estimate | d d P K λ [ g ] ( P ) | by g C 0 [ P * , P * ] .
In the following, we need to use the correct boundaries in the integral defining K λ [ g ] ( P ) (until now, we simply said that g becomes zero if the argument is outside [ P * , P * ] ). For a given P, we compute the τ -value, where the upper (lower) boundary of Ω hits P ψ B . For the upper bound, we have
P ψ B = P * e γ τ τ = 1 γ ln P * P ψ B .
The intersection with the lower boundary is the same expression where we use P * instead of P * . Accordingly, we define
τ ¯ ( P ) = 1 γ ln ( P * / ( P ψ B ) ) for P > ψ B + P * e γ τ ^ τ ^ else , τ ̲ ( P ) = 1 γ ln ( P * / ( P ψ B ) ) for P < ψ B + P * τ ˇ else .
Note that τ ¯ and τ ̲ are Lipschitz continuous, and the derivative is zero for the constant part, respectively,
τ ¯ ( P ) = τ ̲ ( P ) = 1 γ ( P ψ B )
in the non-constant part. As P P * > ψ B , the derivative is bounded. With these preliminaries, we start to estimate.
| d d P K λ [ g ] ( P ) | = | d d P τ ˇ τ ^ g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ | = | d d P τ ̲ ( P ) τ ¯ ( P ) g ( ( P ψ B ) e γ τ ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ | | g ( ( P ψ B ) e γ τ ¯ ( P ) ) e γ τ ¯ ( P ) ρ ( τ ¯ ( P ) ) e 0 τ ¯ ( P ) λ + ρ ( s ) d s τ ¯ ( P ) | + | g ( ( P ψ B ) e γ τ ̲ ( P ) ) e γ τ ̲ ( P ) ρ ( τ ̲ ( P ) ) e 0 τ ̲ ( P ) λ + ρ ( s ) d s τ ̲ ( P ) | + | τ ̲ ( P ) τ ¯ ( P ) g ( ( P ψ B ) e γ τ ) e 2 γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ | .
We estimate the three terms, one after the other.
| g ( ( P ψ B ) e γ τ ¯ ( P ) ) e γ τ ¯ ( P ) ρ ( τ ¯ ( P ) ) e 0 τ ¯ ( P ) λ + ρ ( s ) d s τ ¯ ( P ) | g C 0 [ P * , P * ] e γ τ ^ m a x { 1 , e ( λ ) τ ^ } ρ ( τ ^ ) e 0 τ ¯ ( P ) ρ ( s ) d s | τ ¯ ( P ) | C ( λ ) g C 0 [ P * , P * ]
as τ ¯ ( P ) is Lipschitz-continuous, and ρ ( τ ) e 0 τ ρ ( s ) d s is bounded (Assumption 1). Similarly,
| g ( ( P ψ B ) e γ τ ̲ ( P ) ) e γ τ ̲ ( P ) ρ ( τ ̲ ( P ) ) e 0 τ ̲ ( P ) λ + ρ ( s ) d s τ ¯ ( P ) | C ( λ ) g C 0 [ P * , P * ] .
For the third term, we integrate by parts and proceed
τ ̲ ( P ) τ ¯ ( P ) g ( ( P ψ B ) e γ τ ) e 2 γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ = 1 γ ( P ψ B ) τ ̲ ( P ) τ ¯ ( P ) e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d d τ g ( ( P ψ B ) e γ τ ) d τ = 1 γ ( P ψ B ) e γ τ ¯ ( P ) ρ ( τ ¯ ( P ) ) e 0 τ ¯ ( P ) λ + ρ ( s ) d s g ( ( P ψ B ) e γ τ ¯ ( P ) ) 1 γ ( P ψ B ) e γ τ ̲ ( P ) ρ ( τ ̲ ( P ) ) e 0 τ ̲ ( P ) λ + ρ ( s ) d s g ( ( P ψ B ) e γ τ ̲ ( P ) ) 1 γ ( P ψ B ) τ ̲ ( P ) τ ¯ ( P ) g ( ( P ψ B ) e γ τ ) d d τ e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ .
We again estimate the three terms, where the first two terms follow the same scheme as above,
| 1 γ ( P ψ B ) e γ τ ¯ ( P ) ρ ( τ ¯ ( P ) ) e 0 τ ¯ ( P ) λ + ρ ( s ) d s g ( ( P ψ B ) e γ τ ¯ ( P ) ) | C ( λ ) g C 0 [ P * , P * ]
(note that P ψ B P * ψ B > 0 ) and similarly for the second term. For the third term, we remark
τ ̲ ( P ) τ ¯ ( P ) | d d τ e ( γ λ ) τ ρ ( τ ) e 0 τ ρ ( s ) d s | d τ ( γ + | λ | ) e ( γ + | ( λ ) | ) τ ^ τ ̲ ( P ) τ ¯ ( P ) ρ ( τ ) e 0 τ ρ ( s ) d s d τ + e ( γ + | ( λ ) | ) τ ^ τ ̲ ( P ) τ ¯ ( P ) | ρ ( τ ) | e 0 τ ρ ( s ) d s d τ + e ( γ + | ( λ ) | ) τ ^ τ ̲ ( P ) τ ¯ ( P ) ρ 2 ( τ ) e 0 τ ρ ( s ) d s d τ .
Due to Assumption 1, all three integrals are finite. Hence,
| 1 γ ( P ψ B ) τ ̲ ( P ) τ ¯ ( P ) g ( P ψ B ) d d τ e γ τ ρ ( τ ) e 0 τ λ + ρ ( s ) d s d τ | C ( λ ) g C 0 [ P * , P * ] .
Therewith, the C 0 , 1 estimate for K λ is established, and the compactness is a consequence of the Arzel a –Ascoli theorem. □

References

  1. Codeço, C.T. Endemic and epidemic dynamics of cholera: The role of the aquatic reservoir. BMC Infect. Dis. 2001, 1, 1. [Google Scholar]
  2. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. Contain. Pap. Math. Phys. Character 1927, 115, 700–721. [Google Scholar]
  3. Kermack, W.O.; McKendrick, A.G. Contributions to the mathematical theory of epidemics. II.—The problem of endemicity. Proc. R. Soc. Lond. Ser. Contain. Pap. Math. Phys. Character 1932, 138, 55–83. [Google Scholar]
  4. Garira, W. A complete categorization of multiscale models of infectious disease systems. J. Biol. Dyn. 2017, 11, 378–435. [Google Scholar] [CrossRef]
  5. Gog, J.R.; Pellis, L.; Wood, J.L.; McLean, A.R.; Arinaminpathy, N.; Lloyd-Smith, J.O. Seven challenges in modeling pathogen dynamics within-host and across scales. Epidemics 2015, 10, 45–48. [Google Scholar] [CrossRef] [PubMed]
  6. Wang, X.; Wang, J. Modeling the within-host dynamics of cholera: Bacterial–viral interaction. J. Biol. Dyn. 2017, 11, 484–501. [Google Scholar] [CrossRef]
  7. Anderson, R.M.; May, R.M. Infectious Diseases of Humans: Dynamics and Control; Oxford University Press: Oxford, UK, 1992. [Google Scholar]
  8. Martcheva, M. An Introduction to Mathematical Epidemiology; Springer: Berlin/Heidelberg, Germany, 2015; Volume 61. [Google Scholar]
  9. Hu, Z.; Ma, W.; Ruan, S. Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math. Biosci. 2012, 238, 12–20. [Google Scholar] [CrossRef]
  10. Feng, Z.; Velasco-Hernandez, J.; Tapia-Santos, B. A mathematical model for coupling within-host and between-host dynamics in an environmentally-driven infectious disease. Math. Biosci. 2013, 241, 49–55. [Google Scholar] [CrossRef]
  11. Feng, Z.; Velasco-Hernandez, J.; Tapia-Santos, B.; Leite, M.C.A. A model for coupling within-host and between-host dynamics in an infectious disease. Nonlinear Dyn. 2012, 68, 401–411. [Google Scholar] [CrossRef]
  12. Angulo, O.; Milner, F.; Sega, L. A SIR epidemic model structured by immunological variables. J. Biol. Syst. 2013, 21, 1340013. [Google Scholar] [CrossRef]
  13. Martcheva, M.; Tuncer, N.; St Mary, C. Coupling within-host and between-host infectious diseases models. Biomath 2015, 4, 1510091. [Google Scholar] [CrossRef]
  14. Gandolfi, A.; Pugliese, A.; Sinisgalli, C. Epidemic dynamics and host immune response: A nested approach. J. Math. Biol. 2015, 70, 399–435. [Google Scholar] [CrossRef] [PubMed]
  15. Metz, J.A.; Diekmann, O. The Dynamics of Physiologically Structured Populations; Springer: Berlin/Heidelberg, Germany, 1986; Volume 68. [Google Scholar]
  16. Thieme, H.R. Mathematics in population biology. In Mathematics in Population Biology; Princeton University Press: Princeton, NJ, USA, 2018. [Google Scholar]
  17. Brauer, F.; Shuai, Z.; Van Den Driessche, P. Dynamics of an age-of-infection cholera model. Math. Biosci. Eng. 2013, 10, 1335–1349. [Google Scholar]
  18. Hartley, D.M.; Morris, J.G., Jr.; Smith, D.L. Hyperinfectivity: A critical element in the ability of V. cholerae to cause epidemics? PLoS Med. 2005, 3, e7. [Google Scholar] [CrossRef]
  19. Mukandavire, Z.; Liao, S.; Wang, J.; Gaff, H.; Smith, D.L.; Morris, J.G. Estimating the reproductive numbers for the 2008–2009 cholera outbreaks in Zimbabwe. Proc. Natl. Acad. Sci. USA 2011, 108, 8767–8772. [Google Scholar] [CrossRef]
  20. Cai, L.M.; Modnak, C.; Wang, J. An age-structured model for cholera control with vaccination. Appl. Math. Comput. 2017, 299, 127–140. [Google Scholar] [CrossRef]
  21. Jiang, X. Global Dynamics for an Age-Structured Cholera Infection Model with General Infection Rates. Mathematics 2021, 9, 2993. [Google Scholar] [CrossRef]
  22. Wang, X.; Wang, J. Disease dynamics in a coupled cholera model linking within-host and between-host interactions. J. Biol. Dyn. 2017, 11, 238–262. [Google Scholar] [CrossRef]
  23. Ratchford, C.; Wang, J. Modeling cholera dynamics at multiple scales: Environmental evolution, between-host transmission, and within-host interaction. Math. Biosci. Eng. 2019, 16, 782–812. [Google Scholar] [CrossRef]
  24. Musundi, B. An Immuno-epidemiological Model Linking Between-host and Within-host Dynamics of Cholera. arXiv 2021, arXiv:2105.12675. [Google Scholar]
  25. Reidl, J.; Klose, K.E. Vibrio cholerae and cholera: Out of the water and into the host. FEMS Microbiol. Rev. 2002, 26, 125–139. [Google Scholar] [CrossRef] [PubMed]
  26. Nelson, E.J.; Harris, J.B.; Morris, J.G.; Calderwood, S.B.; Camilli, A. Cholera transmission: The host, pathogen and bacteriophage dynamic. Nat. Rev. Microbiol. 2009, 7, 693–702. [Google Scholar]
  27. Schmid-Hempel, P.; Frank, S.A. Pathogenesis, virulence, and infective dose. PLoS Pathog. 2007, 3, e147. [Google Scholar] [CrossRef]
  28. Alavi, S.; Mitchell, J.D.; Cho, J.Y.; Liu, R.; Macbeth, J.C.; Hsiao, A. Interpersonal gut microbiome variation drives susceptibility and resistance to cholera infection. Cell 2020, 181, 1533–1546. [Google Scholar] [CrossRef]
  29. Cho, J.Y.; Liu, R.; Macbeth, J.C.; Hsiao, A. The interface of Vibrio cholerae and the gut microbiome. Gut Microbes 2021, 13, 1937015. [Google Scholar]
  30. Gillman, A.N.; Mahmutovic, A.; Abel zur Wiesch, P.; Abel, S. The Infectious Dose Shapes Vibrio cholerae Within-Host Dynamics. Msystems 2021, 6, e00659-21. [Google Scholar] [CrossRef]
  31. Heijmans, H.J.A.M. An eigenvalue problem related to cell growth. J. Math. Anal. Appl. 1985, 111, 253–280. [Google Scholar] [CrossRef] [Green Version]
  32. Heijmans, H.J. The dynamical behaviour of the age-size-distribution of a cell population. In The Dynamics of Physiologically Structured Populations; Springer: Berlin/Heidelberg, Germany, 1986; pp. 185–202. [Google Scholar]
  33. Jauffret, M.D.; Gabriel, P. Eigenelements of a general aggregation-fragmentation model. Math. Model. Methods Appl. Sci. 2010, 20, 757–783. [Google Scholar] [CrossRef]
  34. Campillo, F.; Champagnat, N.; Fritsch, C. Links between deterministic and stochastic approaches for invasion in growth-fragmentation-death models. J. Math. Biol. 2016, 73, 1781–1821. [Google Scholar]
  35. Stadler, E. Eigensolutions and spectral analysis of a model for vertical gene transfer of plasmids. J. Math. Biol. 2019, 78, 1299–1330. [Google Scholar] [CrossRef]
  36. Stadler, E.; Müller, J. Analyzing plasmid segregation: Existence and stability of the eigensolution in a non-compact case. Discret. Contin. Dyn. Syst. B 2020, 25, 4127. [Google Scholar] [CrossRef]
  37. Schilling, K.A.; Cartwright, E.J.; Stamper, J.; Locke, M.; Esposito, D.H.; Balaban, V.; Mintz, E. Diarrheal illness among US residents providing medical services in Haiti during the cholera epidemic, 2010 to 2011. J. Travel Med. 2014, 21, 55–57. [Google Scholar] [CrossRef] [PubMed]
  38. Morris, J.G., Jr. Cholera—Modern pandemic disease of ancient lineage. Emerg. Infect. Dis. 2011, 17, 2099. [Google Scholar] [CrossRef]
  39. Gyllenberg, M.; Heijmans, H.J. An abstract delay-differential equation modelling size dependent cell growth and division. SIAM J. Math. Anal. 1987, 18, 74–88. [Google Scholar] [CrossRef]
  40. Yosida, K. Functional Analysis; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  41. Banasiak, J.; Miekisz, J. Multiscale Problems in the Life Sciences: From Microscopic to Macroscopic; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  42. Magal, P.; Ruan, S. Structured Population Models in Biology and Epidemiology; Springer: Berlin/Heidelberg, Germany, 2008; Volume 1936. [Google Scholar]
  43. Webb, G.F. Theory of Nonlinear Age-Dependent Population Dynamics; CRC Press: Boca Raton, FL, USA, 1985. [Google Scholar]
  44. Inaba, H. Age-Structured Population Dynamics in Demography and Epidemiology; Springer: Berlin/Heidelberg, Germany, 2017. [Google Scholar]
  45. Drnovšek, R. Spectral inequalities for compact integral operators on Banach function spaces. Math. Proc. Camb. Philos. Soc. 1992, 112, 589–598. [Google Scholar] [CrossRef]
  46. Mokhtar-Kharroubi, M. On spectral gaps of growth-fragmentation semigroups with mass loss or death. Commun. Pure Appl. Anal. 2022, 21, 1293. [Google Scholar] [CrossRef]
  47. Webb, G.F. An operator-theoretic formulation of asynchronous exponential growth. Trans. Am. Math. Soc. 1987, 303, 751–763. [Google Scholar] [CrossRef]
  48. Müller, J.; Kuttler, C. Methods and models in mathematical biology. In Lecture Notes on Mathematical Modelling in Life Sciences; Springer: Berlin, Germany, 2015. [Google Scholar]
  49. Kuehn, C. Multiple Time Scale Dynamics; Springer: Berlin/Heidelberg, Germany, 2015. [Google Scholar]
  50. Umoh, J.; Adesiyun, A.; Adekeye, J.; Nadarajah, M. Epidemiological features of an outbreak of gastroenteritis/cholera in Katsina, Northern Nigeria. Epidemiol. Infect. 1983, 91, 101–111. [Google Scholar] [CrossRef] [Green Version]
  51. Tuckwell, H.C. Stochastic Processes in the Neurosciences; Chapter IV; SIAM: Philadelphia, PA, USA, 1989. [Google Scholar]
  52. Lymperopoulos, I.N. #Stayhome to contain COVID-19: Neuro-SIR–Neurodynamical epidemic modeling of infection patterns in social networks. Expert Syst. Appl. 2021, 165, 113970. [Google Scholar]
Figure 1. Time course of the pathogen load within a single susceptible individual: booster events (uptake of pathogens) instantaneously increase the load, in-between two booster events the immune system leads to an exponential decrease in pathogens. The time between two booster events is given by i.i.d. random variable T. ψ B = 1 is taken constant here, the initial time span is dismissed as a burn-in phase s.t. time runs here from 2970 to 3010.
Figure 1. Time course of the pathogen load within a single susceptible individual: booster events (uptake of pathogens) instantaneously increase the load, in-between two booster events the immune system leads to an exponential decrease in pathogens. The time between two booster events is given by i.i.d. random variable T. ψ B = 1 is taken constant here, the initial time span is dismissed as a burn-in phase s.t. time runs here from 2970 to 3010.
Mathematics 10 03114 g001
Figure 2. Shape of Ω : The upper and lower bound of Ω are given by characteristic curves, originating in ( 0 , P * ) and ( 0 , P * ) . The points P * resp. P * are constructed in such a way that the characteristics hit ( τ ˇ , P * ψ B ) resp. ( τ ^ , P * ψ B ) .
Figure 2. Shape of Ω : The upper and lower bound of Ω are given by characteristic curves, originating in ( 0 , P * ) and ( 0 , P * ) . The points P * resp. P * are constructed in such a way that the characteristics hit ( τ ˇ , P * ψ B ) resp. ( τ ^ , P * ψ B ) .
Mathematics 10 03114 g002
Figure 3. Black: equilibrium solution Σ ( P ) = S ( τ , P ) d τ (T is assumed to be uniformly distributed between τ ˇ = 2 , τ ^ = 10 , ψ B = 1 , γ = 0.2 ). The support of that function is in the interval [ P * , P * ] , as expected. Gray: β ( P ) = β 0 max { 0 , 1 π 0 / P } , where we use π 0 = 3 and (for this figure) β 0 = 1 .
Figure 3. Black: equilibrium solution Σ ( P ) = S ( τ , P ) d τ (T is assumed to be uniformly distributed between τ ˇ = 2 , τ ^ = 10 , ψ B = 1 , γ = 0.2 ). The support of that function is in the interval [ P * , P * ] , as expected. Gray: β ( P ) = β 0 max { 0 , 1 π 0 / P } , where we use π 0 = 3 and (for this figure) β 0 = 1 .
Mathematics 10 03114 g003
Figure 4. Left: I = F ( a , I ) , right: three simulations differing in the choice of a. Black: I ( t ) , gray: B ( t ) / 20 (B is scaled to be comparable with the prevalence). Parameter values are given in the Appendix A.
Figure 4. Left: I = F ( a , I ) , right: three simulations differing in the choice of a. Black: I ( t ) , gray: B ( t ) / 20 (B is scaled to be comparable with the prevalence). Parameter values are given in the Appendix A.
Mathematics 10 03114 g004
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Musundi, B.; Müller, J.; Feng, Z. A Multi-Scale Model for Cholera Outbreaks. Mathematics 2022, 10, 3114. https://doi.org/10.3390/math10173114

AMA Style

Musundi B, Müller J, Feng Z. A Multi-Scale Model for Cholera Outbreaks. Mathematics. 2022; 10(17):3114. https://doi.org/10.3390/math10173114

Chicago/Turabian Style

Musundi, Beryl, Johannes Müller, and Zhilan Feng. 2022. "A Multi-Scale Model for Cholera Outbreaks" Mathematics 10, no. 17: 3114. https://doi.org/10.3390/math10173114

APA Style

Musundi, B., Müller, J., & Feng, Z. (2022). A Multi-Scale Model for Cholera Outbreaks. Mathematics, 10(17), 3114. https://doi.org/10.3390/math10173114

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop