Global Dynamics of a New Huanglongbing Transmission Model with Quarantine Measures

Yujiang Liu, Chunmei Zeng, Jing Guo, Zhenzhen Liao

Key Laboratory of Jiangxi Province for Numerical Simulation and Emulation Techniques, Gannan Normal University, Gan-zhou, China.

DOI: 10.4236/am.2022.131001                

Abstract

An epidemic model which describes Huanglongbing transmission is proposed with the goal of investigating the effect of quarantine measures on the spread of diseases. First of all, the analytical formula for the basic reproduction number R₀ is obtained by the means of next generation matrix, and the existence of disease-free equilibrium and endemic equilibrium is discussed. Then, the local stability and the global stability of equilibria are investigated by using Routh-Hurwitz criterion and Lyapunov function, respectively. Numerical simulations indicate that comprehensive quarantine measures can effectively control the spread of Huanglongbing. It provides a reliable tactic basis for preventing the epidemic outbreak.

Keywords

Huanglongbing, Basic Reproduction Number, Sensitivity Analysis, PRCC, Global Stability

1. Introduction

The citrus Huanglongbing (HLB), allied with the bacterium Candidatus Liberibacter asiaticus (CLas), poses a great threat to the industry of citrus worldwide [1] [2]. The most vital method of disease breadth occurs by the citrus psyllids that serve as vectors [1]. The citrus psyllids consume sap from an infected host through their stylets, the viruses in the sap enter the salivary glands, circulate within the psyllid, and then cause infection between citrus trees and psyllids [3] [4] [5] [6] [7]. The citrus tree infected HLB shows yellow shoots, leaves with blotchy mottle, small lopsided fruits, overall stunting of the trees, in server cases, and ultimately the death of the entire tree [8] [9] [10]. In 2020, citrus Huanglongbing was included in the List of Crop pests and diseases by the Ministry of Agriculture and Rural Affairs [11].

There are no economically viable curative methods for citrus trees infected HLB due to absence of resistant or tolerant commercial citrus varieties [12]. However, a large number of researchers have studied this disease biologically and come to describe the dynamics of disease mathematically [1] [13]. Recent works have provided different insights into the dynamic mechanism of HLB propagation through differential equation modeling, statistical analysis and individual-based modeling [10]. Jacobsen et al. established a 6-dimensional dynamic model within a citrus orchard and analyzed how the number of each class changes over time due to bacterial transmission between citrus trees and psyllids [7]. Parry et al. fitted a spatially explicit disease model, using specialized statistical methods and software, and then obtained the effect of tree age on transmission parameters, which could be used to predict the spread of early epidemics [14]. Lee et al. [2] conducted a mathematical model that combined experiments and individual-based. It emphasized the need to keep the psyllid populations as low as possible due to the whole grove to be infected before the first symptoms appear on any tree. In fact, removing infected trees and quarantining them are also important measures to prevent as many trees from becoming infected. Specially, if HLB is not yet present, quarantine measures should be enforced to keep it out [15]. The quarantine measures, which began on June 17, 2010, were designed to prevent the transportation of citrus in Florida [7]. Therefore, motivated by these researches, we propose a dynamic model of HLB with quarantine measures for exposed and infected citrus trees, and then perform some theoretical analysis of its properties.

The organization of this paper is as follows: the model is developed and the feasible region, the basic reproduction number and the equilibria of the model are obtained in Section 2. By the Routh-Hurwitz criterion, the local stability of disease-free equilibrium and endemic equilibrium are proved in Section 3. By constructing a suitable Lyapunov function, the global stability of the equilibria is studied in Section 4. In addition, numerical simulations are illustrated in our theoretical results in Section 5 and a brief conclusion is given in the last section.

2. Model and Preliminaries

The dynamics of the model with a fraction of quarantined citrus trees involve citrus-psyllid interactions, therefore, the citrus tree population given as $N_h\left(t\right)$ is divided into susceptible host $S\left(t\right)$, exposed host $E\left(t\right)$, infected host $I\left(t\right)$, quarantined host $Q\left(t\right)$ and removed host $R\left(t\right)$, at any time t. Thus, $N_h\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+Q\left(t\right)+R\left(t\right)$. The citrus psyllid population denoted by $N_v\left(t\right)$, is divided into susceptible class $X\left(t\right)$ and infected class $Y\left(t\right)$ at any time t. Thus, $N_v\left(t\right)=X\left(t\right)+Y\left(t\right)$. We assume that the citrus population grows at a recruitment rate ${\Lambda }_h$, moreover, $q_1$ and $q_2$ represent the rate of the infectious convert to the quarantined host and the rate of exposed to the quarantined host, respectively. According to the principle of the compartmental model, the model with first order nonlinear ordinary differential equations is constructed as follows:

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}={\Lambda }_h-{\beta }_{1}SY-{\mu }_h{S}_h\mathrm{,}\hfill \\ \frac{\text{d}E\left(t\right)}{\text{d}t}={\beta }_{1}SY-\left(q_1+\sigma +{\gamma }_1+{\mu }_h\right)E\mathrm{,}\hfill \\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\sigma E-\left(q_2+{\gamma }_2+{\mu }_h+d\right)I\mathrm{,}\hfill \\ \frac{\text{d}Q\left(t\right)}{\text{d}t}=q_{1}E+q_{2}I-\left({\gamma }_3+{\mu }_h\right)Q\mathrm{,}\hfill \\ \frac{\text{d}R\left(t\right)}{\text{d}t}\mathrm{<mo>\; =\; </mo>}{\gamma }_{1}E+{\gamma }_{2}I+{\gamma }_{3}Q-{\mu }_{h}R\mathrm{,}\hfill \\ \frac{\text{d}X\left(t\right)}{\text{d}t}={\Lambda }_v-{\beta }_2\left(\theta E+I\right)X-{\mu }_{v}X\mathrm{,}\hfill \\ \frac{\text{d}Y\left(t\right)}{\text{d}t}={\beta }_2\left(\theta E+I\right)X-{\mu }_{v}Y\mathrm{,}\hfill \end{array}$ (2.1)

where the parameters involved in the model (2.1) and their values are showed in Table 1. The schematic diagram for the flow of the dynamics is depicted in Figure 1. In the subsequent section, we will explore the basic qualitative properties of the model (2.1).

2.1. Feasible Region Denote

