Dynamics of an age-structured heroin transmission model with vaccination and treatment

1.   Introduction

Heroin is a highly abused opioid and incurs a significant detriment to society worldwide. Heroin usually appears as a white or brown powder or as a black sticky substance, known as "black tar heroin" ^[1], and its most frequent routes of delivery were intravenous injection (25%) and inhalation ^[12]. It crosses the blood-brain barrier within 15–20 seconds, rapidly achieving a high level syndrome in the brain and the central nervous system which causes both the 'rush' experienced by users and the toxicity ^[25]. Heroin users are at high risk for addiction. It is estimated that about 23% of individuals who use heroin become dependent on it.

More recently, a heroin conjugate vaccine attracted much attentions. It was developed through comprehensive evaluation of hapten structure, carrier protein, adjuvant and dosing, which can generate a significant and sustained antidrug IgG titers in each subject and it is effective in rhesus monkeys ^[3]. Also, it is found that immunization of mice with an optimized heroin-tetanus toxoid (TT) conjugate can reduce heroin potency by $>15\%$ and the vaccine effects proved to be durable and persisting for over eight months. Although it is unknown what will happen if the heroin vaccine is used in clinical setting, the heroin vaccine brings much hope for the defence against and control of heroin abuse.

In fact, the spread of heroin habituation and addiction can be well modeled by epidemic type models as "transmission" occurs in the form of peer pressure where established users recruit susceptible individuals into trying and using the drug ^[4,9,16], that is, mathematical modelling is a means to provide a general insight for how classes of drug takers behave, and as such, could hopefully becomes a useful device to aid specialist teams in devising treatment strategies. Modeling heroin addiction and spread in epidemic fashion is not new ^[21]. Recently, Fang et al. ^[10] proposed a age-structured heroin transmission model and proved its global dynamics behaviors. Usually, the population is divided into three classes, namely the number of susceptibles, $S(t)$, the number of drug users not in treatment, $U_1(t)$ and the number of drug users in treatment, $U_2(t)$, respectively. Naturally, we wonder how the heroin vaccine effects the heroin transmission process.

The study of vaccination has been the subject of intense theoretical analysis ^{[2,7,8,14,17,18,19,20,27,29]}. Based on the study of classical epidemic models, Kribs-Zaleta and Velasco-Hern$\acute{a}$ndez ^[17] added a compartment $V$ into an $SIS$ model and studied the vaccination of disease such as pertussis and tuberculosis; Kribs-Zaleta and Martcheva ^[18] studied the effects of a vaccination campaign upon spread of a non-fatal disease which features both acute and chronic infective stages, as well as variable infectivity and recovery rates in the chronic stage. Interestingly, Xiao and Tang ^[29] developed a simple SIV epidemic model including susceptible, infected and imperfectly vaccinated classes, with a nonlinear incidence rate and there would be backward bifurcations; Arino et al. ^[2] also showed that, if vaccines are imperfect, i.e., vaccinated individuals can be infected, there could be backward bifurcations.

To study the role of the heroin vaccine in the control of heroin abuse, adding a compartment $V(t)$ into a heroin transmission model is necessary. To our best knowledge, there is no related work in this field by now. Hopefully, mathematical modeling can provide a new insight into the interaction mechanism between the vaccinated, the susceptibles, the drug users and the individuals in treatment.

Motivated by the development of the heroin vaccine, in this paper, we present an age structured heroin transmission $SU_1U_2V$ model which incorporates a drug reuse rate $\alpha(a)$ dependant on treat-age (i.e., the time since the host has been in treatment) and a vaccine waning rate dependant on vaccination age (i.e., the time since the host has been vaccinated). We also assume the susceptible population is vaccinated at a constant rate $\psi$, the vaccinated individuals can be infected at reduced rate $\sigma\beta$ with $0\leq\sigma\leq 1$. Obviously, $\sigma = 0$ means the vaccine is completely effective in preventing infection, while $\sigma = 1$ means that the vaccine is utterly ineffective. As a result, the system is as follows:

