Analysis and optimal control of a Huanglongbing mathematical model with resistant vector.

Huanglongbing (HLB) is an incurable disease that affects citrus trees. To better understand the transmission of HLB, the mathematical model is developed to investigate the transmission dynamics of the disease between Asian citrus psyllid (ACP) and citrus trees. Through rigorous mathematical derivations, we derive the expression of the basic reproduction number (R 0) of HLB. The findings show that the disease-free equilibrium is locally asymptotically stable if R 0 < 1, and if R 0 > 1 the system is uniformly persistent. By applying the global sensitivity analysis of R 0, we can obtain some parameters that have the greatest influence on the HLB transmission dynamics. Additionally, the optimal control theory is used to explore the corresponding optimal control problem of the HLB model. Numerical simulations are conducted to reinforce the analytical results. These theoretical and numerical results provide useful insights for understanding the transmission dynamics of HLB and may help policy makers to develop intervention strategies for the disease.

Introduction

Huanglongbing (HLB) is one of the most destructive diseases of citrus, which is caused by the bacterium Candidatus Liberibacter asiaticus (CLas) and vectored by the Asian citrus psyllid (ACP) (Bovè, 2006; Gottwald, 2010; Li et al., 2015). It is regarded as the most devastating citrus disease worldwide, and was first found in southern China in 1919, and now it is transmitted in fifty different countries (Bovè, 2006; Wang & Trivedi, 2013). It causes substantial economic burdens to individual growers, citrus industries and governments (Taylor et al., 2016). Recently, citrus production has endured serious economic losses due to the occurrence of citrus HLB disease. This disease reduces fruit quality and yield, infringes on citrus health and has become one of the biggest challenges for citrus growers all over the world (Arredondo Valdés et al., 2016). The United States, Brazil and China, the world's largest citrus producers, are all threatened by HLB disease. Up to now, there is still no effective method to eradicate HLB, and once orchards are infected with the disease, they are often destroyed due to low yields (McCollum & Baldwin, 2017). In view of the serious impact of HLB on the agricultural economy, it is of great importance to understand HLB transmission dynamics between citrus trees and ACP and to implement some effective intervention strategies to curb its transmission (Chiyaka et al., 2012).

Mathematical models have been playing a vital role in understanding HLB disease transmission dynamics and also in decision making processes regarding intervention mechanisms for the disease control (Chowdhury et al., 2019; Jackson & Chen-Charpentier, 2019; Jeger et al., 2018; Li et al., 2020; Meng & Li, 2010; Zhang et al., 2021; Zhang & Georgescu, 2015; Zhao et al., 2017). However, despite their significant impact on the agricultural economy, very few mathematical models have investigated how HLB disease is transmitted between the insect vector and the host (Taylor et al., 2016; Chiyaka et al., 2012; Lee et al., 2015; Vilamiu et al., 1479; Jacobsen et al., 2013; Gao et al., 2018; Khan et al., 2021). For the goal of highlighting the importance of flush for ACP dynamics, Chiyaka and Halbert (Chiyaka et al., 2012) formulated a mathematical model to study the spread of HLB within a tree because the infection payllid transmits among the different flush patches on the tree. To explore the impact of seasonal fluctuations on HLB disease, Gao and Yu (Gao et al., 2018) built an impulsive switching model for the disease with seasonal fluctuations. Jacobsen and Stupiansky (Jacobsen et al., 2013) investigated a deterministic model for HLB transmission within a single citrus garden, which included the intervention strategy of roguing. Other further studies include Tu and Gao (Tu et al., 2019) developed a vector-borne disease model with stage structure and analyzed the effect of the measure in controlling the spread of HLB. We review those models that have been applied to HLB here because they show the main insights models have provided for this disease system.

