Treatment outcome in an SI model with evolutionary resistance: a Darwinian model for the evolution of resistance

Abstract

We consider a Darwinian (evolutionary game theoretic) version of a standard susceptible-infectious SI model in which the resistance of the disease causing pathogen to a treatment that prevents death to infected individuals is subject to evolutionary adaptation. We determine the existence and stability of all equilibria, both disease-free and endemic, and use the results to determine conditions under which the treatment will succeed or fail. Of particular interest are conditions under which a successful treatment in the absence of resistance adaptation (i.e. one that leads to a stable disease-free equilibrium) will succeed or fail when pathogen resistance is adaptive. These conditions are determined by the relative breadths of treatment effectiveness and infection transmission rate distributions as functions of pathogen resistance.

Keywords:

• SI model
• Darwinian dynamics
• evolutionary adaptation
• pathogen resistance to treatment

AMS Subject Classifications:

• 92D30
• 92D15

1. Introduction

In the twenty-first century , the world had seen a large reduction in the incidence of communicable diseases and their associated mortality, primarily due to improved access to and development of novel treatments and preventive vaccines, as well as improving coverage of interventions targeting disease vectors [20]. However, the evolution of resistance to treatment in pathogens and vectors remains a major hurdle to disease control and elimination (see [19] and the World Health Organization (WHO) Global Malaria Programme at https://www.who.int/teams/global-malaria-programme/prevention/vector-control/insecticide-resistance. Increasing antimicrobial resistance has led to reduced efficacy of treatments, resulting in a higher disease burden and higher mortality from infectious diseases. Similarly, increasing insecticide resistance in disease vectors, especially mosquitoes, has reduced the effectiveness of preventive interventions such as insecticide-treated nets, resulting in higher incidences of vector-borne diseases and associated mortality.

There is a substantial body of literature in modelling resistance; however, most studies have focused on the transmission of resistant pathogens or vectors, not the evolution of resistance; subsequently they consider either the impact of interventions on the transmission of resistance or on the effect of resistance on intervention effectiveness [5,15]. The studies that have modelled the evolution of resistance have focused on discrete genotypes that only consider a limited number of phenotypes [3,14–16]. Our main goal here is to investigate the methodology of Darwinian dynamics (evolutionary game theory) [18] as a tool to study the effectiveness of treatment in the presence of a pathogen's ability to develop resistance through evolutionary principles. Unlike most applications of adaptive dynamics [6], we do not focus here on a limited number of pathogen mutants or strains or genotypes, utilize methods of optimizing $R_0,$ nor assume a separation of ecological and evolutionary time scales [7–9]. Instead we consider a pathogen with a continuously distributed phenotypic trait that confers resistance to a treatment and consider a model in which the host dynamics and the evolutionary dynamics of the mean trait are coupled. The model is based on the (so-called first order) methodology of evolutionary game theory in which the mean trait dynamics are governed by Lande's equation (sometimes referred to as Fisher's equation or the canonical equation of evolution) [18]. In this paper we focus only on the existence and stability properties of disease-free and endemic equilibria and how they depend on the way treatment resistance and pathogen infectiousness depend on the evolving mean trait. To do this we extend a basic compartmental susceptible-infectious-recovered model to better understand the evolution of resistance to treatment and the impact of this evolution on treatment effectiveness. We assume in the model that

1. there is a constant recruitment/birth rate into the population;

2. the force of infection follows a mass action law;

3. infection leads to an additional disease-induced death rate; and

4. recovery from infection leads to lifelong immunity.

Since the recovered population does not influence the dynamics of the susceptible or infectious populations, we do not explicitly model the recovered class, and henceforth refer to the model as an SI model.

We consider a treatment that is fully effective in preventing disease-induced mortality and that further reduces the duration of the infectious period, in the absence of resistance. The evolution of resistance increases the fitness of the pathogen by reducing the effectiveness of the treatment in clearing the infection and thereby increasing the duration of infection. However, it also reduces the fitness of the pathogen since reducing treatment effectiveness leads to an increase in host death due to increased disease-induced mortality (and consequently leads to the loss of the pathogen). We furthermore assume that resistance has a cost to the pathogen in onward transmissibility from one host to another, leading to reduced infectivity of individuals infected with resistant pathogens.

We focus on a situation where the disease causing pathogen is controlled (asymptotically eliminated) by a certain level of a treatment in the absence of the pathogen's ability to evolve resistance to the treatment. Using the methodology of evolutionary game theory [18], we derive a Darwinian SI model and using it to address the following questions. Under what conditions will this treatment level still work in eliminating the pathogen if it is capable of evolving resistance to the treatment? If this treatment level fails due to the pathogen's adaptation, will a higher treatment level succeed? Are there conditions under which no level of treatment will successfully eliminate an evolutionarily adapting pathogen?

We first describe the SI model with treatment in the absence of pathogen resistance, including the definition of the basic reproductive number, $R_0$ and the existence and stability of equilibrium solutions. We then derive a Darwinian version of this model that describes the evolution of the fitness of the pathogen in response to treatment. We analytically investigate the existence and stability of disease-free and endemic equilibria. We discuss the criteria for the existence and stability in terms of two important bifurcation parameters: the proportion of people treated and the ratio of the variance of treatment effectiveness to the variance of the infection transmission rate. We illustrate these results with time series plots, phase portrait plots and bifurcation diagrams. We summarize these results with a parameter map showing the criteria for the existence and stability of the equilibria.

The basic SI model we use in this study is described, as is commonly done, by ordinary differential equations. The method of Darwinian dynamics which we apply to the model is also available for discrete time (difference equation) models, which are also often used for modelling epidemics (for example, by A.-A. Yakubu, to whom this paper is dedicated).

2. An SI model with resistance to treatment

A basic disease model with two classes of individuals, susceptible S and infectious I, is [4] (1) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-cSI\\ I^{\mathrm{\prime }}& =cSI-\rho I.\end{array}$(1) Here $\mathrm{\Lambda }>0$ is the constant rate at which susceptible individuals enter the population and $\mu >0$ is the (per capita) natural death rate. The coefficient c>0 measures the transmission rate at which infectious individuals infect susceptible individuals (assumed to follow a mass action law) and cause them to enter the infectious class. Finally $\rho >0$ is the (per capita) rate at which infectious individuals are removed from the population and, because complete immunity is assumed, do not re-enter the susceptible class. There exist two equilibria, a disease-free equilibrium $\begin{array}{r}DFE:\left(\begin{array}{c}S\\ I\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\\ 0\end{array}\right)\end{array}$and, provided the basic reproduction number $\begin{array}{r}{R}_0=\frac{c\mathrm{\Lambda }}{\mu }\frac{1}{\rho }\end{array}$satisfies $R_0>1$, an endemic equilibrium $\begin{array}{r}EE:\left(\begin{array}{c}S\\ I\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{\mathrm{\Lambda }}{\mu R_0}}\\ {\displaystyle \frac{\mu }{c}(R_0-1)}\end{array}\right).\end{array}$A routine linearization analysis shows DFE is stable if $R_0<1$ and is unstable if $R_0>1$ when EE is stable. Thus, a standard transcritical bifurcation and exchange of stability occurs between these equilibria at $R_0=1$. In fact, it is shown in [2] using Liapunov functions that DFE is globally asymptotically stable on $R_{+}^2$ when $R_0\le 1$ and EE is globally asymptotically stable on $R_{+}^2\mathrm{\setminus }\{(S,0),S\ge 0\}$ when $R_0>1$.

We parse the infectious removal rate ρ into the sum of three rates: the natural death rate (which we assume equal to $\mu >0$, the natural death rate of susceptible individuals), the death rate due to the disease $\kappa >0$, and removal rate due to recovery from the infection $\varphi \ge 0$. We are interested in the effect of a treatment for the disease, which we assume is successfully administered to a fraction τ of the infectious class, $0\le \tau \le 1$. This treatment introduces an additional recovery rate $\gamma >0$ to the $\tau I$ individuals treated and prevents death due to the disease. Adding together all the removals from the infectious class (see Table 1), we get $\begin{array}{r}\rho I=\mu \tau I+(\varphi +\gamma )\tau I+\mu (1-\tau )I+\kappa (1-\tau )I+\varphi (1-\tau )I\end{array}$which gives us the model equations $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-cSI\\ I^{\mathrm{\prime }}& =[cS-(\mu +\kappa +\varphi )-\tau (\gamma -\kappa )]I.\end{array}$

Table 1. Removal rates from the infectious class.

Rates	Treated	Untreated
Natural death	$\mu \tau I$	$\mu (1-\tau )I$
Death by disease	0	$\kappa (1-\tau )I$
Recovery	$(\varphi +\gamma )\tau I$	$\varphi (1-\tau )I$

We consider the reproduction number as a function of the treatment fraction τ : $\begin{array}{r}{R}_0\left(\tau \right)=c\frac{\mathrm{\Lambda }}{\mu }\frac{1}{(\mu +\kappa +\varphi )+\tau (\gamma -\kappa )}.\end{array}$We are interested here in the scenario when the disease is endemic in the absence of the treatment, but is eliminated if the treatment reaches a sufficient fraction of the infectious class. To this end, we assume $\begin{array}{r}{R}_0\left(1\right)<1<R_0\left(0\right),\end{array}$inequalities which we can re-write as (2) $\begin{array}{r}\mu +\kappa +\varphi <c\frac{\mathrm{\Lambda }}{\mu }<\gamma +\mu +\varphi .\end{array}$(2) Note that these constraints imply $\gamma >\kappa ,$ namely that the additional recovery rate provided by the treatment exceeds the death rate due to the disease. The treatment is successful in preventing an epidemic if $R_0(\tau )<1$, i.e. if (3) $\begin{array}{r}\tau >{\tau }_0:=\frac{c\frac{\mathrm{\Lambda }}{\mu }-(\mu +\kappa +\varphi )}{\gamma -\kappa }\end{array}$(3) and fails if $\tau <{\tau }_0$.

3. A Darwinian SI model

To derive the fitness of the pathogen for the Darwinian SI model, we re-interpret I as the number of infections in the population, where each individual can have exactly zero or one infection. The individuals with no infection are susceptible, while the individuals with one infection are infectious. The derivation of the model in Section 2 remains exactly the same with this new interpretation, but c now needs to have the appropriate units to convert between individuals and infections. The number of infections increases when susceptible individuals get infected and decreases when infected individuals recover or die (because the death of the host also implies the death of the infection).

We suppose that the disease is caused by a pathogen that can develop resistance to the treatment by Darwinian processes. We hypothesize that resistance is conferred by a (continuous) phenotypic trait v that is subject to natural selection and that is normally distributed throughout the population at all time, with variance ${\sigma }_v^2$ and mean u. Let $\begin{array}{r}p(v-u)=\frac{1}{{\sigma }_v\sqrt{2\pi }}\exp (-\frac{{(v-u)}^2}{2{\sigma }_v^2})\end{array}$denote the probability distribution for v so that there are $p(v-u)I$ infections with trait v when the mean population trait is u. The fraction of these infections that are successfully treated depends on the resistance of the pathogen and $\tau (v)$ is a function v.

We assume there is a trait v at which the pathogen has no resistance to the treatment. We denote this reference point by v=0 and assume that the treatment $\tau (v)$ is a decreasing function of increasing pathogen resistance $v>0.$

In the model considered here, it is assumed that increased resistance to the treatment comes at a cost to the pathogen, namely at the cost of decreased infectiousness. Thus, the infection transmission rate c in model (1(1) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-cSI\\ I^{\mathrm{\prime }}& =cSI-\rho I.\end{array}$(1) ) becomes a decreasing function $c(v)$ of $v>0$.

Consider the sub-class $p(v-u)I$ of infections with resistance v. The fraction of these infections successfully treated is $\tau (v).$ The removal rates from this sub-class are shown in Table 2. The total removal rate is the sum of the entries in Table 2, which equals $\begin{array}{r}[(\mu +\kappa +\varphi )+(\gamma -\kappa )\tau \left(v\right)]Ip(v-u).\end{array}$The removal rate from the susceptible class due to infections with trait v and the entry rate into this infectious sub-class is $c(v)SIp(v-u)$. Thus, the change in this infectious class is $r(v,u,S)I$ where (4) $\begin{array}{r}r(v,u,S):=c\left(v\right)Sp(v-u)-[(\mu +\kappa +\varphi )+(\gamma -\kappa )\tau \left(v\right)]p(v-u).\end{array}$(4) We take this per capita rate as fitness of the pathogen with resistance v.

Table 2. Removal rates of infections with trait v.

Rates	Treated	Untreated
Natural death	$\mu \tau (v)Ip(v-u)$	$\mu (1-\tau (v))Ip(v-u)$
Death by disease	0	$\kappa (1-\tau (v))Ip(v-u)$
Recovery	$(\varphi +\gamma )\tau (v)Ip(v-u)$	$\varphi (1-\tau (v))Ip(v-u)$

The methodology of Darwinian dynamics provides a model for the dynamics of S, I, and the mean trait u [1,10,11,13,18] : $\begin{array}{rl}{S}^{\mathrm{\prime }}& ={\mathrm{\Lambda }-\mu S-c(v)SIp(v-u)|}_{v=u}\\ I^{\mathrm{\prime }}& ={r(v,u,S)|}_{v=u}I\\ u^{\mathrm{\prime }}& =\theta {\frac{\mathrm{\partial }r(v,u,S)}{\mathrm{\partial }v}|}_{v=u}\end{array}$where $\theta >0$ is called the speed of evolution. Since $p(0)=1$ and $p^{\mathrm{\prime }}(0)=0$ the expression (4(4) $\begin{array}{r}r(v,u,S):=c\left(v\right)Sp(v-u)-[(\mu +\kappa +\varphi )+(\gamma -\kappa )\tau \left(v\right)]p(v-u).\end{array}$(4) ) yields the model equations $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\left(u\right)SI\\ I^{\mathrm{\prime }}& =(c\left(u\right)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \left(u\right))I\\ u^{\mathrm{\prime }}& =\theta (c^{\mathrm{\prime }}\left(u\right)S-(\gamma -\kappa ){\tau }^{\mathrm{\prime }}\left(u\right)).\end{array}$(The primes denote differentiation with respect to a function's argument.)

In this paper, we assume both the treatment effectiveness and the infection transmission rates have Gaussian distributions as functions of the trait v. Thus, we replace $\tau (v)$ and $c(v)$, respectively, by (5) $\begin{array}{r}\tau \exp (-\frac{1}{2{\sigma }^2}{v}^2),0\le \tau \le 1\\ c\exp (-\frac{1}{2}{v}^2),c>0.\end{array}$(5) where τ and c now denote positive constants that denote the treatment effectiveness and infection transmission rates for an pathogen with no resistance to the treatment (v=0). The coefficient ${\sigma }^2$ is the variance of the treatment effectiveness relative to the variance of the infection transmission rate (which we assume, by re-scaling the units for v if necessary, to be equal to 1). The resulting Darwinian SI model equations are (6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) The two parameters τ and ${\sigma }^2$ that define the distributions appearing in (5(5) $\begin{array}{r}\tau \exp (-\frac{1}{2{\sigma }^2}{v}^2),0\le \tau \le 1\\ c\exp (-\frac{1}{2}{v}^2),c>0.\end{array}$(5) ) play a central role in our analysis and findings below. τ is the maximal effectiveness of the treatment (which occurs for those infections with trait v=0, i.e. with no resistance to the treatment). The variance ${\sigma }^2$ measures how widely spread the treatment is effective; a larger value of ${\sigma }^2$ indicates that a larger distribution of infections are successfully treated.

Since the equations for I and u in (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) have factors I and u on their right sides and because $\mathrm{\Lambda }>0$ (so that $S^{\mathrm{\prime }}>0$ when $S=0)$ we obtain the following lemma.

Lemma 3.1

For $\begin{array}{r}\left(\begin{array}{c}S\left(0\right)\\ I\left(0\right)\\ u\left(0\right)\end{array}\right)\in R_{+}^3,\end{array}$the solution of (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) satisfies $\begin{array}{r}\left(\begin{array}{c}S\left(t\right)\\ I\left(t\right)\\ u\left(t\right)\end{array}\right)\in R_{+}^3\end{array}$for all $t\ge 0$. Moreover for all t>0

(a)	$S(0)\ge 0$ implies $S(t)>0;$
(b)	$I(0)>0$ implies $I(t)>0$ and $I(0)=0$ implies $I(t)=0;$
(c)	$u(0)>0$ implies $u(t)>0$ and $u(0)=0$ implies $u(t)=0.$

As a result of this lemma, we restrict our analysis of the dynamics of (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) (with a focus on the existence and stability of equilibria) to $R_{+}^3$. In particular, we consider only solutions with non-negative trait components $u(t)$. Note that a change in variables from u to $-u$ leaves the model equations (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) unchanged and results for $u(t)\ge 0$ yield symmetric results for $u(t)\le 0$.

From the right sides of (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ), with u fixed, we calculate the reproduction number $\begin{array}{r}{R}_0(u,\tau )=c\exp (-\frac{1}{2}{u}^2)\frac{\mathrm{\Lambda }}{\mu }\frac{1}{(\mu +\kappa +\varphi )+\tau (\gamma -\kappa )\exp (-\frac{1}{2{\sigma }^2}{u}^2)}.\end{array}$

4. Equilibrium analysis of the Darwinian SI model

The equilibrium equations of (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) are (7a) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\end{array}$(7a) (7b) $\begin{array}{rl}0& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\end{array}$(7b) (7c) $\begin{array}{rl}0& =u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(7c) In the next section we see that these equations are analytically tractable and that there exist two types of disease-free and two types of endemic equilibria.

4.1. Existence

A disease-free equilibrium $\begin{array}{r}\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right)\end{array}$is a solution of the equilibrium equations (7(7a) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\end{array}$(7a) ) with $I_e=0$ and $S_e>0$ and $u_e$ that satisfy the equations (8) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S_e\end{array}$(8) (9) $\begin{array}{rl}0& =u_e(-c\exp (-\frac{1}{2}{u}_e^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}_e^2)).\end{array}$(9) Equation (8(8) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S_e\end{array}$(8) ) implies $S_e=\mathrm{\Lambda }/\mu$ and equation (9(9) $\begin{array}{rl}0& =u_e(-c\exp (-\frac{1}{2}{u}_e^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}_e^2)).\end{array}$(9) ) implies that either $u_e=0$ or $u_e>0$ is a root of the equation (10) $\begin{array}{r}-c\exp (-\frac{1}{2}{u}_e^2)\frac{\mathrm{\Lambda }}{\mu }+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}_e^2)=0.\end{array}$(10) This results in two disease-free equilibria in $R_{+}^3:$ the non-adaptive disease-free equilibrium $\begin{array}{r}DF{E}_0:\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right)=\left(\begin{array}{c}\frac{\mathrm{\Lambda }}{\mu }\\ 0\\ 0\end{array}\right)\end{array}$and the equilibrium $\begin{array}{r}DF{E}_{+}:\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\\ 0\\ {\displaystyle \sqrt{2\frac{{\sigma }^2}{1-{\sigma }^2}\ln \left(\tau \frac{{\sigma }_1^2}{{\sigma }^2}\right)}}\end{array}\right),{\sigma }_1^2:=\frac{\mu }{\mathrm{\Lambda }c}(\gamma -\kappa )\end{array}$with a positive equilibrium trait component $u_e$ in either of two cases depending on the constraints $\begin{array}{rl}(\mathrm{a})0& <{\sigma }^2<1\mathrm{and}\tau >\frac{{\sigma }^2}{{\sigma }_1^2}\\ (\mathrm{b})1& <{\sigma }^2\mathrm{and}\tau <\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$for the two distributional characteristics ${\sigma }^2$ and τ of the treatment effectiveness (5(5) $\begin{array}{r}\tau \exp (-\frac{1}{2{\sigma }^2}{v}^2),0\le \tau \le 1\\ c\exp (-\frac{1}{2}{v}^2),c>0.\end{array}$(5) ). Note that $DF{E}_{+}=DF{E}_0$ when $\tau ={\sigma }^2/{\sigma }_1^2$ where these two disease-free equilibria mathematically undergo a pitchfork bifurcation (taking into consideration the additional $DF{E}_{+}$ with $u_e$ replaced by $-u_e$).