${\Omega }_0=\left\{\left(S\mathrm{,}E\mathrm{,}I\mathrm{,}Q\mathrm{,}R\mathrm{,}X\mathrm{,}Y\right)\in {ℝ}_{+}^7\mathrm{<mo>\; |\; </mo>}S+E+I+Q+R\le \frac{{\Lambda }_h}{{\mu }_h}\mathrm{,}X+Y\le \frac{{\Lambda }_v}{{\mu }_v}\right\}\mathrm{.}$ (2.2)

Lemma 1. The feasible region is given by ${\Omega }_0$ is a positively invariant and globally attractive set for system (2.1).

Proof. Clearly, the solution set $\left(S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}Q\left(t\right)\mathrm{,}R\left(t\right)\mathrm{,}X\left(t\right)\mathrm{,}Y\left(t\right)\right)$ of the system (2.1) is positive.

Since $I\left(t\right)$ is non-negative, by adding the first to fifth equations and the last two equations of the system (2.1), respectively, the total number of citrus tree $N_h\left(t\right)$ and psyllid vector $N_v\left(t\right)$ satisfy

Figure 1. The schematic diagram of the model showing transitions to different categories for trees and psyllids.

Table 1. Parameters in the model (2.1).

${\stackrel{\dot{}}{N}}_h\left(t\right)\le {\Lambda }_h-{\mu }_h{N}_h\left(t\right),{\stackrel{\dot{}}{N}}_v\left(t\right)={\Lambda }_v-{\mu }_v{N}_v\left(t\right)$ (2.3)

From Equation (2.3), we get

$\begin{array}{l}{N}_h\left(t\right)\le \frac{{\Lambda }_h}{{\mu }_h}-\left(\frac{{\Lambda }_h}{{\mu }_h}-N_h\left(0\right)\right)\exp \left\{-{\mu }_{h}t\right\},\\ N_v\left(t\right)=\frac{{\Lambda }_v}{{\mu }_v}-\left(\frac{{\Lambda }_v}{{\mu }_v}-N_v\left(0\right)\right)\exp \left\{-{\mu }_{v}t\right\}.\end{array}$ (2.4)

This implies that

$\underset{t\to +\infty }{\lim \sup }{N}_h\le \frac{{\Lambda }_h}{{\mu }_h},\underset{t\to +\infty }{\lim \sup }{N}_v=\frac{{\Lambda }_v}{{\mu }_v}.$ (2.5)

Thus, all solutions in $\Omega$ are uniformly bounded. This completes our proof.

2.2. The Basic Reproduction Number

Since the fourth and fifth equation is decoupled from the rest in system (2.1) and to better organize the analysis, we denote ${\mu }_1=q_1+\sigma +{\gamma }_1+{\mu }_h$, ${\mu }_2=q_2+{\gamma }_2+{\mu }_h+d$, and consider the following equivalent system of (2.1):

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}={\Lambda }_h-{\beta }_{1}SY-{\mu }_h{S}_h\mathrm{,}\hfill \\ \frac{\text{d}E\left(t\right)}{\text{d}t}={\beta }_{1}SY-{\mu }_{1}E\mathrm{,}\hfill \\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\sigma E-{\mu }_{2}I\mathrm{,}\hfill \\ \frac{\text{d}X\left(t\right)}{\text{d}t}={\Lambda }_v-{\beta }_2\left(\theta E+I\right)X-{\mu }_{v}X\mathrm{,}\hfill \\ \frac{\text{d}Y\left(t\right)}{\text{d}t}={\beta }_2\left(\theta E+I\right)X-{\mu }_{v}Y\mathrm{.}\hfill \end{array}$ (2.6)

It is easy to see

${\Omega }_1=\left\{\left(S\mathrm{,}E\mathrm{,}I\mathrm{,}X\mathrm{,}Y\right)\in {ℝ}_{+}^5\mathrm{<mo>\; |\; </mo>}S+E+I\le \frac{{\Lambda }_h}{{\mu }_h}\mathrm{,}X+Y\le \frac{{\Lambda }_v}{{\mu }_v}\right\}$ (2.7)

is a positively invariant and globally attractive set for system (2.6). From now on, we restrict the analysis of system (2.6) to the positive invariant set ${\Omega }_1$.

Clearly, the disease free equilibrium (DFE) of system (2.6) exists and is given by $P_0=\left(\frac{{\Lambda }_h}{{\mu }_h}\mathrm{,0,0,}\frac{{\Lambda }_v}{{\mu }_v}\mathrm{,0}\right)$.

Now, we calculate the basic reproductive number $R_0$ of system (2.6) by using the next generation operator approach [19]. Considering the infectious compartments $x={\left(E\mathrm{,}I\mathrm{,}Y\right)}^{\text{T}}$, let

$F\left(x\right)={\left(F_1\left(x\right)\mathrm{,}{F}_2\left(x\right)\mathrm{,}{F}_3\left(x\right)\right)}^{\text{T}}=\left(\begin{array}{c}{\beta }_{1}SY\\ 0\\ {\beta }_2\left(\theta E+I\right)X\end{array}\right)\mathrm{,}$

and

$V\left(x\right)={\left(V_1\left(x\right)\mathrm{,}{V}_2\left(x\right)\mathrm{,}{V}_3\left(x\right)\right)}^{\text{T}}=\left(\begin{array}{c}{\mu }_{1}E\\ {\mu }_{2}I-\sigma E\\ {\mu }_{v}Y\end{array}\right)$

represent the rate of appearance of new infection and the rates of transfer into and out of the class of infected states, respectively. Then the Jacobian matrixes $\mathcal{F}\left(x\right)$ and $\mathcal{V}\left(x\right)$ of system (2.6) with respect to $P_0$ are given below:

$\mathcal{F}\left(x\right)=\left(\begin{array}{ccc}0& 0& \frac{{\Lambda }_h{\beta }_1}{{\mu }_h}\\ 0& 0& 0\\ \frac{{\Lambda }_v{\beta }_2\theta }{{\mu }_{v}h}& \frac{{\Lambda }_v{\beta }_2}{{\mu }_{v}h}& 0\end{array}\right)\mathrm{,}\mathcal{V}\left(x\right)=\left(\begin{array}{ccc}{\mu }_1& 0& 0\\ -\sigma & {\mu }_2& 0\\ 0& 0& {\mu }_v\end{array}\right)$ (2.8)

Using the next generation matrix theory, the expression of the basic reproduction number ${\mathcal{R}}_0$ can be obtained as

