Transmission dynamics and optimal control of a Huanglongbing model with time delay

Abstract: In this paper, a mathematical model has been formulated for the transmission dynamics of citrus Huanglongbing considering latent period as the time delay factor. Existence of the equilibria and their stability have been studied on the basis of basic reproduction number in two cases τ = 0 and τ > 0. The results show that stability changes occur through Hopf bifurcation in the delayed system. Optimal control theory is then applied to investigate the optimal strategy for curtailing the spread of the disease using three time-dependent control variables determined from sensitivity analysis. By using Pontryagin’s Maximum Principle, we obtain the optimal integrated strategy and prove the uniqueness of optimal control solution. Analytical and numerical findings suggest that it is feasible to implement control techniques while minimizing the cost of implementation of optimal control strategies.


Introduction
Citrus Huanglongbing (HLB), currently the most economically destructive disease of citrus worldwide, is caused by phloem-restricted bacteria of the Candidatus Liberibacter group with symptoms including stunted growth and poor flowering of citrus trees as well as blotchy mottling and yellowing of their leaves [1][2][3]. It was first reported from southern China in 1919, with a history of more than a century in China [1]. Now it is distributed in nearly 50 countries and regions in Asia, Africa, North and South America, Europe, etc., and its occurrence scope is still expanding at present, which makes it true that over 80% of the planted area can not escape from being infected [4][5][6]. The infected trees are usually destroyed or become unproductive in 2 to 5 years [1,7]. To date, there is no good source of genetic resistance to HLB in the genus citrus, and the disease can not be cured once the trees are infected [8]. Therefore, HLB still pose a significant threat to the development of the citrus industry.
Experimental studies have shown that early stages of HLB, the bacteria can remain undetectable in the trees after initial infection. It might take a 2 to 6 months latent period for the tree to become infectious before new psyllid generations can become infected [9]. HLB latency is highly variable and apparently greatly affected by tree age, horticultural health, and other factors. In order to be more precise description of HLB transmission in a citrus orchard, it is reasonable to incorporate time latency into the delayed system. Delay differential equations (DDEs) are one of the most powerful mathematical tools for studying the effect of delay in the dynamics of the disease [10][11][12][13][14][15]. For instance, AI Basir et al. [12] derived a mathematical model for the transmission dynamics of vector-borne plant disease considering incubation period as the time delay factor. Fan et al. [14] formulated a delay differential equation model for the transmission of West Nile virus between vector mosquitoes and avian hosts that incorporates maturation delay for mosquitoes. They found that a combination of the maturation delay and the vertical transmission of the virus in mosquitoes has substantial influence on the abundance and number of infection peaks of the infectious mosquitoes. Vilamiu et al. [15] presented a delayed HLB model and assess the influence of a long incubation period of HLB in the plant. Motivated by the preceding discussion, our first aim of this paper is to formulate a mathematical model that is able to describe the dynamics of citrus HLB in citrus orchards, considering latent period as the time delay factor, and then to analyze the dynamic properties of the model. Disease management is complicated, because incubation periods are long and diagnosis is difficult. For the devastating citrus disease, the current state of preventive strategy involves insecticide spraying to reduce the abundance of psyllid vector population, removal of infected trees to reduce inoculum sources and nutrient solution injection to reduce the infection of the bacteria. There is an urgent need to find a sustainable and economical solution to the HLB menace through psyllid control by insecticides, removal of infected trees and nutrient solution injection. In recent years, using the Pontryagin's Maximum Principle to determine optimal control strategies for various diseases has flourished, see [16][17][18][19][20][21][22] and the references therein. For instance, the work by Lee et al. [18] introduced a dynamic model of pine nematode disease, in which nutrient injection, tree cutting, insecticide spraying and inhibition of mating of wood louse were used as control variables. In 2016, Pei et al. [22] established a delayed SIS model of infectious diseases with stage structure and introduced three treatment strategies into the model to study the existence and uniqueness of optimal control solution. But for the optimal control research of citrus HLB, there is few literatures (see [19,20]). In [19] and [20], although the authors considered different interventions for disease control, they did not take into account the time latency, which may lead to inappropriate policy making for preventing disease. Therefore, the second aim of this paper is to explore an optimal integrated strategy for a delayed HLB model considering the disease latency.
The remaining of this paper is organized as follows. A model formulation is presented and the basic dynamics behavior of the model is analyzed in Section 2. Then sensitivity analysis is performed in Section 3. In Section 4, by using the results of the sensitivity analysis, an optimal control problem is developed, the optimal strategy is obtained and the uniqueness of the optimal control solution is proved. Numerical results with detailed discussions are presented in Section 5. Finally, this paper ends with a brief conclusion and the potential outlook for future work in Section 6.