To find endemic equilibria we need to solve the equilibrium equations (7a(7a) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\end{array}$(7a) ) for $I_e>0$ and the equilibrium equations become (after a cancellation of I) (11) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ 0& =c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2)\\ 0& =(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2))u.\end{array}$(11) The third equation implies either $u_e=0$ or $u_e>0$ is a root of the parenthetical factor on the right side of the equation.

In the first case we can solve the second equation for S after which we can solve the first equation for I. The result is the non-adaptive endemic equilibrium $\begin{array}{r}E{E}_0:\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{\mathrm{\Lambda }}{\mu }\frac{1}{R_0(0,\tau )}}\\ {\displaystyle \frac{\mu }{c}(R_0(0,\tau )-1)}\\ 0\end{array}\right)\end{array}$in $R_{+}^3$ provided $\begin{array}{r}{R}_0(0,\tau )=c\frac{\mathrm{\Lambda }}{\mu }\frac{1}{(\mu +\kappa +\varphi )+\tau (\gamma -\kappa )}>1\end{array}$i.e. provided $\begin{array}{r}\tau <{\tau }_0.\end{array}$In the second case the equilibrium equations reduce (after a cancellation of u) to the system of equations $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ 0& =c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2)\\ 0& =-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2).\end{array}$By adding the last two equations we get an equation for u alone (with S and I absent) which we can solve for (12) $\begin{array}{r}{u}_e=\sqrt{2{\sigma }^2\ln \left(\frac{1-{\sigma }^2}{{\sigma }^2}\frac{\gamma -\kappa }{\mu +\kappa +\varphi }\tau \right)}\end{array}$(12) after which we can solve the third and first equations (in that order) for (13) $\begin{array}{r}{S}_e=\frac{1}{c}(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1-{\sigma }^2}{2{\sigma }^2}{u}_e^2)\end{array}$(13) and (14) $\begin{array}{r}{I}_e=\frac{1}{c}\frac{\mathrm{\Lambda }-\mu S_e}{S_e}\exp \left(\frac{1}{2}{u}_e^2\right)\end{array}$(14) to obtain the adaptive endemic equilibrium $\begin{array}{r}E{E}_{+}:\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right),\end{array}$provided this vector lies in $R_{+}^3.$ While it is clear that $S_e>0,$ it takes some manipulations with inequalities to show that $I_e>0$ and $E{E}_{+}$ lies in $R_{+}^3$ if and only if τ lies in the interval $\begin{array}{r}\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}<\tau <\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}{\left((1-{\sigma }^2)\frac{\mathrm{\Lambda }c}{\mu }\frac{1}{\mu +\kappa +\varphi }\right)}^{\frac{1}{{\sigma }^2}}.\end{array}$Fixing all other parameters, we can view the various constraints required for the existence of the two disease-free equilibria $DF{E}_0,$ $DF{E}_{+}$ and the two endemic equilibria $E{E}_0,$ and $E{E}_{+}$ in terms of the two parameters ${\sigma }^2$ and τ characterizing the effectiveness of the treatment as a function of pathogen resistance, as given by the assumed distribution (5(5) $\begin{array}{r}\tau \exp (-\frac{1}{2{\sigma }^2}{v}^2),0\le \tau \le 1\\ c\exp (-\frac{1}{2}{v}^2),c>0.\end{array}$(5) ). The relevant region in the $({\sigma }^2,\tau )$-plane, namely the infinite rectangle ${\sigma }^2\ge 0$ and $0\le \tau \le 1$, is divided in to a number of sub-regions in each of which there exist a specific set of equilibria. These regions are shown in the parameter maps in Figure 1, which also include information about the (local asymptotic) stability of the equilibria, a topic we take up in the next section.