${\mathcal{R}}_0=\sqrt{{\mathcal{R}}_0^1+{\mathcal{R}}_0^2},{\mathcal{R}}_0^1=\frac{{\beta }_1{\beta }_2{\Lambda }_h{\Lambda }_v\theta }{{\mu }_h{\mu }_1{\mu }_v^2},{\mathcal{R}}_0^2=\frac{{\beta }_1{\beta }_2{\Lambda }_h{\Lambda }_v\sigma }{{\mu }_h{\mu }_1{\mu }_2{\mu }_v^2},$ (2.9)

Remark 1: ${\mathcal{R}}_0^1$ represents the number of secondary infections by the circular transmission path from E to Y and back to E, and ${\mathcal{R}}_0^2$ indicates the number of secondary infections by the other circular transmission path from I to Y and back to I.

2.3. The Existence of Endemic Equilibrium

Next, we consider the existence of endemic equilibrium. Let the right side of equations of system (2.6) be equal to 0, then we obtain algebraic equations as follows:

$\begin{array}{l}{\Lambda }_h-{\beta }_1{S}^{*}{Y}^{*}-{\mu }_h{S}_h^{*}=0,\\ {\beta }_1{S}^{*}{Y}^{*}-{\mu }_1{E}^{*}=0,\\ \sigma E^{*}-{\mu }_2{I}^{*}=0,\\ {\Lambda }_v-{\beta }_2\left(\theta E^{*}+I^{*}\right)X^{*}-{\mu }_v{X}^{*}=0,\\ {\beta }_2\left(\theta E^{*}+I^{*}\right)X^{*}-{\mu }_v{Y}^{*}=0.\end{array}$ (2.10)

From the third equation of system (2.10), we get

$I^{*}=\frac{\sigma E^{*}}{{\mu }_2}.$ (2.11)

By adding the first second equations and the last two equations of system (2.10) respectively, we find

$S^{*}=\frac{{\Lambda }_h-{\mu }_1{E}^{*}}{{\mu }_h},X^{*}=\frac{{\Lambda }_v}{{\mu }_v}-Y^{*}$ (2.12)

Substituting Equations (2.11) and (2.12) into the second and fifth equation of (2.10), we have

$E^{*}=\frac{{\mu }_h{\mu }_2{\mu }_v^2\left({\mathcal{R}}_0^2-1\right)}{{\beta }_2\left({\mu }_2\theta +\sigma \right)\left({\Lambda }_v{\beta }_1+{\mu }_h{\mu }_v\right)},Y^{*}=\frac{{\mu }_h{\mu }_1{\mu }_2{\mu }_v\left({\mathcal{R}}_0^2-1\right)}{{\beta }_1\left({\Lambda }_h{\beta }_2\left({\mu }_2\theta +\sigma \right)+{\mu }_1{\mu }_2{\mu }_v\right)}.$ (2.13)

If ${\mathcal{R}}_0>1$, for system (2.6) there exists a unique endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}$ in the interior of ${\Omega }_1$. Therefore, from the above calculation, we get the following theorem.

Theorem 2. There is always a disease-free equilibrium $P_0\left(\frac{{\Lambda }_h}{{\mu }_h}\mathrm{,0,0,}\frac{{\Lambda }_v}{{\mu }_v}\mathrm{,0}\right)$, and the unique endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{E}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ of system (2.6) emerges in ${\Omega }_1$ when ${\mathcal{R}}_0>1$ and no endemic equilibrium when ${\mathcal{R}}_0\le 1$.

3. Local Stability of Equilibria

In this section, we will investigate the local stability of the disease-free equilibrim $P_0$ and the endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}$.

Theorem 3. For system (2.6), the disease-free equilibrium $P_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof. The Jacobian matrix of system (2.6) at disease-free equilibrium $P_0\left(\frac{{\Lambda }_h}{{\mu }_h}\mathrm{,0,0,}\frac{{\Lambda }_v}{{\mu }_v}\mathrm{,0}\right)$ is

${\mathcal{J}}_0=\left(\begin{array}{ccccc}-{\mu }_h& 0& 0& 0& -\frac{{\beta }_1{\Lambda }_h{\beta }_1}{{\mu }_h}\\ 0& -{\mu }_1& 0& 0& \frac{{\beta }_1{\Lambda }_h{\beta }_1}{{\mu }_h}\\ 0& \sigma & -{\mu }_2& 0& 0\\ 0& -\frac{{\Lambda }_v{\beta }_2\theta }{{\mu }_{v}h}& -\frac{{\Lambda }_v{\beta }_2}{{\mu }_{v}h}& -{\mu }_v& 0\\ 0& \frac{{\Lambda }_v{\beta }_2\theta }{{\mu }_{v}h}& \frac{{\Lambda }_v{\beta }_2}{{\mu }_{v}h}& 0& -{\mu }_v\end{array}\right)\mathrm{.}$ (3.1)

The characteristic equation of matrix ${\mathcal{J}}_0$ is the following polynomial equation

$\left(\lambda +{\mu }_h\right)\left(\lambda +{\mu }_v\right)\left(\frac{{\beta }_1{\beta }_2{\Lambda }_h{\Lambda }_v\left(\left(\lambda +{\mu }_2\right)\theta +\sigma \right)}{{\mu }_h{\mu }_v}-\left(\lambda +{\mu }_1\right)\left(\lambda +{\mu }_2\right)\left(\lambda +{\mu }_v\right)\right)=0.$ (3.2)

It is obviously that two eigenvalues of ${\mathcal{J}}_0$ are ${\lambda }_1=-{\mu }_h$, ${\lambda }_2=-{\mu }_v$, which are negative. The other three roots are determined by the following cubic equation

