Coevolutionary Dynamics of Host Immune and Parasite Virulence Based on an Age-Structured Epidemic Model

Abstract

Hosts can activate a defensive response to clear the parasite once being infected. To explore how host survival and fecundity are affected by host-parasite coevolution for chronic parasitic diseases, in this paper, we proposed an age-structured epidemic model with infection age, in which the parasite transmission rate and parasite-induced mortality rate are structured by the infection age. By use of critical function analysis method, we obtained the existence of the host immune evolutionary singular strategy which is a continuous singular strategy (CSS). Assume that parasite-induced mortality begins at infection age $\tau$ and is constant v thereafter. We got that the value of the CSS, $c^*$, monotonically decreases with respect to infection age $\tau$ (see Case (I)), while it is non-monotone if the constant v positively depends on the immune trait c (see Case (II)). This non-monotonicity is verified by numerical simulations and implies that the direction of immune evolution depends on the initial value of immune trait. Besides that, we adopted two special forms of the parasite transmission rate to study the parasite’s virulence evolution, by maximizing the basic reproduction ratio ${\mathcal {R}}_0$. The values of the convergence stable parasite’s virulence evolutionary singular strategies $v^*$ and $k^*$ increase monotonically with respect to time lag L (i.e., the time lag between the onset of transmission and mortality). At the singular strategy $v^*$ and $k^*$, we further obtained the expressions of the case mortalities $\chi ^*$ and how they are affected by the time lag L. Finally, we only presented some preliminary results about host and parasite coevolution dynamics, including a general condition under which the coevolutionary singular strategy $(c^*,v^*)$ is evolutionarily stable.

Appendices

Appendix A. Local stability of the positive steady state $E^*$

Proof

Assuming that $\beta (\theta )$ is taken as (80) and $v(\theta )$ is adopted as (24) or (37). If $b_S>\mu$ and $\frac{1}{2}<{\mathcal {R}}_0<{\mathcal {P}}_1(c)<1$,the positive steady state $E^*$ of model (2) is locally stable. To prove this local stability, we make the following perturbation variables around $E^*$,

$$\begin{aligned} {S(t)} = S^* + x(t),\ \ {I}(\theta ,t) = I^*(\theta ) + y(\theta ,t). \end{aligned}$$

(73)

Substituting (73) into model (2), we have the linearized system

