Huanglongbing Model under the Control Strategy of Discontinuous Removal of Infected Trees

: By using differential equations with discontinuous right-hand sides, a dynamic model for vector-borne infectious disease under the discontinuous removal of infected trees was established after understanding the transmission mechanism of Huanglongbing (HLB) disease in citrus trees. Through calculation, the basic reproductive number of the model can be attained and the properties of the model are discussed. On this basis, the existence and global stability of the calculated equilibria are veriﬁed. Moreover, it was found that different I 0 in the control strategy cannot change the dynamic properties of HLB disease. However, the lower the value of I 0 , the fewer HLB-infected citrus trees, which provides a theoretical basis for controlling HLB disease and reducing expenditure.


Introduction
Vector-borne disease is an infectious disease spread by media. These diseases do not only occur in the population, but also affect tree populations. Huanglongbing disease is a kind of vector-borne disease in plant populations. As an infectious disease affecting plants, citrus Huanglongbing (HLB) disease results in inestimable profit loss in terms of the production of citrus fruit. Therefore, the disease is regarded as a cancer in the citrus orchard. The transmission of citrus HLB disease has a history of nearly 100 years in China. However, the transmission of the disease was not reported until the early 20th century, when first identified in south China (in Guangdong and Fujian provinces as well as the Guangxi Zhuang autonomous region). During the 1970s and 1980s, citrus HLB disease spread to multiple provinces (including Jiangxi, Yunnan, and Sichuan provinces) in China, affecting the production of citrus fruit and damaging the agricultural economy in local areas. Moreover, the disease tended to spread to surrounding provinces [1]. As expected, among 19 provinces of China where citrus trees were planted, more than two thirds suffered HLB disease after several years. To date, the HLB disease mainly occurs in more than 50 countries and regions in five continents (except for Europe and Antarctica) [2]. Through investigation, Psyllidae is the main transmitting vector of HLB disease, which leads to an extremely rapid rate of transmission. Moreover, with fruit pulp of poor quality, the yield from HLB-infected citrus trees decreases and the longevity of such trees declines. As a result, the citrus industry fails to develop, which leads to a huge economic loss for local agricultural production. To guarantee the normal growth of citrus trees and the increase in profits from agricultural production, it is essential to explore the development of HLB disease and effectively prevent and control it. To date, citrus species have no real resistance to HLB disease, and there is no effective treatment for this disease. Mathematical models play an important role in analyzing the epidemiology of vector-borne plant pathogens. Van der Plank [3] and Kranz [4] first used the mathematical dynamics model to study the dynamics of plant diseases. In 1998, J. Jeger et al. established in the literature [5] the kinetic model of vector transmission of the African cassava mosaic disease with constant coefficients. The HLB model can be established with similar methods. There are currently several HLB mathematical models to analyze how HLB spreads among individual trees, citrus groves, or forests [6]. Chiyaka et al. [7] first proposed a mathematical model to describe and analyze the transmission of HLB in trees, and studied the effect of flushing on the dynamics of psyllids. Vilamiu et al. [8] proposed a model of HLB transmission between citrus plants, including delay time and continuous human intervention.
In recent decades, an integrated management strategy for diseases has become widely used in treating plant diseases: this includes replanting, removing inferior plants, introducing natural enemies, genetic modification, and the spraying of pesticides. At present, numerous scholars have applied an integrated management strategy for diseases to an epidemic dynamics model. In addition, obtaining the solution of a differential equation is an important means to investigate its dynamic behavior. The common methods of approximate analytical solution are: homotopy perturbation method [9][10][11][12][13]; homotopy analysis method [14,15]; variational iteration method [16][17][18]; Taylor series method [19]; etc. [20][21][22] for discriminating the development of the disease was proposed and it showed that the infection-free periodic solution of the system is stable on condition that R 1 < 1; however, for R 1 > 1, the disease will be permanent. Nevertheless, the strategy of implementing planting and fixed-time removal is baseless. For the management of disease, it is necessary to utilize multi-stage control measures according to the degree of infection. The multi-stage control strategy is mathematically discontinuous. Discontinuous power systems existed in large numbers in real life. Due to the limitation of analytical tools, the development of the theory of differential equations with a discontinuous right-hand side is relatively lagging and rather slow. Scholars of the former Soviet Union such as F. Filippov performed a lot of excellent work on this, in the meaning of the solution defined by it, by systematically analyzing the basic properties of the understanding and the stability of the system. At present, the theory of differential equations with discontinuous right-hand side is mainly used in artificial neural network dynamics, fishery dynamics, and infectious disease dynamics. For example, M. Forti and P. Nisti investigated the global stability of neural network models with discontinuous activation functions for the first time [27]. Guo Zhenyuan performed research on autonomous fishery models with discontinuous catches and infectious disease models with discontinuous treatment strategies [28,29]. Therefore, differential equations with discontinuous right-hand sides were employed in the model for HLB disease spread to establish a dynamic model for HLB disease under the control strategy of the discontinuous removal of infected trees.
This paper is broadly organized as follows. First, we establish the model, and then use the right discontinuous differential equation theory to analyze the dynamics of the modelmainly the existence and stability analysis of the disease-free equilibrium point and the endemic equilibrium point. Finally, the model is substituted into reasonable parameters, and the correctness of the above conclusions was verified through numerical simulation.