for $t>0$, with the initial and boundary conditions:

where $U_{20}(\theta), V_{0}(a)\in L^1_+(0, \infty)$.

In system (1)-(2), $\theta$ is the treat-age, that is the time that has elapsed since a drug user is in treatment; $U_2(\theta, t)$ is the density of drug users in treatment with age $\theta$ at time $t$; $S(t)$ is the density of the susceptibles at time $t$; $U_1(t)$ is the density of drug users not in treatment, initial and relapsed drug users. The positive constant $\Lambda$ is the recruitment of susceptible, $\mu$ the natural death rate of the general population, $\beta$ is the force of drug use per contact with the susceptible per unit time, $\sigma\beta$ is the force of drug use per contact with the vaccinated individuals per unit time, $\gamma$ is the rate of drug users who enter treatment, $\delta_1$ a removal rate that includes drug-related deaths of users not in treatment and a spontaneous recovery rate, individuals not in treatment who stop using drugs but are no longer susceptible, $\delta_2$ a removal rate that includes the drug-related deaths of users in treatment and a rate of successful "cure" that corresponds to recovery to a drug free life and immunity to drug addiction for the duration of the modelling time period. The function $p(\theta)$ is the probability of a drug user in treatment with the treatment-age $\theta$ relapsing to the untreated users.

Throughout the paper, we make the following assumptions: (A1) there is a positive number of drug users not in treatment, i.e., $U_{10}>0$; (A2) the initial conditions $U_{20}(\theta)$ and $V_{0}(a)$ are uniformly bounded respectively for $\theta, a\in(0, +\infty)$; (A3) the maps $\theta\rightarrow p(\theta)$ and $a\rightarrow \alpha(a)$ are almost everywhere bounded and belong to $L^{\infty}_+((0, +\infty), \mathbb{R})\setminus\{0_{L^{\infty}}\}$.

The paper is organized as follows. In the next section, we present some preliminary results of system (1)-(2). In Section 3, we prove the local and global stability of the drug-free steady state of system (1). In Section 4, we present the existence and the stability results of the drug spread steady state of system (1) when the basic reproduction number is larger than one. Finally, in Section 5, a brief discussion and some numerical examples are presented.

2.   Preliminary results

In this section, we give some basic results prepared for the further study of system (1).

For the sake of convenience, we let

By simple calculations, we have that $\mathscr{K} = 1-\mu \mathscr{K}_1$ and $\mathscr{\phi} = \gamma-\gamma(\mu+\delta_2) \mathscr{\phi}_1$.

Naturally, system (1)-(2) has a unique drug-free steady state $E^0(S^0, 0, 0, V^0(a))$ which satisfies that

Solving the last two equations, we have that

Substituting Equation (4) into the first equation of (3), we have that

Thus, we have that

According to the definition of the basic reproduction number in existing literatures ^[5,6,28], we define the basic reproduction number $\mathcal R_0$ as:

According to ^[24], any positive equilibrium ($S^*$, $U_1^*$, $U_2^*(\theta)$, $V^*(a)$) of system (1), if it exists, must be a constant solution of the following equations

Let

By simple calculations, we have that

From the third and the forth equation of (8), we have that

From the fifth and the sixth equation of (8), we have that

Substituting (11) and (12) into the second and the first equation of (8) respectively, we have that the following equations

It follows from the first equation of (13), we have that

Substituting (14) into the second equation of (13) and eliminating $U_1^*$, we have that

Define a function $\mathcal F(U_1^*)$ to be the right hand side of Equation (15). Obviously, $\mathcal F(0) = \mathcal R_0$. It follows from (9) we have that $\mathcal F(U_1^*)$ is a decreasing function of $U_1^*$ and $\mathcal F(U_1^*)\rightarrow0$ as $U_1^*\rightarrow\infty$. Thus, there exists an unique positive root of Equation (15) only if $\mathcal R_0>1$.