Figure 1. Summary map in the bifurcation parameter space showing the criteria for the existence and stability of the disease-free and endemic equilibria. The shaded regions correspond to disease-free asymptotic dynamics and the unshaded regions correspond to endemic asymptotic dynamics.

4.2. Stability

The Jacobian associated with the equations (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) is (15) $\begin{array}{r}\left(\begin{array}{ccc}{\displaystyle -(\mu +cI\exp (-\frac{u^2}{2}))}& {\displaystyle -cS\exp (-\frac{u^2}{2})}& {\displaystyle cuSI\exp (-\frac{u^2}{2})}\\ {\displaystyle cI\exp (-\frac{u^2}{2})}& J_{22}(S,u)& J_{23}(S,I,u)\\ {\displaystyle -c\theta u\exp (-\frac{u^2}{2})}& 0& J_{33}(S,u)\end{array}\right)\end{array}$(15) where $\begin{array}{rl}{J}_{22}(S,u)& :=-(\mu +\kappa +\varphi )+cS\exp (-\frac{u^2}{2})-(\gamma -\kappa )\tau \exp (-\frac{u^2}{2{\sigma }^2})\\ J_{23}(S,I,u)& :=\frac{1}{{\sigma }^2}u[(\gamma -\kappa )\tau \exp (-\frac{u^2}{2{\sigma }^2})-c{\sigma }^{2}S\exp (-\frac{u^2}{2})]I\\ J_{33}(S,u)& :=\frac{\theta }{{\sigma }^4}[\tau (\gamma -\kappa )({\sigma }^2-u^2)\exp (-\frac{u^2}{2{\sigma }^2})+c{\sigma }^{4}S(u^2-1)\exp (-\frac{u^2}{2})]\end{array}$To apply the Linearization Principle to the disease-free equilibrium $DF{E}_0$ we evaluate the Jacobian (15(15) $\begin{array}{r}\left(\begin{array}{ccc}{\displaystyle -(\mu +cI\exp (-\frac{u^2}{2}))}& {\displaystyle -cS\exp (-\frac{u^2}{2})}& {\displaystyle cuSI\exp (-\frac{u^2}{2})}\\ {\displaystyle cI\exp (-\frac{u^2}{2})}& J_{22}(S,u)& J_{23}(S,I,u)\\ {\displaystyle -c\theta u\exp (-\frac{u^2}{2})}& 0& J_{33}(S,u)\end{array}\right)\end{array}$(15) ) at $DF{E}_0$ and obtain $\begin{array}{r}\left(\begin{array}{ccc}-\mu & {\displaystyle -c\frac{\mathrm{\Lambda }}{\mu }}& 0\\ 0& (\gamma -\kappa )({\tau }_0-\tau )& 0\\ 0& 0& {\displaystyle \theta \frac{\gamma -\kappa }{{\sigma }^2}(\tau -\frac{{\sigma }^2}{{\sigma }_1^2})}\end{array}\right)\end{array}$whose eigenvalues are the diagonal entries $\begin{array}{r}{\lambda }_1=-\mu <0,{\lambda }_2=(\gamma -\kappa )({\tau }_0-\tau ),{\lambda }_3=\theta \frac{\gamma -\kappa }{{\sigma }^2}(\tau -\frac{{\sigma }^2}{{\sigma }_1^2}).\end{array}$Recalling that the assumptions (2(2) $\begin{array}{r}\mu +\kappa +\varphi <c\frac{\mathrm{\Lambda }}{\mu }<\gamma +\mu +\varphi .\end{array}$(2) ) imply $\gamma >\kappa ,$ we explore the signs of ${\lambda }_2$ and ${\lambda }_3$ and obtain the following result.