$$\begin{aligned} \left\{ \begin{aligned} \frac{dx(t)}{dt} =&({b_S} - \mu )x(t) +({b_I}(c) + c)\int _0^\infty y(\theta ,t)d\theta \\&-x(t)\int _0^\infty \beta (\theta )I^*(\theta )d\theta -S^*\int _0^\infty \beta (\theta )y(\theta ,t)d\theta ,\\ \frac{\partial y(\theta ,t)}{\partial t} +&\frac{{\partial y(\theta ,t)}}{{\partial \theta }} = -(\mu + v(\theta ) + c)y(\theta ,t),\\ y(0,t) =&S^*\int _0^\infty \beta (\theta )y(\theta ,t)d\theta . \end{aligned}\right. \end{aligned}$$

To analyze the asymptotic behavior near $E^*$, we assume

$$\begin{aligned} x(t) = {{\overline{x}}}{e^{\lambda t}},\ \ y(\theta ,t) = \overline{y}(\theta ){e^{\lambda t}}. \end{aligned}$$

(74)

By use of (74), we get the following equation

$$\begin{aligned} \left\{ \begin{aligned} \lambda {{\overline{x}}} =&(b_S -\mu ){{\overline{x}}} +({b_I}(c) + c)\int _0^{+\infty } {{{\overline{y}}}(\theta )d\theta } \\&-{{\overline{x}}}\int _0^{+\infty }\beta (\theta )I^*(\theta )d\theta -S^*\int _0^{+\infty }\beta (\theta ){{\overline{y}}}(\theta )d\theta ,\\ \frac{d{\overline{y}}(\theta )}{d\theta } =&-(\lambda +\mu + v(\theta )+ c){{\overline{y}}}(\theta ),\\ {{\overline{y}}}(0) =&x\int _0^{+\infty } {\beta (\theta )I^*(\theta )d\theta } -S^*\int _0^{+\infty } {\beta (\theta ){{\overline{y}}}(\theta )d\theta }. \end{aligned} \right. \end{aligned}$$

(75)

It follows from the second and the third equations of Eq. (75), we have

$$\begin{aligned} {{\overline{y}}}(\theta ) ={{\overline{y}}}(0)e^{ - \lambda \theta }w(\theta ,c). \end{aligned}$$

(76)

Substituting (76) into the third equation of (75), we have

$$\begin{aligned} {{\overline{y}}}(0)={{\overline{x}}}\int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta }+ S^*\overline{y}(0)\int _0^{+\infty }{\beta (\theta ){e^{ - \lambda \theta }}w(\theta ,c)d\theta } \end{aligned}$$

It then follows that

$$\begin{aligned} {{\overline{y}}}(\theta ) = \frac{\displaystyle {{\overline{x}}}\int _0^{ + \infty }\beta (\theta )I^*(\theta )d\theta e^{ - \lambda \theta }w(\theta ,c)}{ \displaystyle 1 - S^*\int _0^{ + \infty } \beta (\theta )e^{ - \lambda \theta }w(\theta ,c)d\theta }. \end{aligned}$$

(77)

Note that

$$\begin{aligned} \int _0^{+\infty }\beta (\theta )I^*(\theta )d\theta =\frac{b_S-\mu }{1-\mathcal {P}_1(c)}. \end{aligned}$$

Substituting (77) into the first equation of (75), by some direct calculations, we have

$$\begin{aligned} (\lambda -b_S+\mu ) + \frac{\left[ 1-(b_I(c)+c)\int _0^{+\infty }e^{-\lambda \theta }w(\theta ,c)d\theta \right] }{\left( 1-S^*\int _0^{+\infty }\beta (\theta )e^{-\lambda \theta }w(\theta ,c)d\theta \right) }\frac{(b_S-\mu )}{(1 -{\mathcal {P}}_1(c))}=0. \end{aligned}$$

(78)

Here, we have used the expressions of $S^*$ and $I^*$. It then follows that

$$\begin{aligned} \left| \lambda + \frac{\left[ 1-(b_I(c)+c)\int _0^{+\infty }e^{-\lambda \theta }w(\theta ,c)d\theta \right] }{\left( 1-S^*\int _0^{+\infty }\beta (\theta )e^{-\lambda \theta }w(\theta ,c)d\theta \right) }\frac{(b_S-\mu )}{(1 -{\mathcal {P}}_1(c))}\right| =|b_S-\mu |, \end{aligned}$$

(79)

If $\lambda =a_1+\text {i}a_2$ with $a_1>0$ is a root of (78), then we have the real part

$$\begin{aligned} \begin{aligned}&\Re \left\{ \frac{1-(b_I(c)+c)\int _0^{+\infty }e^{-\lambda \theta }w(\theta ,c)d\theta }{1-S^*\int _0^{+\infty }\beta (\theta )e^{-\lambda \theta }w(\theta ,c)d\theta }\right\} \\&=\Re \left\{ \frac{1-(b_I(c)+c)\int _0^{+\infty }e^{-\lambda \theta }w(\theta ,c)d\theta }{1-\frac{\int _0^{+\infty }\beta (\theta )e^{-\lambda \theta }w(\theta ,c)d\theta }{\int _0^{+\infty }\beta (\theta )w(\theta ,c)d\theta }}\right\} \\&=\Re \left\{ \frac{C_1+\text {i}D_1}{C_2+\text {i}D_2} \right\} \end{aligned} \end{aligned}$$

with

$$\begin{aligned} \begin{aligned} C_1&=1-(b_I(c)+c)\int _0^{+\infty }e^{-a_1\theta }\cos (a_2\theta )w(\theta ,c)d\theta ,\\ D_1&=(b_I(c)+c)\int _0^{+\infty }e^{-a_1\theta }\sin (a_2\theta )w(\theta ,c)d\theta ,\\ C_2&=1-\frac{1}{{\mathcal {R}}_0}\int _0^{+\infty }\beta (\theta )e^{-a_1\theta }\cos (a_2\theta )w(\theta ,c)d\theta ,\\ D_2&=\frac{1}{\mathcal {R}_0}\int _0^{+\infty }\beta (\theta )e^{-a_1\theta }\sin (a_2\theta )w(\theta ,c)d\theta . \end{aligned} \end{aligned}$$

Obviously, $0<C_1<1$ if ${\mathcal {P}}_1(c)<1$.

To determine the signs of $C_2$ and $D_2$, we need the following assumption. Assume that

$$\begin{aligned} \beta (\theta )={{\overline{\beta }}}(1-e^{-\alpha \theta })\ \ \text{ with }\ \ 0<{{\overline{\beta }}}\le b_I(c)+c. \end{aligned}$$

(80)

Then it follows that $0<C_2<1$ if ${\mathcal {P}}_1(c)<1$, and $0\le D_1<1$ and $0\le D_2<1$ if ${\mathcal {P}}_1(c)<1$ and $a_2\ge 0$, and $-1<D_1<0$ and $-1<D_2<0$ if ${\mathcal {P}}_1(c)<1$ and $a_2<0$. If $\frac{1}{2}<{\mathcal {R}}_0<{\mathcal {P}}_1(c)<1$, we have $0<1-\mathcal P_1(c)<1-{\mathcal {R}}_0<{\mathcal {R}}_0$. By use of (80), we have that $D_1>{\mathcal {R}}_0 D_2>0$ if $a_2>0$ and $D_1<{\mathcal {R}}_0 D_2<0$ if $a_2<0$. Then we have

$$\begin{aligned}&\frac{D_1D_2}{(1 -{\mathcal {P}}_1(c))}>\frac{D_1D_2}{1-{\mathcal {R}}_0}>\frac{D_1D_2}{{\mathcal {R}}_0}>D_2^2,\ \ \ \text{ if }\ \ a_2>0; \\&\frac{D_1D_2}{(1 -{\mathcal {P}}_1(c))}>\frac{D_1D_2}{\mathcal {R}_0}=\frac{(-D_1)(-D_2)}{{\mathcal {R}}_0}>D_2^2,\ \ \text{ if }\ \ a_2<0. \end{aligned}$$

It then follows from $0<C_2<1$ and $\frac{C_1}{(1 -\mathcal P_1(c))}>1$ that

$$\begin{aligned} \begin{aligned} \Re \left\{ \frac{C_1+\text {i}D_1}{C_2+\text {i}D_2} \right\} \frac{1}{(1 -{\mathcal {P}}_1(c))}&=\frac{C_1C_2+D_1D_2}{C_2^2+D_2^2}\frac{1}{(1 -{\mathcal {P}}_1(c))}\\&=\frac{\displaystyle \frac{C_1}{(1 -{\mathcal {P}}_1(c))}C_2+\frac{D_1D_2}{(1 -{\mathcal {P}}_1(c))}}{C_2^2+D_2^2}\\&>\frac{C_2+D_2^2}{C_2^2+D_2^2}>1. \end{aligned} \end{aligned}$$

Thus, the left-hand side of (79) (LHS) satisfies

$$\begin{aligned} |b_S-\mu |{} & {} <\Re \left\{ \frac{\left[ 1-(b_I(c)+c)\int _0^{+\infty }e^{-\lambda \theta }w(\theta ,c)d\theta \right] }{\left( 1-S^*\int _0^{+\infty }\beta (\theta )e^{-\lambda \theta }w(\theta ,c)d\theta \right) }\right\} \frac{|b_S-\mu |}{(1 -{\mathcal {P}}_1(c))}\\{} & {} <LHS=|b_S-\mu | \end{aligned}$$

if $b_S>\mu$ and $\frac{1}{2}<{\mathcal {R}}_0<{\mathcal {P}}_1(c)<1$. This contradiction implies that Eq. (78) only has roots with negative real parts and the steady state $E^*=(S^*,I^*)$ of model (2) is locally asymptotically stable if $b_S>\mu$ and $\frac{1}{2}<{\mathcal {R}}_0<{\mathcal {P}}_1(c)<1$. $\square$

Appendix B. Fitness function of hosts immune evolution

To get the the fitness function of hosts immune evolution described by model (2), we need to study the dynamics of the following system

$$\begin{aligned} \left\{ \begin{aligned} \frac{d S_{1}(t)}{d t}=&\, b_S S_{1}(t)-\mu S_{1}(t)+\int _{0}^{+\infty } b_{I}(c) I_{1}(\theta , t) d \theta +\int _{0}^{+\infty } c I_{1}(\theta ,t) d \theta \\ {}&\, -\left[ \int _{0}^{+\infty } \beta (\theta ) I_{1}(\theta ,t) d \theta +\int _{0}^{+\infty } \beta (a) I_{2}(a,t) d a\right] S_{1}(t) \\ \frac{{\partial {I_1}(\theta ,t)}}{{\partial t}} +&\frac{{\partial {I_1}(\theta ,t)}}{{\partial \theta }} = -(\mu + v(\theta ) + c){I_1}(\theta ,t)\\ {I_1}(0,t) =&\, \left[ \int _0^{ + \infty } {\beta (\theta ){I_1}(\theta ,t)d\theta } + \int _0^{ + \infty } {\beta (a){I_2}(a,t)da} \right] {S_1}(t)\\ \frac{d S_{2}(t)}{d t}=&b_S S_{2}(t)-\mu S_{2}(t)+\int _{0}^{+\infty } b_{I}({\hat{c}}) I_{2}(a,t) d a+\int _{0}^{+\infty } {\hat{c}} I_{2}(a,t) da \\ {}&\, -\left[ \int _{0}^{+\infty } \beta (\theta ) I_{1}(\theta ,t) d \theta +\int _{0}^{+\infty } \beta (a) I_{2}(a,t) d a\right] S_{2}(t) \\ \frac{{\partial {I_2}(a,t)}}{{\partial t}} +&\frac{{\partial {I_2}(a,t)}}{{\partial a}} = -(\mu + v(a) + {{\hat{c}}}){I_2}(a,t)\\ {I_2}(0,t) =&\, \left[ \int _0^{ + \infty } {\beta (\theta ){I_1}(\theta ,t)d\theta } + \int _0^{ + \infty } {\beta (a){I_2}(a,t)da}\right] {S_2}(t). \end{aligned} \right. \nonumber \\ \end{aligned}$$

(81)

in which $S_1$ and $I_1$ represent the resident hosts, $S_2$ and $I_2$ represent the mutant hosts. $v(\theta )$ is the parasite-induced mortality rate for the resident hosts, and v(a) is the parasite-induced mortality rate for the mutant hosts. Other parameters are same as that in model (2). It is easy to obtain that, if $b_S>\mu$ and ${\mathcal {P}}_1(c)<1$, there is a boundary steady state $E_1^*=(S^*_1,I^*_1(\theta ),0,0_{L^1})$ with

$$\begin{aligned} S^*_1=\frac{1}{\displaystyle \int _0^{+\infty }\beta (\theta )w(\theta ,c)d\theta }=S^*\ \ \ \text{ and }\ \ \ I^*_1(\theta )=\frac{(b_S-\mu )S^*_1w(\theta ,c)}{1-\mathcal {P}_1(c)}=I^*(\theta ). \end{aligned}$$

In which, $(S^*,I^*(\theta ))$ is the positive steady state of model (2). By use of the immune regulation reproduction ratio ${\mathcal {R}}_H$ in (9), we have the following theorem.

Theorem B. When ${\mathcal {R}}_H<1$, the boundary steady state $E_1^*$ of model (81) is locally asymptotically stable if it exists, otherwise, it is unstable.

Proof

To prove Theorem B, we make the following perturbation variables around $E_1^*$,

$$\begin{aligned} {S_1(t)}= & {} S^*+ x_1(t),\ {I_1}(\theta ,t)= I^*(\theta ) + y_1(\theta ,t),\ {S_2}(t)\nonumber \\= & {} {x_2}(t),\ {I_2}(a,t)= {y_2}(a,t). \end{aligned}$$

(82)

Substituting (82) into model (81), we have the linearized system

$$\begin{aligned} \left\{ \begin{aligned} \frac{dx_1(t)}{dt} =&({b_S} - \mu )x_1(t)+({b_I}(c) + c)\int _0^\infty y_1(\theta ,t)d\theta -x_1(t)\int _0^\infty \beta (\theta )I^*(\theta )d\theta \\&-S^*\left[ \int _0^\infty {\beta (\theta )y_1(\theta ,t)d\theta } + \int _0^\infty {\beta (a){y_2}(a,t)da}\right] ,\\ \frac{\partial y_1(\theta ,t)}{\partial t} +&\frac{{\partial y_1(\theta ,t)}}{{\partial \theta }} = -(\mu + v(\theta ) + c)y_1(\theta ,t),\\ y_1(0,t) =&S^*\left[ \int _0^\infty {\beta (\theta )y_1(\theta ,t)d\theta } + \int _0^\infty {\beta (a){y_2}(a,t)da}\right] + x_1(t)\int _0^\infty {\beta (\theta )I^*(\theta )d\theta }, \\ \frac{d{x_2(t)}}{dt} =&(b_S-\mu ){x_2(t)}+({b_I}({{\hat{c}}})+{{\hat{c}}})\int _0^\infty {{y_2}(a,t)da}- {x_2(t)}\int _0^\infty {\beta (\theta )I^*(\theta )d\theta }, \\ \frac{{\partial {y_2}(a,t)}}{{\partial t}} +&\frac{{\partial {y_2}(a,t)}}{{\partial a}} = -(\mu + v(a) + {{\hat{c}}}){y_2}(a,t),\\ {y_2}(0,t) =&{x_2(t)}\int _0^\infty \beta (\theta )I^*(\theta )d\theta . \end{aligned}\right. \end{aligned}$$

To analyze the asymptotic behavior near $E_1^*$, we assume

$$\begin{aligned} x_1(t) ={{\overline{x}}}_1{e^{\lambda t}},\ y_1(\theta ,t) =\overline{y}_1(\theta ){e^{\lambda t}},\ {x_2}(t) ={{\overline{x}}}_2{e^{\lambda t}},\ {y_2}(a,t) = {{\overline{y}}}_2(a){e^{\lambda t}}. \end{aligned}$$

(83)

By use of (83), we get the following equation

$$\begin{aligned} \left\{ \begin{aligned} \lambda {{\overline{x}}}_1 =&(b_S -\mu ){{\overline{x}}}_1 +({b_I}(c) + c)\int _0^{+\infty } {{{\overline{y}}}_1(\theta )d\theta }- {{\overline{x}}}_1\int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta } \\&-S^*\left[ \int _0^{ + \infty } {\beta (\theta ){{\overline{y}}}_1(\theta )d\theta } + \int _0^{ + \infty } {\beta (a){{{\overline{y}}}_2}(a)da}\right] ,\\ \frac{d{{\overline{y}}}_1(\theta )}{d \theta } =&-(\lambda +\mu + v(\theta ) + c){{\overline{y}}}_1(\theta ),\\ {{\overline{y}}}_1(0) =&{{\overline{x}}}_1\int _0^{+\infty }\beta (\theta )I^*(\theta )d\theta + S^*\left[ \int _0^{+\infty }\beta (\theta ){{\overline{y}}}_1(\theta )d\theta +\int _0^{+\infty }\beta (a){{{\overline{y}}}_2}(a)da\right] ,\\ \lambda {{{\overline{x}}}_2}=&({b_S} - \mu ){{{\overline{x}}}_2} +({b_I}({{\hat{c}}}) + {{\hat{c}}})\int _0^{ + \infty }{{{\overline{y}}}_2}(a)da - {{{\overline{x}}}_2}\int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta },\\ \frac{d{{{\overline{y}}}_2}(a)}{{d a}} =&-(\lambda +\mu + v(a) + {{\hat{c}}}){{{\overline{y}}}_2}(a),\\ {{{\overline{y}}}_2}(0) =&{{{\overline{x}}}_2}\int _0^{ + \infty }\beta (\theta )I^*(\theta )d\theta . \end{aligned} \right. \nonumber \\ \end{aligned}$$

(84)

It follows from the second and the third equations of Eq. (84) that we have

$$\begin{aligned} {{\overline{y}}}_1(\theta ) = {{\overline{y}}}_1(0)e^{ - \lambda \theta }w(\theta ,c). \end{aligned}$$

(85)

It follows from the fifth and the sixth equations of Eq. (84) that we have

$$\begin{aligned} {{{\overline{y}}}_2}(a)={{{\overline{y}}}_2}(0)e^{ - \int _0^a(\mu + v(\tau ) + {{\hat{c}}} + \lambda )d\tau }={{\overline{x}}}_2e^{-\lambda a}w(a,{{\hat{c}}})\int _0^{ + \infty }\beta (\theta )I^*(\theta )d\theta . \end{aligned}$$

(86)

Substituting (85) into the third equation of (84), we have

$$\begin{aligned} {{\overline{y}}}_1(0)&={{\overline{x}}}_1\int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta }+ S^*\left[ \overline{y}_1(0)\int _0^{ + \infty } {\beta (\theta ){e^{ - \lambda \theta }}w(\theta ,c)d\theta }\right. \\ {}&\quad \left. + \int _0^{ + \infty } {\beta (a){\overline{y}_2}(a)da}\right] \end{aligned}$$

It then follows from (86), we have

$$\begin{aligned} {{\overline{y}}}_1(0) = \frac{\displaystyle \int _0^{ + \infty }\beta (\theta )I^*(\theta )d\theta \left( {{\overline{x}}}_1 + S^*{{\overline{x}}}_2\int _0^{ + \infty }\beta (a)e^{-\lambda a}w(a,{{\hat{c}}})da\right) }{ \displaystyle 1 - S^*\int _0^{ + \infty } \beta (\theta )e^{ - \lambda \theta }w(\theta ,c)d\theta }. \end{aligned}$$

(87)

Substituting (87) into Eq. (85), we get

$$\begin{aligned} {{\overline{y}}}_1(\theta ) = \frac{\displaystyle \int _0^{ + \infty }\beta (\theta )I^*(\theta )d\theta \left( {{\overline{x}}}_1+ S^*{{\overline{x}}}_2\int _0^{ + \infty }\beta (a)e^{-\lambda a}w(a,{{\hat{c}}})da\right) }{ \displaystyle 1 - S^*\int _0^{ + \infty } \beta (\theta )e^{ - \lambda \theta }w(\theta ,c)d\theta } e^{ - \lambda \theta }w(\theta ,c).\nonumber \\ \end{aligned}$$

(88)

Substituting (86) and (88) into the first equation of (84), substituting (86) into the forth equation of (84), we have

$$\begin{aligned} \left\{ \begin{aligned}&{{\overline{x}}}_1 H_1(\lambda ) + {{{\overline{x}}}_2}H_2(\lambda )S^*\int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta } \int _0^{ + \infty }\beta (a){e^{ - \lambda a}}w(a,{{\hat{c}}})da = 0,\\&\lambda {{{\overline{x}}}_2} -(b_S- \mu ){{{\overline{x}}}_2} + {{{\overline{x}}}_2}\int _0^{ + \infty } \beta (\theta )I^*(\theta )d\theta \\&\ \ \ \ \ \ -(b_I({{\hat{c}}})+{{\hat{c}}}){{{\overline{x}}}_2}\int _0^{+\infty } \beta (\theta )I^*(\theta )d\theta \int _0^{ + \infty }e^{-\lambda a}w(a,{{\hat{c}}})da=0, \end{aligned} \right. \end{aligned}$$

(89)

where

$$\begin{aligned} H_1(\lambda ) =&\lambda -(b_S- \mu ) - \frac{\left[ 1-\displaystyle ({b_I}(c) + c)\int _0^{ + \infty } e^{ - \lambda \theta }w(\theta ,c)d\theta \right] }{\displaystyle 1-S^*\int _0^{ + \infty } \beta (\theta )e^{ - \lambda \theta }w(\theta ,c)d\theta }\int _0^{ + \infty }\beta (\theta )I^*(\theta )d\theta ,\\ H_2(\lambda )&= 1 -\frac{\displaystyle (b_I(c) + c)\int _0^{ + \infty }{e^{ - \lambda \theta }}w(\theta ,c)d\theta }{\displaystyle 1-S^*\int _0^{+\infty }\beta (\theta )e^{ - \lambda \theta }w(\theta ,c)d\theta } + \frac{\displaystyle S^*\int _0^{ + \infty }{e^{ - \lambda \theta }}w(\theta ,c)d\theta }{\displaystyle 1-S^*\int _0^{ + \infty } \beta (\theta )e^{ - \lambda \theta }w(\theta ,c)d\theta }. \end{aligned}$$

If ${{{\overline{x}}}_2}=0$, from the second equation of (89), we can obtain the characteristic equation

$$\begin{aligned} H_1(\lambda ) =0, \end{aligned}$$

which is same as the equation (78). Then $H_1(\lambda ) =0$ has no roots with positive real parts and the boundary steady state $E_1^*$ is locally asymptotically stable if $b_S>\mu$ and $\frac{1}{2}<{\mathcal {R}}_0<{\mathcal {P}}_1(c)<1$ under the condition that $\beta (\theta )$ is taken as (80) and $v(\theta )$ is adopted as (24) or (37). If ${{{\overline{x}}}_2}\ne 0$, from the second equation of (89), we can obtain the characteristic equation

$$\begin{aligned}{} & {} \lambda - \left( {{b_S} - \mu } \right) + \int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta } - \left( {{b_I}({{\hat{c}}}) + {{\hat{c}}}} \right) \nonumber \\{} & {} \int _0^{ + \infty } {\beta (\theta )I^*(\theta )d\theta } \int _0^{ + \infty } {{e^{ - \lambda a}}w(a,{{\hat{c}}})da} = 0. \end{aligned}$$