Summarizing the above discussion, we have the following result.

Theorem 2.1. System (1) always has a drug-free steady state $E^0$, besides that, it has an unique drug spread steady state $E^*$ only if $\mathcal R_0>1$.

Now, by setting

we deduce from (1) that $N(t)$ satisfies the following ordinary differential equation:

and therefore $\limsup_{t\rightarrow \infty}\limits N(t)\leq\frac{\Lambda}{\mu}$. Denote

Then $\Omega$ is the maximum positively invariant set of system (1) that attracts all positive solutions of of (1). Therefore, we restrict our attention to solutions of (1) with initial conditions in $\Omega$.

In the following, we use the approach introduced by Thieme ^[26]. Consider

and

Let the linear operator $A:Dom(A)\subset\mathcal {X}\rightarrow\mathcal {X}$ defined by

with

where $W^{1, 1}$ is a Sobolev space. Then $\overline{Dom(A)} = \mathcal {X}_0$ is not dense in $\mathcal {X}$. We consider a nonlinear map $\mathcal {F}:\overline{Dom(A)}\rightarrow \mathcal {X}$ which is defined by

and let

Then, we can reformulate system (2.3) as the following abstract Cauchy problem:

By applying the results in Hale ^[11], Magal ^[22], and Magal and Thieme ^[23], we obtain the following theorem.

Theorem 2.2. System (1) generates a unique continuous semiflow $\{\mathcal U(t)\}_{t\geq0}$ on $\mathcal X_{0+}$ that is bounded dissipative and asymptotically smooth. Furthermore, the semiflow $\{\mathcal U(t)\}_{t\geq0}$ has a global compact attractor $\mathcal A$ in $\mathcal X_{0+}$, which attracts the bounded sets of $\mathcal X_{0+}$.

3.   The stability of the drug-free steady state

In this section, by use of characteristic equation, we will prove the local and global stability of the drug-free steady state $E^0$.

For the sake of convenience, we give the following Laplace transforms of the corresponding functions

By direct calculations, we have the following relationships

Theorem 3.1. The drug-free steady state $E^0$ is locally asymptotically stable if $\mathcal{R}_0<1$, and unstable if $\mathcal{R}_0>1$.

Proof. By linearization of system (1) at the drug-free steady state $E^0$, we can obtain the corresponding linearized system. We let

By linearization of system (1) at the drug free steady state $E^0$, we obtain the following system

To analyze the asymptotic behaviors around $E^0$, we let

where $\overline{x}$, $\overline{y}$, $\overline{z}(\theta)$ and $\overline{w}(a)$ can be determined. Thus, we consider the following eigenvalue problem

Solving the third equation of (21), we have

Solving the fifth equation of (21), we have

Substituting (23) and (22) into the first and the second equation of (21), we have the characteristic equation

where

Then the roots of Equation (24) are determined by the following equations

and

Obviously, $\lambda = -\mu$ is the root of Equation (25). Then we need only to consider the root of Equation (26) which can be rewritten as

We also have that

Define a function $\mathcal{H}(\lambda)$ to be the right-hand side of Equation (28). It follows from the definition of the basic reproduction number $\mathcal{R}_0$ (see (7)), we have that $\mathcal{H}(0) = \mathcal{R}_0$. By direct computing, it is easy to show that $\mathcal{H}'(\lambda)<0$, that is, $\mathcal{H}(\lambda)$ is a decreasing function of $\lambda$ with $\lim_{t\rightarrow+\infty}\limits\mathcal{H}(\lambda) = 0$.

Assume that $\lambda = x+iy$ is a root of Equation (28). Then it follows from (28) that

Thus, we can have that $\Re(\lambda)$ is negative if $\mathcal{R}_0<1$, and therefore the steady state $E^0$ is locally asymptotically stable if $\mathcal{R}_0<1$ and it is unstable if $\mathcal{R}_0>1$.

