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

1.   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.

2.   HLB model with latent period

2.1.   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 $\beta$, a susceptible psyllid becomes infected from contact with an infected tree, which happens at a probability $\alpha$.

$(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, $\Lambda$ 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 $\tau$ as the time delay factor can be formulated as follows:

where $\mu_{h}$ and $\mu_{v}$ are the natural mortality of citrus population and psyllid population, respectively. $\gamma$ 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 $\phi$ : $[-\tau, 0]\rightarrow R_{+}^{4}$ equipped with the sup-norm

here $\phi = (\phi_{1}, \phi_{2}, \phi_{3}, \phi_{4})\in C([-\tau, 0], R)$. Therefore, the initial condition for model (2.1) is taken as below:

with

Note that the citrus tree population $N_{h}(t) = S_{h}(t)+I_{h}(t)$ and psyllid population $N_{v}(t) = S_{v}(t)+I_{v}(t)$ satisfies

respectively, which implies that

Define

Obviously, $\Omega$ is a positively invariant set for system (2.1).

2.2.   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)\equiv0$, $I_{v}(t)\equiv0$, we can easily obtain the disease-free equilibrium, given by $E_{0} = (S^{0}_{h}, \; 0, \; S^{0}_{v}, \; 0)$, where $S^{0}_{h} = \frac{rK}{r+\mu_{h}}, \; S^{0}_{v} = \frac{\Lambda}{\mu_{v}}$.

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

Denote

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

where

2.3.   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 "$\tau$" as a bifurcation parameter in this subsection.

When $\tau\geq0$, set the right-hand side of system (2.1) to zero, namely

So we have $E^{\ast} = (S_{h}^{\ast}, \; I_{h}^{\ast}, \; S_{v}^{\ast}, \; I_{v}^{\ast})$, where

Obviously $\frac{\Lambda}{\mu_{v}}-S_{v}^{\ast} > 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^{\ast}$ by using the method in ^[24] when $\tau = 0$.

Theorem 2.1. When $\tau = 0$, the disease-free equilibrium $E_{0}$ of (2.1) is globally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0} > 1$.

Theorem 2.2. When $\tau = 0$, the endemic equilibrium $E^{\ast}$ of (2.1) is locally asymptotically stable if $R_{0} > 1$.

Now we turn to discussing the stability of the endemic equilibrium $E^{\ast}$ when $R_{0} > 1$ in case of $\tau > 0$. In addition, we derive the stability conditions for $E^{\ast}$ as well as the conditions for Hopf bifurcation.

Linearizing system (2.1) at $E^{\ast}$, we get the following system

where $X = (S_{h}, \; I_{h}, \; S_{v}, \; I_{v})^{T}$, $A$ and $B$ are $4\times4$ matrices, given by

In view of $|\lambda I-A-Be^{-\lambda\tau}| = 0$, we have

It follows that

where

To show Hopf bifurcation, we must have a pair of purely imaginary roots of characteristic Eq (2.3). Suppose $i\omega\; (\omega > 0)$ is a purely imaginary root associated with Eq (2.3), thus

Separating real and imaginary parts, we have following equations

Squaring and adding (2.4) and (2.5) we obtain

Let $z = \omega^{2}$, then Eq (2.6) yields

where the coefficients are given by

Now, if the following conditions $q_{1}\geq0$, $q_{2}\geq0$, $q_{3} > 0$, $q_{4}\geq0$ 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 $\tau\geq0$ 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}\geq0$, $q_{2}\geq0$, $q_{3} > 0$ and $q_{4}\geq0$ are satisfied, then the endemic equilibrium $E^{\ast}$ of system (2.1) is local asymptotically stable for all $\tau\geq0$.

However, if the coefficients of Eq (2.7) satisfy $q_{1}\geq0$, $q_{2}\geq0$, $q_{3}\geq0$, $q_{4} < 0$, then Eq (2.7) has at least one positive real root. It means that Eq (2.6) has at least one positive real root for $z = \omega^{2}$.

Let $\omega_{0}$ be a positive real root of Eq (2.6). Then there is a purely imaginary root $i\omega_{0}$ of Eq (2.3). From Eqs (2.4) and (2.5) we obtain the values of $\cos\omega\tau$ as