Model formulation
In this subsection, we will improve a HLB model taking the disease latency into account. Suppose that the citrus trees population totalling N h consists of two categories that include the susceptible S h and the infectious I h . Additionally, the psyllids population, denoted by N v , is divided into two classes, namely, susceptible S v and infectious I v .
In order to better understand our model, the assumptions we made according to the transmission mechanism of HLB should be mentioned first of all: (A1): A susceptible citrus tree is assumed to become infected from contact with an infected psyllid with a probability of successful infection β, a susceptible psyllid becomes infected from contact with an infected tree, which happens at a probability α.
(A2): All newly planted citrus trees are susceptible, moreover, new replanting has taken place at rate r with a proportion, K − N h , of the orchard, where K represents the maximum capacity of citrus trees in the orchard.
(A3): All newborn psyllids are susceptible, Λ is used to denote the constant recruitment rate of psyllid population.
(A4): There is a 2 to 6 months latent period for susceptible citrus tree to have the capacity of contagiousness.
Thus, a mathematical model considering latent period denoted by τ as the time delay factor can be formulated as follows: where µ h and µ v are the natural mortality of citrus population and psyllid population, respectively. γ is the the disease-induced death rate of the infected trees. All parameters of model (2.1) are positive constants. Let C denote the Banach space of continuous functions φ : . Therefore, the initial condition for model (2.1) is taken as below: Note that the citrus tree population N h (t) = S h (t) + I h (t) and psyllid population respectively, which implies that Define Obviously, Ω is a positively invariant set for system (2.1).

The basic reproductive number
As we all know, the basic reproductive number R 0 is often considered as the threshold quantity that determines the dynamic behavior of the model. Following, we will calculate the basic reproductive number R 0 of system (2.1) following the idea of [23].
Suppose I h (t) ≡ 0, I v (t) ≡ 0, we can easily obtain the disease-free equilibrium, given by Next, linearizing system (2.1) at the disease-free equilibrium E 0 (here we only write down the equations of the diseased classes for simplification), this leads to If we set I h0 , I v0 be the numbers of each diseased class at time t = 0, and I h (t), I v (t) be the remaining populations of each class at time t, respectively, then we have Thus, the total numbers of newly infected in each diseased class can be expressed as It then follows We can see that M is the next infection operator, and the basic reproductive number, denoted by R 0 , can be described by the spectral radius of the matrix M. Therefore, R 0 for system (2.1) is given by