Although some of the studies mentioned above took into account different mathematical models of HLB and disease intervention strategies, they did not consider the optimality of these intervention strategies and the resistance of the vector ACP, which may sometimes be limited by resource availability. In particular, a comparative analysis to understand the costs of different intervention strategies is important for decision makers who are often faced with the resource allocation challenges. In view of this, the application of optimal control theory can serve as a useful tool for assessing the effectiveness of various policies and interventions relative to the cost of implementing them (Alzahrani et al., 2021; Khan & Fatmawatic, 2021; Okosun et al., 2013). Based on these considerations, in this paper, we formulate a compartmental model to investigate HLB transmission dynamics with resistant ACP in a single orchard of citrus trees. Through extensive calculations, we obtain the basic reproduction number, i.e., the HLB disease below this threshold disappears but outbreaks above this threshold. The persistence of the system is further explored. Subsequently, optimal control theory is used to study the effectiveness and optimality of all possible combinations of two interventions for HLB disease, i.e., removing trees with HLB symptoms and spraying insecticides. We derive the optimal control conditions by using the Pontryagin s Maximum Principle (Pontryagin, 2018). Furthermore, we obtain through numerical experiments that the optimal control strategy outperforms the constant intervention strategy in reducing the prevalence of the diseased citrus trees, and the cost of implementing the optimal control is far lower than the constant intervention strategy.

The remaining of this paper is organized as follows: the deterministic mathematical model is proposed to study the HLB transmission dynamics which incorporates ACP and citrus trees in Section 2. The basic reproduction number and local stability of the eight-dimensional HLB system are obtained in Section 3. Furthermore, the uniformly persistent of the nonlinear system are explored. The optimal control problem of the compartment model is conducted in Section 4. The global sensitivity and uncertainty analysis are performed, and some suggestions for HLB intervention strategies from the proposed model are summarized in Section 5. Finally, this paper ends with a brief conclusion and the potential outlook for future work in Section 6.

Model formulation

According to the work of (Hu et al., 2018), it was reported that trunk injection of penicillin, streptomycin and oxytetracycline hydrochloride, were effective in reducing the concentration of CLas and slowing down progress of HLB disease. In this paper, the model analyses the transmission dynamics of HLB with fraction of susceptible vaccinated citrus tress. Therefore, the citrus tree population given as N(t) is divided into susceptible individuals (S(t)), vaccinated individuals (V(t)), exposed individuals (E(t)), and infected individuals (I(t)), at any time t. Thus, N(t) = S(t) + V(t) + E(t) + I(t). The ACP population, denoted by N(t), is divided into susceptible sensitive ACP (S(t)), susceptible resistant ACP (S(t)), infected sensitive ACP (I(t)), and infected resistant ACP (I(t)), at any time t. Thus, N(t) = S(t) + S(t) + I(t) + I(t).

We assume that the recruitment to the citrus population is by replanting at a rate proportional r to the difference between the actual number of citrus trees present N and maximum population size K. Resistant classes S and I(t) are increased by sensitive vectors S and I(t) who develop resistance at the rates φ1 and φ2, respectively. Moreover, let η the force of infection from ACP to trees and η the force of infection from trees to ACP. The HLB model incorporating resistance for ACP population is represented as a system of first order nonlinear ordinary differential equations as follows:where

A schematic diagram of the model is depicted in Fig. 1, and the state variables and the parameters are described in Table 1.

Fig. 1

Schematic Diagram of the HLB model (2.1).

Table 1

Description of The variables and their parameters for the HLB model (2.1).

Variable	Description
Sh	Population of susceptible citrus trees
Vh	Population of vaccinated citrus trees
Eh	Population of exposed citrus trees
Ih	Population of infected citrus trees
Sv1	Population of susceptible sensitive ACP
Sv2	Population of susceptible resistant ACP
Iv1	Population of infected sensitive ACP
Iv2	Population of infected resistant ACP
Parameter	Description
K	Environmental carrying capacity of citrus trees
r	Replanting rate of citrus trees
ωh	Vaccine waning rate
bv	ACP biting rate
αh	Transmission probability from Iv1 to Sh
βh	Transmission probability from Iv2 to Sh
vh	Vaccination rate
μh	Natural mortality rate in citrus trees
θ	Vaccine efficacy
δh	Disease progression rate of infectious of exposed citrus trees
γh	Disease related death rate
Λ1	Sensitive ACP recruitment rate
Λ2	Resistant ACP recruitment rate
αv	Transmission probability from Eh to Sv1 and Sv2
βv	Transmission probability from Ih to Sv1 and Sv2
μv	Natural mortality rate of ACP
φ1	Mutation rate from susceptible sensitive ACP to susceptible resistant ACP
φ2	Mutation rate from infected sensitive ACP to infected resistant ACP

The basic qualitative properties of the model (2.1) will be explored in the subsequent section.

Model analysis