${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0,$ (3.3)

where

$\begin{array}{l}{a}_1={\mu }_1+{\mu }_2+{\mu }_v,\\ a_2={\mu }_2\left({\mu }_1+{\mu }_v\right)+{\mu }_1{\mu }_v\left(1-\frac{{\mu }_2\theta }{{\mu }_2\theta +\sigma }{\mathcal{R}}_0^2\right),\\ a_3={\mu }_1{\mu }_2{\mu }_v\left(1-{\mathcal{R}}_0^2\right).\end{array}$

It is obviously that if ${\mathcal{R}}_0<1$, we have $a_1>0,a_2>0,a_3>0$, and

$a_1{a}_2-a_3=\left({\mu }_2+{\mu }_v\right)a_2+{\mu }_1^2{\mu }_2+{\mu }_1^2{\mu }_v\left(1-\frac{{\mu }_2\theta }{{\mu }_2\theta +\sigma }{\mathcal{R}}_0^2\right)+{\mu }_1{\mu }_2{\mu }_v{\mathcal{R}}_0^2>0,$

however $a_3<0$ if ${\mathcal{R}}_0>1$. According to the Routh-Hurwitz criterion, we declare that all eigenvalues of ${\mathcal{J}}_0$ have negative real parts if and only if ${\mathcal{R}}_0<1$. This proof is complete.

Theorem 4. If ${\mathcal{R}}_0>1$, the endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}$ of the system (2.6) is locally asymptotically stable.

Proof. The Jacobian matrix of system (2.6) at endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{E}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is

${\mathcal{J}}^{\mathrm{<mo>\; *\; </mo>}}=\left(\begin{array}{ccccc}-{\mu }_{h1}& 0& 0& 0& -{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}\\ {\beta }_1{Y}^{\mathrm{<mo>\; *\; </mo>}}& -{\mu }_1& 0& 0& {\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}\\ 0& \sigma & -{\mu }_2& 0& 0\\ 0& -{\beta }_2\theta X^{\mathrm{<mo>\; *\; </mo>}}& -{\beta }_2{X}^{\mathrm{<mo>\; *\; </mo>}}& -{\mu }_{v1}& 0\\ 0& {\beta }_2\theta X^{\mathrm{<mo>\; *\; </mo>}}& {\beta }_2{X}^{\mathrm{<mo>\; *\; </mo>}}& {\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}& -{\mu }_v\end{array}\right)\mathrm{,}$ (3.4)

where ${\mu }_{h1}={\beta }_1{Y}^{*}+{\mu }_h$ and ${\mu }_{v1}=-{\beta }_2\theta E^{*}+{\beta }_2{I}^{*}+{\mu }_v$. The characteristic equation of matrix ${\mathcal{J}}^{\mathrm{<mo>\; *\; </mo>}}$ is given by

$\begin{array}{l}\left(\lambda +{\mu }_v\right)(\left(\lambda +{\mu }_{h1}\right)\left(\lambda +{\mu }_1\right)\left(\lambda +{\mu }_2\right)\left(\lambda +{\mu }_{v1}\right)\begin{array}{c}\\ \end{array}\\ -\frac{{\mu }_1{\mu }_2{\mu }_v\left(\lambda +{\mu }_h\right)\left(\theta \lambda +{\mu }_2\theta +\sigma \right)}{{\mu }_2\theta +\sigma })=0.\end{array}$ (3.5)

Clearly, there is a negative eigenvalue ${\lambda }_1=-{\mu }_h$. By ${\beta }_1{S}^{*}=\frac{{\mu }_1{E}^{*}}{Y^{*}}$, ${\beta }_2{X}^{*}=\frac{{\mu }_2{\mu }_v{Y}^{*}}{\left({\mu }_2\theta +\sigma \right)E^{*}}$, the other four eigenvalues satisfy the following quartic polynomial

${\lambda }^4+b_1{\lambda }^3+b_2{\lambda }^2+b_3\lambda +b_4=0,$ (3.6)

where

$\begin{array}{l}{b}_1={\mu }_{h1}+{\mu }_1+{\mu }_2+{\mu }_{v1},\\ b_2={\mu }_{h1}\left({\mu }_1+{\mu }_2+{\mu }_{v1}\right)+{\mu }_2\left({\mu }_1+{\mu }_{v1}\right)+{\mu }_1\left(\frac{{\mu }_v\sigma }{{\mu }_2\theta +\sigma }+{\beta }_2\theta E^{*}+{\beta }_2{I}^{*}\right),\\ b_3={\mu }_{h1}{\mu }_2\left({\mu }_1+{\mu }_{v1}\right)+{\mu }_1{\mu }_2\left({\beta }_2\theta E^{*}+{\beta }_2{I}^{*}\right)+{\mu }_1{\mu }_v\left(\frac{{\mu }_h\sigma }{{\mu }_2\theta }+{\beta }_1{Y}^{*}\right),\\ b_4={\mu }_1{\mu }_2\left({\mu }_{h1}{\mu }_{v1}-{\mu }_h{\mu }_v\right).\end{array}$

It is obviously that if ${\mathcal{R}}_0>1$, we have $b_1>0,b_2>0,b_3>0,b_4>0$, and we can verify that $b_3\left(b_1{b}_2-b_3\right)-b_1^2{b}_4>0$. According to the Routh-Hurwitz criterion, we obtain that all eigenvalues of ${\mathcal{J}}^{\mathrm{<mo>\; *\; </mo>}}$ have negative real parts if and only if ${\mathcal{R}}_0>1$. This proof is complete.

4. Global Stability of Equilibria

In this section, we will discuss the global stability of the DFE $P_0$ and the endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}$.

Theorem 5. The disease-free equilibrium $P_0$ of system (2.6) is globally asymptotically stable (GAS) in the interior of the set ${\Omega }_1$ provided that ${\mathcal{R}}_0\le 1$.

Proof. Constructing Lyapunov function

$L\left(E,I,Y\right)=E+\frac{{\mu }_1}{{\mu }_2\theta +\sigma }I+\frac{{\mu }_1{\mu }_2{\mu }_v}{{\beta }_2{\Lambda }_v\left({\mu }_2\theta +\sigma \right)}Y,$ (4.1)

then the derivative of (4.1) with respect to t along the solution curves of (2.6) is given by