In the following, we will use the Fluctuation Lemma to establish the global stability of the drug-free steady state $E^0$. To this end, we first introduce the notation

Then the Fluctuation Lemma is given as follows.

Lemma 3.2. (Fluctuation Lemma ^[13]) Let $g:\ \mathbb{R}_+\rightarrow\mathbb{R}$ be a bounded and continuously differentiable function. Then there exist sequences $\{s_n\}$ and $\{t_n\}$ such that $s_n\rightarrow\infty$, $t_n\rightarrow\infty$, $g(s_n)\rightarrow g_\infty$, $g^\prime(s_n)\rightarrow 0$, $g(t_n)\rightarrow g^\infty$ and $g^\prime(t_n)\rightarrow 0$ as $n\rightarrow\infty$.

Lemma 3.3. ^[15] Suppose $f:\ \mathbb{R}_+\rightarrow\mathbb{R}$ be a bounded function. Then

where $\|h\|_1 = \int_0^\infty h(s)ds$.

Using integration, $U_2(\theta, t)$ and $V(a, t)$ satisfy the following Volterra formulation:

Theorem 3.4. If $\mathcal{R}_0<1$, then the drug-free steady state $E^0$ is the unique steady state of system (1), and it is globally stable.

Proof. Theorem 3.1 shows that the drug-free steady state $E^0$ of system (1)) is locally stable if $\mathcal{R}_0<1$. To use the Fluctuation Lemma, substituting the expressions of $V(a, t)$ and $U_2(\theta, t)$ into the first two equations of system (1), we have that

where

with $\lim_{t\rightarrow \infty}\limits F_{V}(t) = 0$, $\lim_{t\rightarrow \infty}\limits F_{U}(t) = 0$ and $\lim_{t\rightarrow \infty}\limits F_{UV}(t) = 0$.

Choose the sequences $t_n^1\rightarrow \infty$ such that $S(t_n^1)\rightarrow S^\infty$ and $S'(t_n^1)\rightarrow 0$. Then $F_V(t)\rightarrow 0$ as $n\rightarrow \infty$. With the assistance of the Fluctuation Lemma, it follows from the first equation of (31) we have that

and

Choose the sequences $t_n^2\rightarrow \infty$ such that $U_1(t_n^2)\rightarrow U_1^\infty$ and $U_1'(t_n^2)\rightarrow 0$. Then $F_{U}(t)\rightarrow 0$ and $F_{UV}(t)\rightarrow 0$ as $n\rightarrow \infty$. With the assistance of the Fluctuation Lemma, it follows from the second equation of (31) we have that

Thus we obtain that $U_1^\infty\rightarrow0$ if $\mathcal{R}_0<1$. It then follows from (29) we have that $\lim_{t\rightarrow \infty}\limits U_2(\theta, t) = 0$.

Choose the sequences $t_n^3\rightarrow \infty$ such that $S(t_n^3)\rightarrow S_\infty$ and $S'(t_n^3)\rightarrow 0$. Note that $\lim_{n\rightarrow\infty}\limits U_1(t_n^3) = 0$ and $\lim_{n\rightarrow\infty}\limits F_V(t_n^3) = 0$. It then follows from the first equation of (31) we have that

It follows that

which implies that $\lim_{t\rightarrow \infty}\limits S(t) = \frac{\Lambda}{\mu+\psi-\psi\mathscr{K}}$. It follows from (30) we have that

Therefore, $(S, U_1, U_2, V)\rightarrow E^0$ in $\mathbb{R}_+\times\mathbb{R}_+ \times L_+^1\times L_+^1$ as $t \rightarrow \infty$. This completes the proof of Theorem 3.4.

4.   Stability of the drug spread steady state

This section aims to establish the stability of the drug spread steady state of system (1) in terms of the basic reproduction number $\mathcal R_0$.