To better organize the analysis, we simplify some terms in the model by setting d1 = r + μ, d2 = ω + μ, d3 = v + μ, d4 = δ + μ, d5 = γ + μ, d6 = φ1 + μ, d7 = φ2 + μ, d8 = r + v + μ and d9 = v + ω + μ. We denote

Denote . Let the initial data Γ(0) ≫ 0. Then the solutions Γ(t) of the model (2.1) are non-negative for all t > 0. Furthermore

The region is positively-invariant for the model (2.1) with non-negative initial conditions.

Proof. Let the initial data of the model (2.1) Γ(0) ≫ 0. It is obvious that , so we have S(t) ≥ 0, for all t > 0. Similarly, we can get V(t) ≥ 0, E(t) ≥ 0, I(t) ≥ 0, S(t) > 0, S(t) > 0, I(t) ≥ 0 and I(t) ≥ 0, for all t > 0.

Adding the first four equations of model (2.1), the total number of citrus trees N satisfiesthus,

This implies that

Adding the fifth to eighth equations of model (2.1), the total number of ACP N satisfiesthus,

Further, we know that , if and , if . Therefore, the region G is positively-invariant.

Clearly, the model (2.1) has a disease free equilibrium (DFE) point, , whereDenote

Then the Jacobian matrix of system (2.1) with respect to E0 iswhere , , , ,

, , , .

Denote

To better organize the analysis, we denote p0 = μd4d5d7, p1 = k5d5μ(k1 + k3), p2 = k5d5φ2(k2 + k4), p3 = k6δμ(k1 + k3), p4 = k6δφ2(k2 + k4), p5 = k7d5(k2 + k4)d7, p6 = k8δ(k2 + k4)d7, and

All eigenvalues of the matrix J(E0) (3.1) have negative real parts if and only if R0 < 1.

Proof. The characteristic equation of J11 is the following polynomial equation

According to Vieta theorem, the two roots of equation (3.3) are negative. It is easy to see that two eigenvalue of J33 are − μ and − d6, which are negative. Denote

We can see that all eigenvalues of the matrix M are the roots of the following quartic polynomial equation

Obviously, if R0 < 1, then B > 0, C > 0 and D > 0. The Hertz determinants of the first to fourth order polynomial are as follows

It is obvious that Δ2 > 0, Δ3 > 0 and Δ4 > 0 if R0 < 1.By the criterion of Routh-Hurwitz, we obtain that all eigenvalues of M have negative real parts if and only if R0 < 1. Therefore all eigenvalues of J(E0) have negative real parts if R0 < 1.

Lemma 3.2 implies the following result holds.

For system (2.1), the disease-free equilibrium E0 is locally asymptotically stable if R0 < 1 and unstable if R0 > 1.

In order to analyze the permanence of system (2.1), we first give the following lemma.

(See (Zhao, 2003)) Assume that

f(X0) ⊂ X0 and f has a global attractor Ω.

The maximal compact invariant set Ω∂ = Ω ∩ M∂ off in ∂X0, possibly empty, has an acyclic covering with the following properties:

Mi is isolated in .

Ws(Mi) ∩ X0 = Φ for each 1 ≤ i ≤ k.

Then f is uniformly persistent with respect to (X0, ∂X0), i.e., there is η > 0 such that for any compact internally chain transitive set L with L⊄{Mi, for all 1 ≤ i ≤ k}, .

If R0 > 1, then system (2.1) is uniformly persistent.

Proof. We define X = {(S, V, E, I, S, S, I, I) ∈ R8: S ≥ 0, V ≥ 0, E ≥ 0, I ≥ 0, S ≥ 0, S ≥ 0, I ≥ 0, I ≥ 0}, X0 = {(S, V, E, I, S, S, I, I) ∈ X: E > 0, I > 0, I > 0, I > 0}, and ∂X0 = X∖X0.

In order to show that system (2.1) is uniformly persistent, we only need to prove that ∂X0 repels uniformly the solutions of X0. Obviously, both X and X0 are positively invariant and ∂X0 is relatively closed in X and (2.1) is point dissipative. Let

We now show that