Theorem 4.1

The disease-free equilibrium $DF{E}_0$ is stable if $\begin{array}{r}\tau >{\tau }_{0}and\tau <\frac{{\sigma }^2}{{\sigma }_1^2}.\end{array}$If either inequality is reversed, $DF{E}_0$ is unstable.

To apply the Linearization Principle to the disease-free equilibrium $DF{E}_{+}$ we calculate the Jacobian (15(15) $\begin{array}{r}\left(\begin{array}{ccc}{\displaystyle -(\mu +cI\exp (-\frac{u^2}{2}))}& {\displaystyle -cS\exp (-\frac{u^2}{2})}& {\displaystyle cuSI\exp (-\frac{u^2}{2})}\\ {\displaystyle cI\exp (-\frac{u^2}{2})}& J_{22}(S,u)& J_{23}(S,I,u)\\ {\displaystyle -c\theta u\exp (-\frac{u^2}{2})}& 0& J_{33}(S,u)\end{array}\right)\end{array}$(15) ) at $DF{E}_{+}$ and obtain $\begin{array}{r}\left(\begin{array}{ccc}-\mu & {\displaystyle -\frac{\mathrm{\Lambda }}{\mu }c{e}^{-\frac{1}{2}{u}_e^2}}& 0\\ 0& J_{22}{\displaystyle (\frac{\mathrm{\Lambda }}{\mu },u_e)}& 0\\ -cu\theta {\displaystyle e^{-\frac{1}{2}{u}_e^2}}& 0& J_{33}{\displaystyle (\frac{\mathrm{\Lambda }}{\mu },u_e)}\end{array}\right)\end{array}$where $\begin{array}{r}{u}_e=\sqrt{2\frac{{\sigma }^2}{1-{\sigma }^2}\ln \left(\tau \frac{{\sigma }_1^2}{{\sigma }^2}\right)}.\end{array}$The eigenvalues of this matrix are the diagonal entries: $\begin{array}{r}{\lambda }_1=-\mu <0,{\lambda }_2=J_{22}(\frac{\mathrm{\Lambda }}{\mu },u_e),{\lambda }_3=J_{33}(\frac{\mathrm{\Lambda }}{\mu },u_e).\end{array}$Since $u_e$ satisfies Equation (10(10) $\begin{array}{r}-c\exp (-\frac{1}{2}{u}_e^2)\frac{\mathrm{\Lambda }}{\mu }+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}_e^2)=0.\end{array}$(10) ) and hence the equation $\begin{array}{r}\frac{1}{{\sigma }^2}(\gamma -\kappa )\tau \exp (-\frac{1}{2}\frac{u_e^2}{{\sigma }^2})=c\frac{\mathrm{\Lambda }}{\mu }\exp (-\frac{1}{2}{u}_e^2)\end{array}$we obtain $\begin{array}{rl}{\lambda }_2& =-(\mu +\kappa +\varphi )-c({\sigma }^2-1)\frac{\mathrm{\Lambda }}{\mu }\exp (-\frac{1}{2}{u}_e^2)\\ {\lambda }_3& =\theta \frac{{\sigma }^2-1}{{\sigma }^2}c\frac{\mathrm{\Lambda }}{\mu }{u}_e^2\exp (-\frac{1}{2}{u}_e^2).\end{array}$If ${\sigma }^2>1$ then ${\lambda }_3>0$ and $DF{E}_{+}$ is unstable. If ${\sigma }^2<1$ then stability by linearization is determined by the sign of ${\lambda }_2$ which we treat as a function of τ and obtain the following result.

Theorem 4.2

The disease-free equilibrium $DF{E}_{+}$ is stable if $\begin{array}{r}{\sigma }^2<1\mathrm{and}\tau >\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}{\left((1-{\sigma }^2)\frac{\mathrm{\Lambda }c}{\mu }\frac{1}{\mu +\kappa +\varphi }\right)}^{\frac{1}{{\sigma }^2}}\end{array}$and unstable if either inequality is reverse.