Furthermore, substituting $\omega = \omega_{0}$ into (2.8), then delay factor $\tau$ can be presented

Hence $(\omega_{0}, \tau_{k})$ is the solution of Eq (2.3), that is to say, there exists a pair of purely imaginary roots $\lambda = \pm i\omega_{0}$ for Eq (2.3) when $\tau = \tau_{k}$.

Note that $\tau = \tau_{0}$ is the minimum value of $\tau$ when the purely imaginary root $\lambda = \pm i\omega_{0}$ of Eq (2.3) occurs. Based on above discussion, we state following lemmas which will be essential to our claims.

Lemma 2.1. Eq (2.3) has a pair of purely imaginary roots $\lambda = \pm i\omega_{0}$ if $q_{1}\geq0$, $q_{2}\geq0$, $q_{3}\geq0$, $q_{4} < 0$, $\tau = \tau_{0}$.

Set the characteristic value of Eq (2.3) to be $\lambda(\tau) = \alpha(\tau)+i\omega(\tau)$ satifying the condition of $\alpha(\tau_{k}) = 0$ and $\omega(\tau_{k}) = \omega_{0}$, the transversality condition is verified below.

Lemma 2.2. If the condition $f'(\omega_{0}^{2}) = 4\omega_{0}^{6}+3q_{1}\omega_{0}^{4}+2q_{2}\omega_{0}^{2}+q_{3}\neq0$ holds, then $\frac{d Re\lambda(\tau)}{d\tau}\mid_{\tau = \tau_{k}}\neq0$.

Proof. Differentiating both sides of Eq (2.3) with respect to $\tau$, we have

Obviously,

It follows from $\lambda(\tau_{k}) = i\omega_{0}$ and (2.10) that

After substituting $\lambda = i\omega_{0}$ into Eq (2.3) we can see that

Considering $|e^{-i\omega_{0}\tau}| = 1$ for $e^{-i\omega_{0}\tau} = \cos\omega_{0}\tau-i\sin\omega_{0}\tau$, we derive from (2.12) that

Consequently,

Finally, from (2.11) and (2.13) we get

This completes the proof.

Therefore, based on Hopf bifurcation theory and the above discussions we have the following theorem.

Theorem 2.4. If the conditions $q_{1}\geq0$, $q_{2}\geq0$, $q_{3}\geq0$, $q_{4} < 0$, $f'(\omega_{0}^{2})\neq0$ and $R_{0} > 1$ hold, then the endemic equilibrium $E^{\ast}$ of system (2.1) is asymptotically stable for $\tau < \tau_{0}$; the endemic equilibrium $E^{\ast}$ of system (2.1) is unstable for $\tau > \tau_{0}$; Hopf bifurcation occurs at $E^{\ast}$ for $\tau = \tau_{k}\; \; (k = 0, 1, 2, \cdots)$.

3.   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 $\beta$,

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, $\mu_{v}$. Other important parameters include infection rate from infected citrus trees to susceptible citrus psyllid, $\alpha$, the infection rate from infected citrus psyllid to susceptible citrus trees, $\beta$, and the disease-induced mortality of the infected trees, $\gamma$. Since $\mathfrak{\gamma}^{R_{0}}_\beta = 0.5$, decreasing (or increasing) $\beta$ by 10% increases (or decreases) $R_{0}$ by 5%. Similarly, as $\mathfrak{\gamma}^{R_{0}}_\gamma = -0.4417$, increasing (or decreasing) $\gamma$ 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 $\mu_{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 $\alpha$, $\beta$, $\Lambda$, $r$ and $K$ have the positive impact on $R_{0}$. For example, a decrease in $\beta$ decreases $R_{0}$, whereas other parameters (i.e., $\mu_{h}$, $\mu_{v}$, $\gamma$, $\tau$) 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.

Therefore, it follows from sensitivity analysis that the disease in the grove can be effectively reduced using control strategies:

$\bullet$ Nutrient solution injecting $u_{1}(t)$ (i.e., decrease the value of the parameter $\beta$);

$\bullet$ Infected trees removal $u_{2}(t)$ (i.e., increase the value of the parameter $\gamma$);

$\bullet$ Insecticides spraying $u_{3}(t)$ (i.e., increase the value of the parameter $\mu_{v}$).

4.   Optimal control problem

The results of the sensitivity analysis suggest that control strategies that reduce infection rate ($\beta$), increase disease-induced mortality of the infected trees ($\gamma$) and psyllid natural death rate ($\mu_{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 $\varsigma_{1}$, $\varsigma_{2}$ and $\varsigma_{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}^{\ast}, u_{2}^{\ast}, u_{3}^{\ast})$ 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.

4.1.   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}(\theta), I_{h}(\theta), S_{v}(\theta), I_{v}(\theta)$ are non-negative continuous functions for $\theta\in[-\tau, 0]$ and

Then system (4.1) can be rewritten in matrix form as follows:

where

and $V_{\tau}(t) = V(t-\tau)$. It is easy to see that system (4.1) is a nonlinear system with bounded coefficients.

Denote

and

here,

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\big\{L_{1}, L_{2}, \|A\|\big\} < \infty$, 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)\geq0$, 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}^{\ast}, u_{2}^{\ast}, u_{3}^{\ast})\in 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 $\big(u_{1}(t), u_{2}(t), u_{3}(t)\big)\in 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}+\frac{B_{1}}{2}u_{1}^{2}+\frac{B_{2}}{2}u_{2}^{2}+\frac{B_{3}}{2}u_{3}^{2}$, is convex on the control set $U$. Besides, we can see that, there exists constants $\eta_{1}$, $\eta_{2} > 0$ and $\rho > 1$ such that