Stability of equilibria and Hopf bifurcation
Qualitative analysis, including existence and stability are important topics in epidemic dynamics, we shall provide the stability analysis of equilibria for system (2.1) and look for the possibility of Hopf bifurcation by taking time delay "τ" as a bifurcation parameter in this subsection.
When τ ≥ 0, set the right-hand side of system (2.1) to zero, namely Obviously Λ µ v − S * v > 0 when R 0 > 1, which leads to the existence of the unique endemic equilibrium. We have the following results on the stability of equilibria E 0 and E * by using the method in [24] when τ = 0.
Theorem 2.2. When τ = 0, the endemic equilibrium E * of (2.1) is locally asymptotically stable if R 0 > 1. Now we turn to discussing the stability of the endemic equilibrium E * when R 0 > 1 in case of τ > 0. In addition, we derive the stability conditions for E * as well as the conditions for Hopf bifurcation.
Linearizing system (2.1) at E * , we get the following systeṁ where X = (S h , I h , S v , I v ) T , A and B are 4 × 4 matrices, given by It follows that To show Hopf bifurcation, we must have a pair of purely imaginary roots of characteristic Eq (2.3). Suppose iω (ω > 0) is a purely imaginary root associated with Eq (2.3), thus = − c 0 ω 2 cos ωτ + ic 1 ω cos ωτ + ic 0 ω 2 sin ωτ + c 1 ω sin ωτ.
Separating real and imaginary parts, we have following equations Squaring and adding (2.4) and (2.5) we obtain Let z = ω 2 , then Eq (2.6) yields where the coefficients are given by Now, if the following conditions q 1 ≥ 0, q 2 ≥ 0, q 3 > 0, q 4 ≥ 0 are satisfied, then there will be no positive roots of Eq (2.7), which implies there is no any purely imaginary roots of characteristic Eq (2.3). In the case for τ ≥ 0 all characteristic roots of Eq (2.3) have negative real part by Rouchet theorem [25]. Hence, the following result is claimed. Theorem 2.3. If the conditions R 0 > 1, q 1 ≥ 0, q 2 ≥ 0, q 3 > 0 and q 4 ≥ 0 are satisfied, then the endemic equilibrium E * of system (2.1) is local asymptotically stable for all τ ≥ 0.
Let ω 0 be a positive real root of Eq (2.6). Then there is a purely imaginary root iω 0 of Eq (2.3). From Eqs (2.4) and (2.5) we obtain the values of cos ωτ as Furthermore, substituting ω = ω 0 into (2.8), then delay factor τ can be presented Hence (ω 0 , τ k ) is the solution of Eq (2.3), that is to say, there exists a pair of purely imaginary roots λ = ±iω 0 for Eq (2.3) when τ = τ k .
Note that τ = τ 0 is the minimum value of τ when the purely imaginary root λ = ±iω 0 of Eq (2.3) occurs. Based on above discussion, we state following lemmas which will be essential to our claims.

Sensitivity analysis
Uncertainty analyses and sensitivity analyses are necessary to explore the behavior of many complex models, because the structural complexity of the models are coupled with a high degree of uncertainty in estimating the values of many of the input parameter [26]. In this section, we employed a global sensitivity analysis to assess the impact of uncertainty and the sensitivity of the outcomes of the numerical simulations to variations in each parameter of model (2.1) using Latin Hypercube Sampling (LHS) and partial rank correlation coefficient (PRCC). LHS is a stratified sampling without replacement technique which allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [26][27][28][29]. Further, PRCC is used to measure the strength of the relationship between the model outcome and the parameters, state the degree of the effect that each parameter has on the outcome [26][27][28][29]. Therefore, by sensitivity analyses we can assess the parameters with the most significant impact on the outcome of the numerical simulations of the model. Initial disease transmission is directly related to the reproductive number R 0 , it is a combination of parameters related to both the trees and the psyllids. Therefore, in determining the control measure which is the most effective in controlling HLB transmission, we use R 0 as the response function. Firstly, we calculate the sensitivity indices of R 0 to the parameters in the model. These indices tell us how crucial each parameter is to disease transmission and prevalence. All numerical parameter values of model (2.1) are presented in Table 1.
Definition 3.1. [30] The normalized forward sensitivity index of a variable, u, that depends differentially on a parameter, p, is defined as: Then we derive an analytical expression for the sensitivity of R 0 , to each of the nine different parameters described in Table 2. For example, the sensitivity index of R 0 with respect to β, does not depend on any parameter values. It follows from Table 2 that the most sensitive parameter is the natural death rate of citrus psyllid, µ v . Other important parameters include infection rate from infected citrus trees to susceptible citrus psyllid, α, the infection rate from infected citrus psyllid to susceptible citrus trees, β, and the disease-induced mortality of the infected trees, γ. Since γ R 0 β = 0.5, decreasing (or increasing) β by 10% increases (or decreases) R 0 by 5%. Similarly, as γ R 0 γ = −0.4417, increasing (or decreasing) γ by 10% decreases (or increases) R 0 by 4.417%.
Next, we use sensitivity analysis to discover parameters that have a high impact on the basic reproductive number R 0 , and should be targeted by intervention strategies. Figure 1 presents our sensitivity analysis on R 0 , which involved computing the PRCC of R 0 using the LHS method. We see from Figure 1 that the natural mortality of psyllid population µ v has the most effect on R 0 of all the constant parameters. This occurs because the parameter is involved in both directions of transmission: from a tree to a psyllid and vice versa. In addition, we can see that the parameters α, β, Λ, r and K have the positive impact on R 0 . For example, a decrease in β decreases R 0 , whereas other parameters (i.e., µ h , µ v , γ, τ) are negatively correlated with R 0 .
Identification of these key parameters is important to the formulation of effective control strategies necessary for reducing the prevalence of the disease in the grove. In particular, the results of this sensitivity analysis suggest that a strategy that reduces the parameters with positive impact on R 0 will adequately reduce the spread of HLB in the grove. Furthermore, a strategy that increases the parameters with negative impact on R 0 will be effective in curtailing the spread of HLB in the grove.  Figure 1. Sensitivity analysis of the basic reproduction number R 0 .
Therefore, it follows from sensitivity analysis that the disease in the grove can be effectively reduced using control strategies: • Nutrient solution injecting u 1 (t) (i.e., decrease the value of the parameter β); • Infected trees removal u 2 (t) (i.e., increase the value of the parameter γ); • Insecticides spraying u 3 (t) (i.e., increase the value of the parameter µ v ).