By a similar discussion as Theorem 3.5 in ^[10], we have the uniform persistence result as following.

Theorem 4.1. Suppose the heroin spread is initially present, i.e., $U_{10}>0$. If $\mathcal R_0>1$, then the semiflow generated by system (1) is uniformly persistent, i.e., there exists $\varepsilon>0$ which is independent of initial values such that

and $\liminf_{t\rightarrow\infty}\|V(\cdot, t)\|_{L_+^1}\geq\varepsilon.$

For the sake of convenience, we let

It then follows that

In the following, we try to study the stability of the drug spread steady state $E^*$ ($S^*$, $U_1^*$, $U_2^*(\theta)$, $V^*(a)$) of system (1). We let

and $V(a, t) = \widetilde{V}(a, t)+V^*(a)$. By linearization of system (1) at the steady state $E^*$, we obtain the following system

To analyze the asymptotic behaviors around $E^*$, we let

where $\overline{x}$, $\overline{y}$, $\overline{z}(\theta)$ and $\overline{w}(a)$ can be determined. Then, we consider the following eigenvalue problem

Solving the third equation of (35), we have

Solving the fifth equation of (35), we have

where

Substituting (37) into the first equation of (35), we have

where

It follows from (33) we have that

Substituting (36) and (37) into the second equation of (35), by use of $\sigma\beta\psi S^*\mathscr{K}_\pi = (\mu+\delta_1+\gamma)-\beta S^*-\phi$ (i.e., the second equation of (13)), we have that

Substituting (38) into the above equation and dividing both sides by $\overline{y}$ leads to the characteristic equation

It follows from Equation (39) we have that

Due to the fact that the characteristic equation (39) is too complex, it is very difficult to determine the distribution of the eigenvalues. In the following, we will study the stability of the drug spread steady state $E^*$ in three special cases of system (1) respectively.

Case (ⅰ) We assume that $\alpha(a) = 0$, that is, the vaccine does not wane once the host is vaccinated. which is reasonable since the heroin does not change the toxic (pathogenic) substance. Mathematically, letting $\int_0^\infty V(a, t)da = V(t)$ in this case, we can rewrite system (1) as

Then we have the following result.

Lemma 4.2. If $\alpha(a)\equiv0$, the drug spread steady state $E^*$ of system (41) is locally asymptotically stable if it exists.

Proof. Linearizing system (41) at its positive steady state $(S^*, U_1^*, U_2^*(\cdot), V^*)$, we have the following characteristic equation

By simple calculations, we have the following

where

Obviously, if $\sigma\neq0$, $\lambda = -(\mu+\psi+\beta U_1^*)$ and $\lambda = -(\mu+\sigma\beta U_1^*)$ are not the roots of Equation (43). Dividing both side of Equation (43) by $B_1(\lambda)B_2(\lambda)$ leads to

where

Now, assume that $\lambda$ is a root of Equation (44) with $\Re(\lambda)\geq 0$. If we can prove the real part of $Z$ is positive (see Appendix A for its detailed proof), then by using $|\widehat{\Phi}(\lambda)|\leq\phi$, we have

This contradiction implies that all roots of Equation (44) have negative real parts. Hence, the positive steady state $(S^*, U_1^*, U_2^*(\cdot), V^*)$ of system (41) is locally stable if it exists.

Case (ⅱ) We assume that $\sigma = 0$, i.e., the heroin vaccine can provide a prefect protection for the vaccinated individuals to avoid the heroin drug addiction. It then follows that $\pi(a) = k_1(a)$. Mathematically, in this case, system (1) can be rewritten as

Lemma 4.3. If $\sigma = 0$, the drug spread steady state $E^*$ of system (45) is globally asymptotically stable if it exists.

Proof. It then follows from (18) and (19) we have that Equation (40) can be modified as