$\begin{array}{c}\frac{\text{d}L\left(E\mathrm{,}I\mathrm{,}Y\right)}{\text{d}t}=\frac{\text{d}E}{\text{d}t}+\frac{{\mu }_1}{{\mu }_2\theta +\sigma }\frac{\text{d}I}{\text{d}t}+\frac{{\mu }_1{\mu }_2{\mu }_v}{{\beta }_2{\Lambda }_v\left({\mu }_2\theta +\sigma \right)}\frac{\text{d}Y}{\text{d}t}\\ ={\beta }_{1}SY-{\mu }_{1}E+\frac{{\mu }_1}{{\mu }_2\theta +\sigma }\left(\sigma E-{\mu }_{2}I\right)\\ +\frac{{\mu }_1{\mu }_2{\mu }_v}{{\beta }_2{\Lambda }_v\left({\mu }_2\theta +\sigma \right)}\left({\beta }_2\left(\theta E+I\right)X-{\mu }_{v}Y\right)\\ \le \frac{{\beta }_1{\Lambda }_{h}Y}{{\mu }_h}-{\mu }_{1}E+\frac{{\mu }_1}{{\mu }_2\theta +\sigma }\left(\sigma E-{\mu }_{2}I\right)\\ +\frac{{\mu }_1{\mu }_2{\mu }_v}{{\beta }_2{\Lambda }_v\left({\mu }_2\theta +\sigma \right)}\left(\frac{{\beta }_2\left(\theta E+I\right){\Lambda }_v}{{\mu }_v}-{\mu }_{v}Y\right)\\ =\frac{{\mu }_1{\mu }_2{\mu }_v^2\left({\mathcal{R}}_0^2-1\right)}{{\beta }_2{\Lambda }_v\left({\mu }_2\theta +\sigma \right)}Y\mathrm{.}\end{array}$ (4.2)

Thus, ${\mathcal{R}}_0\le 1$ ensures that $L{\left(E\mathrm{,}I\mathrm{,}Y\right)}^{\prime }\le 0$ for all $E,I,Y>0$. It is easy to verify that the disease-free equilibrim $P_0\left(\frac{{\Lambda }_h}{{\mu }_h}\mathrm{,0,0,}\frac{{\Lambda }_v}{{\mu }_v}\mathrm{,0}\right)$ is the largest invariant set in $\left\{\left(S,E,I,X,Y\right)\in {\Omega }_1:L^{\prime }\left(E,I,Y\right)=0\right\}$, and hence by the LaSalle’s invariance principle [20], we conclude that all trajectories starting in ${\Omega }_1$ approach $P_0$ for ${\mathcal{R}}_0\le 1$. That is to say, $P_0$ is globally asymptotically stable in ${\Omega }_1$ if ${\mathcal{R}}_0\le 1$.

Next, we will investigate the global stability of the endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}$ in the positively invariant set ${\Omega }_1$ by constructing Lyapunov function.

Theorem 6. The endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}$ of system (2.6) is globally asymptotically stable in the interior of the set ${\Omega }_1$ if ${\mathcal{R}}_0>1$.

Proof. From (2.10), we know that $S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{E}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}$ satisfy

$\{\begin{array}{l}{\mu }_h=\frac{{\Lambda }_h-{\beta }_1{S}^{*}{Y}^{*}}{S^{*}},\hfill \\ {\mu }_1=\frac{{\beta }_1{S}^{*}{Y}^{*}}{E^{*}},\hfill \\ {\mu }_2=\frac{\sigma E^{*}}{I^{*}},\hfill \\ {\mu }_v=\frac{{\beta }_2\left(\theta E^{*}+I^{*}\right)X^{*}}{Y^{*}},\hfill \end{array}$ (4.3)

then the system (2.6) can be rewritten as the following form

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}={\Lambda }_h\left(1-\frac{S}{S^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{,}\hfill \\ \frac{\text{d}E\left(t\right)}{\text{d}t}={\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}-\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{,}\hfill \\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\sigma E^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}-\frac{I}{I^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{,}\hfill \\ \frac{\text{d}X}{\text{d}t}={\Lambda }_v\left(1-\frac{X}{X^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{,}\hfill \\ \frac{\text{d}Y}{\text{d}t}={\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{.}\hfill \end{array}$ (4.4)

To illustrate the uniqueness and global stability of $P^{\mathrm{<mo>\; *\; </mo>}}$, we set $V_1\left(t\right)=\frac{1}{{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\left(S-S^{\mathrm{<mo>\; *\; </mo>}}-S^{\mathrm{<mo>\; *\; </mo>}}\ln \frac{S}{S^{\mathrm{<mo>\; *\; </mo>}}}\right)$, $V_2\left(t\right)=\frac{1}{{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\left(E-E^{\mathrm{<mo>\; *\; </mo>}}-E^{\mathrm{<mo>\; *\; </mo>}}\ln \frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\right)$, $V_3\left(t\right)=\frac{1}{\sigma E^{\mathrm{<mo>\; *\; </mo>}}}\left(I-I^{\mathrm{<mo>\; *\; </mo>}}-I^{\mathrm{<mo>\; *\; </mo>}}\ln \frac{I}{I^{\mathrm{<mo>\; *\; </mo>}}}\right)$, $V_4\left(t\right)=X-X^{\mathrm{<mo>\; *\; </mo>}}-X^{\mathrm{<mo>\; *\; </mo>}}\ln \frac{X}{X^{\mathrm{<mo>\; *\; </mo>}}}$ and $V_5\left(t\right)=Y-Y^{\mathrm{<mo>\; *\; </mo>}}-Y^{\mathrm{<mo>\; *\; </mo>}}\ln \frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}$. Then the derivative of functions $V_1\left(t\right)\mathrm{,}{V}_2\left(t\right)\mathrm{,}{V}_3\left(t\right)\mathrm{,}{V}_4\left(t\right)\mathrm{,}{V}_5\left(t\right)$ along solutions of system (2.6) yield

$\begin{array}{c}{{V^{\prime }}_1\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}=\frac{1}{{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\left(1-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}\right)\left({\Lambda }_h\left(1-\frac{S}{S^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\\ \le \left(1-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}\right)\left(1-\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ =1-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}+\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\end{array}$ (4.5)

$\begin{array}{c}{{V^{\prime }}_2\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}=\frac{1}{{\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\left(1-\frac{E^{\mathrm{<mo>\; *\; </mo>}}}{E}\right)\left({\beta }_1{S}^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}-\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\\ =1-\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}-\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}+\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\end{array}$ (4.6)

$\begin{array}{c}{{V^{\prime }}_3\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}=\frac{1}{\sigma E^{\mathrm{<mo>\; *\; </mo>}}}\left(1-\frac{I^{\mathrm{<mo>\; *\; </mo>}}}{I}\right)\left(\sigma E^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}-\frac{I}{I^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\\ =1-\frac{I^{\mathrm{<mo>\; *\; </mo>}}}{I}-\frac{E{I}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}I}+\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\end{array}$ (4.7)

$\begin{array}{c}{{V^{\prime }}_4\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}=\left(1-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}\right)({\Lambda }_v\left(1-\frac{X}{X^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ +{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right))\\ \le \left(1-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}\right)\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\\ ={\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}+\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ +{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}+\frac{I}{I^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{,}\end{array}$ (4.8)

and

$\begin{array}{c}{{V^{\prime }}_5\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}=\left(1-\frac{Y^{\mathrm{<mo>\; *\; </mo>}}}{Y}\right)\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}\right)+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\\ ={\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}-\frac{EX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}+\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ +{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}-\frac{IX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}+\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\mathrm{.}\end{array}$ (4.9)

Define Lyapunov function as follows:

$V\left(t\right)=c_1{V}_1\left(t\right)+c_2{V}_2\left(t\right)+c_3{V}_3\left(t\right)+c_4{V}_4\left(t\right)+V_5\left(t\right),$ (4.10)

where $c_i>\mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}\left(i=1,2,3,4\right)$ are left unspecified. Then the derivative of function $V\left(t\right)$ with respect to t along solutions of system (2.6) is given by

$\begin{array}{c}{V^{\prime }\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}={c_1{V^{\prime }}_1\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}+{c_2{V^{\prime }}_2\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}+{c_3{V^{\prime }}_3\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}+{c_4{V^{\prime }}_4\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}+{{V^{\prime }}_5\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}\\ \le c_1\left(1-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}+\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}\right)+c_2\left(1-\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}-\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}+\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ +c_3\left(1-\frac{I^{\mathrm{<mo>\; *\; </mo>}}}{I}-\frac{E{I}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}I}+\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ +{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}-\frac{EX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}+\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ +c_4\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}+\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\end{array}$

$\begin{array}{c}+c_4\left({\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}+\frac{I}{I^{\mathrm{<mo>\; *\; </mo>}}}\right)\right)\\ +{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}-\frac{IX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}+\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ =c_1+c_2+c_3+\left(c_4+1\right)\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\right)-c_1\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-c_1\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\\ -c_2\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}-c_3\frac{E{I}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}I}-c_4\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\right)\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}\end{array}$