which completes the existence of an optimal control.

4.2.   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

and the Hamiltonian function $H$ is

where $\lambda_{1}(t)$, $\lambda_{2}(t)$, $\lambda_{3}(t)$ and $\lambda_{4}(t)$ are the adjoint variables.

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_{\tau}(t)$ defined above, then there exists a continuous function $\lambda(t)$ on $[0, T]$ satisfying equations given below, namely, the state equations

the optimal condition

and the adjoint equation

where $H_{\lambda}$, $H_{u}$, $H_{X}$ and $H_{X_{\tau}}$ denote the derivatives of the Hamiltonian $H$ with respect to $\lambda$, $u$, $X$ and $X_{\tau}$, respectively. Now, applying the necessary condition to Hamiltonian $H$ in (4.8) with $X = X^{\ast}$, following result hence obtain.

Theorem 4.2. Let $X^{\ast} = (S_{h}^{\ast}$, $I_{h}^{\ast}$, $S_{v}^{\ast}$, $I_{v}^{\ast})$ be an optimal state solution associated with the optimal control variable $u^{\ast} = (u_{1}^{\ast}$, $u_{2}^{\ast}$, $u_{3}^{\ast}$) for the optimal control problem (4.1)-(4.4). Then there exist four adjoint variables $\lambda_{i}(t), (i = 1, \; 2, \; 3, \; 4)$ satisfying the adjoint equations as follows

with transversality conditions

Moreover, the optimal controls $u_{1}^{\ast}$, $u_{2}^{\ast}$ and $u_{3}^{\ast}$ 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^{\ast}(t)$ and $X_{\tau}(t) = X_{\tau}^{\ast}(t)$, we obtain

Then substituting the corresponding derivatives $H_{X}(t)$ and $H_{X_\tau}(t)$ (here $X = S_{h}^{\ast}, I_{h}^{\ast}, S_{v}^{\ast}, I_{v}^{\ast}$) into (4.16), we can obtain the adjoint equations (4.13). Further, from the optimal conditions (4.11), we have