To apply the Linearization Principle to the disease-free equilibrium $E{E}_0$ (which, as we saw above, is positive if and only if $\tau <{\tau }_0$) we calculate the Jacobian (15(15) $\begin{array}{r}\left(\begin{array}{ccc}{\displaystyle -(\mu +cI\exp (-\frac{u^2}{2}))}& {\displaystyle -cS\exp (-\frac{u^2}{2})}& {\displaystyle cuSI\exp (-\frac{u^2}{2})}\\ {\displaystyle cI\exp (-\frac{u^2}{2})}& J_{22}(S,u)& J_{23}(S,I,u)\\ {\displaystyle -c\theta u\exp (-\frac{u^2}{2})}& 0& J_{33}(S,u)\end{array}\right)\end{array}$(15) ) at $E{E}_0$ and obtain $\begin{array}{r}\left(\begin{array}{ccc}{\displaystyle -\frac{\mathrm{\Lambda }c}{(\mu +\kappa +\varphi )+(\gamma -\kappa )\tau }}& -(\mu +\varphi +\kappa +(\gamma -\kappa )\tau )& 0\\ {\displaystyle -\mu \frac{\mu +\kappa +\varphi +(\gamma -\kappa )\tau -\frac{c\mathrm{\Lambda }}{\mu }}{\mu +\kappa +\varphi +(\gamma -\kappa )\tau }}& 0& 0\\ 0& 0& {\displaystyle \theta (\gamma -\kappa )\frac{1-{\sigma }^2}{{\sigma }^2}}\\ & & (\tau -\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2})\end{array}\right)\end{array}$which is block diagonal. The eigenvalues are the two eigenvalues ${\lambda }_1$ and ${\lambda }_2$ of the $2\times 2$ submatrix $\begin{array}{r}\left(\begin{array}{cc}{\displaystyle -\frac{\mathrm{\Lambda }c}{(\mu +\kappa +\varphi )+(\gamma -\kappa )\tau }}& -(\mu +\varphi +\kappa +(\gamma -\kappa )\tau )\\ {\displaystyle -\mu \frac{\mu +\kappa +\varphi +(\gamma -\kappa )\tau -\frac{c\mathrm{\Lambda }}{\mu }}{\mu +\kappa +\varphi +(\gamma -\kappa )\tau }}& 0\end{array}\right)\end{array}$and the diagonal term $\begin{array}{r}{\lambda }_3=\theta (\gamma -\kappa )\frac{1-{\sigma }^2}{{\sigma }^2}(\tau -\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}).\end{array}$Once again recall that (2(2) $\begin{array}{r}\mu +\kappa +\varphi <c\frac{\mathrm{\Lambda }}{\mu }<\gamma +\mu +\varphi .\end{array}$(2) ) imply $\gamma >\kappa .$ The trace-determinant criteria sufficient to guarantee that the eigenvalues of a $2\times 2$ matrix lie in the left half complex plane are that the trace be negative, which it clearly is, and that the determinant $\begin{array}{r}-\mu (\mu +\kappa +\varphi +(\gamma -\kappa )\tau -c\frac{\mathrm{\Lambda }}{\mu })=\mu (\gamma -\kappa )({\tau }_0-\tau )\end{array}$be positive, which it is since $\tau <{\tau }_0$. Thus, the stability by linearization of $E{E}_0$ is determined by the sign of ${\lambda }_3.$ This leads to the the following result.

Theorem 4.3

The endemic equilibrium $E{E}_0$ is stable if ${\sigma }^2>1$ or if $\begin{array}{r}{\sigma }^2<1\mathrm{and}\tau <\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}.\end{array}$It is unstable if $\begin{array}{r}{\sigma }^2<1\mathrm{and}\tau >\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}.\end{array}$

Finally we consider the endemic equilibrium $E{E}_{+}.$ In the second equilibrium equation (11(11) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ 0& =c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2)\\ 0& =(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2))u.\end{array}$(11) ), we see that $\begin{array}{r}(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}_e^2)=c\exp (-\frac{1}{2}{u}_e^2)S-(\mu +\kappa +\varphi )\end{array}$and hence $J_{22}(S_e,u_e)=0$. Using this fact, together with $\begin{array}{r}c\exp (-\frac{u_e^2}{2})=\frac{\mathrm{\Lambda }-\mu S_e}{S_e{I}_e}\end{array}$which follows from the first equilibrium equation (11(11) $\begin{array}{rl}0& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ 0& =c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2)\\ 0& =(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2))u.\end{array}$(11) ), we find that the Jacobian at $E{E}_{+}$ is (16) $\begin{array}{r}\left(\begin{array}{ccc}-(\mu +\frac{\mathrm{\Lambda }-\mu S_e}{S_e})& {\displaystyle -\frac{\mathrm{\Lambda }-\mu S_e}{I_e}}& u_e(\mathrm{\Lambda }-\mu S_e)\\ {\displaystyle \frac{\mathrm{\Lambda }-\mu S_e}{S_e}}& 0& J_{23}(S_e,I_e,u_e)\\ {\displaystyle -\theta u\frac{\mathrm{\Lambda }-\mu S_e}{S_e{I}_e}}& 0& J_{33}(S_e,u_e)\end{array}\right).\end{array}$(16) As shown by the parameter maps in Figure 1, the endemic equilibrium $E{E}_{+}$ exists only in restricted regions of the $({\sigma }^2,\tau )$-plane that require, among other constraints, that ${\sigma }^2$ be not too large, namely, that ${\sigma }^2<{\sigma }_0^2$. Although we do not have a general analysis of this Jacobian for all parameter values for which $E{E}_{+}$ exists, we do have the following result for small ${\sigma }^2.$

Theorem 4.4

For ${\sigma }^2$ sufficiently small, the endemic equilibrium $E{E}_{+}$ is locally asymptotically stable.

Proof.

From $\begin{array}{r}\underset{{\sigma }^2\to 0}{lim}{\sigma }^2\ln \left(\frac{1-{\sigma }^2}{{\sigma }^2}\right)=0\end{array}$and $\begin{array}{r}{u}_e^2=2{\sigma }^2\ln \left(\frac{\frac{1-{\sigma }^2}{{\sigma }^2}(\gamma -\kappa )\tau }{\mu +\kappa +\varphi }\right)=2{\sigma }^2\ln \left(\frac{1-{\sigma }^2}{{\sigma }^2}\right)+2{\sigma }^2\ln \left(\frac{(\gamma -\kappa )\tau }{\mu +\kappa +\varphi }\right)\end{array}$we get that $\begin{array}{r}\underset{{\sigma }^2\to 0}{lim}{u}_e=0.\end{array}$From the formulas (12(12) $\begin{array}{r}{u}_e=\sqrt{2{\sigma }^2\ln \left(\frac{1-{\sigma }^2}{{\sigma }^2}\frac{\gamma -\kappa }{\mu +\kappa +\varphi }\tau \right)}\end{array}$(12) ) and (13(13) $\begin{array}{r}{S}_e=\frac{1}{c}(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1-{\sigma }^2}{2{\sigma }^2}{u}_e^2)\end{array}$(13) ) we have, after some algebraic manipulations, that $\begin{array}{r}{S}_e=\frac{1}{c}(\gamma -\kappa )\tau {\left(\frac{\mu +\kappa +\varphi }{(\gamma -\kappa )\tau }\right)}^{1-{\sigma }^2}\frac{1}{{\sigma }^2}{\left(\frac{{\sigma }^2}{1-{\sigma }^2}\right)}^{1-{\sigma }^2}\end{array}$and since $\begin{array}{r}\underset{\sigma \to 0}{lim}\frac{1}{{\sigma }^2}{\left(\frac{{\sigma }^2}{1-{\sigma }^2}\right)}^{1-{\sigma }^2}=1\end{array}$we find that $\begin{array}{r}\underset{{\sigma }^2\to 0}{lim}{S}_e=\frac{1}{c}(\mu +\kappa +\varphi ).\end{array}$Finally from these two limits and the formula (14(14) $\begin{array}{r}{I}_e=\frac{1}{c}\frac{\mathrm{\Lambda }-\mu S_e}{S_e}\exp \left(\frac{1}{2}{u}_e^2\right)\end{array}$(14) ) we get $\begin{array}{rl}\underset{{\sigma }^2\to 0}{lim}{I}_e& =\underset{{\sigma }^2\to 0}{lim}(\frac{1}{c}\frac{\mathrm{\Lambda }-\mu S_e}{S_e}\exp \left(\frac{1}{2}{u}_e^2\right))\\ & =\frac{1}{c}\frac{\mu }{\mu +\kappa +\varphi }(\frac{c\mathrm{\Lambda }}{\mu }-(\mu +\kappa +\varphi )).\end{array}$Thus, (17) $\begin{array}{r}\underset{{\sigma }^2\to 0}{lim}\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{c}(\mu +\kappa +\varphi )}\\ {\displaystyle \frac{1}{c}\frac{\mu }{\mu +\kappa +\varphi }(\frac{c\mathrm{\Lambda }}{\mu }-(\mu +\kappa +\varphi ))}\\ 0\end{array}\right).\end{array}$(17) The first column, third row entry in the Jacobian (16(16) $\begin{array}{r}\left(\begin{array}{ccc}-(\mu +\frac{\mathrm{\Lambda }-\mu S_e}{S_e})& {\displaystyle -\frac{\mathrm{\Lambda }-\mu S_e}{I_e}}& u_e(\mathrm{\Lambda }-\mu S_e)\\ {\displaystyle \frac{\mathrm{\Lambda }-\mu S_e}{S_e}}& 0& J_{23}(S_e,I_e,u_e)\\ {\displaystyle -\theta u\frac{\mathrm{\Lambda }-\mu S_e}{S_e{I}_e}}& 0& J_{33}(S_e,u_e)\end{array}\right).\end{array}$(16) ) tends to 0 as ${\sigma }^2\to 0$. Therefore, its three eigenvalues are, for small ${\sigma }^2$, approximately equal to the two eigenvalues ${\lambda }_1({\sigma }^2)$ and ${\lambda }_2({\sigma }^2)$ of the upper left $2\times 2$ matrix (18) $\begin{array}{r}\left(\begin{array}{cc}{\displaystyle -(\mu +\frac{\mathrm{\Lambda }-\mu S_e}{S_e})}& -S_{e}c{e}^{-\frac{1}{2}{u}_e^2}\\ {\displaystyle \frac{\mathrm{\Lambda }-\mu S_e}{S_e}}& 0\end{array}\right)\end{array}$(18) and the real eigenvalue ${\lambda }_3({\sigma }^2)=J_{33}(S_e,u_e)$.