$\begin{array}{c}-{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\frac{EX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}-{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\frac{IX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}\\ +\left(c_1-{\beta }_2\left(\theta E^{\mathrm{<mo>\; *\; </mo>}}+I^{\mathrm{<mo>\; *\; </mo>}}\right)X^{\mathrm{<mo>\; *\; </mo>}}\right)\frac{Y}{Y^{\mathrm{<mo>\; *\; </mo>}}}+\left(c_2-c_1\right)\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\\ +\left(c_3-c_2+c_4{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\right)\frac{E}{E^{\mathrm{<mo>\; *\; </mo>}}}+\left(c_4{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}-c_3\right)\frac{I}{I^{\mathrm{<mo>\; *\; </mo>}}}\\ +{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-c_4\right)\frac{EX}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(1-c_4\right)\frac{IX}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{.}\end{array}$ (4.11)

To choose the suitable constants $c_i>\mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}\left(i=1,2,3,4\right)$ such that ${V^{\prime }\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}$ is negative definite or semidefinite, we eliminate the terms $Y\mathrm{,}SY\mathrm{,}E\mathrm{,}I\mathrm{,}EX\mathrm{,}IX$ by taking

$c_1={\beta }_2\left(\theta E^{*}+I^{*}\right)X^{*},c_2={\beta }_2\left(\theta E^{*}+I^{*}\right)X^{*},c_3={\beta }_2{X}^{*}{I}^{*},c_4=1.$ (4.12)

Consequently, the function (4.10) is specified. Then (4.11) is in turn given as

$\begin{array}{c}{V^{\prime }\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}\le c_1+c_2+c_3+\left(c_4+1\right)\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\right)-c_1\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-c_1\frac{SY}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}}\\ -c_2\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}-c_3\frac{E{I}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}I}-c_4\left({\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}+{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\right)\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}\\ -{\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\frac{EX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}-{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\frac{IX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}\\ ={\beta }_2\theta E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(4-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{EX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}\right)\\ +{\beta }_2{I}^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}\left(5-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{IX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}-\frac{E{I}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}I}\right)\end{array}$ (4.13)

By the relationship between the arithmetic and the associated geometric means, we have

$\begin{array}{l}4-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{EX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}\le \mathrm{0,}\\ 5-\frac{S^{\mathrm{<mo>\; *\; </mo>}}}{S}-\frac{SY{E}^{\mathrm{<mo>\; *\; </mo>}}}{S^{\mathrm{<mo>\; *\; </mo>}}{Y}^{\mathrm{<mo>\; *\; </mo>}}E}-\frac{X^{\mathrm{<mo>\; *\; </mo>}}}{X}-\frac{IX{Y}^{\mathrm{<mo>\; *\; </mo>}}}{I^{\mathrm{<mo>\; *\; </mo>}}{X}^{\mathrm{<mo>\; *\; </mo>}}Y}-\frac{E{I}^{\mathrm{<mo>\; *\; </mo>}}}{E^{\mathrm{<mo>\; *\; </mo>}}I}\le 0.\end{array}$ (4.14)

That is, ${V^{\prime }\left(t\right)|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.6\; <mo\; stretchy="false">\; )\; </mo>}}\le 0$ and the equality holds if and only if $S=S^{*}$, $X=X^{*}$, $\frac{E}{E^{*}}=\frac{Y}{Y^{*}}=\frac{I}{I^{*}}$. It can be easily verified that the largest invariant set of system (2.6) on the set

$\left\{\left(S,E,I,X,Y\right):S=S^{*},X=X^{*},\frac{E}{E^{*}}=\frac{Y}{Y^{*}}=\frac{I}{I^{*}}\right\}$

is the singleton $\left\{P^{\mathrm{<mo>\; *\; </mo>}}\right\}$. Therefore, by the LaSalle's invariance principle [20], it follows that the endemic equilibrium $P^{\mathrm{<mo>\; *\; </mo>}}\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{E}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is globally asymptotically stable in ${\Omega }_1$ when ${\mathcal{R}}_0>1$. This completes the proof.

Remark 2: The basic reproduction number ${\mathcal{R}}_0$ gives a sharp threshold that completely determines their global dynamics [21].

5. Numerical Simulation

In this section, we first provide results from numerical simulations of model (2.1) that illustrate and support our theoretical results. In the model (2.1), all parameters are in months and their values are shown in Table 1. According to the above parameters of the model (2.1), we conduct a global sensitivity analysis on the basic reproductive number ${\mathcal{R}}_0$ by employing Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCCs) [22] [23]. Figure 2 and Figure 3 depict our sensitivity and uncertainty analysis, which involved computing the PRCCs of ${\mathcal{R}}_0$ using the LHS method [22]. From Figure 2, we