Optimal control problem
The results of the sensitivity analysis suggest that control strategies that reduce infection rate (β), increase disease-induced mortality of the infected trees (γ) and psyllid natural death rate (µ v ) will be effective in curtailing the spread of HLB in the grove. However, large quantities use of insecticides and removal of infected trees in large numbers have been associated with the increased cost, which is an important problem in sustainable control of HLB. For this, there is an urgent need to find an optimal control strategy in terms of possible combination of strategies to prevent the spread of citrus HLB while minimizing the implementation cost. In this section, we suggest the following delayed control system with three-time-dependent control variables u 1 (t), u 2 (t), and u 3 (t): with the initial conditions where ς 1 , ς 2 and ς 3 denote the maximum injecting rate of nutrient solution, removing rate and killing rate, respectively, and the control functions u(t) = (u 1 (t), u 2 (t), u 3 (t)) belong to the control set U defined by Here T is the control period. For the optimal control problem of system (4.1), we now consider the following objective functional Generally, cost function for each control is assumed to be quadratic. The parameters A 1 , A 2 and B 1 , B 2 , B 3 represent the desired weights on the benefit and cost. The objective functional J formulates the optimization problem of identifying the most effective strategies, which involves the minimization of the number of the infected citrus trees and psyllids and total implementation cost over a finite time interval [0, T ]. The goal of the optimal control problem is to seek an optimal control (u * 1 , u * 2 , u * 3 ) such that where the control set U is defined in (4.3) and subject to control system (4.1) with initial conditions (4.2).
In what follows, we use Pontryagin's Maximum Principle [35] to solve this optimal control problem with time delay.

Existence of optimal control
Note that non-negative bounded solutions to the state system exists for the bounded Lebesgue measurable controls and non-negative initial conditions. We assume not only the control variables u 1 (t), u 2 (t) and u 3 (t) satisfy constraints (4.3) but also S h (θ), I h (θ), S v (θ), I v (θ) are non-negative continuous functions for θ ∈ [−τ, 0] and Then system (4.1) can be rewritten in matrix form as follows: . It is easy to see that system (4.1) is a nonlinear system with bounded coefficients. Denote and It then follows from (4.6) that there exist positive constants L 1 and L 2 such that where L 1 and L 2 are independent of state variables S h (t), I h (t), S v (t) and I v (t). Consequently, we obtain where L = max L 1 , L 2 , A < ∞, which shows that G is uniformly Lipschitz continuous, hence by definition of U with the restriction on S h (t), I h (t), S v (t), I v (t) ≥ 0, we can see that the solution of system (4.1) exists following the Theorem 4.1 and its corresponding corollary by Fleming and Rishel [36]. Theorem 4.1. There exists an optimal control (u * 1 , u * 2 , u * 3 ) ∈ U such that subject to the control system (4.1) with the initial conditions (4.2).
Proof. The existence of an optimal control problem can be verified by using the result in [36]. It is obvious to see that both the control variables and the state variables are non-negative values, which satisfies the necessary convexity of the objective functional. Note that the set of all the control variables u 1 (t), u 2 (t), u 3 (t) ∈ U is convex and closed by definition. The optimal system is bounded which determines the compactness needed for the existence of the optimal control. Also, the integrand in the functional (4.4), i.e., A 1 I h + A 2 I v + B 1 2 u 2 1 + B 2 2 u 2 2 + B 3 2 u 2 3 , is convex on the control set U. Besides, we can see that, there exists constants η 1 , η 2 > 0 and ρ > 1 such that which completes the existence of an optimal control.