which indicates

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}^{\ast}(t)$, $u_{2}^{\ast}(t)$ and $u_{3}^{\ast}(t)$ is positive, which shows that the optimal problem is minimum at controls $u_{1}^{\ast}(t)$, $u_{2}^{\ast}(t)$ and $u_{3}^{\ast}(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.

with initial conditions (4.2) and transversality conditions (4.14).

4.3.   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 $\mid\lambda_{i}\mid < MT$ on $[0, T]$.

Theorem 4.3. For $T$ sufficiently small, the solution to the optimality system (4.17) is unique.

Proof. Suppose that $(S_{h}, I_{h}, S_{v}, I_{v})$ and $(\overline{S_{h}}, \overline{I_{h}}, \overline{S_{v}}, \overline{I_{v}})$ are two distinct solutions to the optimality system (4.17). Let $m_{0} > 0$ be chosen such that

It follows from (4.15), (4.18) and (4.19) that

Denote

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

holds, where $x = \omega, \; v, \; j, \; k$, and $a_\omega = m_{0}+r+\mu_{h}$, $a_v = m_{0}+\mu_{h}+\gamma$, $a_j = m_{0}+\mu_{v}$, $a_k = m_{0}+\mu_{v}$.

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

and

Subtracting the above two equations, we have

Multiplying (4.22) by $v-\bar{v}$ and then integrating from 0 to $T$, we can derive the following equality by transversality conditions (4.14).

From (4.22) and (4.23), we get

Note that $e^{m_{0}t-\mu_{h}\tau} < e^{m_{0}T}$, $0 < \chi_{[0, T-\tau]}(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\leq (a^{2}+b^{2})/2$. It is obvious to show that

Besides, we can see that there exists a positive constant $M_{1}$ such that

Next we will specifically analyze formula $I_{3}$. In view of $\bar{f}(t+\tau)\bar{j}(t+\tau)-f(t+\tau)j(t+\tau) = \bar{f}(t+\tau)\bar{j}(t+\tau)-\bar{f}(t+\tau)j(t+\tau)+\bar{f}(t+\tau)j(t+\tau)-f(t+\tau)j(t+\tau)$, we have

where $M_{2}$, $M_{3}$ are the upper bounds for $\bar{f}(t+\tau)$ and $j(t+\tau)$, respectively.

Further, by using the inequality $ab\leq (a^{2}+b^{2})/2$, we have

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

From (4.26), (4.27), (4.31) and (4.32), we can conclude

Claim II . There exist positive constants $C_y$ and $D_y$ such that

where $y = h, \; q, \; f, \; l$, and $a_h = m_{0}+r+\mu_{h}$, $a_q = m_{0}+\mu_{h}+\gamma$, $a_f = m_{0}+\mu_{v}$, $a_l = m_{0}+\mu_{v}$.

The proof of Claim II is similar to Claim I , we omit it.

By Claim I and Claim II , we can see that

which implies 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}\geq0$ holds true, that is,

For the choice of $m_{0}$, we have $T < \ln\frac{m_{0}-C}{D}$. Therefore, the solution to the optimality system (4.17) is unique provided that $T$ sufficiently small. This completes the proof.

5.   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}, u_{3}) = (u_{10}, u_{20}, u_{30})$ over the interval $[-\tau, T]$. Then, the initial conditions in (4.2) are used to solve the state system forward in time on interval $[0, T]$. The results obtained for the state variables are plugged into the adjoint equations (4.13). These adjoint equations are then solved backward in time on intervals $[T-\tau, T]$ and $[0, T-\tau]$, 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 :

$\bullet$ Strategy $I$: Combined nutrient solution and infected trees removal strategy (i.e., $u_{1}\geq0$, $u_{2}\geq0$, $u_{3} = 0$).

$\bullet$ Strategy $II$: Combined nutrient solution and insecticide strategy (i.e., $u_{1}\geq0$, $u_{3}\geq0$, $u_{2} = 0$).

$\bullet$ Strategy $III$: Combined infected trees removal and insecticide strategy (i.e., $u_{2}\geq0$, $u_{3}\geq0$, $u_{1} = 0$).

$\bullet$ Strategy $IV$: Combined nutrient solution, infected trees removal and insecticide strategy (i.e., $u_{1}\geq0$, $u_{2}\geq0$, $u_{3}\geq0$).

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 $\varsigma_{1} = 0.3$, $\varsigma_{2} = 0.6$ and $\varsigma_{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 $\tau$. 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.

6.   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 $\tau = 0$ and $\tau > 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^{\ast}$ is locally asymptotically stable if $R_{0} > 1$ in the case of $\tau = 0$; whereas stability changes occur through Hopf bifurcation in our delayed system provided that $\tau > 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.

Acknowledgments

The research has been supported by the Natural Science Foundation of China (11961003, 11901110), The Natural Science Foundation of Jiangxi Province (No.20192ACBL20004), The Science and Technology Project of Education Department of Jiangxi Province (GJJ201406, GJJ1907, GJJ209906), and The Graduate Innovation Project of Jiangxi Province (YC2020-S600).

Conflict of interest

All authors declare no conflicts of interest in this paper.