If we re-write $\begin{array}{r}{\lambda }_3\left({\sigma }^2\right)=\frac{\theta }{{\sigma }^4}[(({\sigma }^2-u_e^2)+{\sigma }^4(u_e^2-1))\frac{\mathrm{\Lambda }-\mu S_e}{I_e}-(\mu +\kappa +\varphi )({\sigma }^2-u_e^2)]\end{array}$and observe from (13(13) $\begin{array}{r}{S}_e=\frac{1}{c}(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1-{\sigma }^2}{2{\sigma }^2}{u}_e^2)\end{array}$(13) ) and (14(14) $\begin{array}{r}{I}_e=\frac{1}{c}\frac{\mathrm{\Lambda }-\mu S_e}{S_e}\exp \left(\frac{1}{2}{u}_e^2\right)\end{array}$(14) ) that $\begin{array}{r}\frac{\mathrm{\Lambda }-\mu S_e}{I_e}=c{S}_e\exp (-\frac{1}{2}{u}_e^2)=\frac{\mu +\kappa +\varphi }{1-{\sigma }^2}\end{array}$we find that $\begin{array}{r}{\lambda }_3\left({\sigma }^2\right)=-\frac{\theta }{{\sigma }^2}{u}_e^2(\mu +\kappa +\varphi )<0.\end{array}$As ${\sigma }^2\to 0$ we see, by using (17(17) $\begin{array}{r}\underset{{\sigma }^2\to 0}{lim}\left(\begin{array}{c}{S}_e\\ I_e\\ u_e\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{c}(\mu +\kappa +\varphi )}\\ {\displaystyle \frac{1}{c}\frac{\mu }{\mu +\kappa +\varphi }(\frac{c\mathrm{\Lambda }}{\mu }-(\mu +\kappa +\varphi ))}\\ 0\end{array}\right).\end{array}$(17) ), that the matrix (18(18) $\begin{array}{r}\left(\begin{array}{cc}{\displaystyle -(\mu +\frac{\mathrm{\Lambda }-\mu S_e}{S_e})}& -S_{e}c{e}^{-\frac{1}{2}{u}_e^2}\\ {\displaystyle \frac{\mathrm{\Lambda }-\mu S_e}{S_e}}& 0\end{array}\right)\end{array}$(18) ) approaches the matrix $\begin{array}{r}\left(\begin{array}{cc}{\displaystyle -c\frac{\mathrm{\Lambda }}{\mu +\kappa +\varphi }}& -(\mu +\kappa +\varphi )\\ {\displaystyle \frac{\mu }{\mu +\kappa +\varphi }(\frac{c\mathrm{\Lambda }}{\mu }-(\mu +\kappa +\varphi ))}& 0\end{array}\right).\end{array}$Since the trace of this matrix is negative and its determinant is positive, it follows that its two eigenvalues lie in the left half complex plane. Thus, ${\lambda }_1({\sigma }^2)$ and ${\lambda }_2({\sigma }^2)$ lie in the left half complex plane for ${\sigma }^2$ small.

4.3. Summary of equilibrium analysis

The existence and stability results for the equilibria of the Darwinian SI model (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) that we obtained in the preceding sections are geometrically summarized by the two maps in the $({\sigma }^2,\tau )$-plane shown in Figure 1. The relevant infinite rectangle $\begin{array}{r}R:0\le \tau \le 1,{\sigma }^2>0\end{array}$in the $({\sigma }^2,\tau )$-plane is partitioned into regions that result from the constraints determining the existence and stability of the four types of equilibria. These regions are determined by the following curves and lines:

• the vertical straight line at ${\sigma }^2=1;$

• the horizontal line at $\tau ={\tau }_0;$

• the straight line L given by the equation $\begin{array}{r}\tau =\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$where $\begin{array}{r}{\sigma }_1^2:=\frac{\mu }{\mathrm{\Lambda }c}(\gamma -\kappa );\end{array}$

• the curves given respectively by the equations $\begin{array}{r}{C}_0:\tau =\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2},0<{\sigma }^2<{\sigma }_0^2\end{array}$and $\begin{array}{r}{C}_1:\tau =\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}{\left((1-{\sigma }^2)\frac{\mathrm{\Lambda }c}{\mu }\frac{1}{\mu +\kappa +\varphi }\right)}^{\frac{1}{{\sigma }^2}},{\sigma }_2^2<{\sigma }^2<{\sigma }_0^2\end{array}$that describe τ as a function of ${\sigma }^2$.

The horizontal line at $\tau ={\tau }_0$, the line L, and the curves $C_0$ and $C_1$ all intersect at the point $({\sigma }^2,\tau )=({\sigma }_0^2,{\tau }_0)$ where $\begin{array}{r}{\sigma }_0^2:=\frac{\mu }{\mathrm{\Lambda }c}(\frac{\mathrm{\Lambda }c}{\mu }-(\mu +\kappa +\varphi )).\end{array}$The line L and the curve $C_1$ intersect the horizontal line at $\tau =1$ at $({\sigma }_1^2,1)$ and $({\sigma }_2^2,1)$, respectively, where ${\sigma }^2={\sigma }_2^2$ is the solution of the equation $\begin{array}{r}1=\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}{\left((1-{\sigma }^2)\frac{\mathrm{\Lambda }c}{\mu }\frac{1}{\mu +\kappa +\varphi }\right)}^{\frac{1}{{\sigma }^2}}.\end{array}$Within each region appearing in the rectangle R there exists a certain number of equilibria, as indicated by the tables in Figure 1. In each region there is exactly one stable equilibrium, with the exception of region IV. For parameter points $({\sigma }^2,\tau )$ lying in region IV, we show in the Appendix that the asymptotic dynamics are disease-free in the sense that $\underset{t\to \mathrm{\infty }}{lim}I(t)=0$, even though orbits do not equilibrate. As indicated in Figure 1 the stability of $E{E}_{+}$ in regions A, B, and C is qualified in that it has been rigorously established only for ${\sigma }^2$ small and is conjectured to be true throughout these regions.

As one moves around on the parameter maps in Figure 1, one passes through various regions and by so doing loses or gains equilibria and/or equilibrium stability. These occurrences result from equilibrium bifurcations which can be displayed in bifurcation diagrams. For example, fixing the variance ${\sigma }^2<{\sigma }_2^2$ and varying τ from 0 to 1 one follows a vertical path on the left side of the parameter map in Figure 1, passing through regions D, C, B and A, along which the equilibrium counts change from 2 to 3 to 4 to 3. A sample bifurcation diagram for such a path is shown in Figure 2 along with some other sample bifurcation diagrams along other vertical paths.