Characterization of optimal control
Based on the well-known technique Pontryagin's Maximum Principle, we manage to obtain the optimal solution of the control system (4.1). For convenience, we set: To determine the precise formulation of our optimal control, we consider the optimal control problem (4.1)-(4.4) to find the Lagrangian function as well as Hamiltonian function. The Lagrangian function L is given by
For simplicity, we introduce a indicative function Notice that a necessary condition for optimality problems can be found in [35]. If we consider X(t) and X τ (t) defined above, then there exists a continuous function λ(t) on [0, T ] satisfying equations given below, namely, the state equationsẊ the optimal condition 0 = H u (X, X τ , u, λ)(t), (4.11) and the adjoint equation where H λ , H u , H X and H X τ denote the derivatives of the Hamiltonian H with respect to λ, u, X and X τ , respectively. Now, applying the necessary condition to Hamiltonian H in (4.8) with X = X * , following result hence obtain.
be an optimal state solution associated with the optimal control variable u * = (u * 1 , u * 2 , u * 3 ) for the optimal control problem (4.1)-(4.4). Then there exist four adjoint variables λ i (t), (i = 1, 2, 3, 4) satisfying the adjoint equations as follows (4.14) Moreover, the optimal controls u * 1 , u * 2 and u * 3 are characterized in the following Proof. The proof follows from Pontryagin's Maximum Principle. By using the adjoint equation (4.12) and differentiating the Hamiltonian (4.8) with respect to X(t) = X * (t) and X τ (t) = X * τ (t), we obtain (4.16) Then substituting the corresponding derivatives H X (t) and H X τ (t) (here X = S * h , I * h , S * v , I * v ) into (4.16), we can obtain the adjoint equations (4.13). Further, from the optimal conditions (4.11), we have Taking into account the property of control set in (4.3), we obtain This can be rewritten in compact notations: This completes our proof.
Formulae represented in (4.15) are the characterization of the optimal controls. In addition, the second derivative of the Lagrangian L with respect to u * 1 (t), u * 2 (t) and u * 3 (t) is positive, which shows that the optimal problem is minimum at controls u * 1 (t), u * 2 (t) and u * 3 (t). Further, the optimality system, which consists of the adjoint equations (4.13) and the state equations (4.1), the transversality conditions (4.14) with the explicit representations (4.15) for the controls is given below.

Uniqueness of the optimality system
The optimal control depends on adjoint variables and state variables, thus the uniqueness of the optimal control derives from the uniqueness of the optimality system. In this subsection, we attempt to prove that the optimal system has a unique solution. Note that S h (t), I h (t), S v (t), I v (t) are bounded, then the adjoint equation (4.13) has a bounded right-hand side which is dependent on the final time T . Therefore, there exists a positive constant M, which is dependent on the coefficients of the state equation and the uniform bound for S h (t), I h (t), S v (t) and I v (t) such that | λ i |< MT on [0, T ]. It follows from (4.15), (4.18) and (4.19) that To prove our main result, we firstly show the following two claims are valid. Claim I. There exist positive constants C x and D x such that Without loss of generality, we take x = v to analyze. Other cases can be discussed similarly. Substituting (4.18)-(4.21) into the sixth equation of system (4.17), we have Subtracting the above two equations, we have Multiplying (4.22) by v −v and then integrating from 0 to T , we can derive the following equality by transversality conditions (4.14) .