Assume (S(0), V(0), E(0), I(0), S(0), S(0), I(0), I(0)) ∈ M. It suffices to show that E(t) = 0, I(t) = 0, I(t) = 0 and I(t) = 0 for all t ≥ 0. Suppose not. Then there exists a t0 ≥ 0 such that E(t0), I(t0), I(t0), I(t0) at least one is greater than 0. Without loss of generality, we only discuss the case E(t0) > 0, S(t0) = 0, V(t0) = 0, I(t0) = 0, S(t0) = 0, S(t0) = 0, I(t0) = 0 and I(t0) = 0. Sinceit follows that there is an ϵ0 > 0 small enough such that S(t) > 0, S(t) > 0, E(t) > 0 and I(t) > 0 for all t0 < t < t0 + ϵ0. Furthermore, let , then we have S(t1) > 0, S(t1) > 0, E(t1) > 0 and I(t1) > 0. If I(t1) > 0, thenthis implies that I(t) > 0 for all t ≥ t1. If I(t1) = 0, we get

It then follows that there is an ϵ1 () such that I(t) > 0 for all t1 < t < t1 + ϵ1. Similarly, we can show that there exists an ϵ2 () such that I(t) > 0 for all , respectively. Thus, for all we have E(t) > 0, I(t) > 0, I(t) > 0 and I(t) > 0. This contradicts the assumption that (S(0), V(0), E(0), I(0), S(0), S(0), I(0), I(0))∉M. Thus (3.5) holds.

It is obvious that is the unique equilibrium in M. We now demonstrate that E0 repels the solutions in X0.

By the proof of Lemma 3.2, we know that det(M) < 0 provided that R0 > 1, here M is defined in (3.4). Therefore, if R0 > 1, then there exist ζ1 > 0 and ζ2 > 0 small enough such thatandhold, where , , , , , , , , , , , and

Suppose (S(t), V(t), E(t), I(t), S(t), S(t), I(t), I(t)) is a solution of system (2.1) with initial value (S(0), V(0), E(0), I(0), S(0), S(0), I(0), I(0)) ∈ X0. We now claim that

Suppose that the claim is not valid. Then there is a T0 > 0 such that

It follows from system (2.1) and (3.10) that

Thus, there exists a T1 > T0 such that

From the first and second equations of the model (2.1) and (3.11), we havefor t ≥ T1. Consider the following comparison system

We can restrict ζ2 to be small enough such that system (3.12) admits a unique positive equilibrium , which is globally asymptotically stable (the proof is stated in Appendix A). By the comparison principle in differential equations and (3.8), there is a T2 > 0 such that

From the fifth equation of system (2.1) and (3.11), we getfor t ≥ T1. According to the comparison principle and (3.8), there exists a T3 > T2, such that

It follows from the sixth equation of system (2.1), (3.11) and (3.14) thatfor t ≥ T3. Then, there exists a T4 > T3, such that

Consequently, from system (2.1) and (3.13)–(3.15), we obtain that for t > T4,

Consider the following auxiliary system:

Obviously, system (3.16) has zero equilibrium point. By simple calculation, we obtain that the Jacobian matrix of system (3.16) at zero equilibrium point is (defined in (3.7)). Since admits positive off-diagonal elements, by Perron-Frobenius Theorem, we can see that there exists a positive eigenvector for the maximin eigenvalue of . Extensive calculations yield that λ > 0 since (3.6) and (3.7) hold. Consequently, we can get , , , . Furthermore, by the comparison principle, we obtain , , , . This contradicts E(t) ≤ ζ2, I(t) ≤ ζ2, I(t) ≤ ζ2, I(t) ≤ ζ2, for all t ≥ T0. Hence, (3.9) holds and the claim is proved. This indicates W(E0) ∩ X0 = Φ. Obviously, every forward orbit in M converges to E0. It follows from Lemma 3.3 that system (2.1) is uniformly persistent with respect to (X0, ∂X0). This completes the proof.

Optimal control