(90)

Note that

$$\begin{aligned} \int _0^{ + \infty } \beta (\theta )I^*(\theta )d\theta =\frac{(b_S-\mu )}{1-{\mathcal {P}}_1(c)}. \end{aligned}$$

Then Eq. (90) is equivalent to

$$\begin{aligned} \lambda + \frac{(b_S-\mu )}{1-\mathcal {P}_1(c)}=(b_S-\mu )\left( 1+\frac{(b_I({{\hat{c}}}) + {{\hat{c}}})\int _0^{+\infty } e^{-\lambda a}w(a,{{\hat{c}}})da}{1-\mathcal P_1(c)}\right) . \end{aligned}$$

(91)

If $\lambda$ with $\Re \lambda >0$ is a root of Eq. (91) when $b_S>\mu$ and $0<{\mathcal {P}}_2({{\hat{c}}})<{\mathcal {P}}_1(c)<1$, i.e., ${\mathcal {R}}_H<1$, it follows that

$$\begin{aligned} \begin{aligned} \left| \frac{(b_S-\mu )}{1-{\mathcal {P}}_1(c)}\right|<&\left| \lambda + \frac{(b_S-\mu )}{1-{\mathcal {P}}_1(c)}\right| \\ =&\left| (b_S-\mu )\left( 1+\frac{(b_I({{\hat{c}}}) + {{\hat{c}}})\int _0^{+\infty } e^{-\lambda a}w(a,{{\hat{c}}})da}{1-{\mathcal {P}}_1(c)}\right) \right| \\ \le&\left| (b_S-\mu )\right| \left( 1+\frac{\left| (b_I({{\hat{c}}}) + {{\hat{c}}})\int _0^{+\infty } e^{-\lambda a}w(a,{{\hat{c}}})da\right| }{1-{\mathcal {P}}_1(c)}\right) \\<&|(b_S-\mu )|\left( 1+\frac{{\mathcal {P}}_2({{\hat{c}}})}{1-{\mathcal {P}}_1(c)}\right) \\ =&|(b_S-\mu )|\frac{\left( 1-{\mathcal {P}}_1(c)+{\mathcal {P}}_2({{\hat{c}}})\right) }{1-{\mathcal {P}}_1(c)}<\left| \frac{(b_S-\mu )}{1-\mathcal P_1(c)}\right| . \end{aligned} \end{aligned}$$

This contradiction shows that the characteristic equation (90) has no roots with positive real parts and the boundary steady state $E_1^*$ is locally asymptotically stable if ${\mathcal {R}}_H<1$. Otherwise, the boundary steady state $E_1^*$ is unstable. $\square$