(4.24)
Note that e m 0 t−µ h τ < e m 0 T , 0 < χ [0,T −τ] (t) < 1 for 0 < t < T , (4.24) yields where In the following, we will calculate formulae I 1 -I 4 and enlarge them into the form of sum of squares by using the inequality ab ≤ (a 2 + b 2 )/2. It is obvious to show that (4.26) Besides, we can see that there exists a positive constant M 1 such that (4.27) Next we will specifically analyze formula I 3 . In view off (t + τ) where M 2 , M 3 are the upper bounds forf (t + τ) and j(t + τ), respectively. Further, by using the inequality ab ≤ (a 2 + b 2 )/2, we have (4.29) By similar above discussion, we can yield It follows from (4.28), (4.29) and (4.30) that Similarly, we can obtain that there exist positive constants M 4 and M 5 such that where C and D depend on all coefficients and bounds on all solution variables. It is easy to see that there exist m 0 > 0 and T > 0 sufficiently small such that m 0 − C − De m 0 T ≥ 0 holds true, that is, For the choice of m 0 , we have T < ln m 0 −C D . Therefore, the solution to the optimality system (4.17) is unique provided that T sufficiently small. This completes the proof.

Numerical simulation
The optimality system (4.17), consisting of the state and the adjoint equations, is solved numerically by using the forward-backward sweep numerical method [37] in the software Matlab. The process begins with an initial estimate for the control pair (u 1 , u 2 , , respectively, using the transversality conditions (4.14) and the approximated solution of the state system. Both the state and adjoint values are then used to update the control, and the process is repeated until the current state, adjoint, and controls values converge sufficiently [38].
Hence, to illustrate the optimal control strategy, different control intervention schemes are compared. The parameter values tabulated in Table 1 and the estimated values of weight constants in the objective functional, are used in the numerical simulations. Based on the actual application, we do not consider single control strategies but focus on evaluating the population-level effectiveness of the following coupled and threefold control combination strategies on the transmission of HLB : • Strategy I: Combined nutrient solution and infected trees removal strategy (i.e., u 1 ≥ 0, u 2 ≥ 0, u 3 = 0).
The system (4.1) is now simulated to assess the impact of different control strategies on HLB transmission between citrus tree population and psyllid population. For the aforementioned strategies, the optimal solutions for system (4.1) are solved over a time period of 60 months (i.e., T = 60 months) and the initial conditions are taken as follows: S h (0) = 1650, I h (0) = 50, S v (0) = 250 and I v (0) = 5. In the simulations, the maximum injecting rate, removing rate and killing rate are taken to be ς 1 = 0.3, ς 2 = 0.6 and ς 3 = 0.8. For the weight factors we choose A 1 = 1000, A 2 = 1, B 1 = 200, B 2 = 1000 and B 3 = 100. It should be pointed out that the weights in the simulations here are only of theoretical sense to illustrate the control strategies proposed in this section. The basic reproductive number in the absence of intervention is given as R 0 = 2.2505. This indicates that if not controlled, HLB will spread rapidly and be epidemic in the grove.
The results of the optimal control simulations of the HLB model (4.1) are presented in Figures 2, 3, 4 and 5. Firstly, we observe from Figure 2 that the infected citrus trees I h and infected psyllids I v increase rapidly without control, while with control, they decrease continuously and tend to be zero eventually. In other words, HLB will die out based on the time-varying optimal control strategy whereas it would break out without control.
Secondly, we find that the numbers of susceptible citrus trees S h , infected citrus trees I h , susceptible psyllids S v and infected psyllids I v under different control strategies are different (see Figure 3). We can see from Figure 3 that the control Strategies III and IV have a very desirable effect upon the population of the citrus trees and the citrus psyllids. To be specific, it is clear that the number of infected individuals reduces significantly in presence of the controls u 1 , u 2 and u 3 , or u 2 and u 3 . This indicates that nutrient solution has little impact on the prevention of HLB. Moreover, after computing total optimal cost and the total infection averted under different strategies (the values are listed in Table 3), we can conclude that the total infection averted is the largest and the total cost is the lowest when implementing Strategy IV. Therefore, for our HLB model, minimizing the total cases of the disease with minimum total cost can be achieved by optimal control strategy IV. We may conclude that Strategy IV is the optimal strategy among the four strategies and it appears to offer a promising measure for HLB control. Additionally, we can see that the total cost of implementing Strategy IV is dramatically less than that of implementing Strategy I or II, and the total infection averted of implementing Strategy IV is more than that of implementing Strategy I or II. This implies application of insecticides and removal of infected trees play an important role in managing HLB transmission. The corresponding time-dependent controls u 1 (t), u 2 (t) and u 3 (t) are depicted in Figure 4. We can see that the optimal control variable u 1 starts at 0.15, which then increases to 0.96 rapidly in the first 6 months before gradually decreasing to the lower bound 0 at the end of the simulation period. While the optimal control variable u 2 starts at the upper bound 1 for over 10 months and slowly decreases to the lower bound 0 at the end of the simulation period. The optimal control variable u 3 is similar to u 2 , it starts at the upper bound 1 for over 20 months and slowly decreases to the lower bound 0 at the end of the simulation period. These results suggest that, in order to economically and effectively control the spread of HLB, nutrient solution control should be only implemented in the start for a short time, while the infected trees removal control and the insecticide control are maintain at high levels for a long time.  Finally, we present Figure 5 to investigate the effect of delay. Apparently, the longer time delay results in the higher total cost for the stationary terminal time. And we can also see from Figure 5 that there is an increase in the number of infected trees I h and infected psyllids I v with the increase of time delay τ. Thus, the longer latent delay leads to HLB control be more complicated, which coincides with the reality. This also implies that ignoring time delay would underestimate the transmission risk of HLB.

Conclusions
In this paper, we present a new deterministic delayed vector-borne plant model considering disease latency to analyze the spread of HLB, which healthy, infected individuals in the citrus trees and healthy, infected individuals vector populations are taken into consideration. We derive the basic reproductive number R 0 and analyze dynamics behavior such as the existence of the equilibria and their stability on the basis of basic reproductive number in two cases τ = 0 and τ > 0. It shows that the disease-free equilibrium E 0 of model (2.1) is globally asymptotically stable if R 0 < 1 and unstable if R 0 > 1, and the endemic equilibrium E * is locally asymptotically stable if R 0 > 1 in the case of τ = 0; whereas stability changes occur through Hopf bifurcation in our delayed system provided that τ > 0.
Optimal control theory is then applied to investigate the optimal strategy for HLB transmission using three time dependent control variables determined from sensitivity analysis. To the best of our knowledge, no HLB modeling studies have formulated models that incorporate time delay and integrated control. The novel model developed in the current study to explore an optimal integrated strategy for a delayed HLB model considering the disease latency. By using Pontryagin's Maximum Principle, we obtain the exact optimal control formula and prove the uniqueness of optimal control solution. Numerical simulations illustrate the effectiveness of the proposed control problem. The following results are observed from our numerical simulations: (i) effective management of psyllid population will be helpful to curtail the spread of citrus HLB; (ii) the application of three time-dependent controls u 1 (t), u 2 (t) and u 3 (t) obtained from the sensitivity analysis appears to offer a promising measure for HLB control; (iii) threefold control strategy is superior to the coupled control strategies in reducing the spread of HLB; (iv) ignoring time delay would underestimate the transmission risk of HLB.
As we all know, application of insecticides is the most widely followed option for reducing citrus psyllid population at present. However, some existing literature indicates that the over-reliance on insecticides, their injudicious use and intense selection pressure has resulted in the development of varying levels of resistance in the citrus psyllid population [39,40]. Therefore, there is an urgent need to assess the impact of intensity of insecticides resistance on the transmission dynamics of HLB and explore resistance management strategy that maintains the effective use of all classes of insecticides for citrus psyllids control. We leave this interesting but challenging problem for future investigation.
Additionally, in [41], the authors illustrate that the seasonal population dynamics of citrus psyllid could be described by trimodal curve. This means the growth rate of psyllid population is strongly impacted by environmental conditions, especially temperature. Therefore, some important features, including the role of temperature-dependent infection rates as well as recruitment of psyllids should be considered in our future work.