Since $\lambda = 0$ is not the roots of Equation (46). Both sides of Equation (46) are divided by $\left(\lambda+\phi-\widehat{\Phi}(\lambda)\right)$, we obtain

where

Assuming $\lambda$ is a root of Equation (47), we can prove that $\Re(\lambda)$ (i.e., the real part of $\lambda$) is negative. Supposed $\Re(\lambda)$ is nonnegative, then the real part of $Z$ is nonnegative (see Appendix B). It follows from (10) and (47) we have that

which is a contradiction. So, all roots of Equation (47) have negative real part, and therefore all roots of Equation (46) have negative real part and the drug spread steady state $E^*$ of system (45) is locally asymptotically stable if it exists.

Based on the persistence results in Theorem 4.1 and the local stability results of $E^*$ of system (45), we will use a suitable Lyapunov functional to prove the global stability of $E^*$.

Set $g(x) = x-1-\ln x$, for $x\in\mathbb R_+$. The function $g(x)$ has a global minimum at $x = 1$ with $g(1) = 0$. For our presentation here, we define

Note that $\varepsilon_1(\theta)$, $\varepsilon_2(a)>$ for all $\theta>0$ and $a>0$ respectively. We can easily check that $\varepsilon_1(0) = \phi/\gamma$, $\varepsilon_2(0) = \int_0^\infty\alpha(a)V^*(a)da$ and

Now, we define the following Lyapunov functional

where

Then $W$ is bounded. Then we calculate the time derivative of $W(t)$ along with the solutions of system (45). Following the results in the proof of Theorem 2.2 in ^[10] and the proof of Theorem 3.11 in ^[30], we have that

By adding $\frac{d W_{S}(t)}{dt}$ and $\frac{d W_{V}(t)}{dt}$ together, though some simple calculations, we have that

Combining the above four compartments of the Lyapunov functionals, through some simple calculations, we obtain

Notice that equality holds only if $S(t) = S^*$, $U_1(t) = U_1^*$ and $U_2(\theta, t) = U_2^*(\theta)$. Thus we conclude that the largest positive invariant is the singleton $\{E^*\}$. By Lyapunov-LaSalle invariance principle, we conclude that the drug spread steady state $E^*$ is globally asymptotically stable when it exists.

Case (ⅲ) We assume that $\sigma = 1$, i.e., the heroin vaccine is noneffective and cannot provide any protection from heroin drug addicted, the vaccination makes no sense. In this case, mathematically, the compartments $S$ and $V$ can be combined, system (1) can be rewritten as the system in ^[10] (see Appendix C for more details) and the drug spread steady state $E^*$ is locally and globally stable.

5.   Numerical simulation and discussion

In this paper, we have studied an age structured heroin transmission model with treatment and vaccination, in which the vaccination can only provide an imperfect protection and the vaccinated wanes the protection as vaccination age goes.

The basic reproduction number $\mathcal R_0$ of our system (1) has been found by the definition. When $\mathcal R_0<1$, system (1) has only the drug free steady state $E^0$ and it is globally asymptotically stable, which implies that the heroin drug will die out eventually. Meanwhile, from the expression of $\mathcal R_0$, we find that the vaccination, although it is imperfect, plays an important role in the control of heroin spread. When $\mathcal R_0>1$, system (1) has only the drug spread steady state $E^*$ and it is uniformly persistent provided that the heroin spread is initially present. Due to the fact that the characteristic equation of system (1) at the drug spread steady state is very complex, it is difficult to discuss the distribution of its eigenvalues. From the biology angle, we have recast system (1) into three special cases and obtained the local and global stability of the drug spread steady state $E^*$ if it exists (see Lemmas 4.2-4.3).

Recalling that the expression of $\mathcal R_0$ in (7), it follows from $0\leq\sigma\leq1$ we have that

It implies that $\mathcal R_0(\psi)\leq \mathcal R_0(0)$, i.e., the vaccination plays an important role in the basic reproduction number which can reduce the reproduction number, although the vaccine provides an imperfect protection. Thus, heroin vaccine will definitely benefit the people.