Assumptions and Establishment of the Model
In this paper, a mathematical model of vector infectious disease under the control strategy of discontinuous disease tree removal was established on the basis of existing infectious disease models [30,31]. Citrus trees are divided into susceptible and infected trees, which are recorded as S h (t), I h (t), respectively. The vector Psyllidae is partitioned into susceptible and infected Psyllidae, which are separately marked as S v (t) and I v (t), respectively. Moreover, the following assumptions are made:

1.
The total number of citrus trees remains unchanged, i.e., S h (t) + I h (t) = K. If plants are dead, the same number of susceptible citrus trees are re-planted; 2.
Only Psyllidae which sucks the tender bud of infected citrus trees is regarded as infectious Psyllidae. The susceptible citrus trees which are sucked by infected Psyllidae are considered infected; 3.
All newly planted citrus trees and new Psyllidae are not infected with HLB disease; 4.
HLB virus cannot be detected in susceptible citrus trees immediately after infection by Psyllidae but requires a certain time τ after sucking. The citrus trees are only infectious in this case; 5.
The reproduction of Psyllidae is controlled by spraying pesticide to control the transmission of citrus HLB disease.
According to the aforementioned assumptions, the following dynamic models are established: where Λ represents the constant input rate of susceptible Psyllidae; β 1 and β 2 denote the infection rates of susceptible citrus trees and Psyllidae; µ 1 and µ 1 represent the natural mortality rate and mortality rate due to disease in citrus trees, respectively; d, p, and CI h denote the natural mortality rate of Psyllidae, the delousing rate arising from the use of pesticide, and the rate of removal of infected trees, respectively; CI h is a piecewise function of I h , such that: where c 1 < c 2 . According to practical biological meaning, it is supposed that all of the aforementioned parameters are positive.
It is noted that In this context, it can be verified that Ω is a positive invariant set in model (1). Given that N v (t) = S v (t) + I v (t), by adding the third and the fourth equations in system (1), we have: The model is mainly to explore the changing trend of diseased citrus and diseased psyllids.
the changes of diseased citrus and diseased psyllids in the system (1) will become more and more similar to the following model over time: Based on the published literature [32], the solution to model (2) under Filippov conditions can be defined as thus: the vector function (I h (t), I v (t)) satisfies the following conditions: For almost all t ∈ [0, +∞), (I h (t), I v (t)) shows the following differential inclusions: where Thus, (I h (t), I v (t)) represents the solution to model (2) meeting the initial conditions I h (0) = I h0 and I v (0) = I v0 under Filippov conditions.
It can be shown that (I h , is an upper semi-continuous set-valued mapping with a non-void compact convex set. Before conducting the dynamic analysis of the model, it is shown that model (2) satisfies a positive initial condition I h (0) = I h0 ≥ 0. The solution to I v (0) = I v0 ≥ 0 shows positivity and boundedness, satisfying biological significance criteria.

Existence of Equilibria
Due to the piecewise linearity of control function C(I h ), the phase plane of (I h (t), I v (t)) is divided into three parts: Therefore, model (2) is expressed as follows within the range of G − : Within the range of G + , it is as follows: Within the range of Σ, it is written as follows: Let the right-hand side of the equation in system I be equal to zero: through calculation, the solution I h (t) ≡ 0 and I v (t) ≡ 0 is obtained, i.e., the disease-free equilibrium E 0 = (0, 0, ) of system I. Afterwards, system I is subjected to linearization at the diseasefree equilibrium E 0 = (0, 0, ). In this case: According to other work [33], the basic reproductive number R 1 of the system I can be calculated, with the following specific steps: It is defined that I h0 , I v0 separately refer to the initial values of two infected classes I h (t), I v (t) at = 0. On this condition, I h (t) = I h0 e −(µ 1 +µ 2 +c 1 )t and I v (t) = I v0 e −(µ 1 +µ 2 +c 1 )t . Therefore, the total numbers among the population infected ab initio are separately: Thus: Based on the spectral radius of the matrix M 0 with basic reproductive number, i.e., the basic reproductive number R 1 of the model: Above all, in system I, the basic reproductive number is R 1 = β 1 β 2 N * K (d+p)(µ 1 +µ 2 +c 1 ) and the disease-free equilibrium is E 0 = (0, 0); when R 1 > 1, there is an endemic equilibrium in system I, and E − * can be expressed as Similarly, for system I I, it can be found that: the basic reproductive number is R 2 = β 1 β 2 N * K (d+p)(µ 1 +µ 2 +c 2 ) and the disease-free equilibrium is E 0 = (0, 0). When R 2 > 1, an endemic equilibrium develops: By analyzing systems I and I I, the endemic equilibria of model (2) are calculated according to the following theorems. Proof.
When c 1 < c 2 , it can be seen that I − h * > I + h * : 1.
, the endemic equilibrium E − * of system I always exists.
In the case that R 1 = Then, E − * of system I is the true equilibrium state of model (2) and E + * of system I I is the pseudo-equilibrium state of model (2). Therefore, the endemic equilibrium E − * = E + * of model (2) exists;

2.
When R 2 > K(d+p)+β 2 KI 0 K(d+p)−(d+p)I 0 , endemic equilibrium E + * always exists in system I I. If Then, E + * of system I I is the true equilibrium state of model (2) and E − * of system I is the pseudo-equilibrium state of model (2). Therefore, there is an endemic equilibrium E * = E + * in model (2); In this case, E − * ∈ G + , E + * ∈ G − are both pseudo-equilibrium states of model (2), i.e., they are not endemic equilibria of model (2).
In system I I I, according to the definitions of equilibria of differential inclusion, the endemic equilibrium of system III satisfies: where co[C( According to the second equation in (4), we have I v = β 2 N * I h β 2 I h +d+p . By substituting I v = β 2 N * I h β 2 I h +d+p into the first equation, Equation (5) is as follows: Given Therefore, I 0 satisfies (5) so E * = I 0 , β 2 N * I 0 β 2 I 0 +d+p is the endemic equilibrium of system I I I. Moreover, due to E * ∈ ∑, E * is the true equilibrium state of model (2), i.e., the endemic equilibrium of model (2).

Stability of Equilibria
The stability of the equilibrium of model (2) is analyzed by utilizing the generalized Lyapunov theory in the system with discontinuous right-hand sides, as follows: Theorem 2. In model (2), the disease-free equilibrium E 0 = (0, 0) shows global asymptotic stability when R 1 ≤ 1.

Proof.
For model (2), according to the measurable selection theorem [9], it can be seen that there is a measurable function c(t) ∈ co[C(I h )], which can make most of t satisfy the following equations: The Lyapunov function is constructed on the positive invariant set Ω: Through direct verification, it can be determined that V(t) is a positively definite local Lipschitz regular function on t. Moreover, for random l > 0, the set L 0l = {t ∈ R|V 1 (t) ≤ l} is bounded.
According to the generalized chain rule, it can be seen that nearly all t satisfies: When R 1 < 1, there is dV dt ≤ 0, in which dV dt = 0 exists only at I h (t) = 0. When R 0 = 1, there is dV dt = 0, so it can be inferred that I v (t) = 0. Therefore, based on the generalized LaSalle invariance principle [10], the disease-free equilibrium E 0 = (0, 0) of model (2) is globally asymptotically stable.
Proof. The existence and uniqueness of the endemic equilibrium E * of model (2) were derived above (refer to Theorem 1). The stability of the endemic equilibrium E * is verified as follows.
For model (2), according to the measurable selection theorem [34], it can be found that there is a measurable function c(t) ∈ co[C(I h )] such that nearly all t satisfies: Additionally, c(t) is differentiable for nearly all t and the derivative is equal to zero. According to: it can be found that: where the measurable function c(t) ∈ co[C(I h * )].
The following Lyapunov functional is constructed: Through direct verification, it can be found that V(t) is a positively definite local Lipschitz regular function of t. Moreover, for random l > 0, the set L 0l = {t ∈ R + |V 1 (t) ≤ l} is bounded.
Based on the generalized chain rule, it can be seen that for nearly all t, there exists: Therefore, according to the generalized LaSalle invariance principle [35], no matter which range the endemic equilibrium E * of the model (2) is within and no matter what the value of I 0 is, E * shows globally asymptotic stability.

Numerical Simulation
For the condition with only disease-free equilibrium, the parameter selection of model (2) is expressed as follows: In this case, the basic reproductive numbers are R 1 = 0.97 and R 2 = 0.75. According to Theorems 1 and 2, it can be seen that there is no endemic equilibrium and the only a disease-free equilibrium E 0 = (0, 0) exhibits globally asymptotic stability.
It can be seen from Figure 1 that the solutions (I h (t), I v (t)) under different initial conditions finally tend to a same value E 0 = (0, 0). Therefore, the disease-free equilibrium is globally asymptotically stable.
In the case that an endemic equilibrium exists, the parameter selection for model (2) is as follows: When the basic reproductive numbers are R 1 = 1.55, R 2 = 1.31, E − * = (0.64, 0.78) and E + * = (0.38, 0.55). According to Theorems 1 and 3, the following conclusions can be drawn:
As shown in Figures 2-4, when I 0 = 1 > 0.64, the solutions all tend to the endemic equilibrium E − * = (0.64, 0.78); in the case of I 0 = 0.3 < 0.38, the solutions all stabilize at the endemic equilibrium E + * = (0.38, 0.55); when I 0 = 0.5 ∈ [0.38, 0.64], the solutions all approximate to the endemic equilibrium E * = (0.5, 0.67). Above all, when I 0 varies, the model will show different endemic equilibria; however, regardless of the value of I 0 , the solutions (I h (t), I v (t)) under different initial values eventually approximate to the endemic equilibrium. Therefore, the endemic equilibrium exhibits globally asymptotic stability.

Conclusions
According to the characteristics of HLB disease and the characteristics of control measures for removing infected trees, reasonable assumptions are proposed. Moreover, by selecting proper parameters and reasonable functions, the influence of the control strategy of discontinuously removing infected trees on HLB disease is described. Furthermore, a corresponding mathematical model was established to explore the dynamic properties of the system. Through numerical simulation, the existence and global stability of the resulting equilibria were validated. Additionally, it was found that different values of I 0 do not change the dynamic properties of HLB disease. However, the lower the value of I 0 is, the fewer HLB-infected citrus trees there are. When R 1 is less than or equal to 1, the model only has a disease-free equilibrium point that is globally asymptotically stable. This means that fruit growers can artificially increase the removal rate of diseased trees so that R 1 is less than or equal to 1. Huanglongbing disease will gradually disappear at this time. Of course, the cost of control will be more. When R 1 is greater than 1, the model has a globally asymptotically stable endemic equilibrium point, which means that Huanglongbing disease will continue to develop. The value of I 0 will not change the dynamic properties of Huanglongbing disease; however, the smaller the value of I 0 is, the less citrus will be infected. A small value of I 0 means that the removal of diseased trees is increased when there are few infected citrus. At this time, the cost of removing diseased trees will also increase, and after the final stabilization, the less diseased citrus will be, and the corresponding citrus sales revenue will be increased. When R 1 is less than or equal to 1, the model only has a disease-free equilibrium point that is globally asymptotically stable. This means that fruit growers can artificially increase the removal rate of diseased trees so that R 1 is less than or equal to 1. Huanglongbing disease will gradually disappear at this time. Of course, the cost of control will be increased. When R 1 is greater than 1, the model has a globally asymptotically stable endemic equilibrium point, which means that Huanglongbing disease will continue to develop. The value of I 0 will not change the dynamic properties of Huanglongbing disease, but the smaller the value of I 0 , the less citrus will be infected. A small value of I 0 means that the removal of diseased trees is increased when there are few infected citrus. At this time, the cost of removing diseased trees will also increase, and after the final stabilization, the less diseased citrus there will be, and the corresponding citrus sales revenue will be increased. This research provides a theoretical basis for fruit growers trying to control HLB disease and reducing expenditure thereon.
Author Contributions: Conceptualization, P.W. and W.L.; methodology, P.W.; writing-original draft preparation, X.L. and L.X.; writing-review and editing, P.W. and W.L.; software, P.W.; funding acquisition, W.L. All authors have read and agreed to the published version of the manuscript.