Figure 2. PRCC performed on the system (1). Parameters are listed in Table 1. * indicates that PRCCs are very significantly different from zero ( $p<0.05$ ), and ** denotes that PRCCs are very significantly different from zero ( $p<0.01$ ).

Figure 3. The distribution of ${\mathcal{R}}_0$ under LHS parameters.

can easily observe that ${\mathcal{R}}_0$ is sensitive to the natural death rate of citrus population ( ${\mu }_h$ ), infection intensity relative to infected ( $\theta$ ), psyllid recruitment rate ( ${\Lambda }_v$ ), citrus recruitment rate ( ${\Lambda }_h$ ), virus transmission probability from psyllid to plant ( ${\beta }_1$ ), virus transmission probability from plant to psyllid ( ${\beta }_2$ ), natural death rate of psyllid population ( ${\mu }_v$ ), where ${\mathcal{R}}_0$ increases with the increase of $\theta$, ${\Lambda }_v$, ${\Lambda }_h$, ${\beta }_1$ and ${\beta }_2$, but decreases with the increase of ${\mu }_h$ and ${\mu }_v$. This shows that reducing $\theta$, ${\Lambda }_v$, ${\Lambda }_h$, ${\beta }_1$, ${\beta }_2$, or increasing ${\mu }_h$, ${\mu }_v$ will help control the spread of the disease. For uncertainty analysis, it can be seen from Figure 3 that the probability of ${\mathcal{R}}_0>1$ is significantly higher than that of ${\mathcal{R}}_0<1$ (about 2:1), which indicates that disease control should be strengthened for preventing the occurrence of endemic diseases in the current situation.

To better understand the spread of the disease, we further do numerical simulations. Select parameter values as: ${\Lambda }_h=6.6$, ${\mu }_h=0.0033$, ${\beta }_1=0.003125$, $\sigma =0.25$, $d=0.025$, ${\gamma }_1=0.1$, ${\gamma }_2=0.05$, ${\Lambda }_v=2000$, ${\mu }_v=5$, ${\beta }_2=0.0025$, $\theta =0.1$, and $q_1=0$, $q_2=0$, ${\gamma }_3=0$ (i.e., do not take quarantine measures). We get ${\mathcal{R}}_0=3.4133$, $I^{*}=54.3824$, the disease is persistent (Figure 4). Take certain isolation measures $q_1=0.3$, $q_2=0.3$, ${\gamma }_3=0.2$, we have ${\mathcal{R}}_0=1.2066$, $I^{*}=2.0847$, there are still endemic disease (Figure 5). One interesting thing is that for ${\mathcal{R}}_0>1$, the number of the infected populations has a distinct wavy pattern, with multiple peaks and troughs, which means that at some time the infected individual goes to zero and that doesn't mean that the epidemic is under control, it's going to erupt in the near future. Further, if $q_1=0.45$, $q_2=0.45$, ${\gamma }_3=0.2$, we have ${\mathcal{R}}_0=0.9444$, the endemic disease disappears and the disease is well controlled (Figure 6). This shows that the quarantine measure is effective in controlling the spread of HLB.

In addition, the rate of the exposed to the quarantined ( $q_1$ ) and rate of the infectious convert to the quarantined ( $q_2$ ) also have a certain influence on ${\mathcal{R}}_0$. Therefore, in order to further analyze the influence of $q_1$ and $q_2$ on ${\mathcal{R}}_0$, the sensitivity of parameters $q_1\mathrm{,}{q}_2$ for ${\mathcal{R}}_0$ is shown in Figure 7. It can be seen that controlling $q_2$ alone cannot make ${\mathcal{R}}_0$ drop below 1, and controlling $q_1$

Figure 4. Time series plot of infectious citrus tree without quarantine ( $q_1=0$, $q_2=0$, ${\gamma }_3=0$ ).

Figure 5. Time series plot of infectious citrus tree with quarantine ( $q_1=0.3$, $q_2=0.3$, ${\gamma }_3=0.2$ ).

Figure 6. Time series plot of infectious citrus tree with quarantine ( $q_1=0.45$, $q_2=0.45$, ${\gamma }_3=0.2$ ).

alone can only make ${\mathcal{R}}_0$ close to 1, but the joint control of $q_1$, $q_2$ can make ${\mathcal{R}}_0<1$.

According to PRCC analysis (see Figure 2), virus transmission probability from psyllid to plant ( ${\beta }_1$ ), virus transmission probability from plant to psyllid ( ${\beta }_2$ ), natural death rate of citrus population ( ${\mu }_h$ ) and natural death rate of psyllid population ( ${\mu }_v$ ) have great influence on ${\mathcal{R}}_0$. We conduct sensitivity analysis of these four parameters to ${\mathcal{R}}_0$, as shown in Figure 8 and Figure 9. It implies that disease control can be achieved by inhibiting ${\beta }_1$ and ${\beta }_2$ or increasing ${\mu }_h$ and ${\mu }_v$ simultaneously.

Figure 7. The sensitivity analysis of for ${\mathcal{R}}_0$ with the parameters $q_1\mathrm{,}{q}_2$.

Figure 8. The sensitivity analysis of for ${\mathcal{R}}_0$ with the parameters ${\beta }_1\mathrm{,}{\beta }_2$.

Figure 9. The sensitivity analysis of for ${\mathcal{R}}_0$ with the parameters ${\mu }_h\mathrm{,}{\mu }_v$.

6. Conclusions

In this work, a citrus-psyllid dynamic model with quarantine measure is formulated. Based on the method of next-generation matrix, we obtain the expression of basic reproductive number ${\mathcal{R}}_0$. The global stability of disease-free equilibrium and endemic equilibrium are demonstrated by constructing Lyapunov functions.

Numerical simulation shows that when ${\mathcal{R}}_0>1$, the number of infected hosts has obvious wavy, with multiple peaks and troughs, which implies that at a time infected individuals tending to zero does not mean that the disease is under control, instead, it may break out soon. Additionally, our investigations show that a certain degree quarantine measure is effective, and it is more effective to isolate the exposed than the infected host. Moreover, comprehensive quarantine can more effectively control the outbreak of disease. These results can provide a reference for fruit industry to conduct comprehensive management of HLB.