Insecticide applications and removing infected trees are two tactics being employed for the control of HLB (Bovè, 2006; Gottwald, 2010). There is an urgent need to explore 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 paper, we introduce into the transmission model (2.1) a time dependent control variable u(t), representing removing effort of infected citrus trees, and also consider a time dependent control variable u(t), representing killing effort of sensitive ACP and resistant ACP. Thus, model (2.1) with the time dependent control variables u(t) and u(t) becomes:where ρ1, ρ2, ρ3 is the control intensity coefficient. Our goal is to minimize the cost function defined assubject to system (4.1), where t is the control period. This performance specification involves minimizing the numbers of exposed and infected citrus tree, along with the cost of applying the controls u(t) = (u(t), u(t)). The quadratic costs have been frequently used (Agusto, 2013; Pei et al., 2016; Zhang et al., 2020). The coefficients, A, i = 1, …, 5 represent the desired weights on the benefit and cost the controls u(t) = (u(t), u(t)) is a bounded Lebesgue integrable functions (Agusto & Lenhart, 2013; Pei et al., 2016). And we need to find an optimal control , such thatwhere the control set,where u(t) is Lebesgue measurable. The necessary conditions that an optimal control quintuple must satisfy derive from the Pontryagin's Minimum Principle (Pontryagin, 2018). This principle converts (4.1) and (4.2) into a problem of minimizing pointwise a Hamiltonian H, with respect to the controls u(t) and u(t). In order to obtain the optimality conditions, we first formulate the following Hamilton functional from the objective functional (4.2) and the governing dynamics (4.1):where , , , , , , , are the associated adjoints for the states S, V, E, I, S, S, I, I. The system of adjoint equations is found by taking the appropriate partial derivatives of the Hamiltonian (4.4) with respect to the associated state and control variables.

Given an optimal control variables , and solutions , , , , , , , of the corresponding state system (4.1) that minimizes over U . Then there exists adjoint variables , , , , , , , satisfyingand with transversality conditions

The optimality conditions is given as

Furthermore, the time dependent control variables are given as

Proof. According to the result in (Fleming & Rishel, 2012), we can easily obtain the existence of an optimal control. Consequently, the differential equations governing the adjoint variables are obtained by the differentiation of the Hamiltonian function, evaluated at the optimal controls. Then the adjoint system yields:

Evaluated at the optimal controls and corresponding state variables, results in the stated adjoint system (4.4) and (4.5) given by

In addition, differentiating the Hamiltonian functional with respect to the control variables in the interior of the control set U, and then solving for controls result in the optimality conditions given as follows:

By simple calculation, we have

Taking into account the property of control set in (4.3), the characterization (4.6) can be obtained. This completes the proof.

In the next part, we will simulate numerically the solutions of the optimality system and the corresponding optimal control, and then give the interpretations from various cases.

Numerical simulation

In this section, we mainly present some numerical simulation results to verify or extend our results, which explore the impact of insecticides resistance on HLB transmission between citrus tree population and ACP population, and assess the effects of different control strategies against HLB. By using the Latin Hypercube Sampling (LHS) method (Blower & Dowlatabadi, 1994) and the forward-backward sweep method (Lenhart & Workman, 2007). Table 2 contains the value of the parameters that will be used in numerical simulations.

Table 2

Numerical values of the parameters for HLB model (2.1).

Parameter	Baseline value	Unit	Reference
K	1000	–	Zhang et al. (2020)
r	0.6	year−1	Zhang et al. (2020)
ωh	0.05	year−1	Assumed
bv	600	year−1	Assumed
αh	4.8830 × 10−4	–	Taylor et al. (2016)
βh	0.8 × 4.8830 × 10−4	–	Assumed
vh	0.7	year−1	Assumed
μh	0.04	year−1	Mingxue (2009)
θ	0.8	–	Assumed
δh	12.97	year−1	Chiyaka et al. (2012)
γh	0.1	year−1	Assumed
Λ1	3.3253 × 105	year−1	Taylor et al. (2016)
Λ2	0.8 × 3.3253 × 105	year−1	Assumed
αv	0.9 × 3.9064 × 10−4	–	Assumed
βv	3.9064 × 10−4	–	Taylor et al. (2016)
μv	5.9441	year−1	Taylor et al. (2016)
φ1	0.3	year−1	Assumed
φ2	0.24	year−1	Assumed