Figure 2. Four sample bifurcation diagrams showing the equilibrium components I and u plotted against the bifurcation parameter τ with model parameter values $\mathrm{\Lambda }=1,$ $\mu =0.1,$$c=0.1,$ $\kappa =0.1,$ $\varphi =0.1,$ and $\gamma =1.$ The fours diagrams correspond to vertical paths in the parameter map in Figure 1(A) located at the displayed value of ${\sigma }^2$. For the upper left diagram ${\sigma }^2=0.2<{\sigma }_2^2=0.47$ and for the upper right diagram, ${\sigma }^2=0.8$ is between ${\sigma }_0^2=0.7$ and ${\sigma }_1^2=0.9$ in Figure 1 For the lower left diagram ${\sigma }^2=0.9$ is between ${\sigma }_1^2$ and 1 and for the lower right diagram ${\sigma }^2=1.1>1$ in Figure 1.

5. Concluding remarks

As shown in Section 2, under the conditions (2(2) $\begin{array}{r}\mu +\kappa +\varphi <c\frac{\mathrm{\Lambda }}{\mu }<\gamma +\mu +\varphi .\end{array}$(2) ) on the parameters in the non-evolutionary SI model (1(1) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-cSI\\ I^{\mathrm{\prime }}& =cSI-\rho I.\end{array}$(1) ), any treatment level $\tau >{\tau }_0$ will eliminate the disease pathogen (in the sense that the disease-free equilibrium is stable), but the treatment will fail (in the sense that there is a stable endemic equilibrium) if $\tau <{\tau }_0$. We can use the results of the equilibrium analysis of the Darwinian SI model (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ), summarized in Figure 1, to determine the success or failure of the treatment when the disease pathogen can develop resistance to the treatment by means of evolutionary adaptation.

Note that the shaded regions in both maps correspond to disease-free asymptotic dynamics while the unshaded regions correspond to endemic asymptotic dynamics. This implies that the treatment is successful if and only if τ and ${\sigma }^2$ are appropriately large; otherwise the disease becomes endemic.

Recall that τ in the Darwinian SI model (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) is the treatment effectiveness against the disease pathogens who have no resistance to the treatment (i.e. trait v=0) and that the threshold treatment level ${\tau }_0$ is defined by (3(3) $\begin{array}{r}\tau >{\tau }_0:=\frac{c\frac{\mathrm{\Lambda }}{\mu }-(\mu +\kappa +\varphi )}{\gamma -\kappa }\end{array}$(3) ). It is therefore no surprise that a treatment $\tau <{\tau }_0$ will also fail to eliminate the pathogen in the Darwinian model (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ), as the parameter maps in Figure 1 show (in which regions B, C, D and E correspond to stable endemic equilibria). That is to say, if the treatment fails to eradicate a non-evolving pathogen, then it will fail to eradicate an evolving pathogen.

It is also clear from the parameter maps in Figure 1 that a treatment level $\tau >{\tau }_0$ that eliminates the pathogen in the non-evolutionary model will not necessary eliminate the pathogen in the Darwinian model. In order to successfully defeat an evolving pathogen (i.e. for there to be a stable disease-free equilibrium), not only τ but also the ratio of the variance of the treatment effectiveness to the ratio of the variance of the infection transmission rate, ${\sigma }^2$, must be sufficiently large. More precisely, the pair $({\sigma }^2,\tau )$ must lie in the shaded region shown in the parameter maps. This leads to three distinct possibilities, depending on the variance ${\sigma }^2$ of the treatment effectiveness:

1. ${\sigma }^2>{\sigma }_0^2$ implies any treatment level $\tau >{\tau }_0$ leads asymptotically to a disease-free state;

2. ${\sigma }_2^2<{\sigma }^2<{\sigma }_0^2$ implies a treatment level τ leads asymptotically to a disease-free state if and only if it meets a higher threshold (placing the point $({\sigma }^2,\tau )$ above the curve $C_1$), namely $\begin{array}{r}\tau >{\tau }_1:=\frac{\mu +\kappa +\varphi }{\gamma -\kappa }\frac{{\sigma }^2}{1-{\sigma }^2}{\left((1-{\sigma }^2)\frac{\mathrm{\Lambda }c}{\mu }\frac{1}{\mu +\kappa +\varphi }\right)}^{\frac{1}{{\sigma }^2}};\end{array}$

3. $0<{\sigma }^2<{\sigma }_2^2$ implies no treatment level will eradicate the disease causing pathogen and it will become endemic by evolutionary adaptation.

In the case (i), the distribution of treatment effectiveness as a function of pathogen resistance v is so wide that the treatment succeeds despite the evolutionary adaptation of the pathogen. In case (ii) in which the distribution is narrower, the treatment succeeds only if it meets a higher threshold than what it would need in the absence of evolutionary adaptation. Numerically calculated time series illustrating the adjusted threshold ${\tau }_1$ in case (ii) are shown in Figure 3. The orbits associated with those times series are shown in the three dimensional phase space, along with other sample orbits, in Figure 4. Finally, in case (iii) the distribution of treatment effectiveness is too narrowly constrained around those pathogens with little resistance to be successful in eradicating the disease.

Figure 3. Shown are the time series of $I(t)$ and $u(t)$ for the solutions of the Darwinian SI model (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) with parameter values $\mathrm{\Lambda }=5,$ $\mu =1,$ $c=1,$ $\kappa =0.5,$ $\varphi =0.2,$ $\gamma =10,$ ${\sigma }^2=0.35,$ $\theta =1$ and initial conditions $S(0)=1,$ $I(0)=1$ $u(0)=0.05$ when (a) $\tau =0.4$ and (b) $\tau =0.8$. The corresponding point $({\sigma }^2,\tau )$ is shown as a solid dot in the (numerically calculated) parameter map. In case (a) the point is below the threshold curve $C_1$ and, as predicted the endemic equilibrium $E{E}_{+}$ is stable. In this case the treatment is insufficient to eradicate the adapting disease agent which asymptotically attains a resistance level of approximately u=1. In case (b) the point is above the threshold curve $C_1$ and as predicted the disease-free equilibrium $DF{E}_{+}$ is stable. In this case the treatment eradicates the disease, even though it develops a level of resistance. In this example ${\tau }_0=0.347$, ${\sigma }_0^2=0.660$, and ${\sigma }_2^2=0.281$.

Figure 4. The phase portrait orbits of the time series solutions in Figure 3 are displayed in three dimensional phase space, together with a selection of other orbits (with initial conditions denoted by open circles). The graphs (a) and (b) correspond to the graphs (a) and (b) in Figure 3. In (a) all orbits tend to an equilibrium point indicated by the solid black circle which, because it located in the interior of the positive octant, is an endemic equilibrium. In (b) all orbits tend to an equilibrium point indicated by the solid black circle which, because it located on the vertical coordinate plane I=0, is a disease-free equilibrium.

We also note that a decreased treatment level τ always results in a higher disease prevalence $I_e$ in endemic equilibria. (It is an elementary calculus exercise to show that $d{I}_e/d\tau <0$ for the adaptive endemic equilibria $E{E}_{+}$ and $E{E}_0$). Therefore, in this model a decreased treatment level τ always results in a higher disease prevalence $I_e$.

Although we applied this novel methodology to a relatively simple model to allow for tractability, the derivation of Darwinian dynamics described in Section 3 could in principle be applied to many other compartmental infectious disease models (in continuous or discrete time). The key requirement would be to identify a class that represents the number of infections, from which the fitness of the pathogen can be derived. This methodology could similarly be applied to mosquito population dynamics to model the evolution of insecticide resistance, that can then be combined with models of vector-borne disease dynamics. However, applying Darwinian dynamics to more complicated disease transmission models would likely result in systems of equations where analytical solutions of equilibria and their Jacobians may not be tractable, requiring a greater reliance on numerical simulations. Nevertheless, such models could provide greater understanding of the evolution of resistance in specific pathogens, and the role of the evolution of insecticide resistance in vector-borne disease control.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Junpyo Park was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) [No. NRF-2021R1A4A1032924, NRF-2023R1A2C1004893]. Nakul Chitnis was supported by the Bill and Melinda Gates Foundation [INV025569].

Appendix

In the triangle in Figure 1(A) defined by the inequalities (A1) $\begin{array}{r}{\sigma }^2>1,{\tau }_0<\tau <1,\tau >\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$(A1) in the case when $\begin{array}{r}{\sigma }_1^2:=\frac{\mu }{\mathrm{\Lambda }c}(\gamma -\kappa )>1\end{array}$the only equilibrium $DF{E}_0$ is unstable.