Acknowledgements

This research was supported by the Education Department of Jiangxi Province (GJJ190740, GJJ201406).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Kumar, P., Suat Ertürk, V. and Nisar, K.S. (2021) Fractional Dynamics of Huanglongbing Transmission within a Citrus Tree. Mathematical Methods in the Applied Sciences, 44, 11404-11424. https://doi.org/10.1002/mma.7499
[2]	Lee, J.A., Halbert, S.E., Dawson, W.O., Robertson, C.J., Keesling, J.E. and Singer, B.H. (2015) Asymptomatic Spread of Huanglongbing and Implications for Disease Control. Proceedings of the National Academy of Sciences, 112, 7605-7610. https://doi.org/10.1073/pnas.1508253112
[3]	Jackson, M. and Chen-Charpentier, B.M. (2017) Modeling Plant Virus Propagation with Delays. Journal of Computational and Applied Mathematics, 309, 611-621. https://doi.org/10.1016/j.cam.2016.04.024
[4]	Shi, R., Zhao, H. and Tang, S. (2014) Global Dynamic Analysis of a Vector-Borne Plant Disease Model. Advances in Difference Equations, 2014, Article No. 59. https://doi.org/10.1186/1687-1847-2014-59
[5]	Loebenstein, G., Berger, P.H. and Brunt, A.A. (2001) Virus and Viruslike Diseases of Potatoes and Production of Seed-Potatoes. Springer Science & Business Media, Berlin. https://doi.org/10.1007/978-94-007-0842-6
[6]	Whitfield, A.E., Falk, B.W. and Rotenberg, D. (2015) Insect Vector Mediated Transmission of Plant Viruses. Virology, 479, 278-289. https://doi.org/10.1016/j.virol.2015.03.026
[7]	Jacobsen, K., Stupiansky, J. and Pilyugin, S.S. (2013) Mathematical Modeling of Citrus Groves Infected by Huanglongbing. Mathematical Biosciences & Engineering, 10, 705-728. https://doi.org/10.3934/mbe.2013.10.705
[8]	Chiyaka, C., Singer, B.H., Halbert, S.E., Morris, J.G. and van Bruggen, A.H. (2012) Modeling Huanglongbing Transmission within a Citrus Tree. Proceedings of the National Academy of Sciences, 109, 12213-12218. https://doi.org/10.1073/pnas.1208326109
[9]	Sanchez, L., Pant, S., Xing, Z., Mandadi, K. and Kurouski, D. (2019) Rapid and Noninvasive Diagnostics of Huanglongbing and Nutrient Deficits on Citrus Trees with a Handheld Raman Spectrometer. Analytical and Bioanalytical Chemistry, 411, 3125-3133. https://doi.org/10.1007/s00216-019-01776-4
[10]	Taylor, R.A., Mordecai, E.A., Gilligan, C.A., Rohr, J.R. and Johnson, L.R. (2016) Mathematical Models Are a Powerful Method to Understand and Control the Spread of Huanglongbing. PeerJ, 4, e2642. https://doi.org/10.7717/peerj.2642
[11]	Of Agriculture, M. and of the People’s Republic of China, R.A. (2020) Notice of Public Comment on the First Class of Crop Diseases and Insect Pests. http://www.moa.gov.cn/xw/bmdt/202006/t20200604_6345940.htm
[12]	Bassanezi, R., Belasque Jr., J. and Montesino, L. (2013) Frequency of Symptomatic Trees Removal in Small Citrus Blocks on Citrus Huanglongbing Epidemics. Crop Protection, 52, 72-77. https://doi.org/10.1016/j.cropro.2013.05.012
[13]	Wang, J., Feng, F., Guo, Z., Lv, H. and Wang, J. (2019) Threshold Dynamics of a Vector-Borne Epidemic Model for Huanglongbing with Impulsive Control. Applied Mathematics, 10, 196-211. https://doi.org/10.4236/am.2019.104015
[14]	Parry, M., Gibson, G.J., Parnell, S., Gottwald, T.R., Irey, M.S., Gast, T.C. and Gilligan, C.A. (2014) Bayesian Inference for an Emerging Arboreal Epidemic in the Presence of Control. Proceedings of the National Academy of Sciences, 111, 6258-6262. https://doi.org/10.1073/pnas.1310997111
[15]	Bové, J.M. (2006) Huanglongbing: A Destructive, Newly-Emerging, Century-Old Disease of Citrus. Journal of Plant Pathology, 88, 7-37.
[16]	Luo, L., Gao, S., Ge, Y. and Luo, Y. (2017) Transmission Dynamics of a Huanglongbing Model with Cross Protection. Advances in Difference Equations, 2017, Article No. 355. https://doi.org/10.1186/s13662-017-1392-y
[17]	Pelz-Stelinski, K., Brlansky, R., Ebert, T. and Rogers, M. (2010) Transmission Parameters for Candidatus Liberibacter Asiaticus by Asian Citrus Psyllid (Hemiptera: Psyllidae). Journal of Economic Entomology, 103, 1531-1541. https://doi.org/10.1603/EC10123
[18]	Gao, S., Yu, D., Meng, X. and Zhang, F. (2018) Global Dynamics of a Stage-Structured Huanglongbing Model with Time Delay. Chaos, Solitons & Fractals, 117, 60-67. https://doi.org/10.1016/j.chaos.2018.10.008
[19]	Van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[20]	La Salle, J.P. (1976) The Stability of Dynamical Systems. SIAM, Philadelphia. https://doi.org/10.1137/1.9781611970432
[21]	Shuai, Z. and van den Driessche, P. (2013) Global Stability of Infectious Disease Models Using Lyapunov Functions. SIAM Journal on Applied Mathematics, 73, 1513-1532. https://doi.org/10.1137/120876642
[22]	Marino, S., Hogue, I.B., Ray, C.J. and Kirschner, D.E. (2008) A Methodology for Performing Global Uncertainty and Sensitivity Analysis in Systems Biology. Journal of Theoretical Biology, 254, 178-196. https://doi.org/10.1016/j.jtbi.2008.04.011
[23]	Zhang, F., Qiu, Z., Huang, A. and Zhao, X. (2021) Optimal Control and Cost-Effectiveness Analysis of a Huanglongbing Model with Comprehensive Interventions. Applied Mathematical Modelling, 90, 719-741. https://doi.org/10.1016/j.apm.2020.09.033