To examine the sensitivity of model results to the uncertainty of parameters, we do sensitivity and uncertainty analysis. Based on the parameters of model (2.1), we perform a global sensitivity analysis on the basic reproduction number R0. By using the Latin Hypercube Sampling (LHS) method (Blower & Dowlatabadi, 1994), we compute the Partial Rank Correlation Coefficients (PRCCs) of R0. The sensitivity and uncertainty analysis are presented in Fig. 2. It shows how uncertainty in model parameters may influence R0. We can observe from Fig. 2(a) that R0 is very sensitive to environmental carrying capacity of citrus trees K, ACP biting rate b, Vaccination rate v, transmission probability from ACP to citrus trees α and β, vaccine efficacy θ, disease related death rate γ, ACP recruitment rate Λ1 and Λ2, transmission probability from I to ACP β, natural mortality rate of ACP μ, but not sensitive to replanting rate of citrus trees r, natural mortality rate in citrus trees μ, disease progression rate of infectious of exposed citrus trees δ, transmission probability from E to ACP α, mutation rate from sensitive ACP to infected resistant ACP φ1 and φ2. For both transmission probabilities, β has a greater impact on R0. Thus, decreasing β is very effective in reducing R0, that is to say, decreasing the transmission probability from I to ACP plays a key role in the control of HLB. Fig. 2(b) (uncertainty analysis), we can see that about 88.7% of the distribution of R0 is greater than 1. This indicates that persistent HLB bacterial infection is likely to occur. Furthermore, by numerical calculation, we obtain the mean and standard deviation of R0 are 1.5611 and 0.4973, respectively.

Fig. 2

Sensitivity analysis and uncertainty analysis of the basic reproduction number R0. (a) shows the sensitivity indices of R0, and (b) shows histogram obtained from LHS using a sample size of 1000 for R0.

We consider the influence of two model parameter γ and μ on the basic reproduction number R0 in Fig. 3. We can observe from Fig. 3 that R0 is very sensitive to γ and μ, it is seen that as γ and μ increase, the basic reproduction number R0 decreases quickly.

Fig. 3

The basic reproduction number of model (2.1) as a function of disease related death rate γ and natural mortality rate of ACP μ. The red curve indicates that the basic regeneration number R0 is equal to one.

The dynamics of the basic reproduction number of our model as a function of insecticide resistance intensities φ1 and φ2 is depicted in Fig. 4. It is seen that as the resistance coefficients φ1 and φ2 increase, the basic reproduction number R0 decreases slowly. This indicates that the intensity of insecticide resistance present in ACP population may slightly reduce the effectiveness of control measures, which is in agreement with results of sensitivity analysis (see Fig. 2(a)).

Fig. 4

The basic reproduction number of model (2.1) as a function of resistance intensities φ1 and φ2.

In order to reduce the number of trees infected by HLB disease, two possible control measures have been considered in this paper, including removal of diseased trees and application of insecticides. The objective of this paper is to analyze the effect of three controls on the transmission of HLB, and explore the optimal solution of the optimality system (4.1) and the corresponding optimal control by using the forward-backward sweep method. The initial numbers of susceptible citrus trees and ACP are assumed that S(0) = 500, V(0) = 492, S(0) = 30000, S(0) = 20000. Then we assume that the initial numbers of infectious citrus trees and ACP are E(0) = 3, I(0) = 5, I(0) = 20, I(0) = 10, respectively. In addition, the control intensity coefficients ρ1 = 10μ, ρ2 = 20μ and ρ1 = 5μ. In the simulations, the values of other parameters are taken from Table 2.

In (Yan & Zou, 2008), the authors point out fixing the right weights in practical problems is a very difficult task which requires an amount of work on data analysis and fitting. For the weight factors in the objective function, we choose A1 = 1000, A2 = 1000, A3 = 0.01, A4 = 100 and A5 = 100. Then the weights in the simulations here are only of theoretical sense to illustrate the control strategies proposed.

The following algorithm is used to compute the optimal controls and state values using a Runge-Kutta method of the fourth order. As illustrated in (Agusto & Khan, 2018), we should first give an initial estimate for the control quintuple, then solve the state variables forward in time by using the dynamics (4.1). The results obtained for the state variables are substituted into the adjoint equation (4.6). These adjoint equations with given final conditions (4.8) are then solved backward in time. 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.

The numbers of infectious citrus trees and infectious ACP under different control levels are depicted in Fig. 5. In view of reducing the total number of the diseased citrus trees (or the diseased ACP), we can observe that the optimal control strategy is consistent with the upper bound control strategy and is superior to the constant control. The enhancement of the level of control can achieve significant effects on both the number of hosts and vectors. Moreover, it follows from Table 3 that the cost of the optimal control is less than that of the upper bound control. We may conclude that the optimal control strategy appears to offer a promising measure for HLB control.

Fig. 5

Number of infectious citrus trees and infectious ACP under different control levels, (a) exposed citrus trees, (b) infected citrus trees, (c) infected sensitive ACP, (d) infected resistant ACP.