In the following, we will present some numerical simulations to study the dynamic behaviors of system (1) under the condition that the basic reproduction number is lager than one, i.e., $\mathcal R_0>1$. Before that, for simplicity, we take one month as the unit time. Note that the function $\theta\rightarrow p(\theta)$ and $a\rightarrow\alpha(a)$ are both almost everywhere bounded and belong to $L^{\infty}_+((0, +\infty), \mathbb{R})\setminus\{0_{L^{\infty}}\}$. In this section we assume that the vaccine waning rate of the vaccinated individuals is

and the drug reuse rate of the individuals in treatment is

for $\theta\in[0, \overline{\theta}]$, where $\overline{a}$, $\overline{\theta}$ are the maximum value of vaccination age and treat-age respectively. For simplify, we adopt that $\overline{a} = \overline{\theta} = 50$ (months).

To study the stability of the drug spread steady state of system (1), we adopt the other parameters in system (1) as follows

It follows from the expression of the basic reproduction number $\mathcal{R}_0$ that $\mathcal{R}_0 = 1.0117>1$. By use of the parameter values adopted in (54)-(56) and appropriate initial conditions, we will perform some numerical simulations with the help of Matlab. The numerical simulations show that the drug spread steady state $E^*$ is asymptotically stable if it exists (see Figure 1).

Furthermore, we want to illustrate that the stability of the drug spread steady state is not dependent on the initial conditions by numerical simulations. Let

We obtain that $\mathcal{R}_0 = 7.6102>1$ and simulate the solutions of system (1) under four pairs of initial values (see Figure 2). The numerical simulation results show that the stability of the drug spread steady state is not dependent on the initial conditions. In this case, we may conjecture that the drug spread steady state is globally asymptotically stable whenever it exists.

Appendix A

In the course of the proof of Lemma 4.2, we let $\Re(\lambda)\geq0$ and want to prove that $\Re(Z)>0$, where

To study the sign of $\Re(Z)$, we let $\lambda = x + iy$ with $x\geq0$ which is assumed nonnegative in the proof of Lemma 4.2. Substituting $\lambda = x + iy$ into $(A.1)$, by some calculations, we can have that

where

Summing the above three terms, we have that

where

It then follows that the real part of $Z$ is positive, i.e., $\Re(Z)>0$.

Appendix B

To consider the roots of Equation (47), we let $\Re(\lambda)\geq0$ and want to prove that $\Re(Z_1)>0$, where

Let $\lambda = x + iy$ with $x>0$. By substituting $\lambda = x + iy$ into $(B.1)$, we have that

and

where $E_1 = x+\phi-\int_0^\infty\Phi(\theta)e^{-x\theta}\cos (y\theta) d\theta$ and $E_2 = y+\int_0^\infty\Phi(\theta)e^{-x\theta}\sin (y\theta) d\theta$. It follows from

we have that $\Re(Z_1)$ is nonnegative.

Appendix C

If $\sigma = 1$. Set $V(t) = \int_0^\infty V(a, t)da$. It follows from the last equation of system (1) that we have

Denote that $\widehat{S}(t): = S(t)+V(t)$. By dropping the hat, we have that

Then system (1) can be rewritten as the main system in ^[10] and the drug spread steady state $E^*$ is locally stable.

Acknowledgments

We would be very grateful to anonymous referees for their comments and suggestions that helped to improve this paper. This work is supported partially by China Postdoctoral Science Foundation 2017M621523; X. Li is supported partially by the National Natural Science Foundation of China (11771017); M. Martcheva is supported partially through grant DMS-1515661. Part of this work was done when XD was a visiting scholar at the Department of Mathematics, University of Florida. XD would like to thank the Department for kind hospitality he received there.

Conflict of interest

The authors declare there is no conflict of interest.