Abstract
Understanding both the epidemiological and evolutionary dynamics of antimicrobial resistance is a major public health concern. In this paper, we propose a nested model, explicitly linking the within- and between-host scales, in which the level of resistance of the bacterial population is viewed as a continuous quantitative trait. The within-host dynamics is based on integro-differential equations structured by the resistance level, while the between-host scale is additionally structured by the time since infection. This model simultaneously captures the dynamics of the bacteria population, the evolutionary transient dynamics which lead to the emergence of resistance, and the epidemic dynamics of the host population. Moreover, we precisely analyze the model proposed by particularly performing the uniform persistence and global asymptotic results. Finally, we discuss the impact of the treatment rate of the host population in controlling both the epidemic outbreak and the average level of resistance, either if the within-host scale therapy is a success or failure. We also explore how transitions between infected populations (treated and untreated) can impact the average level of resistance, particularly in a scenario where the treatment is successful at the within-host scale.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
Antimicrobial resistance (AMR) is one of the major challenges we face in the modern area (Larsson and Flach 2022). An antimicrobial substance is a chemical agent interacting with the physiology of a bacterial cell. The antimicrobial activity on a given bacterium’s (\({\mathcal {S}}\)) is an increasing function of its concentration in the medium (\({\mathcal {C}}\)), such that \({\mathcal {S}}\left( 0\right) =0\) and \({\mathcal {S}}\left( {\mathcal {C}}\right) \rightarrow {\mathcal {S}}_{\mathrm{{sat}}}\) as \({\mathcal {C}}\rightarrow {\mathcal {C}}_{\mathrm{{sat}}}\), where \({\mathcal {S}}_{\mathrm{{sat}}}\) and \({\mathcal {C}}_{\mathrm{{sat}}}\) are saturating constants. This intuitive approach implies that there exists \({\mathcal {C}}^{\star }\) in \(\left( 0,{\mathcal {C}}_{\mathrm{{sat}}}\right) \) such that \({\mathcal {S}}\left( {\mathcal {C}}^{\star }\right) \) is equal to the intrinsic rate of increase and reverses the growth of a bacterial population. Such a threshold concentration at which a bacterial population does not grow (at least in in vitro) is called the Minimum Inhibitory Concentration (MIC). The level of resistance to a given antimicrobial is then a continuous trait by nature referred to as antimicrobial quantitative resistance (qAMR), at least at the population level, and qAMR is key to better understanding the evolutionary dynamics of AMR (Djidjou-Demasse et al. 2023). Here, we introduce a quantitative descriptor \(x \in {\mathbb {R}}\)—a label of the bacterial strain with resistance level x– describing the level of resistance. Most of the modelling approaches devoted to AMR tackling the case of qualitative (or “binary”) resistance are generally based on the dynamical interaction between two parasite strains resulting in a discrete and finite formulation of MICs (Blanquart 2019). This analysis ignores the evolutionary short-term transient dynamics which lead to the emergence of resistance (e.g., Lipsitch and Levin 1997; Kepler and Perelson 1998; Day and Read 2016; Djidjou-Demasse et al. 2021; Tazzyman and Bonhoeffer 2014; Millan et al. 2014; D’Agata et al. 2008).
Here, we proposed a nested (or embedded) model explicitly linking the within- and between-host evolutionary dynamics. Such a nested structure is particularly important because, over the past few decades, it is clear that ecological and evolutionary dynamics are influenced by processes operating across scales (Elderd et al. 2022). Very few studies considered the continuous nature of AMR in the context of this work (e.g., Djidjou-Demasse et al. 2023), and few studies have implemented a nested model in this context so far (e.g., Beardmore et al. 2017; Shen et al. 2019). The bacterial population is assumed to be phenotypically (and genetically) diverse through the level of antimicrobial resistance x. This quantitative trait affects different components of the bacterial population life cycle, such as growth and death rates. In addition to those effects on the death and birth rates, bacterial population resistance level also mitigates the antimicrobial efficiency with respect to that population. From a theoretical point of view, properties of the within-host model proposed here are based on previous analytical quantitative genetics results developed in Djidjou-Demasse et al. (2017) and Burie et al. (2020).
An integro-differential equation is used to model the within-host dynamics of the bacterial population. Such a within model formulation is previously proposed in Djidjou-Demasse et al. (2023). Each host individual is classified as either a treated host, labeled as T, or an untreated host, labeled as U. The model tracks the dynamics of a bacterial population within a treated host (\(b^T\)) or an untreated host (\(b^U\)). At time \(\tau \), the bacteria density with resistance level \(y\in {\mathbb {R}}\) within a treated and untreated host is quantified by \(b_i^T(\tau ,y)\) and \(b_i^U(\tau ,y)\) respectively. The subscript “i” (with \(i\in {\mathcal {I}}=\{1,2,\cdots ,n\}\)) represents an individual immune system and then allows taking into account the immune system heterogeneity in the host population. A bacteria with resistance level y generate offspring with resistance level x at a per-capita rate \(J(x-y) p(y) b_i^{\vartheta }(\tau ,y)\), where p(y) is the bacterial intrinsic growth rate, and \(J(x-y)\) is the probability for a bacterial population with resistance level y to mutate towards a level x during the reproduction process. Therefore, the total number of bacteria produced at time \(\tau \) with a resistance level x is quantified by \( \left( 1+\int _{{\mathbb {R}}} b_i^{\vartheta }(\tau ,x)dx\right) ^{-\kappa } \int _{{\mathbb {R}}}J(x-y)p(y) b_i^{\vartheta }(\tau ,y) \textrm{d}y\), where \(\kappa \) is a positive parameter. The parameter \(\kappa >0\) is introduced to impose the bacterial population homeostasis. The within-host model reads as
where \(B_i^{\vartheta }(\tau )= \int _{{\mathbb {R}}} b_i^{\vartheta }(\tau ,x)\textrm{d}x,\) is the total bacteria load. The term \(\xi _i^{\vartheta }\) accounts for the individual clearance of bacterial cells with resistance level x, either by the immune system (\(\mu _i\)) or by the efficiency of antimicrobial pressure (k). Thus, \(\xi _i^U(x)= \mu _i(x)\), for untreated host, and \(\xi _i^T(x)=\mu _i(x)+k(x)\), for treated host. Here, it is assumed that bacteria are subject to a biocidal antimicrobial pressure, .ie. killing and not diminishing the birth rate of bacteria. Note that, the within-host model (1.1) allows to follow evolutionary parameters such as the average level of resistance for treated (\(\bar{x}_i^T(\tau )\)) and untreated (\(\bar{x}_i^U(\tau )\)) individuals \(\tau \)-time post infection given by
At the between-host scale, the host population is subdivided into three states. At any time t, an individual—with the immune system’s response level \(i \in {\mathcal {I}}\)—can be susceptible to the infection \(S_i(t)\) or infected \(I_i^\vartheta (t,\tau ,\bar{x}_i^{\vartheta }(\tau ))\), \(\vartheta =\{T,U\}\). The variables \(\tau \) and \(\bar{x}_i^{\vartheta }(\tau )\) respectively represent the time post-infection and the average resistance level of the infected host. It is important to clearly understand the meaning of infected individuals \(I^\vartheta (t,\tau ,{\bar{x}}_i^\vartheta (\tau ))\). Indeed, each infected individuals is potentially infected with multiple bacteria strains with variable frequencies and resistance levels. Therefore, \(\bar{x}_i^\vartheta (\tau )\) represent the individual resistance level quantified by the within-host dynamics through the above formula. However, for simplicity, and without loss of generality, we will note \(I_i^\vartheta (t,\tau ,{\bar{x}}_i^\vartheta (\tau ))\equiv I_i^\vartheta (t,\tau )\) for \(\vartheta \in \{T,U\}\) and \(i\in {\mathcal {I}}\).
Individual transmission and loss rates at the between-host scale, \(\tau \)-time since infection, \(\beta _i^{\vartheta }(\tau )\) and \(\alpha _i^{\vartheta }(\tau )\), are linked to the within-host dynamics at time \(\tau \). As an example, these parameters can be represented as Holling functions of type II (or similarly the Beddington-DeAngelis functional response) such that, for all \(\vartheta \in \{T,U\}\),
where \(\beta _0\) and \(\alpha _0\) are scaling constants, and \(r_0\) is the half-saturation constant for the total bacterial load \(B_i^{\vartheta }\). Note that, for the loss rate of infected individuals \(\alpha _i^{\vartheta }\), the term \(\alpha _0B_i^{\vartheta }/(r_0+ B_i^{\vartheta })\) represents the loss due to the disease induced mortality while \(\gamma ^{\vartheta }\) is the loss due to recovery. We can assume that the function \(\gamma ^{\vartheta }\) is of the form
where \(B_{min}\) is the threshold below which the infection becomes undetectable such that the infected individuals is considered as recovered.
The force of infection induced by infected individuals at time t is then given by
The nested model proposed here then makes it possible to simultaneously track the epidemiological dynamics of the host population as well as evolutionary quantities such as the average level of resistance at both the within- and between-host scales. Such an approach is original and to our knowledge, no study has considered nested models for the evolutionary dynamics of AMR, viewed as a continuous quantitative trait.
The dynamics of newly infected individuals (i.e. \(\tau =0\)) in each group (treated or untreated) is thus defined by (for \(\vartheta \in \{T,U\}\))
where \(q_i^T\in (0,1)\) is the treatment rate in the host population and \(q_i^U=1-q_i^T\). During their infection, treated individuals can stop the treatment at rate \(\omega _U^T(\tau )\), and untreated infections can join the treated group at rate \(\omega _T^U(\tau )\). The loss rate of infected individuals \(\tau \)-time post infection occurs at rate \(\alpha _i^{\vartheta }(\tau )\). Susceptible individuals are recruited at a constant rate \(\Lambda _i\) and the natural death rate of the host population is \(\mu _h\). The between-host model then reads
Finally, the nested within-host (1.1) and between-host model (1.4)–(1.5) are summarised by Fig. 1. The main variables and parameters are listed in Table 1. We emphasize that the dynamical properties (that we will recall later) of within-host model (1.1) are precisely analyzed in Djidjou-Demasse et al. (2023). Therefore, our main objective here is devoted to the analysis of the nested model (1.1)–(1.5).
The rest of this work is organized as follows. In Sect. 2, we state the main results of the nested model that are obtained in this work. These include the existence of the globally defined non-negative semiflow and the existence of the unique positive equilibrium for the within-host model, and global threshold analysis results for the between-host model. The model’s typical dynamics are provided in Sect. 3. This includes the within- and between-host models parameterization and the characterization of the evolutionary parameters such as the average levels of resistance. In Sect. 4, we delve into the effects of various parameters on the equilibrium structure of the host population, along with addressing the parameterization issue within nested models. Section 5 focuses on providing preliminary results. Specifically, it addresses the existence and uniqueness of solutions, derivation of the basic reproduction number, and the existence of a unique endemic equilibrium for System (1.4)–(1.5). Finally, Sect. 6 is devoted to the proof of the global asymptotic results.
2 Main results
This section is devoted to the main results of the nested model (1.1)–(1.5). Such results include the existence of the unique maximal bounded semiflow, and a precise description of the unique positive equilibrium of Model (1.4)–(1.5). By providing global stability results, we will also conduct a precise threshold analysis of the between-host model (1.4)–(1.5).
First of all, for biological feasibility of the nested model (1.1)–(1.5), we make use of the following assumptions. More precisely, the within-host model (1.1) is formulated based on the following assumption
Assumption 2.1
-
1.
Functions \(\mu _i\), k, \(\xi _i^{\vartheta }\), and p are always positive on \({\mathbb {R}}\), with \(\vartheta \in \{T,U\} \). Furthermore, p is a bounded function on \({\mathbb {R}}\) and \(\kappa >0\). Finally, the function \(\frac{p}{\xi _i^{\vartheta }}\) is continuous on \({\mathbb {R}}\) and satisfies \(\frac{p}{\xi _i^{\vartheta }} > 0\) and \(\lim \nolimits _{|x|\rightarrow \infty }\frac{ p}{\xi _i^{\vartheta }}(x)=0\).
-
2.
The mutation kernel J is bounded and integrable on \({\mathbb {R}}^+\), positive almost everywhere, and satisfies \(\int _{{\mathbb {R}}^+}J(x)dx>0\), \(J(-x)=J(x)\), for all x.
-
3.
The mutation kernel J decays rather rapidly towards infinity in the sense that \(J(x)=O\left( \frac{1}{|x|^{\infty }}\right) \) as \(|x|\rightarrow \infty \). In other words, \(\lim \nolimits _{|x|\rightarrow \infty }|x|^nJ(x)=0\), for all \(n\in {\mathbb {N}}\).
Furthermore, the between-host model’s parameters satisfy the following assumption.
Assumption 2.2
-
1.
Recruitment rate \(\Lambda _i\) (\(i\in {\mathcal {I}}\)) and natural death rate \(\mu _h\) are positive constants.
-
2.
The treatment rates \(q_i^{\vartheta }\) (\(i\in {\mathcal {I}}\), \(\vartheta \in \{T,U\}\)) are positive constants.
-
3.
The rates \(\omega _U^T\), \(\omega ^U_T\) belongs to \(L^{\infty }({\mathbb {R}}_+)\), with respective essential upper bounds \({\overline{\omega }}_U^T\), \({\overline{\omega }}^U_T\) and positive essential lower bounds \(\underline{\omega }_U^T\), \(\underline{\omega }^U_T\).
-
4.
Parameters \(\beta _i^\vartheta \) and \(\alpha _i^\vartheta \) (\(i\in {\mathcal {I}}\), \(\vartheta \in \{T,U\}\)) are such that \(\beta _i^\vartheta , \alpha _i^\vartheta \in L^\infty ({\mathbb {R}}_+)\).
-
5.
The transmission rates \(\beta _i^\vartheta (\cdot )\),s are Lipschitz continuous almost everywhere on \({\mathbb {R}}_+\).
2.1 Summary key findings on the within-host dynamics
The dynamical properties of the within-host model (1.1) have been precisely investigated in Djidjou-Demasse et al. (2023). The first result of Model (1.1) is about the existence of the unique maximal bounded semiflow. Such a result reads as,
Theorem 2.3
Let Assumption 2.1 be satisfied. Let \(b_{i0}^{\vartheta }\in L^1_+\). Then,
-
1.
There exists a unique global solution \(v(\cdot ,b_{i0}^{\vartheta }):[0,\infty )\rightarrow L_+^1({\mathbb {R}})\) of (1.1) with \(v(0,b_{i0}^{\vartheta })=b_{i0}^{\vartheta }\) and \(v(\tau ,b_{i0}^{\vartheta })=b_i^{\vartheta }(\tau ,\cdot )\) for all \(\tau >0\).
-
2.
The semi-flow defined by \(\{v(\tau ,b_{i0}^{\vartheta })\}_{\tau }\) is bounded dissipative and asymptotically smooth, and hence, its admits a global attractor in \(L_+^1({\mathbb {R}})\).
-
3.
The semi-flow \(\{v(\tau ,b_{i0}^{\vartheta })\}_{\tau }\) is such that for any \(b_{i0}^{\vartheta }\in L^1_+({\mathbb {R}})\setminus \{0\}\), \(b_i^{\vartheta }(\tau ,x)>0\), for all \(\tau >0\), \(x\in {\mathbb {R}}\).
The basic reproduction number \({\mathcal {N}}_{i0}^{\vartheta }\)—defined as the expected number of bacteria arising from one bacterium in a bacteria-free environment—of the bacteria population with resistance level x, within an individual with immune system level i, is calculated as
Next, a non-trivial equilibrium of Model (1.1) is strongly related to the principal eigenpair of the below linear integral operator \(H^{\vartheta }_{i}\) defined on \(L^p({\mathbb {R}})\) (for any \(p\ge 1\)), by
We then have the following result.
Theorem 2.4
Let \(r(H^{\vartheta }_{i})\), the spectral radius of the operator \(H^{\vartheta }_{i}\), and \(\phi >0\) the associated eigenfunction normalized such that \(||\phi ||_{L^1}=1\).
-
1.
When \(r(H^{\vartheta }_{i})\le 1\), the bacteria-free equilibrium \(F_{i0}^{\vartheta }\) is the unique equilibrium of Model (1.1).
-
2.
When \(r(H^{\vartheta }_{i})> 1\), in addition to \(F_{i0}^{\vartheta }\), Model (1.1) has a unique equilibrium \({\overline{F}}_{i}^{\vartheta }>0\) such that
$$\begin{aligned} {\overline{F}}_{i}^{\vartheta }(x)= \left( \frac{\left( r(H^{\vartheta }_{i})\right) ^{\frac{1}{\kappa }}-1}{\int _{{\mathbb {R}}}\frac{\phi }{\sqrt{p\; \xi _i^{\vartheta }}}\textrm{d}y} \right) \frac{\phi (x)}{\sqrt{p(x)\xi _i^{\vartheta }(x)}}. \end{aligned}$$(2.3)Furthermore, the semi-flow \(\{v(\tau ,b_{i0}^{\vartheta })\}_{\tau }\) is uniformly persistent, that is, there exists a constant \(\eta \) such that for any \(b_{i0}^{\vartheta }\in L^1_+({\mathbb {R}}){\setminus }\{0\}\), the unique solution \(v(\tau ,b_{i0}^{\vartheta })=b_i^{\vartheta }(\tau ,\cdot )\) of Model (1.1) with initial data \(b_{i0}^{\vartheta }\) satisfies \(\lim \limits _{\tau \rightarrow \infty }\inf \Vert b_i^{\vartheta }(\tau ,\cdot )\Vert _{L^1}>\eta .\)
-
3.
The bacteria-free equilibrium \(F^{\vartheta }_{i0}\) of Model (1.1) is asymptotically stable if \(r(H^{\vartheta }_{i})<1\) and unstable if \(r(H^{\vartheta }_{i})>1\).
-
4.
When \(r(H^{\vartheta }_{i})<1\), the bacteria-free equilibrium \(F^{\vartheta }_{i0}\) is globally asymptotically stable in \(L^1_+({\mathbb {R}})\), that is, for any solution \(b_i^{\vartheta }(\tau ,\cdot )\) with initial \(b_{i0}^{\vartheta }\in L^1_+({\mathbb {R}}){\setminus }\{0\}\), we have \(b_i^{\vartheta }(\tau ,\cdot )\rightarrow 0\) in \(L^1_+({\mathbb {R}})\), as \(\tau \rightarrow \infty \).
We recall that the within-host model (1.1) is precisely analyzed in Djidjou-Demasse et al. (2023). We then refer to Appendices F–I in Djidjou-Demasse et al. (2023) for the detailed proof of Theorems 2.3 and 2.4. Note that the linear operator \(H^{\vartheta }_{i}\) naturally emerges when characterizing the positive equilibrium of the within-host model (1.1) (Djidjou-Demasse et al. 2023).
Furthermore, the estimate (2.3) gives that the endemic equilibrium \({\overline{F}}_{i}^{\vartheta }\) of the within-host model (1.1) basically relied to the principal eigenfunction of the linear operator \(H^{\vartheta }_{i}\) for any given probability kernel J satisfying Assumption 2.1. However, the profile of the endemic equilibrium \({\overline{F}}_{i}^{\vartheta }\) with respect to \(x\in {\mathbb {R}}\) can be precisely described when the mutation kernel J depends on a small positive parameter (let say \(\varepsilon<\!\!<1\)) with the scaling form
where \(\varepsilon >0\) represents the mutation variance in the phenotypic space. More precisely, when \(\varepsilon >0\) is small, then the endemic equilibrium \({\overline{F}}_{i}^{\vartheta }\) concentrates on the set \({\mathcal {S}}_i^{\vartheta }\) defined by
The set \({\mathcal {S}}_i^{\vartheta }\) is referred to as the set of Evolutionary Attractors (or dominant strains) of the within-host model in the classical adaptive dynamics theory ( e.g., Geritz et al. 1997; Metz et al. 1996). Furthermore, when the function \({\mathcal {N}}_{i0}^{\vartheta }\) is at least of class \({\mathcal {C}}^1\), with a finite number of maximum, it is shown in Djidjou-Demasse et al. (2017) that these dominant strains coincide with the set \({\mathcal {S}}_i^{\vartheta }\). Denoting by \(H^{\vartheta }_{i,\varepsilon }\), the operator \(H^{\vartheta }_i\)—by replacing the kernel J by \(J_\varepsilon \)—by results in Djidjou-Demasse et al. (2017) (Theorem 2.2), the spectral radius \(r(H^{\vartheta }_{i,\varepsilon })\) of \(H^{\vartheta }_{i,\varepsilon }\) satisfied, for \(\varepsilon \) sufficiently small
By the above estimate, \(\textrm{sign}\left[ r\left( H^{\vartheta }_{i,\varepsilon }\right) -1 \right] = \textrm{sign}\left[ {\mathcal {N}}_{i0}^{\vartheta }(x^*)-1 \right] \), for all \(x^*\in \) and \(\varepsilon \) sufficiently small. Furthermore, if \(\varepsilon \ll 1\), \({\mathcal {S}}_i^{\vartheta }=\{x_i^*\}\) and \(N_{i0}^{\vartheta }(x_i^*)>1\), then the unique positive stationary state \({\overline{F}}_{i}^{\vartheta }\equiv {\overline{F}}_{i,\varepsilon }^{\vartheta } \), given by (2.3), of the within-host model (1.1) is concentrated around the evolutionary attractor \(x_i^*\) in the space of resistance level \({\mathbb {R}}\). In other words, \(x_i^*\) is the average bacterial resistance level at the within-host scale equilibrium and we have \(\lim _{\varepsilon \rightarrow 0}\int _{{\mathbb {R}}}u(x) {\overline{F}}_{i,\varepsilon }^{\vartheta }(x)\textrm{x}=u\left( x_i^*\right) \) for any continuous function \(u \in {\mathcal {C}}\left( {\mathbb {R}}\right) \). We refer to Theorem 2.3 in Djidjou-Demasse et al. (2017) for such a concentration phenomenon.
2.2 Key findings of the nested within- and between-host dynamics
At the between-host scale, by setting \({{\textbf {S}}}(t)=(S_i(t))_{i\in {\mathcal {I}}}\), \({{\textbf {I}}}_i(t,\tau )=(I_i^T(t,\tau ),I_i^U(t,\tau ))\), \({{\textbf {I}}}(t,\tau )=({{\textbf {I}}}_i(t,\tau ))_{i\in {\mathcal {I}}}\), \(\varvec{\alpha }_i(\tau )= \text {diag}(\alpha _i^T(\tau ),\alpha _i^U(\tau ))\), \(\varvec{\beta }_i(\tau )=(\beta _i^T(\tau ), \beta _i^U(\tau ))\), \({{\textbf {q}}}_i=(q_i^T,q_i^U)\), \(\varvec{\omega }(\tau )=\left( \begin{array}{cc} 0 &{} \omega _T^U(\tau ) \\ \omega _U^T(\tau ) &{} 0 \\ \end{array}\right) \), and \({\varvec{e}}= \begin{pmatrix} 0 &{} 1\\ 1 &{} 0 \end{pmatrix} \), System (1.4)–(1.5) rewrites into the following compact form,
where \(\lambda (t)=\sum _i\int _0^{\infty }\left\langle \varvec{\beta }_i(\tau ),{{\textbf {I}}}_i(t,\tau )\right\rangle \textrm{d}\tau \), \(\varvec{\Lambda }=(\Lambda _i)_{i\in {\mathcal {I}}}\), \({{\textbf {q}}}=({{\textbf {q}}}_i)_{i\in {\mathcal {I}}}\), \(\varvec{\Phi }(\tau )=(\varvec{\Phi }_i(\tau ))_{i\in {\mathcal {I}}}\), with \(\varvec{\Phi }_i(\tau )= {\varvec{e}}\varvec{\omega }(\tau ) +\varvec{\alpha }_i(\tau )+\mu _h\).
Using the next-generation operator approach (e.g., Diekmann et al. 1990; Inaba 2012), the basic reproduction number \({\mathcal {R}}_0^i\) of the whole infected individuals of group i, is given by
where
and where \(\varvec{\Pi }_i(\tau _2,\tau _1)\), \(0\le \tau _1\le \tau _2<\infty \), is the evolutionary system generated by the linear operator \(\left[ -\varvec{\Phi }_i(\tau )+ \varvec{\omega }(\tau ) \right] \); see Remark 2.5 for some details on \(\varvec{\Pi }_i\). Moreover, the basic reproduction number \({\mathcal {R}}_0\) at the whole between-host scale is such that
We refer to Sect. 5.2 for details of the computation of \({\mathcal {R}}_0^i\),s and \({\mathcal {R}}_0\).
Note that the parameter \(\chi _k\) quantifies the overall infectiousness of the whole infected individuals of group \(k\in {\mathcal {I}}\). A more explicit expression of the infectiousness \(\chi _k\) is difficult to obtain in general. However, one can go further steps in some particular configurations of the treatment status transition rates \(\varvec{\omega }(\tau ).\) Indeed, assume that we can find \(\tau _0>0\) and \(\tau _1>0\) such that
In the above scenario, the regimen \((0,\tau _0)\)-post infection may corresponds to the initial phase where each infections, either treated or untreated, remain to their initial treatment status. The second regimen \((\tau _0,\tau _0+\tau _1)\)-post infection may corresponds to the phase during which previously untreated infections becomes treated while treated infections remain to their initial status. In such a configuration, we have (see Sect. 5.2 for details)
with \(c_k=\mu _h + \inf _\tau \alpha _k^U(\tau ) + \inf _\tau \alpha _k^T(\tau )\),
and
Note that parameters \(\Gamma _{0}^{k,\vartheta }\),s and \(\Gamma _{1}^{k,\vartheta }\),s are survival probabilities during phases \((0,\tau _0)\) and \((\tau _0,\tau _0+\tau _1)\)-post infection of infected individuals of group \(k\in {\mathcal {I}}\), treated (\(\vartheta =T\)) or untreated (\(\vartheta =U\)).
Remark 2.5
Let \(\varvec{\Pi }_i(\tau _2,\tau _1)\), \(0\le \tau _1\le \tau _2<\infty \), the evolutionary system generated by the linear operator \({\varvec{A}}_i(\tau ):= -\varvec{\Phi }_i(\tau )+ \varvec{\omega }(\tau ) \). It means that \(\varvec{\Pi }_i\) is generated from the following evolutionary system
If, for example, the linear operator \({\varvec{A}}_i\) is diagonal, we have
In such a configuration we explicitly have \(\varvec{\Pi }_i(\tau _2,\tau _1)= e^{\int _{\tau _1}^{\tau _2} {\varvec{A}}_i(\eta ) \textrm{d}\eta }.\) However, obtaining an explicit expression for \(\varvec{\Pi }_i\) may not always be straightforward or possible in general. A naive approach would be to solve problem (2.8) as above, but it is well known that such an exponent formula does not give a solution to the problem at hand.
In addition to the disease-free equilibrium—the DFE—\({{\textbf {E}}}^0=({{\textbf {S}}}^0,{{\textbf {0}}}_{L^1((0,\infty ),{\mathbb {R}}^{2n})})\), with \({{\textbf {S}}}^0=\left( \Lambda _i/\mu _h \right) _{i\in {\mathcal {I}}}\), which is always an equilibrium of Model (1.4)–(1.5), this model also exhibits an endemic equilibrium given by the following result
Theorem 2.6
Let Assumptions 2.1 and 2.2 hold. If \({\mathcal {R}}_0>1\), then system (1.4)–(1.5) has a unique endemic equilibrium \({{\textbf {E}}}^*=({{\textbf {S}}}^*,{{\textbf {I}}}^*(\tau ))\), such that \(\forall i\in {\mathcal {I}}\),
where \(\lambda ^*=\mu _h({\mathcal {R}}_0-1)\).
Therefore, the threshold dynamics of Model (1.4)–(1.5) is summarized as follows
Theorem 2.7
Let Assumptions 2.1 and 2.2 hold. Then,
-
(i)
If \({\mathcal {R}}_0\le 1\) or \(\sum _{i\in {\mathcal {I}}}\int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ),{{\textbf {I}}}_{i0}(\tau ) \right\rangle \textrm{d}\tau =0\), then the disease-free equilibrium \({{\textbf {E}}}^0=({{\textbf {S}}}^0,\,{{\textbf {0}}}_{L^1((0,\infty ),{\mathbb {R}}^{2n})})^t\) of system (1.4)-(1.5) is globally asymptotically stable in the sens that
$$\begin{aligned} \lim _{t\rightarrow \infty } \left( S_i(t), I_i^T(t,\cdot ), I_i^U(t,\cdot ) \right) _{i \in {\mathcal {I}}} = {{\textbf {E}}}^0, \end{aligned}$$where the above convergence holds for the topology of \({\mathbb {R}}^n\times L^1((0,\infty ), {\mathbb {R}}^{2n})\).
-
(ii)
If \({\mathcal {R}}_0>1\) and \(\sum _{i\in {\mathcal {I}}}\int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ),{{\textbf {I}}}_{i0}(\tau ) \right\rangle \textrm{d}\tau >0\), then the endemic equilibrium \({{\textbf {E}}}^*\) of system (1.4)-(1.5) is globally asymptotically stable, that is,
$$\begin{aligned} \lim _{t\rightarrow \infty } \left( S_i(t), I_i^T(t,\cdot ), I_i^U(t,\cdot ) \right) _{i \in {\mathcal {I}}} = {{\textbf {E}}}^*, \end{aligned}$$for the topology of \({\mathbb {R}}^n\times L^1((0,\infty ), {\mathbb {R}}^{2n})\).
3 Numerical illustrations
Here, we present a series of numerical simulations employing semi-explicit finite difference numerical schemes. We refer to Djidjou-Demasse (2021) for an example of a code repository within the context of the model proposed here. We illustrate an example of typical dynamics that can be simulated by the nested model (1.1)–(1.5). The model simultaneously captures the outbreak dynamics as well as the evolutionary dynamics of the average resistance level within the host population. The within-host model parameters are basically the same as in Djidjou-Demasse et al. (2023). Intuitively there exist two threshold levels, assumed here 0 and 1 (called reference “sensitive” and “resistant” strains) such that, a strain with resistance level x can be classically referred to as “sensitive”, “intermediate”, or “resistant” depending on whether \(x < 0\), \(0< x < 1\), or \(x > 1\). For sake of simplicity, we assume that the host population is homogeneous in terms of immune system level, i.e., \(\text {card}({\mathcal {I}})=1\). For all illustrative scenarios, we will have \({\mathcal {R}}_0>1\) such that the disease is persistent at the between-host scale (Theorem 2.7). The probability density function at the within-host scale (\(J\equiv J_\varepsilon \)) is assumed of type (2.4). Specifically, we define \(J_\varepsilon \) as a Gaussian distribution \(J_\varepsilon (x)=\frac{1}{\varepsilon \sqrt{2\pi }}e^{-\frac{1}{2}\left( \frac{x}{\varepsilon }\right) ^2}\), where \(\varepsilon >0\) represents a small parameter that signifies the mutation variance within the phenotypic space.
Within-host parameterization The antimicrobial killing rate function \(k(\cdot )\) is a decreasing function with respect to the resistance level x such that, \(k(x)=k_0\left( \frac{k_1}{k_0}\right) ^x,\) where \(k_0\) and \(k_1\) are the antimicrobial activity undergone by the reference sensitive and resistant strains. Moreover, knowing \(p_0\) and \(p_1\), respectively the intrinsic growth rate of reference strains 0 and 1, a suitable expression for function of p is \(p(x)= p_m \left[ 1+\left( \frac{p_m-p_0}{p_0}\right) \left( \frac{p_0}{p_1}\cdot \frac{p_m-p_1}{p_m-p_0} \right) ^x \right] ^{-1},\) where \(p_m\) is the upper bound of the intrinsic growth rate p and \(0<p_1<p_0<p_m\). The qualitative behaviour of functions k and p can be found in (Djidjou-Demasse et al. 2023, Fig. 2). We assume that the clearance rate of the bacteria cell due to the immune response, \(\mu (\cdot )\), is a constant function given by \(\mu (x)=\mu \). Furthermore, the average fitness cost-benefit ratio of resistance within a bacterial population can be expressed as \(c_b=\frac{\log (\Delta )}{\log (1+\delta )}\). Here, \(\Delta =\frac{(p_m-p_1)/p_1}{(p_m-p_0)/p_0}>1\) quantifies the relative cost of resistance, while \(\delta =\frac{k_0-k_1}{k_1}>0\) measures the fitness advantage of the reference resistant strain (see Djidjou-Demasse et al. 2023 for details).
Between-host parameterization Parameters \(\beta _i^{\vartheta }\) and \(\alpha _i^{\vartheta }\) are defined using Holling type functional responses introduce by (1.2). For all simulations, the threshold \(B_{min}\), introduced by (1.3), below which the infection becomes undetectable such that the infected individuals is considered as recovered is fixed as \(B_{min}=10^{-3}B_{0}\), with \(B_{0}=B(0)\) the initial total bacteria load. The total bacteria load (\(B_i^{\vartheta }(\tau )\)) and the recovery probability (\(1-\exp (-\int _0^\tau \gamma _{i}^\vartheta (s) \textrm{d}s\)), \(\tau \)-time post infection, are illustrated in Fig. 2.
Furthermore, an untreated infected individual joins the treated compartment when her total bacteria load is above a threshold \((1+\theta )\, B_{0}\), with \(\theta \ge 0\). Therefore, the influx rate from untreated to treated is assumed to be a function with respect to time \(\tau \) and is defined as follows
Similarly, we assume that an infected individual under treatment can drop down such a treatment when the bacteria load reach the same range as before the treatment. Therefore, the influx rate from treated to untreated is given by
Initial conditions and model outputs The initial bacterial population \(b_0^{\vartheta }(x)\) is assumed to be composed by a sensitive bacterial population with average resistance level \(x=0\). Hence, we set \(b_0^{\vartheta }(x)=m_0\times {\mathcal {N}}(0,\sigma _0,x),\) where \({\mathcal {N}}(0,\sigma _0,x)\) stands for the normalized density function of the Gaussian distribution at x with mean 0 and variance \(\sigma _0^2\). This means that the initial bacterial population is mostly composed of the reference “sensitive” strain. At the between-host scale, the initial condition of the epidemiological model is taken such that the susceptible population starts close to its disease-free equilibrium. More precisely, assuming an initial infection prevalence denoted as \(\mathrm{P_{rev}}=10\%\), we derive the initial susceptible population as \(S_{0}=(1-\mathrm{P_{rev}}) \frac{\Lambda }{\mu _h}\), along with the initial distribution of infectives which consists of \(I_{0}^T(\tau )=0\) and \(I_{0}^U(\tau )=\mathrm{P_{rev}} \frac{\Lambda }{\mu _h} \times L(\tau )\) for all \(\tau \ge 0\). Here, \(L(\tau )=10\ln (10)\times 10^{-10\tau }\), and it is important to note that L symbolizes the arbitrary initial distribution of individuals who have been infected since time \(\tau \). This distribution is scaled so that \(\int _{\mathbb {R}}L(\tau ) \textrm{d}\tau =1\).
The average level of resistance at within-host scale (\(\eta (t)\)) of the host population at time t is such that
where \(\bar{x}^{\vartheta }\),s are the individual average level of resistance and \(I(t)=\int _0^{\infty } \big (I^T(t,\tau ) +I^U(t,\tau ) \big ) \textrm{d}\tau .\)
Simulated scenarios Two simulated scenarios are considered, the first when the treatment is successful at the within-host level, and the second when the treatment failed at the within-host level. For all our simulated scenarios, the infection is assumed here to be always successful for untreated individuals, i.e., the immune system alone is no more enough to control the infection such that \(\max _{x\in {\mathbb {R}}} {\mathcal {N}}_{0}^{U}(x) >1\), leading to the bacterial persistence for untreated infections.
Our first scenario is for the case where the treatment is successful at the within-host level, i.e., the basic reproduction number of treated individuals \({\mathcal {N}}_{0}^{T}\) is such that \(\max _{x\in {\mathbb {R}}} {\mathcal {N}}_{0}^{T}(x) <1\) (Fig. 3C). In such a situation, the bacterial load is under control in the relatively short term for treated individuals (Fig. 3A), while it remains persistent for untreated individuals (Fig. 3B). At the between-host scale, the treatment rate have a strong effect on the epidemic outbreak (Fig. 3F–H). More precisely, increasing the treatment rate \(q^T\) in the host population strongly reduce the overall epidemic size (Fig. 3F–H), with \({\mathcal {R}}_0=4.7815\), 2.6750 and 0.5685, respectively. Furthermore, the average resistance level in the host population rapidly reach an equilibrium for which the level of resistance is moderately high compared to the initial resistance level of the host population (Fig. 3E).
In the second scenario, the treatment is assumed unsuccessful at the within-host level, i.e., the basic reproduction number of treated individuals \({\mathcal {N}}_{0}^{T}\) is such that \(\max _{x\in {\mathbb {R}}} {\mathcal {N}}_{0}^{T}(x) >1\) (Fig. 4C). In such a situation, the bacterial load remains persistent for both treated and untreated infections (Fig. 4A,B). Indeed, while we can observe an apparent decreasing of the bacteria load for some period of time for treated infections (Fig. 4A), at the end, we have the re-emergence of the bacteria population at within-host scale. Such a transient dynamics is explain by the fact that the initiation of treatment modifies the fitness landscape by shifting the maximum point of the within-host basic reproduction number \({\mathcal {N}}_{0}^{T}\) to the point \(x = x^* > 0\) (Fig. 4C). In contrast to the treatment success scenario (Fig. 3), increasing the treatment rate \(q^T\) in the host population have marginal effect in controlling the epidemic outbreak (Fig. 4F–H). In fact, with \(q^T=0.1\), 0.5 and 0.9, the outbreak remains persistent with \({\mathcal {R}}_0=10.9451\), 9.5855 and 8.2259, respectively. Significantly, it is worth noting that although the average resistance level in the host population continues to rise with the treatment rate \(q^T\) at equilibrium (Fig. 4E), there is a substantial increase in the range of resistance levels compared to the initial resistance level. This stands in contrast to the treatment success scenario (Fig. 3E), where the range of resistance levels remained relatively low. Additionally, in the treatment success scenario (Fig. 3E), the average resistance level in the host population quickly reaches equilibrium. However, in the treatment failure scenario (Fig. 4E), there is a comparatively longer transient period before the average resistance level in the host population reaches equilibrium. Overall, during the transient regimen, there is an initial subsequent increase in the average level of resistance to significantly higher levels, followed by a small decrease (Fig. 4E). This behavior is mostly attributed to the alteration of the fitness landscape caused by the treatment (Fig. 4C).
4 Discussion
Optimizing the treatment rate in the host population is key to controlling both the epidemic outbreak and the average level of resistance Increasing the treatment rate in the host population contributes to reducing the epidemic size at the between-host scale, although the effect is quite marginal in the treatment failure scenario at the within-host scale (Figs. 3F–H, 4F–H). This can be primarily attributed to the fact that, even in cases of treatment failure, the within-host infection remains controlled for a certain period (Fig. 4A). However, the subsequent phase is characterized by an increase in bacterial population density (Fig. 4A). Conversely, raising the treatment rate within the host population leads to an elevation in the average resistance level of that population, regardless of the treatment scenario at the within-host scale (Figs. 3E, 4E). However, this increase in the average resistance level is particularly pronounced in the case of treatment failure (Fig. 4E), in comparison to the case of treatment success (Fig. 3E). In the case of treatment success, the resistance level remains highly similar to that of the initial bacterial population. Consequently, depending on the treatment regimen implemented within the host population, it becomes essential and intriguing to determine an optimal treatment rate to effectively manage both the outbreak and the average level of antimicrobial resistance.
The delay for treating infections can impact the epidemic outbreak as well as average level of resistance Assume the scenario where the treatment is successful at the within-host scale (Fig. 3A) and at least 50% of infected individuals are under treatment (\(q^T\ge 0.5\)). Two configurations are introduced. In the first configuration, untreated individuals begin treatment at a rate \(\omega ^U_T\) [defined by (3.1)] when their total bacterial load reaches a threshold value of \(1.5\times B_0\), i.e., \(\theta =0.5\). In the second configuration, untreated individuals start treatment at a rate \(\omega ^U_T\) when their total bacterial load reaches a threshold value of \(2\times B_0\), i.e., \(\theta =1\). The case of \(\theta =0.5\) indicates a situation where the delay before initiating the treatment is very short. On the other hand, in the case of \(\theta =1\), the delay before starting the treatment is relatively more significant (Fig. 2A). In general, the average resistance level in the host population decreases as the delay before initiating treatment increases (Fig. 5A, B). One possible explanation is that untreated infected individuals do not significantly contribute to the increase in the average resistance level within the host population. Nevertheless, although early treatment effectively controls the epidemic outbreak for both moderate and high treatment rates (Fig. 5C, D)—with \({\mathcal {R}}_0=0.3922\) and 0.1120, respectively—delaying the treatment of infected individuals leads to an epidemic that is out of control, except in cases where the treatment rate is exceptionally high (Fig. 5E, F)—with \({\mathcal {R}}_0=2.6750\) and 0.5685, respectively.
Nested models parameterization issue Explicitly connecting the within- to and between-host scales are crucial to gain a more realistic picture aiming to integrate into the same modelling framework the epidemic dynamics and the evolutionary dynamics of antimicrobial resistance. However, such an approach requires making assumptions about the parameters at the within-host scale that are equally unknown as the parameters at the between-host scale and hence leading to uncertainty about the appropriate parametrization (Uecker and Bonhoeffer 2021). Such uncertainty is amplified by our approach where the level of resistance is considered a continuous quantitative trait, compared to the classical qualitative (or “binary”) approaches (Djidjou-Demasse et al. 2023).
Time-scale separation hypothesis Within the context of nesting within- and between-host scales, for the modelling of the epidemiology and evolution of pathogens, some studies assumed that the epidemiological and evolutionary time scales are distinct, i.e., the within-host dynamic is fast relative to the between-host dynamic such that the within-host model remains at equilibrium, e.g. (Gilchrist and Coombs 2006; Xue and Bloom 2020; Almocera et al. 2018; Boldin and Diekmann 2008; Coombs et al. 2007; André and Gandon 2006). Although such an assumption on the within-host equilibrium dynamic’ might be appropriate for chronic infections, it leads to a population-scale model that does not explicitly account for the individual time-dependent infectiousness dynamics (Hart et al. 2020). Furthermore, our illustrative examples (Figs. 3, 4) strongly highlight the infectiousness’ time-dependency of infected individuals and show that the between-host dynamic is not necessarily faster compared to the within-host dynamic.
5 Preliminaries and technical materials
We will go through details on the proof of our main results, namely Theorem 2.7. We will first discuss the existence of a positive global solution of the nested model. Next, we will give details on the derivation of the basic reproduction number of System (1.4)–(1.5) for individuals with an immune system of level i. We will also derive the existence of a unique endemic equilibrium of System (1.4)–(1.5) when \({\mathcal {R}}_0>1\), as well as the long-term persistence of the epidemic in such a case.
5.1 Existence of the semiflow
We establish the existence of a positive global solution of the system (2.5). We first formulate system (2.5) in an abstract Cauchy problem. For that, we introduce the Banach space \({\mathcal {X}}={\mathbb {R}}^n\times {\mathbb {R}}^{2n} \times L^1((0,\infty ), {\mathbb {R}}^{2n})\), endowed with the usual product norm \(\Vert \cdot \Vert _{{\mathcal {X}}}\) as well as its positive cone \({\mathcal {X}}_+\). Let \(A\,:\, D(A)\subset {\mathcal {X}}\rightarrow {\mathcal {X}}\) be the linear operator defined by \(D(A)={\mathbb {R}}^n\times \{{{\textbf {0}}}_{{\mathbb {R}}^{2n}}\}\times W^{1,1}((0,\infty ),{\mathbb {R}}^{2n})\) and
Let us introduce the non-linear map \(F\;:\; \overline{D(A)}\rightarrow {\mathcal {X}}\) defined by
By identifying \(\varvec{\varphi }(t)\) together with \(({{\textbf {S}}}(t),{{\textbf {0}}}_{L^1},{{\textbf {I}}}(t,\cdot ))^t\) and by setting \(\varvec{\varphi }_{0}=({{\textbf {S}}}_0,{{\textbf {0}}}_{L^1},{{\textbf {I}}}_0(\cdot ))^t\) the associated initial condition, system (2.5) becomes
By setting \({\mathcal {X}}_0 = D(A)\) and \({\mathcal {X}}_{0+} = {\mathcal {X}}_0 \cap {\mathcal {X}}_+\), the positivity and boundedness of the solutions of system (2.5) are provided by the following result.
Theorem 5.1
There exists a unique strongly continuous semiflow \(\big \{ \varvec{\Psi }(t,\cdot ): {\mathcal {X}}_0 \rightarrow {\mathcal {X}}_0\big \}_{t\ge 0}\) such that, for each \(\varvec{\varphi }_{0} \in {\mathcal {X}}_{0+}\), the map \(\varvec{\varphi } \in {\mathcal {C}} \left( [0,\infty ), {\mathcal {X}}_{0+} \right) \) defined by \(\varvec{\varphi }=\varvec{\Psi }(\cdot ,\varvec{\varphi }_{0})\) is a mild solution of (5.2). That is, \(\int _0^t \varvec{\varphi }(s) \textrm{d} s \in D(A)\) and \(\varvec{\varphi }(t)= \varvec{\varphi }_{0} + A \int _0^t \varvec{\varphi }(s)d s + \int _0^t F\left( \varvec{\varphi }(s)\right) \textrm{d} s\) for all \(t\ge 0\). Moreover, \(\left\{ \varvec{\Psi }(t,\cdot ) \right\} _t\) satisfies the following properties:
-
1.
Let \(\varvec{\Psi }(t,\varvec{\varphi }_{0})=\left( {{\textbf {S}}}(t),{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}(t,\cdot ) \right) ^t\), then the following Volterra formulation holds true for all \(i\in {\mathcal {I}}\)
$$\begin{aligned} {{\textbf {I}}}_i(t,\tau )=\left\{ \begin{array}{ll} \displaystyle \varvec{\Pi }_i(\tau ,\tau -t)\, {{\textbf {I}}}_{i0}(\tau -t), &{} \hbox {if} \quad t\le \tau , \\ \displaystyle \lambda (t-\tau )\;S_i(t-\tau )\;\varvec{\Pi }_i(\tau ,0)\, {{\textbf {q}}}_i, &{} \hbox {if}\quad t>\tau , \end{array} \right. \end{aligned}$$(5.3)coupled with the \(S_i(t)\) equation of (2.5), and where \(\varvec{\Pi }_i(\tau _2,\tau _1)\), \(0\le \tau _1\le \tau _2<\infty \), is the evolutionary system generated by the linear operator \(\left[ -\varvec{\Phi }_i(\tau )+ \varvec{\omega }(\tau ) \right] \).
-
2.
For all \(\varvec{\varphi }_{0} \in {\mathcal {X}}_{0+}\), and for all \(t\ge 0\), one has
$$\begin{aligned} \sum _{i\in {\mathcal {I}}}\left( S_i(t) + \int _0^\infty \Big (I_i^T(t,\tau )+ I_i^U(t,\tau ) \Big )\textrm{d}\tau \right) \le \max \left\{ \frac{{\overline{\Lambda }}}{\mu _h},N_0 \right\} , \end{aligned}$$(5.4)where \({\overline{\Lambda }}=\sum _{i\in {\mathcal {I}}}\Lambda _i\) and \(N_0= \sum _{i\in {\mathcal {I}}}\left( S_{i0} + \int _0^\infty \Big (I_{i0}^T(\tau )+ I_{i0}^U(\tau ) \Big )\textrm{d}\tau \right) \). Furthermore, the subset of the phase space
$$\begin{aligned} \left\{ ({{\textbf {S}}},{{\textbf {I}}})\in {\mathbb {R}}^n\times L^1{((0,\infty ),{\mathbb {R}}^{2n})}\, \Big |\, \displaystyle \sum _{i\in {\mathcal {I}}}\left( S_i(t) + \int _0^\infty \Big (I_i^T(t,\tau )+ I_i^U(t,\tau ) \Big )\textrm{d}\tau \right) \le \frac{{\overline{\Lambda }}}{\mu _h} \right\} , \end{aligned}$$is positively invariant and attracts all nonnegative solutions.
-
3.
The semiflow \(\left\{ \varvec{\Psi }(t,\cdot ) \right\} _t\) generated by (2.5) is bounded dissipative, that is, there exists a bounded set \({\mathcal {B}} \subset {\mathcal {X}}_{0}\) such that for any bounded set \(U \subset {\mathcal {X}}_{0}\), we can find \(\sigma = \sigma (U,{\mathcal {B}}) \ge 0\) such that \(\varvec{\Psi }(t,U) \subset {\mathcal {B}}\) for \(t\ge \sigma \).
Proof
It is easy to check that the operator A is a Hille-Yosida operator. Then standard results apply to provide the existence and uniqueness of a mild solution to (2.5) [we refer to Magal and Ruan (2009) and Thieme (2011) for more details]. The Volterra formulation is also standard and we refer to Iannelli (1995), Webb (1985) for more details.
The \(S_i\) equation of (2.5) gives \(\dot{S}_i(t) \le \Lambda _i-\mu _{h} S_i(t)\), that is
Next, for estimate (5.4), let \(\varvec{\varphi }_{0} \in {\mathcal {X}}_{0+}\), then adding up the \(S_i\),s equation together with the \(I_i^T\),s and \(I_i^U\),s equations of (1.5) yields for all \(i\in {\mathcal {I}}\)
It comes
with \({\overline{\Lambda }}= \sum _{i\in {\mathcal {I}}} \Lambda _i\). From where one deduces estimate (5.4) and which ends item 2. of the theorem.
The bounded dissipativity of the semiflow \(\left\{ \varvec{\Psi }(t,\cdot ) \right\} _t\) is a direct consequence of estimate 2. \(\square \)
The following result is straightforward.
Lemma 5.2
Let \(0\le \tau _1\le \tau _2<\infty \). By setting \(\displaystyle \varvec{\Pi }_i(\tau _2,\tau _1)= \Big (\Pi _{i,j}^k(\tau _2,\tau _1)\Big )_{k,j}\), with \(k,j\in \{T,U\}\), we have
where \(\gamma _0\), \(\gamma _1>0\) and \(\alpha _0= \max _i(\sup \alpha _i^T,\sup \alpha _i^U)\).
Proof
Let \(0\le \tau _1\le \tau _2<\infty \). Note that, for all \(\tau \ge 0\), we have \(-\omega _1\le \omega _U^T(\tau )\le \omega _1\) and \(-\omega _2\le \omega ^U_T(\tau )\le \omega _2\), with \(\omega _i>0\). Then, for all \(\tau \ge 0\),
Therefore,
where
Note that
Thus,
from where inequality (5.5) follows, and this ends the proof of the lemma. \(\square \)
5.2 The basic reproduction number
In the absence of infection, that is \({{\textbf {I}}}(t,\tau )={{\textbf {0}}}_{L^1((0,\infty ),{\mathbb {R}}^{2n})}\), the system (2.5) has a disease-free equilibrium (DFE) given by \({{\textbf {E}}}^0=({{\textbf {S}}}^0,{{\textbf {0}}}_{L^1((0,\infty ),{\mathbb {R}}^{2n})})\), with \({{\textbf {S}}}^0=\left( \Lambda _i/\mu _h \right) _{i\in {\mathcal {I}}}\). Let \(\Theta _i(t)\) be the number of new infections in the host population of group i at time t. Then in an initially infection-free population, by (2.5), we have
where \(\left\langle \cdot ; \cdot \right\rangle \) is the usual scalar product.
Linearizing the Volterra formulation (5.3) at the DFE, it comes
From where,
where \(f_i(t)\) is the number of new infections produced by the initial population. Therefore, the basic reproduction number \({\mathcal {R}}_0^i\) of individuals of group i is calculated as
with
The term \(\chi _k(\tau )= \left\langle \varvec{\beta }_k(\tau ), \varvec{\Pi }_k(\tau ,0) {{\textbf {q}}}_k \right\rangle \) quantifies the infectiousness at \(\tau \)-time post infection of the whole infected individuals of group \(k \in {\mathcal {I}}\).
Next, let \(\Theta (t)=(\Theta _i(t))_{i\in {\mathcal {I}}}\), the number of new infections in all groups at time t and, \(f(t)=(f_i(t))_{i\in {\mathcal {I}}}\) the number of new infections produced by the initial population. We have
Due to the above formulation, the basic reproduction number \({\mathcal {R}}_0\) of all individuals is calculated as the spectral radius of the matrix \((a_{i,k})_{i,k\in {\mathcal {I}}}\), where
Some calculations give
A more explicit expression of the infectiousness \(\chi _k\) is difficult to obtain in general. However, one can go further steps in some particular configurations of the treatment status transition rates \(\varvec{\omega }(\tau )=\left( \begin{array}{cc} 0 &{} \omega _T^U(\tau ) \\ \omega _U^T(\tau ) &{} 0 \\ \end{array}\right) .\) Indeed, assume that (2.7) holds. In such a configuration, we have
with \(c_k=\mu _h + \inf _\tau \alpha _k^U(\tau ) + \inf _\tau \alpha _k^T(\tau )\), and
From where, by setting
it comes
where
5.3 Proof of Theorem 2.6
The equilibrium of system (2.5) is obtained by solving the following system for all \(i\in {\mathcal {I}}\)
where
Solving (5.6) for \(S_i^*\) and \({{\textbf {I}}}_i^*\) yields
Replacing (5.8) in (5.7) leads to \(\left( 1+\mu _h^{-1}\lambda ^*\right) \lambda ^*= {\mathcal {R}}_0 \lambda ^*\) and since \(\lambda ^*>0\), we have \(\lambda ^*=\mu _h({\mathcal {R}}_0-1)\).
It follows that system (2.5) has a unique positive endemic equilibrium when \({\mathcal {R}}_0>1\), such that \(\forall i\in {\mathcal {I}}\),
where \(\lambda ^*=\mu _h({\mathcal {R}}_0-1)\).
5.4 Technical materials
Before proceed to the proof Theorem 2.7, we introduce some technical materials including the existence of a global compact attractor for the solution semiflow of Model (2.5), the spectral properties of the linearized semiflow of Model (2.5) at any given equilibrium, and the uniform persistence of Model (2.5) when \({\mathcal {R}}_0>1\).
5.4.1 Global compact attractor
To derive the global properties of the solution dynamics, it is necessary to show that the semiflow generated by system (2.5) has a global compact attractor. Denote by
and endow the set \({\mathcal {Y}}\) with the norm
For any initial condition \(\varvec{\varphi }_{0}\in {\mathcal {Y}}\), the solution semiflow of system (2.5) in \({\mathcal {Y}}_+\) is denoted by \(\displaystyle \varvec{\Psi }^*(t,\varvec{\varphi }_{0})=({{\textbf {S}}}(t),{{\textbf {I}}}(t,\cdot ))^t\). From the Volterra formulation (5.3), we rewrite system (2.5) as follows for all \(i\in {\mathcal {I}}\):
where
We need to prove the following claim.
Claim 5.3
Let Assumption 2.2 be satisfied. Then, function \(\lambda (\cdot )\) is Lipschitz continuous on \({\mathbb {R}}_+\).
Proof of Claim 5.3
Let \(C_0\ge \max \left\{ \frac{{\overline{\Lambda }}}{\mu _h},\Vert \varvec{\varphi }_0\Vert \right\} \), \(\Vert \beta _i\Vert _{\infty }=\max \Big \{\Vert \beta _i^T\Vert _{\infty }, \Vert \beta _i^U\Vert _{\infty }\Big \}\) and \(\Vert \beta \Vert _{\infty }=\max _{i\in {\mathcal {I}}}\Vert \beta _i\Vert _{\infty }\). Then, \(|\lambda (t)|\le C_0\Vert \beta \Vert _{\infty }\). Let \(t>0\) and \(h>0\). It comes that
Recalling (5.3) and combining the integrals, we obtain
We have \(|\gamma _1e^{-\mu _hh}-1|\le |e^{-\mu _hh}-1|\le \mu _hh\). Using the Lipschitzianity of \(\beta _i^{\vartheta }\), we find a positive constant \(C_{\beta }\) such that
where \(C_{\lambda }=C_0^2\Vert \beta \Vert ^2_{\infty }+C_0 \Vert \beta \Vert _{\infty }\mu _h+C_{\beta }C_0\). \(\square \)
Next, we will show that system (2.5) has a global attractor. By using the similar method as in Martcheva and Thieme (2003) and Cheng et al. (2018), we can state the following result.
Lemma 5.4
There exists \({\mathcal {A}}_{0}\), a compact subset of \({\mathcal {Y}}_+\), which is a global attractor for the solution semiflow of system (2.5). Moreover, \({\mathcal {A}}_{0}\) is invariant under the solution semiflow, that is
Proof
We show that \(\varvec{\Psi }^*\) satisfies the assumptions of Lemma 3.2.3 and Theorem 3.4.6 in Hale (2010). To this end, we split the solution semiflow into two parts. For any initial condition \(\varvec{\varphi }_{0}\in {\mathcal {Y}}_+\), we let \(\varvec{\Psi }^*(t,\varvec{\varphi }_{0})=\widehat{\varvec{\Psi }}^*(t,\varvec{\varphi }_{0})+ \widetilde{\varvec{\Psi }}^*(t,\varvec{\varphi }_{0})\), where
In such a way, we need to prove the following claim:
Claim 5.5
-
(1)
\(\widehat{\varvec{\Psi }}^*(t,\varvec{\varphi }_{0})\rightarrow 0\) as \(t\rightarrow \infty \) for every \(\varvec{\varphi }_{0}\) in \({\mathcal {Y}}\).
-
(2)
For a fixed t and any bounded set B in \({\mathcal {Y}}\), the set \(\{\widetilde{\varvec{\Psi }}^*(t,\varvec{\varphi }_{0})\;:\; \varvec{\varphi }_{0}\in B\}\) is precompact.
Proof of Claim 5.5
Now, we show that the first claim holds.
From (5.9) and Lemma 5.2, we have
Note that for any bounded \(\varphi _{0}\), \(2\gamma _1 e^{-\mu _ht} \Vert \varvec{\varphi }_{0}\Vert \rightarrow 0\) as \(t\rightarrow \infty \). This completes the first claim.
To show that the second claim holds, let \(B\subset {\mathcal {Y}}\) be a bounded subset such that \(\varvec{\Psi }^*(t,\cdot )B\subset B\). Choose \(C_0>0\) such that \(\Vert \varvec{\varphi }_{0}\Vert \le C_0\) for all \(\varvec{\varphi }_{0}\in B\). From Theorem 5.1 Item 2, \(\cup _{\varvec{\varphi }_{0}\in B}\{{{\textbf {S}}}(t)\}\) is bounded in \({\mathbb {R}}^n\) and then is precompact in \({\mathbb {R}}^n\). Hence, to show the compactness, it suffices to show that the set \(\widetilde{\varvec{\Psi }}^*(t,\varvec{\varphi }_{0})B\) is precompact for
By Frechet–Kolmogorov theorem [see Theorem B.2 in Smith and Thieme (2011)], it is sufficient to verify the following conditions:
-
(i)
\(\sup \limits _{\varvec{\varphi }\in B}\sum _{i\in {\mathcal {I}}}\int _{0}^{\infty }\Big ({\widetilde{I}}_i^T(t,\tau )+{\widetilde{I}}_i^U(t,\tau )\Big )\textrm{d}\tau <\infty \),
-
(ii)
\(\lim \limits _{h\rightarrow \infty }\sum _{i\in {\mathcal {I}}}\int _{h}^{\infty } \Big (|{\widetilde{I}}_i^T(t,\tau )|+ |{\widetilde{I}}_i^U(t,\tau )|\Big )\textrm{d}\tau =0\) uniformly with respect to \(\varvec{\varphi }_{0}\in B\).
-
(iii)
\(\lim \limits _{h\rightarrow 0}\sum _{i\in {\mathcal {I}}}\int _{0}^{\infty } \Big (|{\widetilde{I}}_i^T(t,\tau )-{\widetilde{I}}_i^T(t,\tau +h)|+ |{\widetilde{I}}_i^U(t,\tau )-{\widetilde{I}}_i^U(t,\tau +h)|\Big )\textrm{d}\tau =0\) uniformly with respect to \(\varvec{\varphi }_{0}\in B\).
-
(iv)
\(\lim \limits _{h\rightarrow 0}\sum _{i\in {\mathcal {I}}}\int _{0}^{h} \Big (|{\widetilde{I}}_i^T(t,\tau )|+ |{\widetilde{I}}_i^U(t,\tau )|\Big )\textrm{d}\tau =0\) uniformly with respect to \(\varvec{\varphi }_{0}\in B\).
By (5.11) we have for all \(i\in {\mathcal {I}}\)
It follows that above conditions (i), (ii) and (iv) are satisfied.
Now, we show that condition (iii) holds. We have for \(i\in {\mathcal {I}}\) and \(h\le t\),
By Lemma 5.2, and the boundedness of the semiflow, we can find a positice constant \(C_0\) such that
Again by Lemma 5.2, we have
with \(C_{S_i}=\Lambda _i+C_{i0}(C_0\Vert \beta \Vert _{\infty }+\mu _h)\). By (5.12) and (5.13) one concludes that the criterion (iii) holds, and then the second claim holds.
This completes the proof of the lemma. \(\square \)
5.4.2 Spectral properties of the linearized semiflow
The next result is concerned with spectral properties of the linearized semiflow \(\varvec{\Psi }\) of Model (2.5) at a given equilibrium point \(\widetilde{\varvec{\varphi }}\in {\mathcal {X}}_{0+}\). The associated linearized system (2.5) at the point \(\widetilde{\varvec{\varphi }}\) reads as
where A is the linear operator defined in (5.1) while \(G[\widetilde{\varvec{\varphi }}]\in {\mathcal {L}}({\mathcal {X}}_{0},{\mathcal {X}})\) is the bounded linear operator defined by:
where \( {\widetilde{\lambda }}=\sum _i\int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ), \widetilde{{{\textbf {I}}}}_i(\tau )\right\rangle \textrm{d}\tau \) and \( \lambda =\sum _i\int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ), {{\textbf {I}}}_i(\tau )\right\rangle \textrm{d}\tau \). We then have the following lemma.
Lemma 5.6
Let us set \(\Omega =\{\nu \in {\mathbb {C}}\,: \, Re(\nu )>-\mu _h\}\). Then, the spectrum \(\sigma (A+G[\widetilde{\varvec{\varphi }}])\cap \Omega \ne \emptyset \) only consists of the point spectrum and one has
where function \(\Delta (\cdot ,\varvec{{\widetilde{\varphi }}}):\Omega \longrightarrow {\mathbb {C}}\) is defined by
with \({\mathcal {R}}_\nu [\widetilde{\varvec{\varphi }}]=\sum _{i\in {\mathcal {I}}} {\widetilde{S}}_i\; \int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ), \varvec{\Pi }_i(\tau ,0)\;{{\textbf {q}}}_i \right\rangle e^{-\nu \tau }\textrm{d}\tau . \)
Proof
Let us denote by \(A_{0}\,:\, D(A_{0})\subset {\mathcal {X}}_0\rightarrow {\mathcal {X}}_0\) the part of A in \({\mathcal {X}}_0=D(A)\), which is defined by
Then, it is the infinitesimal generator of a \(C_0\)-semigroup on \({\mathcal {X}}_0\) denoted by \(\{T_{A_{0}}(t)\}_{t\ge 0}\). Let \(\varvec{\varphi }=({{\textbf {S}}},{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}(\cdot ))^t\). We find that
Then, for \(t\ge \tau _0\), we have \(\Vert T_{A_{0}}(t-\tau _0)\varvec{\varphi }\Vert _{{\mathcal {X}}}\le e^{-\mu _h(t-\tau _0)}\Vert \varvec{\varphi }\Vert _{{\mathcal {X}}}, \forall t\ge \tau _0.\) We deduce that the growth rate \(\omega _0(A_{0})=\lim \limits _{t\rightarrow +\infty } \frac{\ln \Big (\Vert T_{A_{0}}(t)\Vert _{{\mathcal {L}}({\mathcal {X}})}\Big )}{t}\) of this semigroup satisfies \(\omega _0(A_{0})\le -\mu _h.\) Since operator \(G_i[\widetilde{\varvec{\varphi }}]\) is compact, the results in Arino et al. (1998) or Ducrot et al. (2008) apply and provided that the essential growth rate of \(\Big \{ T_{(A+G[\widetilde{\varvec{\varphi }}])_0}(t)\Big \}_{t\ge 0}\)-the \(C_0\)-semigroup generated by the part of \((A+G[\widetilde{\varvec{\varphi }}])\) in \({\mathcal {X}}_0\) satisfies
By results in Engel and Nagel (2001) and Webb (1987), the latter inequality ensures that \(\Omega \cap \sigma (A+G[\widetilde{\varvec{\varphi }}])\ne \emptyset \), and it is only composed of point spectrum of \((A+G[\widetilde{\varvec{\varphi }}])\).
It remains to derive the characteristic equation. Let \(\nu \in \rho (A+G[\widetilde{\varvec{\varphi }}])\), where \(\rho (\cdot )\) stands for the resolvent. For \(\widehat{\varvec{\varphi }}= (\widehat{{{\textbf {S}}}},\widehat{{{\textbf {u}}}}, \widehat{{{\textbf {I}}}}(\cdot ))^t\in {\mathcal {X}}\) and \(\varvec{\varphi }= ({{\textbf {S}}},0_{L^1}, {{\textbf {I}}}(\cdot ))^t\in D(A)\), we have \((\nu I-A-G[\widetilde{\varvec{\varphi }}])\varvec{\varphi }=\widehat{\varvec{\varphi }}\), that is \((\nu I- A)\varvec{\varphi }-G[\widetilde{\varvec{\varphi }}]\varvec{\varphi }=\widehat{\varvec{\varphi }}\), and from where
Since
we find that
Thus, for all \(i\in {\mathcal {I}}\), equality (5.15) rewrites as
Substituting (5.16) into expression for \(\lambda \), it comes
where \(\displaystyle {\mathcal {R}}_\nu [\varvec{{\widetilde{\varphi }}}]=\sum _{i\in {\mathcal {I}}} {\widetilde{S}}_i\; \int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ), \varvec{\Pi }_i(\tau ,0){{\textbf {q}}}_i \right\rangle \; e^{-\nu \tau }\;\textrm{d}\tau \) and \( \widehat{{{\textbf {y}}}}_i(\tau )= e^{-\nu \tau }\varvec{\Pi }_i(\tau ,0)\widehat{{{\textbf {u}}}}_i+ \int _{0}^{\tau } \varvec{\Pi }_i(\tau ,s)\widehat{{{\textbf {I}}}}_i(s) e^{-\nu (\tau -s)}\textrm{d}s \). Therefore, we can isolate \(\lambda \) in system (5.17) if and only if \(\displaystyle \Delta (\nu ,\varvec{{\widetilde{\varphi }}})=1- {\mathcal {R}}_{\nu } [\varvec{{\widetilde{\varphi }}}] + \frac{{\mathcal {R}}_\nu [\widetilde{\varvec{\varphi }}]\,{\widetilde{\lambda }}}{\nu +\mu _h+{\widetilde{\lambda }}}\ne 0\). \(\square \)
5.5 Uniform persistence
Our next technical material concerns the uniform persistence of Model (2.5) when \({\mathcal {R}}_0>1\) by using the method developed in Theorem 5.2 in Smith and Thieme (2011).
For the invariant sets of uniform persistence, we introduce
where
For the unique solution \(\varvec{\varphi }=({{\textbf {S}}},{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}})\) of system (2.5) associated to the initial condition \(\varvec{\varphi }_0=({{\textbf {S}}}_0,{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}_0(\cdot ))\in {\mathcal {M}}\), we define \(\varvec{\Psi }(t,\varvec{\varphi }_0)=({{\textbf {S}}}(t),{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}(t,\cdot ))\) the semiflow of Model (2.5) passing through \(\varvec{\varphi }_{0}\). Next, we first claim that
Claim 5.7
The subsets \({\mathcal {M}}\) and \(\partial {\mathcal {M}}\) are positively invariant with respect to the semiflow \(\varvec{\Psi }(t,\cdot )\) generated by system (2.5). Furthermore, \(\lim \nolimits _{t\rightarrow \infty }\varvec{\Psi }(t,\varvec{\varphi }_{0})=({{\textbf {S}}}^0,{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {0}}}_{L_+^1((0,\infty ),{\mathbb {R}}^{2n})})^T\) for each \(\varvec{\varphi }_{0}\in \partial {\mathcal {M}}\).
Proof of Claim 5.7
Let \(\varvec{\varphi }_{0}=({{\textbf {S}}}_{0},{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}_{0}(\cdot ))^T\in {\mathcal {M}}\) be given and \(\varvec{\Psi }(t,\varvec{\varphi }_{0})=({{\textbf {S}}},{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}(t,\cdot ))^T,\) the orbit passing through \(\varvec{\varphi }_{0}\). Since \(\varvec{\varphi }_{0}\in {\mathcal {M}}\), then \(\lambda (0)>0\). Through a direct calculation, we have
where \(\alpha _0=\max _{i\in {\mathcal {I}}}\{\sup \alpha _i^T,\sup \alpha _i^U\}\) and \(\omega _0=\max \{{\overline{\omega }}_T^U,{\overline{\omega }}_U^T\}\). Thus, one obtains that
for \(t\ge 0\). This complete the fact that \({\mathcal {M}}\) is positively invariant.
Now, let \(\varvec{\varphi }_{0}\in \partial {\mathcal {M}}\). Since for all \(i\in {\mathcal {I}}\), \(S_i(t)\le S_i^0\) as t is large enough, the comparison principle implies that
where \(\widetilde{{{\textbf {I}}}}_{i}(t,\tau )\) is the solution of the following system
By the Volterra formulation, we have from (5.19) that
where \({\widetilde{\lambda }}(t)\) satisfies
The initial condition with \(\displaystyle \sum _{k\in {\mathcal {I}}}\int _0^{\infty } \left\langle \varvec{\beta }_k(\tau ), {{\textbf {I}}}_{k0}(\tau ) \right\rangle \textrm{d}\tau =0\), leads to
Since \({\widetilde{\lambda }}(0)=0\), we have \({\widetilde{\lambda }}(t)=0\) for all \(t\ge 0\), and then \(\widetilde{{{\textbf {I}}}}_{i}(t,\cdot )=0\) for all \(i\in {\mathcal {I}}\) and \(t\ge 0\). The comparison in (5.18) implies that \({{\textbf {I}}}_{i}(t,\cdot )=0\) for all \(i\in {\mathcal {I}}\) and \(t\ge 0\) and then \(\partial {\mathcal {M}}\) is positively invariant under the semiflow \(\varvec{\Psi }(t,\cdot )\). In addition, it is clear for the solution remaining in \(\partial {\mathcal {M}}\), we have for all \(i\in {\mathcal {I}}\), \(S_i\rightarrow S_i^0\). Hence, \(\lim \nolimits _{t\rightarrow \infty } \varvec{\Psi }(t,\varvec{\varphi }_{0})=({{\textbf {S}}}^0,{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {0}}}_{L_+^1((0,\infty ),{\mathbb {R}}^{2n})})^T\) for each \(\varvec{\varphi }_{0}\in \partial {\mathcal {M}}\). This ends the proof of Claim 5.7. \(\square \)
Finally, we end this technical material section by establishing the uniform persistence of system (2.5).
Theorem 5.8
The semiflow \(\{\varvec{\Psi }(t,\cdot )\}_{t\ge 0}\) generated by system (2.5) is uniformly persistent in \({\mathcal {M}}\) with respect to \(({\mathcal {M}},\partial {\mathcal {M}})\), that is, there exists a constant \(\eta >0\) such that for each \(\varvec{\varphi }_{0}\in {\mathcal {M}}\),
and
Furthermore, there exists compact global attractor \({\mathcal {A}}_1\) in \({\mathcal {M}}\) for the semiflow \(\{\varvec{\Psi }(t,\cdot )\}_{t\ge 0}\).
Proof
In the following, we will prove that \(W^S(\{{{\textbf {E}}}^0\})\cap {\mathcal {M}}=\emptyset \), where
Since from Claim 5.7 the disease-free equilibrium \({{\textbf {E}}}^0\) is globally asymptotically stable in \(\partial {\mathcal {M}}\), we need only to study the behavior of the solution starting in \({\mathcal {M}}\) in some neighborhood of \({{\textbf {E}}}^0\). To this end, it is sufficient to show that there exists \(\sigma >0\) satisfying for each \(\varvec{\varphi }\in \{{{\textbf {v}}}\in {\mathcal {M}}\,:\, \Vert {{\textbf {E}}}^0-{{\textbf {v}}}\Vert \le \sigma \}\) there exists \(t_0\ge 0\) such that \(\Vert \varvec{\Psi }(t,\varvec{\varphi }_0)-{{\textbf {E}}}^0\Vert >\sigma \).
By the way of contradiction, suppose that for each integer \(n\ge 0\) there exists a \(\varvec{\varphi }_{0}^n=({{\textbf {S}}}_{0}^n,{{\textbf {0}}}_{L^1},{{\textbf {I}}}_{0}^n)\in \{{{\textbf {v}}}\in {\mathcal {M}}\,:\, \Vert {{\textbf {E}}}^0-{{\textbf {v}}}\Vert \le \sigma \}\) such that
Denote \(\varvec{\Psi }(t,\varvec{\varphi }_{0}^n)=({{\textbf {S}}}^n(t),{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {I}}}^n(t,\cdot ))\), then for all \(t\ge 0\) we have
It follows that for all i, we have \(S_i^n(t)\ge S_i^0-\frac{1}{n+1}\) for all \(t\ge 0\). Consider the following system
where \(\displaystyle \lambda ^n(t)= \sum _{i\in {\mathcal {I}}} \int _0^t \left\langle \varvec{\beta }_i(\tau ), {{\textbf {I}}}^n_i(t,\tau ) \right\rangle \textrm{d}\tau \), \({{\textbf {S}}}_0^n=({{\textbf {S}}}_{i0}^n)_{i\in {\mathcal {I}}}\) and \({{\textbf {I}}}_0^n(\cdot )=({{\textbf {I}}}_{i0}^n(\cdot ))_{i\in {\mathcal {I}}}\). By the comparison principle, we have
where \(\widetilde{{{\textbf {I}}}}^n(t,\cdot )\) is the solution of the following auxiliary system
which gives for all i,
For ease of notation, let us rewrite the system (5.22) as the following form:
\(\widetilde{{{\textbf {v}}}}^n_i(0)\in \overline{D({\widetilde{A}}^n_i)}\), the closure of \(D({\widetilde{A}}^n_i)=\left\{ {{\textbf {0}}}_{{\mathbb {R}}^2}\right\} \times W^{1,1}((0,\infty ),{\mathbb {R}}^2)\), where \(\widetilde{{{\textbf {v}}}}^n_i(t)=({{\textbf {0}}}_{{\mathbb {R}}^2},\widetilde{{{\textbf {I}}}}^n_i(t,\cdot ))^t\) and the operators \({\widetilde{A}}^n_i\) and \(\widetilde{{\mathcal {L}}}^n_i\) are defined as
and
Similarly to the proof of Lemma 5.6, we can derive the characteristic equation \(\Delta (\nu _n)=0\) for system (5.22), where
Since \({\mathcal {R}}_0>1\), there exists \(n_0>0\) large enough such that for \(n\ge n_0\),
The largest eigenvalue \(\nu _n^*\) of system (5.23) satisfies the characteristic equation \(\Delta (\nu _n)=0\). Furthermore, \({\mathcal {R}}_0^n>1\) implies the existence of a dominant eigenvalue \(\nu _n^*>0\) such that \(\Delta (\nu _n^*)=0\). Therefore, \(\nu _{n}^*>0\) is a simple dominant eigenvalue of \(({\widetilde{A}}^n_i+\widetilde{{\mathcal {L}}}^n_i)\). From Lemma 5.6, we have shown that \(\omega _{0,ess}({\widetilde{A}}^n_i+\widetilde{{\mathcal {L}}}^n_i)\le -\mu _h\) and since the semigroup \(\left\{ T_{({\widetilde{A}}^n_i+\widetilde{{\mathcal {L}}}^n_i)}(t)\right\} _{t\ge 0}\) is irreducible, it follows from Corollary 4.6.8 in Magal and Ruan (2018) that \(\left\{ T_{({\widetilde{A}}^n_i+ \widetilde{{\mathcal {L}}}^n_i)}(t)\right\} _{t\ge 0}\) has asynchronous exponential growth with intrinsic growth constant \(\nu _n^*\in {\mathbb {R}}\). Therefore, using Theorem 3.9 in Magal and McCluskey (2013), we have
and there exist constants \(\varepsilon _0>0\) and \(\eta _0>0\) such that
where \({\widetilde{\Pi }}_{\nu _{n}^*}\) is the projector on the generalized eigenspace associated with the largest eigenvalue \(\nu _n^*>0\). We deduce that
Since \(\nu _n^*>0\), it follows that \(\lim \limits _{t\rightarrow \infty }\Vert {\widetilde{\Pi }}_{\nu _{n}^*}\,{\widetilde{v}}_{i}^n(t)\Vert _{L^1}=+\infty \). Therefore, \(\lim \limits _{t\rightarrow \infty }\Vert \widetilde{{{\textbf {I}}}}_i^n(t,\cdot )\Vert _{L^1}=+\infty \) and from (5.21), we have \(\lim \limits _{t\rightarrow \infty }\Vert {{\textbf {I}}}_i^n(t,\cdot )\Vert _{L^1}=+\infty \), which is a contradiction to the boundedness of the solution. Thus, \(W^S(\{{{\textbf {E}}}^0\})\cap {\mathcal {M}}=\emptyset \) and we derive from Theorem 4.2 in Hale and Waltman (1989) that the semiflow \(\{\varvec{\Psi }(t,\cdot )\}_{t\ge 0}\) is uniform persistent with respect to the pair \(({\mathcal {M}},\partial {\mathcal {M}})\). Moreover, by Theorem 3.7 in Magal and Zhao (2005), there exists a compact global attractor \({\mathcal {A}}_1\subset {\mathcal {M}}\) for the semiflow \(\{\varvec{\Psi }(t,\cdot )\}_{t\ge 0}\). \(\square \)
6 Proof of Theorem 2.7
The proof of Theorem 2.7 is decomposed into two parts. The first part is devoted to the global stability of the disease-free equilibrium, while the second part is devoted to the global stability of the endemic equilibrium.
6.1 Proof of Theorem 2.7 (i): global stability of the disease-free equilibrium
When the initial condition of System (2.5) satisfies \(\varvec{\varphi }_{0}\in \partial {\mathcal {M}}\), i.e., \(\sum _{i\in {\mathcal {I}}}\int _0^{\infty } \left\langle \varvec{\beta }_i(\tau ),{{\textbf {I}}}_{i0}(\tau ) \right\rangle \textrm{d}\tau =0\), by Claim 5.7, it comes that the semiflow \(\varvec{\Psi }(t,\varvec{\varphi }_{0})\) generated by system (2.5) is such that \(\lim \limits _{t\rightarrow \infty }\varvec{\Psi }(t,\varvec{\varphi }_{0})=({{\textbf {S}}}^0,{{\textbf {0}}}_{{\mathbb {R}}^{2n}},{{\textbf {0}}}_{L_+^1((0,\infty ),{\mathbb {R}}^{2n})})\). It then remain to prove the global stability of the disease-free equilibrium when \({\mathcal {R}}_0\le 1\).
Theorem 6.1
If \({\mathcal {R}}_0\le 1\) then, the disease-free equilibrium \({{\textbf {E}}}^0=({{\textbf {S}}}^0,\,{{\textbf {0}}}_{L^1((0,\infty )\times {\mathbb {R}},{\mathbb {R}}^{2n})})\) of system (2.5) is globally asymptotically stable.
Proof
By Theorem 5.8, we introduce the following well defined Lyapunov functional \(V(t)=V_1(t)+V_2(t)\), with
where \({{\textbf {c}}}_k(\tau )=(c_k^T(\tau ),\; c_k^U(\tau ))^T\) is a vector of positive constants such that
with \(\displaystyle {\overline{S}}^0=\sum _{k\in {\mathcal {I}}} S_k^0\) and \({{\textbf {c}}}_k(\tau )\rightarrow 0\) as \(\tau \rightarrow \infty \). From (6.1), we have
with
Differentiating \(V_1(t)\) and using \(\Lambda _k=\mu _h S_k^0\), we have
Differentiating \(V_2(t)\), we have
By integrating by parts, we have
and
Replacing these expressions in \(\frac{dV_2(t)}{dt}\) and using the fact that \(I_k^{T}(t,0)=q_k^T\lambda (t)S_k(t)\) and \(I_k^{U}(t,0)=q_k^U\lambda (t)S_k(t)\), we have
Finally, combining \(\frac{dV_1(t)}{dt}\) and \(\frac{dV_2(t)}{dt}\), gathering some terms and using (6.2), it follows that
Using (6.1), we have
That implies that for all k, \(\left\langle {{\textbf {c}}}_k(0),{{\textbf {q}}}_k \right\rangle \le 1\), when \({\mathcal {R}}_0\le 1\). Therefore, we have \(\frac{dV(t)}{dt}\le 0\) when \({\mathcal {R}}_0\le 1\). The strict equality holds only if \(S_k(t)=S_k^0\) hold simultaneously with either \({\mathcal {R}}_0=1\) or \({{\textbf {I}}}_k(t,0)=0\). It is easy to verify that largest invariant set in \(\left\{ \frac{dV}{dt}=0\right\} \) is the singleton \(\{{{\textbf {E}}}^0\}\). Thus, all solutions of system (2.5) converge to the disease-free equilibrium \({{\textbf {E}}}^0\). Hence, \({{\textbf {E}}}^0\) is globally asymptotically stable when \({\mathcal {R}}_0\le 1\). \(\square \)
6.2 Proof of Theorem 2.7 (ii): global stability of the endemic equilibrium
Theorem 6.2
Assume \({\mathcal {R}}_0>1\), then the endemic equilibrium \({{\textbf {E}}}^*=({{\textbf {S}}}^*,{{\textbf {I}}}^*)^t\) of system (2.5) is globally asymptotically stable in \({\mathcal {Y}}_+\).
Proof
By Theorem 5.8, we introduce the following well defined Lyapunov functional \(L(t)=L_1(t)+L_2(t)+L_3(t)\), where
with h the function defined by \(h(z)=z-1-\ln z\) \((z\in {\mathbb {R}}_+)\), and \({{\textbf {d}}}_k(\tau )=(d_k^T(\tau ),\; d_k^U(\tau ))^T\) a vector of positive constants given by
where \({\overline{S}}^*=\sum _{k\in {\mathcal {I}}}S_k^*\) and \({{\textbf {d}}}_k(\tau )\rightarrow 0\) as \(\tau \rightarrow \infty \). From (6.3), we have
By using the property of function h, we find that the function L(t) is nonnegative with its global minimum point \({{\textbf {E}}}^*\).
Step 1: Differentiating \(L_1(t)\) along the solution of system (1.5) and using \(\Lambda _k=S_k^*\lambda ^*+\mu _h S_k^*\), we obtain
Step 2: Note that
Using (6.5), we have
Differentiating \(L_2(t)\) and, using (6.6) and integration by parts, we obtain
Since \(I_k^{T*}(0)=q_k^T\lambda ^*S_k^*\) and \(I_k^{T}(t,0)=q_k^T\lambda (t)S_k(t)\), then we have
By a similarly manner, the derivative of \(L_3(t)\) gives
Step 3: Finally, combining \(\frac{dL_1(t)}{dt}\), \(\frac{dL_2(t)}{dt}\) and \(\frac{dL_3(t)}{dt}\), we obtain
We observe that \(d_k^T(\tau )\) and \(d_k^U(\tau )\) satisfy
By using (6.7), (6.5), and (6.4), we obtain
Moreover, by using (5.7) and (5.8), we have
Thus, using (6.3) and (6.9), it can be verified that
which implies that for all k, \(\left\langle {{\textbf {d}}}_k(0),{{\textbf {q}}}_k \right\rangle =d_k^T(0) q_k^T+ d_k^U(0) q_k^U= 1\). In addition, note that
Replacing (6.7), (6.8), (6.10) and (6.11) in \(\frac{dL(t)}{dt}\) and gathering some terms, we obtain
Note that
Hence, we have
Note that
Finally, we have
Thus, \(\dfrac{dL(t)}{dt}\le 0\) with equality if and only if \(S_k(t)=S_k^*\), \(I_k^{T}(t,\tau )=I_k^{T*}(\tau )\) and \(I_k^{U}(t,\tau )=I_k^{U*}(\tau )\). Then, it can be verified that largest invariant set in \(\left\{ \frac{dL}{dt}=0\right\} \) is the singleton \(\{{{\textbf {E}}}^*\}\). It follows that the compact global attractor \({\mathcal {A}}_0\), stated by Lemma 5.4, is such that \({\mathcal {A}}_0=\{{{\textbf {E}}}^*\}\). Therefore, the endemic equilibrium \({{\textbf {E}}}^*\) is globally asymptotically stable in \({\mathcal {Y}}_+\) when \({\mathcal {R}}_0> 1\). \(\square \)
References
Almocera AES, Nguyen VK, Hernandez-Vargas EA (2018) Multiscale model within-host and between-host for viral infectious diseases. J Math Biol 77(4):1035–1057
André J-B, Gandon S (2006) Vaccination, within-host dynamics, and virulence evolution. Evol; Int J Org Evol 60(1):13–23
Arino O, Axelrod D, Kimmel M, Capasso V, Fitzgibbon W, Jagers P, Kirschner D, Mode C, Novak B, Sachs R, Stephan W, Swierniak A, Thieme H, Boussouar A (1998) Advances in mathematical population dynamics? Molecules, cells and man. In: Advances in mathematical population dynamics ? Molecules, cells and man, volume 6 of series in mathematical biology and medicine. World Scientific, pp 1–910
Beardmore RE, Peña-Miller R, Gori F, Iredell J (2017) Antibiotic cycling and antibiotic mixing: Which one best mitigates antibiotic resistance? Mol Biol Evol 34(4):802–817
Blanquart F (2019) Evolutionary epidemiology models to predict the dynamics of antibiotic resistance. Evol Appl 12(3):365–383
Boldin B, Diekmann O (2008) Superinfections can induce evolutionarily stable coexistence of pathogens. J Math Biol 56(5):635–672
Burie J-B, Djidjou-Demasse R, Ducrot A (2020) Asymptotic and transient behaviour for a nonlocal problem arising in population genetics. Eur J Appl Math 31(1):84–110
Cheng C-Y, Dong Y, Takeuchi Y (2018) An age-structured virus model with two routes of infection in heterogeneous environments. Nonlinear Anal Real World Appl 39:464–491
Coombs D, Gilchrist MA, Ball CL (2007) Evaluating the importance of within- and between-host selection pressures on the evolution of chronic pathogens. Theor Popul Biol 72(4):576–591
D’Agata EMC, Dupont-Rouzeyrol M, Magal P, Olivier D, Ruan S (2008) The impact of different antibiotic regimens on the emergence of antimicrobial-resistant bacteria. PLoS One 3(12):e4036
Day T, Read AF (2016) Does high-dose antimicrobial chemotherapy prevent the evolution of resistance? PLoS Comput Biol 12(1):e1004689
Diekmann O, Heesterbeek JAP, Metz JAJ (1990) On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol 28(4):365–382
Djidjou-Demasse R (2021) Antimicrobial-Quantitative-Resistance-v1.0. Zenodo
Djidjou-Demasse R, Ducrot A, Fabre F (2017) Steady state concentration for a phenotypic structured problem modeling the evolutionary epidemiology of spore producing pathogens. Math Models Methods Appl Sci 27(02):385–426
Djidjou-Demasse R, Alizon S, Sofonea MT (2021) Within-host bacterial growth dynamics with both mutation and horizontal gene transfer. J Math Biol 82(3):16
Djidjou-Demasse R, Sofonea MT, Choisy M, Alizon S (2023) Within-host evolutionary dynamics of antimicrobial quantitative resistance. Math Model Nat Phenom 18:24
Ducrot A, Liu Z, Magal P (2008) Essential growth rate for bounded linear perturbation of non-densely defined Cauchy problems. J Math Anal Appl 341(1):501–518
Elderd BD, Mideo N, Duffy MA (2022) Looking across scales in disease ecology and evolution. Am Nat 199(1):51–58
Engel K-J, Nagel R (2001) One-parameter semigroups for linear evolution equations. Semigroup Forum 63(2):278–280
Geritz SAH, Metz JAJ, Kisdi É, Meszéna G (1997) Dynamics of adaptation and evolutionary branching. Phys Rev Lett 78(10):2024–2027
Gilchrist MA, Coombs D (2006) Evolution of virulence: interdependence, constraints, and selection using nested models. Theor Popul Biol 69(2):145–153
Hale JK (2010) Asymptotic behavior of dissipative systems. American Mathematical Society
Hale JK, Waltman P (1989) Persistence in infinite-dimensional systems. SIAM J Math Anal 20(2):388–395
Hart WS, Maini PK, Yates CA, Thompson RN (2020) A theoretical framework for transitioning from patient-level to population-scale epidemiological dynamics: influenza A as a case study. J R Soc Interface 17(166):20200230
Inaba H (2012) On a new perspective of the basic reproduction number in heterogeneous environments. J Math Biol 65(2):309–348
Kepler TB, Perelson AS (1998) Drug concentration heterogeneity facilitates the evolution of drug resistance. Proc Natl Acad Sci USA 95(20):11514–11519
Larsson DGJ, Flach C-F (2022) Antibiotic resistance in the environment. Nat Rev Microbiol 20(5):257–269
Lipsitch M, Levin BR (1997) The population dynamics of antimicrobial chemotherapy. Antimicrob Agents Chemother 41(2):363–373
Iannelli M (1995) Mathematical Theory of Age-Structured Population Dynamics. Giardini Editori e Stampatori in Pisa
Magal P, McCluskey C (2013) Two-group infection age model including an application to nosocomial infection. SIAM J Appl Math 73(2):1058–1095
Magal P, Ruan S (2009) On semilinear Cauchy problems with non-dense domain. Adv Differ Equ 14(11–12):1041–1084
Magal P, Ruan S (2018) Theory and applications of abstract semilinear cauchy problems, applied mathematical sciences, vol 201. Springer, Cham
Magal P, Zhao X-Q (2005) Global attractors and steady states for uniformly persistent dynamical systems. SIAM J Math Anal 37(1):251–275
Martcheva M, Thieme HR (2003) Progression age enhanced backward bifurcation in an epidemic model with super-infection. J Math Biol 46(5):385–424
Metz JAJ, Geritz SAH, Meszena G, Jacobs FJA, van Heerwaarden JS (1996) Adaptive dynamics: a geometrical study of the consequences of nearly faithful reproduction. North-Holland, Amsterdam
Millan AS, Peña-Miller R, Toll-Riera M, Halbert ZV, McLean AR, Cooper BS, MacLean RC (2014) Positive selection and compensatory adaptation interact to stabilize non-transmissible plasmids. Nat Commun 5(1):5208
Shen M, Xiao Y, Rong L, Zhuang G (2019) Global dynamics and cost-effectiveness analysis of HIV pre-exposure prophylaxis and structured treatment interruptions based on a multi-scale model. Appl Math Model 75:162–200
Smith HL, Thieme HR (2011) Dynamical systems and population persistence. American Mathematical Society
Tazzyman SJ, Bonhoeffer S (2014) Plasmids and evolutionary rescue by drug resistance. Evolution 68(7):2066–2078
Thieme HR (2011) Global stability of the endemic equilibrium in infinite dimension: Lyapunov functions and positive operators. J Differ Equ 250(9):3772–3801
Uecker H, Bonhoeffer S (2021) Antibiotic treatment protocols revisited: the challenges of a conclusive assessment by mathematical modelling. J R Soc Interface 18(181):20210308
Webb GF (1985) Theory of nonlinear age-dependent population dynamics. CRC Press
Webb GF (1987) An operator-theoretic formulation of asynchronous exponential growth. Trans Am Math Soc 303(2):751–763
Xue KS, Bloom JD (2020) Linking influenza virus evolution within and between human hosts. Virus Evol 6(1):veaa010
Acknowledgements
This work was supported by the ANR (https://anr.fr/en/) grant QUASAR (grant agreement ANR-21-CE45-0004). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors of this article declare that they have no financial conflict of interest with the content of this article.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Mann-Manyombe, M.L., Mendy, A., Seydi, O. et al. Linking within- and between-host scales for understanding the evolutionary dynamics of quantitative antimicrobial resistance. J. Math. Biol. 87, 78 (2023). https://doi.org/10.1007/s00285-023-02008-1
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00285-023-02008-1