Table 3

The costs of the objective function under different weights and control measures.

[1 pt] (A1, A2, A3, A4, A5)	J(0, 0)	J(0.2, 0.2)	J(0.4, 0.4)	Jmax(1, 1)	J_opt(uh∗,uv∗)
(1000,1000,0.01,100,100)	10524	8064	7502	6716	6703
(800,1000,0.01,1000,1000)	10424	8077	7739	8484	7517
(1000,10000,0.01,10000,1000)	9191	7521	7203	7297	6955
(500,1000,0.01,500,600)	10248	7951	7512	7498	7076
(1000,2000,0.01,1000,2000)	19567	15589	14930	15694	14052

Fig. 6 gives the optimal control profile for u and u. It shows that, the control variable u starts at the upper bound 1 for about 90% time and then gradually decreases to the lower bound 0 at the end of the simulation period. While the control variables u start at the upper bound for a shorter time and slowly decrease to the lower bound at the end of the simulation period.

Fig. 6

Optimal control law of u and u.

Fig. 7 presents the number of infectious citrus trees and infectious ACP under different control strategies. These control strategies are: (I) u and u; (II) u, (u = 0); (III) u, (u = 0). From Fig. 7, we can observe that: (1) the cost of implementing strategy (II) is lower than that of implementing strategy (III); (2) the number of infectious citrus trees is the smallest when implementing strategy (I), while the number of infectious citrus trees is the largest when implementing strategy (II); (3) in addition to strategy (I), the final number of diseased ACP is the least when applying strategy (III), and the final number of diseased ACP is the most when applying strategy (II). These results imply that multiple control strategies should be adopted simultaneously to control the spread of HLB.

Fig. 7

Number of infectious citrus treess and infectious ACP under different optimal control measures, i.e. (I) u and u; (II) u, (u = 0); (III)u, (u = 0).

The optimal control trajectories under different control measures are presented in Fig. 8. It can be seen that in the early phase of HLB outbreak, the upper bound control strategies should be adopted in most cases to suppress the transmission of the disease.

Fig. 8

Optimal control law of: (I) u, (u = 0); (II) u, (u = 0).

Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13 show the trajectories of the infectious citrus trees, the infectious ACP, and the optimal controls with different weights in objective function, respectively. Based on the simulation results, we can see that the weights have a great impact on the control cost (see Table 3), but have little or no effect on the spread of HLB. Consequently, the weights in the objective function have little impact on the optimal control strategy (Fig. 13).

Fig. 9

Comparison of the number of exposed citrus treess E(t) with different weights.

Fig. 10

Comparison of the number of infectious citrus treess I(t) with different weights.

Fig. 11

Comparison of the number of infected sensitive ACP I(t) with different weights.

Fig. 12

Comparison of the number of infected resistant ACP I(t) with different weights.

Fig. 13

Optimal control law of u and u with different weights.

Conclusion

In this paper, a deterministic mathematical model was first used to formulate the transmission dynamics of HLB between the ACP and the citrus tress, which incorporated the insecticide resistance of ACP. The stability of equilibria of the model was analyzed using the basic reproduction number R0, derived by the next generation matrix. Theoretical results have shown that the disease free equilibrium is locally asymptotically stable if R0 < 1, whereas the HLB system (2.1) is persistent if R0 > 1.

Furthermore, we used the optimization theory and the three time-dependent control variables to establish an optimal control strategy for exclusion of HLB infection by minimizing the number of exposed and infected trees. The necessary conditions for the existence of optimal solution of control problem are obtained by using Pontryagin's Maximum Principle. Finally, we verified the analytical results by numerical simulations. Based on analytical and the numerical simulation results, the main conclusions of this paper can be summarized as follows:

The intensity of insecticide resistance present in ACP population may slightly reduce the effectiveness of control measures.

The optimal control strategy is superior to the constant control strategy in decreasing the prevalence of the infected citrus trees, and the cost of implementing optimal control is much lower than that of the constant control strategy.

In the early phase of the transmission of HLB, spraying insecticides to kill ACP is more effective than other control strategies in reducing the number of the infected ACP.

The weights in the objective function have little impact on the optimal control strategy.

Data availability

All data generated or analyzed during this study are included in this article.

Declaration of competing interest

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