Theorem A.1

Assume (A1(A1) $\begin{array}{r}{\sigma }^2>1,{\tau }_0<\tau <1,\tau >\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$(A1) ) and ${\sigma }_1^2>1$. For solutions $(S(t),I(t),u(t))$ of  (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) with initial conditions $(S(0),I(0),u(0))\in R_{+}^3$ the following hold:

(a)	$lim\underset{t\to \mathrm{\infty }}{sup}S(t)\le \mathrm{\Lambda }/\mu$
(b)	$u(0)>0$ implies $\underset{t\to \mathrm{\infty }}{lim}u(t)=+\mathrm{\infty }$
(c)	$\underset{t\to \mathrm{\infty }}{lim}I(t)=0.$

Proof.

From Lemma 3.1 it follows that $S(t)>0$ and $I(t)\ge 0$ for all t>0.

1. From the first equation in (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) we see that $\begin{array}{r}{S}^{\mathrm{\prime }}\le \mathrm{\Lambda }-\mu S.\end{array}$By standard comparison theorems $\begin{array}{r}0\le S\left(t\right)\le x\left(t\right)\end{array}$for all $t\ge 0$ where $x(t)$ is the solution of the initial value problem $\begin{array}{rl}{x}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu x\\ x\left(0\right)& =S\left(0\right)\ge 0.\end{array}$Since $lim\underset{t\to \mathrm{\infty }}{sup}x(t)=\mathrm{\Lambda }/\mu$ the result follows.

2. By (a), for any $\epsilon >0$ there exists a $t^{\ast }(\epsilon )>0$ such that $\begin{array}{r}0\le S\left(t\right)<\frac{\mathrm{\Lambda }}{\mu }+\epsilon \mathrm{for}t\ge t^{\ast }\left(\epsilon \right).\end{array}$

From the third equation in (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) $\begin{array}{rl}{u}^{\mathrm{\prime }}& >\theta u(-c\exp (-\frac{1}{2}{u}^2)(\frac{\mathrm{\Lambda }}{\mu }+\epsilon )+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}exp(-\frac{1}{2{\sigma }^2}{u}^2))\\ & =\theta u\exp (-\frac{1}{2{\sigma }^2}{u}^2)(-c\exp (-\frac{1}{2}\left(\frac{{\sigma }^2-1}{{\sigma }^2}\right)u^2)(\frac{\mathrm{\Lambda }}{\mu }+\epsilon )+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}).\end{array}$Since ${\sigma }^2>1$ in (A1(A1) $\begin{array}{r}{\sigma }^2>1,{\tau }_0<\tau <1,\tau >\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$(A1) ) we get $\begin{array}{rl}{u}^{\mathrm{\prime }}& >\theta u\exp (-\frac{1}{2{\sigma }^2}{u}^2)(-c(\frac{\mathrm{\Lambda }}{\mu }+\epsilon )+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2})\\ & =\theta u\exp (-\frac{1}{2{\sigma }^2}{u}^2)(\frac{\gamma -\kappa }{{\sigma }^2}(-\frac{\mathrm{\Lambda }c}{\mu }\frac{1}{\gamma -\kappa }{\sigma }^2+\tau )-c\epsilon )\\ & =\theta u\exp (-\frac{1}{2{\sigma }^2}{u}^2)(\frac{\gamma -\kappa }{{\sigma }^2}(\tau -\frac{{\sigma }^2}{{\sigma }_1^2})-c\epsilon ).\end{array}$If we choose ε so that $\begin{array}{r}0<\epsilon <\frac{1}{c}\frac{\gamma -\kappa }{{\sigma }^2}(\tau -\frac{{\sigma }^2}{{\sigma }_1^2})\end{array}$(which we can do by (A1(A1) $\begin{array}{r}{\sigma }^2>1,{\tau }_0<\tau <1,\tau >\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$(A1) )), then $u^{\mathrm{\prime }}>0$ for all $t>t^{\ast }(\epsilon )$, i.e. $u(t)>0$ is eventually monotonically increasing and if we define $\begin{array}{r}{u}_{\mathrm{\infty }}:=\underset{t\to \mathrm{\infty }}{lim}u\left(t\right)\end{array}$then either $u_{\mathrm{\infty }}<+\mathrm{\infty }$ or $u_e=+\mathrm{\infty }$. If $u_{\mathrm{\infty }}<+\mathrm{\infty }$, then the first two equations of (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) constitute an asymptotically autonomous system (A2) $\begin{array}{r}\begin{array}{l}{S}^{\mathrm{\prime }}=\mathrm{\Lambda }-\mu S-cexp(-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}=(cexp(-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\end{array}\end{array}$(A2) for S and I with limiting system $\begin{array}{r}\begin{array}{l}{S}^{\mathrm{\prime }}=\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}_{\mathrm{\infty }}^2)SI\\ I^{\mathrm{\prime }}=(c\exp (-\frac{1}{2}{u}_{\mathrm{\infty }}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}_{\mathrm{\infty }}^2))I.\end{array}\end{array}$(Note that right sides of (A2(A2) $\begin{array}{r}\begin{array}{l}{S}^{\mathrm{\prime }}=\mathrm{\Lambda }-\mu S-cexp(-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}=(cexp(-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\end{array}\end{array}$(A2) ) converge uniformly as $t\to \mathrm{\infty }$ for $(S,I)$ in any compact set in $R^2.$) This limiting system is an SI model of the form (1(1) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-cSI\\ I^{\mathrm{\prime }}& =cSI-\rho I.\end{array}$(1) ). As pointed out in Section 2 solutions globally equilibrate on $R_{+}^2$ to an equilibrium $(S_{\mathrm{\infty }},I_{\mathrm{\infty }})$ [2]. Then we can invoke Theorem A.2 to conclude that $(S(t),I(t),u(t))$ approaches the equilibrium $(S_{\mathrm{\infty }},I_{\mathrm{\infty }},u_{\mathrm{\infty }}).$ However, this is a contraction to the fact that on the triangular region (A1(A1) $\begin{array}{r}{\sigma }^2>1,{\tau }_0<\tau <1,\tau >\frac{{\sigma }^2}{{\sigma }_1^2}\end{array}$(A1) ) the only equilibrium present in the system is $DF{E}_0$ for which the u component equals 0. Thus, $u_{\mathrm{\infty }}<+\mathrm{\infty }$ cannot hold and it follows that $u_e=+\mathrm{\infty }$.

(c) Since $u_e=+\mathrm{\infty }$ the first two equations for S and I in the Darwinian model (6(6) $\begin{array}{rl}{S}^{\mathrm{\prime }}& =\mathrm{\Lambda }-\mu S-c\exp (-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}& =(c\exp (-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\\ u^{\mathrm{\prime }}& =\theta u(-c\exp (-\frac{1}{2}{u}^2)S+(\gamma -\kappa )\tau \frac{1}{{\sigma }^2}\exp (-\frac{1}{2{\sigma }^2}{u}^2)).\end{array}$(6) ) are asymptotically autonomous system (A2(A2) $\begin{array}{r}\begin{array}{l}{S}^{\mathrm{\prime }}=\mathrm{\Lambda }-\mu S-cexp(-\frac{1}{2}{u}^2)SI\\ I^{\mathrm{\prime }}=(cexp(-\frac{1}{2}{u}^2)S-(\mu +\kappa +\varphi )-(\gamma -\kappa )\tau \exp (-\frac{1}{2{\sigma }^2}{u}^2))I\end{array}\end{array}$(A2) ) for $S(t)$ and $I(t)$ with the limiting equation $\begin{array}{r}\begin{array}{l}{S}^{\mathrm{\prime }}=\mathrm{\Lambda }-\mu S\\ I^{\mathrm{\prime }}=-(\mu +\kappa +\varphi )I\end{array}\end{array}$all of whose solutions equilibria to $(S,I)=(\mathrm{\Lambda }/\mu ,0).$ Invoking Theorem A.2 below we conclude that that $\underset{t\to \mathrm{\infty }}{lim}I(t)=0$.

Theorem A.2

(Theorem 1.3 in [17] or Theorem 7 in [12]) The omega limit set of any orbit of the asymptotically autonomous system either contains equilibria of the limiting system or is the union of periodic orbits of the limiting system.