Modeling the D. citri-HLB Pathosystem

1. Introduction

Citrus Huanglongbing (HLB) represents one of the most devastating threats to citrus production worldwide. This disease is caused by bacteria of the genus Candidatus Liberibacter spp. (variants asiaticus, americanus and africanus), the asiaticus strain being the one present in Colombia, where it is transmitted by the psyllid Diaphorina citri Kuwayama. It was first detected in the Tolima Department of Colombia in 2008 and it has since become established in at least six departments in the north of the country. HLB causes a variety of symptoms, including foliar discoloration, fruit deformation, and eventual tree decline, drastically reducing citrus crop productivity and quality, rendering fruit inedible, and ultimately destroying trees [1,2].

Globally, inadequate management strategies have caused severe regional economic impacts: annual losses of USD 1 billion in Florida, USA [3], an increase in production costs of USD 153 million in Mexico [4] and a 68% decline in Brazil’s fruit production in five years [5]. Therefore, it is important to reassess production and management strategies. In Colombia, although the disease arrived later, its presence has been confirmed in approximately 20% of the citrus production regions [6].

Control strategies for HLB focus on integrated measures, including the use of resistant or tolerant varieties, cultural management such as removal of infected trees, chemical control of the vector, orchard management practices to reduce disease incidence [3], and biological control through parasitoids specific of D. citri, such as the wasp T. radiata [7]. These practices not only aim to mitigate the spread of pathogens, but also to maintain the long-term viability of citrus production systems.

We propose a mathematical model to describe the dynamics of the D. citri-HLB pathosystem, considering some of the most common control practices. We aim to simulate integrated management strategies and to compare the effectiveness of these control practices —individually or in combination— with the goal of evaluating potential bioeconomic strategies for more effective HLB control in the future. We carry out the analysis of a reduced dynamical system model, including local and global existence and uniqueness of solutions, global preservation of positivity, and the existence and stability of up to three steady states:

• The extinction equilibrium —locally asymptotically stable when the basic reproduction number of D. citri, ${\mathcal{R}}_0^{\mathit{D}}$, is smaller than 1 (in which case it is the only steady state), and unstable when it is larger than 1;
• an equilibrium with D. citri nymphs and adults but no T. radiata that exists only when ${\mathcal{R}}_0^{\mathit{D}}>1$; and
• a coexistence equilibrium of vectors and parasitoids.

2. Materials and Methods

In this section, a nonlinear ODE system is described to model the life cycle of the psyllid summarized in just the two compartments that are important for the dynamics of citrus trees infections with HLB. They are D. citri adults and stages-3-to-5 nymphs, because later nymphal stages are considered the most efficient at acquiring and replicating the bacterium [8] and additionally these are the only stages that can be parasitized by T. radiata [9]. Furthermore, group nymphs in stages 3–5 can be merged into a single compartment because there is no appreciable difference in the interactions of those three stages of nymphs with the citrus trees, with the HLB bacteria, or with the T. radiata wasps. Given the demonstrated low bacterial uptake capacity in adult D. citri [10], it restricts the infectious D. citri adults to psyllids that acquired CLas as nymphs.

In summary, we neglect the egg, nymph-1, and nymph-2 stages of D. citri and simplify the model to collapse oviposition, transition from egg to nymph, and the first two moultings, into a single transition from adult to stage-3 nymph.

To keep the model simple, we do not consider the difference between infected and infectious D. citri, that is, we do not consider an incubation period.

Consequently, the psyllid population is divided into four compartments, with the following sizes at time $t\ge 0$:

• $N\left(t\right)$:= number of D. citri stage 3–4–5 nymphs not carrying bacteria CLas (i.e., uninfected nymphs);
• $N_i\left(t\right)$:= number of D. citri stage 3–4–5 nymphs carrying bacteria CLas (i.e., infected/infectious nymphs);
• $A\left(t\right)$:= number of D. citri adults that do not transmit the HLB.
• $A_i\left(t\right)$:= number of D. citri adults that transmit the HLB.

We also divide the citrus tree population into four compartments, with the following sizes at time $t\ge 0$:

• $T_1\left(t\right)$:= number of citrus trees with high sprouting and fruit production, and free of HLB
• ${T_1}_i\left(t\right)$:= number of citrus trees with high sprouting and fruit production, and infected with HLB
• $T_2\left(t\right)$:= number of citrus trees with low sprouting and fruit production, and free of HLB.
• ${T_2}_i\left(t\right)$:= number of citrus trees with low sprouting and fruit production, and infected with HLB

Throughout this paper, we will use the expression infected trees to mean citrus trees infected with HLB and, a fortiori, also with HLB-infected D. citri. However, uninfected trees will mean they are free of HLB though not necessarily of D. citri. If an uninfected tree has D. citri nymphs on its shoots, they must necessarily be HLB-free nymphs. The T. radiata population considered consists only of adult females because this is the stage in which they can parasitize D. citri nymphs. Thus, the dynamics of the wasp population is modeled with a single compartment, with the following size at time $t\ge 0$:

• $W\left(t\right)$:= number of adult female T. radiata

We make several assumptions in order to simplify the model while still capturing the essential features of interactions among the various species and also allowing for the different types of interventions that we want to consider:

• The per capita “oviposition" rate of adult D. citri, $\rho$, is actually the product of four factors: their natural per capita egg-laying (= oviposition) rate, the probability of survival from egg to nymph-1 stage, the probability of survival from nymph-1 to nymph-2 stage, and the probability of survival from nymph-2 to nymph-3 stage. Hence the quotation marks around ‘oviposition’.
• The per capita oviposition rate of adult D. citri is the same whether they are healthy or infected; thus the generation of new nymphs does not depend on whether they originate from a healthy or an infected adult.
• Nymphs that hatch on a healthy tree are considered healthy, while those that hatch on an HLB-infected tree are considered infected, regardless of the HLB-infection status of the adult parent. This is justified because, as far as we know, HLB is not transmitted vertically from adult female moths to their offspring to daughter.
• Because oviposition and nymph development take place only on citrus tree shoots, we assume that the availability or scarcity of the latter impacts the natural oviposition rate through a factor that is an increasing function of the total number of shoots of the corresponding infection status, with values between 0 and 1.
• Because D. citri nymphs are immobile, they do not have logistic per capita mortality driven by intra-species competition.
• Because the livelihood of adult D. citri is intimately tied to the availability of citrus tree shoots to feed on, logistic mortality from intra-species competition is assumed to occur at a rate that decreases with increasing availability of shoots of any infection status.
• Because adult D. citri who get infected with HLB from feeding on shoots from an infected tree have too short a life-span to develop a pathogen load high enough to further transmit infection to uninfected trees when feeding on their shoots, we do not consider such feeding as a source of infected adults since we only want the compartment $A_i$ to contain infectious adults and not infected-but-not-infectious.
• Because T. radiata deposes its eggs only in D. citri nymphs —only one egg per nymph— and this results in the death of the nymph, the per capita T. radiata-induced mortality rate of D. citri nymphs is assumed to be proportional to the number of T. radiata wasps.
• D. citri adults are considered infected only if they originate from infected nymphs.
• D. citri nymphs may die due to natural causes, to the effect of agrochemical products used for their control, or to parasitism by T. radiata.
• D. citri adults may die due to natural causes or as a result of agrochemical control.
• We do not include the planting of new trees. Therefore, the total tree population can only decrease.
• Trees are classified into high-fitness —producing a commercially acceptable amount of fruit— and low-fitness —with very low and commercially negligible production.
• Irrespective of their fitness, trees can develop HLB only if an infected D. citri adult feeds on them.
• The only way for a high-fitness tree to transition to low fitness is through neglect, irrespective of whether it is healthy or HLB-infected.
• A healthy low-fitness tree can recover and transition to high fitness with proper care.
• An HLB-infected low-fitness tree cannot recover and regain high fitness, because of the disease.
• Irrespective of their fitness, healthy trees produce healthy shoots and HLB-infected trees produce infected shoots.
• HLB-infected trees, whether of high or low fitness, may be eradicated as a control measure.
• T. radiata wasps experience natural mortality and may also die from the effect of agrochemicals applied for the control of D. citri.

The dynamics of the system is described by the following system of nine ordinary differential equations, four for the psyllid, four for the citrus trees, and one for the parasitoid.

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dN}{dt}}& =& \rho (A+A_i){\displaystyle \frac{b}{b+b_i}}-{\mu }_{N}N-\gamma N-\xi WN-{\delta }_N\left(t\right)N\hfill \end{array}$$

(2.1a)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{N}_i}{dt}}& =& \rho (A+A_i){\displaystyle \frac{b_i}{b+b_i}}-{\mu }_N{N}_i-\gamma N_i-\xi W{N}_i-{\delta }_N\left(t\right)N_i\hfill \end{array}$$

(2.1b)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dA}{dt}}& =& \gamma N-\left({\mu }_A+{\displaystyle \frac{C_A}{b+b_i}}(A+A_i)\right)A-{\delta }_A\left(t\right)A\hfill \end{array}$$

(2.1c)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{A}_i}{dt}}& =& \gamma N_i-\left({\mu }_A+{\displaystyle \frac{C_A}{b+b_i}}(A+A_i)\right)A_i-{\delta }_A\left(t\right)A_i\hfill \end{array}$$

(2.1d)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{T}_1}{dt}}& =& \tau T_2-\lambda T_1-{\theta }_T{\displaystyle \frac{A_i}{\eta +A+A_i}}{T}_1\hfill \end{array}$$

(2.1e)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{T_1}_i}{dt}}& =& {\theta }_T{\displaystyle \frac{A_i}{\eta +A+A_i}}{T}_1-{\lambda }_i{T_1}_i-{\alpha }_1{T_1}_i\hfill \end{array}$$

(2.1f)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{T}_2}{dt}}& =& \lambda T_1-{\mu }_{T_2}{T}_2-\tau T_2-{\theta }_T{\displaystyle \frac{A_i}{\eta +A+A_i}}{T}_2\hfill \end{array}$$

(2.1g)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{T_2}_i}{dt}}& =& {\lambda }_i{T_1}_i+{\theta }_T{\displaystyle \frac{A_i}{\eta +A+A_i}}{T}_2-{\mu }_{{T_2}_i}{T_2}_i-{\alpha }_2{T_2}_i\hfill \end{array}$$

(2.1h)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dW}{dt}}& =& \Lambda \left(t\right)+\xi \left({\displaystyle \frac{N+N_i}{\nu +N+N_i}}\right)W-{\mu }_{W}W-{\delta }_W\left(t\right)W\hfill \end{array}$$

(2.1i)

where

$$\begin{array}{ccc}b\hfill & =m_1{T}_1+m_2{T}_2\hfill & =\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l} \mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{s}\mathrm{h}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{s} \mathrm{i}\mathrm{n} \mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s} \mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e} \mathrm{o}\mathrm{f} \mathrm{H}\mathrm{L}\mathrm{B},\hfill \\ b_i\hfill & =m_{1i}{T_1}_i+m_{2i}{T_2}_i\hfill & =\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l} \mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{s}\mathrm{h}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{s} \mathrm{i}\mathrm{n} \mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s} \mathrm{i}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{H}\mathrm{L}\mathrm{B}.\hfill \end{array}$$

(2.2)

It is assumed that the average number of shoots for each of the four types of trees —$m_1$, $m_2$, $m_{1i}$, $m_{2i}$— stays constant in time. Concerning control measures, assume that agrochemicals affect D. citri the same way irrespective of their infection status, but differentially for nymphs and adults. We also assume that agrochemicals induce additional mortality for T. radiata as an undesired side effect.

Table 1. Description of parameters of the D. citri model.

Symbol	Description
$\rho$	maximal per capita production rate of D. citri 3rd-instar nymphs (per adult)
$m_1$	average number of shoots per high-sprouting uninfected tree
${m_1}_i$	average number of shoots per high-sprouting infected tree
$m_2$	average number of shoots per low-sprouting uninfected tree
${m_2}_i$	average number of shoots per low-sprouting infected tree
${\mu }_N$	natural per capita mortality rate of stages 3–5 nymphs
${\mu }_A$	natural per capita mortality rate of adult D. citri
${\mu }_{T_2}$	per capita mortality rate of low-productivity uninfected trees
${\mu }_{{T_2}_i}$	per capita mortality rate of low-productivity infected trees
${\mu }_W$	per capita mortality rate of T. radiata
$\gamma$	per capita transition rate of stage 3–5 nymphs to adulthood
$\xi$	per capita oviposition rate of T. radiata
${\delta }_N\left(t\right)$	per capita mortality rate of D. citri nymphs from agrochemicals
${\delta }_A\left(t\right)$	per capita mortality rate of adult D. citri from agrochemicals
${\delta }_W\left(t\right)$	per capita mortality rate of T. radiata from agrochemicals
$C_A$	intra-specific competition constant in adult D. citri mortality per sprout
${\theta }_T$	per capita infection rate of trees by adult D. citri when all of them are infected
$\tau$	per capita transition rate of low-productivity trees to high productivity
$\lambda$	per capita transition rate of uninfected trees from high to low productivity
${\lambda }_i$	per capita transition rate of infected trees from high to low productivity
${\alpha }_1$	per capita eradication rate of high-productivity trees with HLB
${\alpha }_2$	per capita eradication rate of low-productivity trees with HLB
$\Lambda \left(t\right)$	per capita release rate of T. radiata
$\nu$	number of nymphs in stages 3–5 resulting in 50% reduction of T. radiata
$\eta$	fraction of adult D. citri to keep tree infection rate continuous at $A=A_i=0$

Table 1. Description of parameters of the D. citri model.

Symbol	Description
$\rho$	maximal per capita production rate of D. citri 3rd-instar nymphs (per adult)
$m_1$	average number of shoots per high-sprouting uninfected tree
${m_1}_i$	average number of shoots per high-sprouting infected tree
$m_2$	average number of shoots per low-sprouting uninfected tree
${m_2}_i$	average number of shoots per low-sprouting infected tree
${\mu }_N$	natural per capita mortality rate of stages 3–5 nymphs
${\mu }_A$	natural per capita mortality rate of adult D. citri
${\mu }_{T_2}$	per capita mortality rate of low-productivity uninfected trees
${\mu }_{{T_2}_i}$	per capita mortality rate of low-productivity infected trees
${\mu }_W$	per capita mortality rate of T. radiata
$\gamma$	per capita transition rate of stage 3–5 nymphs to adulthood
$\xi$	per capita oviposition rate of T. radiata
${\delta }_N\left(t\right)$	per capita mortality rate of D. citri nymphs from agrochemicals
${\delta }_A\left(t\right)$	per capita mortality rate of adult D. citri from agrochemicals
${\delta }_W\left(t\right)$	per capita mortality rate of T. radiata from agrochemicals
$C_A$	intra-specific competition constant in adult D. citri mortality per sprout
${\theta }_T$	per capita infection rate of trees by adult D. citri when all of them are infected
$\tau$	per capita transition rate of low-productivity trees to high productivity
$\lambda$	per capita transition rate of uninfected trees from high to low productivity
${\lambda }_i$	per capita transition rate of infected trees from high to low productivity
${\alpha }_1$	per capita eradication rate of high-productivity trees with HLB
${\alpha }_2$	per capita eradication rate of low-productivity trees with HLB
$\Lambda \left(t\right)$	per capita release rate of T. radiata
$\nu$	number of nymphs in stages 3–5 resulting in 50% reduction of T. radiata
$\eta$	fraction of adult D. citri to keep tree infection rate continuous at $A=A_i=0$

3. Results

3.1. Mathematical Analysis of a Reduced Model

Let us denote the total adult D. citri population size as $x\left(t\right)=A\left(t\right)+A_i\left(t\right)$, and the total D. citri nymph population size as $y\left(t\right)=N\left(t\right)+N_i\left(t\right)$. Then, summing equations (2.1c) with (2.1d) and (2.1a) with (2.1b), we see that D. citri adults, D. citri nymphs, and T. radiata satisfy the following system of ODEs:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{dt}}& =& \gamma y-\left({\mu }_A+{\delta }_A\left(t\right)+{\displaystyle \frac{C_A}{B}}x\right)x,\hfill \end{array}$$

(3.1)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dy}{dt}}& =& \rho x-\left({\mu }_N+\gamma +\xi W+{\delta }_N\left(t\right)\right)y,\hfill \end{array}$$

(3.2)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dW}{dt}}& =& \Lambda \left(t\right)+\xi \left({\displaystyle \frac{y}{\nu +y}}\right)W-{\mu }_{W}W-{\delta }_W\left(t\right)W,\hfill \end{array}$$

(3.3)

where $B=B\left(t\right)=b\left(t\right)+b_i\left(t\right)$ is the total number of shoots on the plantation’s trees, which we will assume constant for this analysis. We shall also assume that the only control measure used is the parasitoid T. radiata, so that ${\delta }_W\equiv 0$, and that no wasps are released after an initial release that left their population’s fate up to having enough nymphs to lay their eggs in, i.e. $\Lambda \left(t\right)\equiv 0$.

3.1.1. Well-Posedness

The well-posedness of the reduced system

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{dt}}& =& \gamma y-\left({\mu }_A+{\displaystyle \frac{C_A}{B}}x\right)x,\hfill \\ \hfill {\displaystyle \frac{dy}{dt}}& =& \rho x-\left({\mu }_N+\gamma +\xi W\right)y,\hfill \\ \hfill {\displaystyle \frac{dW}{dt}}& =& \xi \left({\displaystyle \frac{y}{\nu +y}}\right)W-{\mu }_{W}W,\hfill \end{array}$$

(3.4)

follows from the general existence and uniqueness theorem for systems of ODEs. That is,

Theorem 1. 

Let $t_0\in \mathbb{R}$ and $X_0={(x_0,y_0,W_0)}^T\in {\mathbb{R}}^3$ be arbitrary, with $y_0\ne -\nu$. Then, there exists $h>0$ such that the system (3.4) has a unique solution $X\left(t\right)={[x\left(t\right),y\left(t\right),W\left(t\right)]}^T$ for $t\in (t_0-h,t_0+h)$ which satisfies $X\left(t_0\right)=X_0$.

We can also prove that the system also gives ecologically meaningful solutions. That is,

Theorem 2. 

For any initial condition $(x_0,y_0,W_0)\in {\mathbb{R}}_{+}^3$ (the non-negative cone of ${\mathbb{R}}^3$), the system (3.4) has a unique, non-negative solution $X\left(t\right)={[x\left(t\right),y\left(t\right),W\left(t\right)]}^T$ for $t\in [t_0,\infty )$ which satisfies $X\left(t_0\right)=X_0$.

Proof. 

First note that, for any function $y\in L_{loc}^1\left([0,\infty )\right)$ bounded away from $-\nu$, for any $W_0\ge 0$, and for all $t\ge 0$,

$$W\left(t\right)=W_{0}exp\left({\int }_0^t\left[{\displaystyle \frac{\xi y\left(s\right)}{\nu +y\left(s\right)}}-{\mu }_W\right]ds\right)\ge 0,$$

(3.5)

and $W=0$ if, and only if, $W_0=0$. Concerning the system for x and y, we first discard the trivial case when no adults or nymphs are initially present, that is $x_0=y_0=0$. Then, the only solution of system (3.4) is $(0,0,W_0{e}^{-{\mu }_{W}t})$, by (3.5), and the theorem holds trivially.

Let us then assume that not both $x_0=0$ and $y_0=0$ are true. We will prove that, in this case, $x\left(t\right),y\left(t\right)>0$ for all $t>0$, so that a fortioriy is bounded away from $-\nu$. Note that if either $x_0=0$ or $y_0=0$ (but not both), then the derivative of that which is initially zero is positive and, therefore, for all $t>0$small enough, by continuity, both $x\left(t\right)$ and $y\left(t\right)$ will be positive. Hence, we may assume, without loss of generality, that $x_0>0$ and $y_0>0$. Using an integrating factor to solve explicitly for y as a function of x in system (3.4),

$$y\left(t\right)=y_0{e}^{-{\int }_0^t\left({\mu }_N+\gamma +\xi W\left(s\right)\right)ds}+\rho {\int }_0^t{e}^{-{\int }_{\tau }^t\left({\mu }_N+\gamma +\xi W\left(s\right)\right)ds}x\left(\tau \right)d\tau .$$

(3.6)

Suppose now that x and y do not remain positive for all time and let $t_0$ be the first positive time for which either x or y vanish. Equation (3.6) proves that it is not possible that $y\left(t_0\right)=0$ because then it would necessarily exist $\overline{t}\in (0,t_0)$ for which $x\left(\overline{t}\right)=0$, contradicting the definition of $t_0$. Therefore, $x\left(t_0\right)=0$, $y\left(t_0\right)>0$, and $\overline{y}=min\{y\left(t\right):t\in [0,t_0]\}>0$ by continuity and definition of $t_0$. Because x is continuous and $x\left(t_0\right)=0$, there exists $h>0$ such that $x\left(t\right)<min\left\{{\displaystyle \frac{\gamma \overline{y}}{{\mu }_A+\frac{C_A}{B}}},1\right\}$ for all $t\in [t_0-h,t_0]$. But then Equation (3.4) implies that $x^{\prime }\left(t\right)>\gamma \overline{y}-\left({\mu }_A+{\displaystyle \frac{C_A}{B}}\right)x\left(t\right)>0$ for all $t\in [t_0-h,t_0]$ and, because $x\left(t\right)>0$ for all $t\in (t_0-h,t_0)$, it is impossible for $x\left(t_0\right)$ to vanish, contradicting our assumption. Consequently, x and y remain positive por all time.

□

3.1.2. Steady States of the System

The steady states of the system (3.4) are the solutions $(\overline{x},\overline{y},\overline{W})\in {\mathbb{R}}_{+}^3$ of the following quadratic algebraic system:

$$\begin{array}{cc}\hfill \gamma y-\left({\mu }_A+{\displaystyle \frac{C_A}{B}}x\right)x& =0,\end{array}$$

(3.7)

$$\begin{array}{cc}\hfill \rho x-\left({\mu }_N+\gamma +\xi W\right)y& =0,\end{array}$$

(3.8)

$$\begin{array}{cc}\hfill \xi \left({\displaystyle \frac{y}{\nu +y}}\right)W-{\mu }_{W}W& =0.\end{array}$$

(3.9)

It follows from (3.8) that

$$\overline{x}=0\iff \overline{y}=0.$$

(3.10)

We see from (3.9) that

$$\overline{W}=0or\overline{y}={\displaystyle \frac{{\mu }_W\nu }{\xi -{\mu }_W}}.$$

(3.11)

Next, Equations (3.7) and (3.8) imply that

$$\overline{y}=\left({\displaystyle \frac{{\mu }_A}{\gamma }}+{\displaystyle \frac{C_A}{\gamma B}}\overline{x}\right)\overline{x}and\overline{y}=\left({\displaystyle \frac{\rho }{\gamma +{\mu }_N+\xi \overline{W}}}\right)\overline{x},$$

(3.12)

whereby

$$\overline{x}=0or\overline{x}={\displaystyle \frac{\gamma B}{C_A}}\left({\displaystyle \frac{\rho }{\gamma +{\mu }_N+\xi \overline{W}}}-{\displaystyle \frac{{\mu }_A}{\gamma }}\right).$$

(3.13)

We can establish the following results concerning the equilibrium points of the system (3.4).

Theorem 3. 

There are two equilibria without T. radiata: the extinction equilibrium

$$E_0=(0,0,0),$$

(3.14)

and a unique “no -T. radiata" equilibrium

$$E_W=(\widehat{x},\widehat{y},0),\widehat{x}>0,\widehat{y}>0.$$

(3.15)

Proof. 

It follows from (3.9) and (3.10) that the only steady state with $\overline{x}=0$ or $\overline{y}=0$ is $E_0$. Also, (3.12) and (3.13) imply that a steady state of the form (3.15) must have

$$\begin{array}{|c|}\hline \widehat{x}={\displaystyle \frac{\gamma B}{C_A}}\left({\displaystyle \frac{\rho }{\gamma +{\mu }_N}}-{\displaystyle \frac{{\mu }_A}{\gamma }}\right)and\widehat{y}={\displaystyle \frac{\gamma \rho B}{C_A(\gamma +{\mu }_N)}}\left({\displaystyle \frac{\rho }{\gamma +{\mu }_N}}-{\displaystyle \frac{{\mu }_A}{\gamma }}\right)\\ \hline\end{array},$$

(3.16)

where the basic reproduction number for D. citri must be larger than 1:

$${\mathcal{R}}_0^{\mathit{D}}:={\displaystyle \frac{\rho \gamma }{{\mu }_A(\gamma +{\mu }_N)}}>1.$$

(3.17)

□

Note that this basic reproduction number is the product of the expected number of nymphs that a D. citri adult will produce, ${\displaystyle \frac{\rho }{\gamma +{\mu }_N}}$ and the expected number of adults that those nymphs will develop into, ${\displaystyle \frac{\gamma }{{\mu }_A}}$.

Concerning coexistence steady states, we see from (3.7) that, for any steady state $(\overline{x},\overline{y},\overline{W})$ with $\overline{x}\ne 0$ and $\overline{y}\ne 0$,

$$\overline{y}={\displaystyle \frac{1}{\gamma }}\left({\mu }_A+{\displaystyle \frac{C_A}{B}}\overline{x}\right)\overline{x},$$

which is equivalent to

$${\displaystyle \frac{\overline{x}}{\overline{y}}}={\displaystyle \frac{\gamma }{{\mu }_A+\frac{C_A}{B}\overline{x}}},$$

(3.18)

and in turn, equivalent to

$$\rho {\displaystyle \frac{\overline{x}}{\overline{y}}}-(\gamma +{\mu }_N)={\displaystyle \frac{\rho \gamma }{{\mu }_A+\frac{C_A}{B}\overline{x}}}-(\gamma +{\mu }_N)=(\gamma +{\mu }_N)\left({\displaystyle \frac{\rho \gamma }{{\mu }_A(\gamma +{\mu }_N)+\frac{C_A(\gamma +{\mu }_N)}{B}\overline{x}}}-1\right).$$

(3.19)

Finally, we can describe the coexistence equilibrium that requires the T. radiata basic reproduction number to be larger than 1:

$${\mathcal{R}}_0^{\mathit{T}}:={\displaystyle \frac{\xi }{{\mu }_W}}>1.$$

(3.20)

.

Theorem 4. 

Let ${\mathcal{R}}_0^{\mathit{T}}>1$, let $\overline{y}={\displaystyle \frac{{\mu }_W\nu }{\xi -{\mu }_W}}$, and $\overline{x}={\displaystyle \frac{\sqrt{{\left({\mu }_{A}B\right)}^2+4B{C}_A\gamma \overline{y}}-{\mu }_{A}B}{2C_A}}$. If ${\mathcal{R}}_0^{\mathit{D}}>1+{\displaystyle \frac{C_A\overline{x}}{B{\mu }_A}}$, there exists a unique coexistence equilibrium

$$E^{*}=(x^{*},y^{*},W^{*}),with{x}^{*},y^{*},W^{*}>0.$$

(3.21)

Proof. 

It follows respectively from (3.7), (3.11), and (3.8) that

$$\begin{array}{|c|}\hline x^{*}={\displaystyle \frac{\sqrt{{\left({\mu }_{A}B\right)}^2+4B{C}_A\gamma y^{*}}-{\mu }_{A}B}{2C_A}},y^{*}={\displaystyle \frac{{\mu }_W\nu }{\xi -{\mu }_W}},W^{*}={\displaystyle \frac{1}{\eta \xi }}\left({\displaystyle \frac{\rho x^{*}}{y^{*}}}-\gamma -{\mu }_N\right)\\ \hline\end{array},$$

(3.22)

where $x^{*}$ in (3.22) is positive if $y^{*}$ is, while $y^{*}$ in (3.22) is positive if, and only if, ${\mathcal{R}}_0^{\mathit{T}}>1$. Finally, $W^{*}$ in (3.22) is positive if, and only if, $\rho \gamma >{\mu }_A(\gamma +{\mu }_N)+{\displaystyle \frac{C_A(\gamma +{\mu }_N)}{B}}{x}^{*}$ because of (3.19), or equivalently, ${\displaystyle \frac{\rho \gamma }{{\mu }_A(\gamma +{\mu }_N)}}={\mathcal{R}}_0^{\mathit{D}}>1+{\displaystyle \frac{C_A{x}^{*}}{B{\mu }_A}}$.

□

3.1.3. Local Stability of Equilibrium Points

We shall now study the stability properties of these steady states by linearizing system (3.4) about each steady state.

The Jacobian matrix of the flow of system (3.4) is

$$J(x,y,W)=\left(\begin{array}{ccc}-{\mu }_A-2{\displaystyle \frac{C_a}{B}}x& \gamma & 0\\ \rho & -{\mu }_N-\gamma -\xi W& -\xi y\\ 0& {\displaystyle \frac{\xi \nu W}{{(\nu +y)}^2}}& {\displaystyle \frac{\xi y}{\nu +y}}-{\mu }_W\end{array}\right).$$

(3.23)

We begin with the extinction steady state, $E_0$.

Theorem 5. 

The extinction equilibrium is locally asymptotically stable if ${\mathcal{R}}_0^{\mathit{D}}<1$ and unstable if ${\mathcal{R}}_0^{\mathit{D}}>1$.

Proof. 

We see from (3.23) that

$$J\left(E_0\right)=\left(\begin{array}{ccc}-{\mu }_A& \gamma & 0\\ \rho & -{\mu }_N-\gamma & 0\\ 0& 0& -{\mu }_W\end{array}\right)$$

(3.24)

has eigenvalues $-{\mu }_W$ and the two roots of ${\lambda }^2+({\mu }_A+{\mu }_N+\gamma )\lambda +{\mu }_A({\mu }_N+\gamma )-\rho \gamma =0$, which has two roots with negative real part if, and only if, its coefficients are positive. This is equivalent to ${\mu }_A({\mu }_N+\gamma )-\rho \gamma >0$, that is

$${\mathcal{R}}_0^{\mathit{D}}={\displaystyle \frac{\rho \gamma }{{\mu }_A(\gamma +{\mu }_N)}}<1.$$

Similarly, when ${\mathcal{R}}_0^{\mathit{D}}>1$, the characteristic equation has two real roots, one negative and the other positive, making the equilibrium unstable and finishing our proof.

□

Next we turn to the stability of the “no-T. radiata" steady state.

Theorem 6. 

Let ${\mathcal{R}}_0^{\mathit{D}}>1$ so that $E_W$ exists. Then, $E_{W}isLASif{\mathcal{R}}_0^{\mathit{T}}<1+{\displaystyle \frac{\nu }{\widehat{y}}}$ and it is unstable if ${\mathcal{R}}_0^{\mathit{T}}>1+{\displaystyle \frac{\nu }{\widehat{y}}}$.

Proof. 

The Jacobian matrix (3.23) at the equilibrium $E_W$ is

$$J\left(E_W\right)=\left(\begin{array}{ccc}{\mu }_A-2{\displaystyle \frac{\gamma \rho }{\gamma +{\mu }_N}}& \gamma & 0\\ \rho & -{\mu }_N-\gamma & -\xi \widehat{y}\\ 0& 0& {\displaystyle \frac{\xi \widehat{y}}{\nu +\widehat{y}}}-{\mu }_W\end{array}\right),$$

(3.25)

with eigenvalues

$${\lambda }_1={\displaystyle \frac{\xi \widehat{y}}{\nu +\widehat{y}}}-{\mu }_W$$

(3.26)

and the two roots of

$${\lambda }^2+\left(\gamma +{\mu }_N+{\displaystyle \frac{2\gamma \rho }{\gamma +{\mu }_N}}-{\mu }_A\right)\lambda +\gamma \rho -{\mu }_A(\gamma +{\mu }_N)=0.$$

(3.27)

The discriminant of this quadratic is positive:

$${\left(\gamma +{\mu }_N+\frac{2\gamma \rho }{\gamma +{\mu }_N}-{\mu }_A\right)}^2-4\left(\gamma \rho -{\mu }_A(\gamma +{\mu }_N)\right)={\left({\displaystyle \frac{2\gamma \rho }{\gamma +{\mu }_N}}-{\mu }_A\right)}^2+{(\gamma +{\mu }_N)}^2+2{\mu }_A(\gamma +{\mu }_N),$$

and hence both its roots are real. Because ${\mathcal{R}}_0^{\mathit{D}}>1$, all the coefficients of (3.27) are positive, making both its roots negative. Finally, it follows from (3.26) that a necessary and sufficient additional condition for $E_W$ to be locally asymptotically stable is ${\displaystyle \frac{\xi \widehat{y}}{\nu +\widehat{y}}}<{\mu }_W$, which is automatic if ${\mathcal{R}}_0^{\mathit{T}}\le 1$. In any case, this condition is equivalent to ${\displaystyle \frac{\xi }{{\mu }_W}}<1+{\displaystyle \frac{\nu }{\widehat{y}}}$ and the theorem is proved.

□

We finally turn to the stability of $E^{*}$.

Theorem 7. 

Let ${\mathcal{R}}_0^{\mathit{T}}>1$ and ${\mathcal{R}}_0^{\mathit{D}}>1+{\displaystyle \frac{C_A{x}^{*}}{B{\mu }_A}}$. Then, the coexistence equilibrium $E^{*}$ is locally asymptotically stable.

Proof. 

The Jacobian matrix (3.23) at the coexistence equilibrium is

$$J\left(E^{*}\right)=\left(\begin{array}{ccc}-d_1& \gamma & 0\\ \rho & -\rho {\displaystyle \frac{x^{*}}{y^{*}}}& -\xi y^{*}\\ 0& {\displaystyle \frac{\xi \nu W^{*}}{{(\nu +y)}^2}}& 0\end{array}\right),$$

(3.28)

where $d_1=\sqrt{{\mu }_A^2+{\displaystyle \frac{4\gamma C_A{y}^{*}}{B}}}$. The characteristic polynomial is

$${\lambda }^3+p_2{\lambda }^2+p_1\lambda +p_0,$$

where

$$p_2=d_1+\rho {\displaystyle \frac{x^{*}}{y^{*}}}>0,p_1=\rho {\displaystyle \frac{x^{*}}{y^{*}}}{d}_1+{\displaystyle \frac{\xi {\mu }_W\nu W^{*}}{\nu +y^{*}}}-\rho \gamma ,p_0=d_1{\displaystyle \frac{\xi {\mu }_W\nu W^{*}}{\nu +y^{*}}}>0.$$

(3.29)

In order to prove that all eigenvalues have negative real part, we shall use the Routh-Hurwitz criterion: $p_2>0$, and $0<p_0<p_1{p}_2$. We just need to prove that $p_0<p_1{p}_2$. We see that

$$p_0<p_1{p}_2\iff 0<d_1\left(\rho {\displaystyle \frac{x^{*}}{y^{*}}}{d}_1-\rho \gamma \right)+\rho {\displaystyle \frac{x^{*}}{y^{*}}}\left(\rho {\displaystyle \frac{x^{*}}{y^{*}}}{d}_1+{\displaystyle \frac{\xi {\mu }_W\nu W^{*}}{\nu +y^{*}}}-\rho \gamma \right),$$

(3.30)

and it suffices to show that ${\displaystyle \frac{x^{*}}{y^{*}}{d}_1>\gamma }$. Using (3.18), it follows that

$$\begin{array}{ccc}\hfill {\displaystyle \frac{x^{*}}{y^{*}}}{d}_1>\gamma & \iff & d_1^2={\mu }_A^2+{\displaystyle \frac{4\gamma C_A{y}^{*}}{B}}>{\left({\displaystyle \frac{\gamma y^{*}}{x^{*}}}\right)}^2={\left({\mu }_A+{\displaystyle \frac{C_A}{B}}{x}^{*}\right)}^2\hfill \\ & \iff & 4\gamma y^{*}>2{\mu }_A{x}^{*}+{\displaystyle \frac{C_A}{B}}{\left(x^{*}\right)}^2\iff 4{\displaystyle \frac{\gamma y^{*}}{x^{*}}}>2{\mu }_A+{\displaystyle \frac{C_A}{B}}{x}^{*}\hfill \\ & \iff & 4\left({\mu }_A+{\displaystyle \frac{C_A}{B}}{x}^{*}\right)>2{\mu }_A+{\displaystyle \frac{C_A}{B}}{x}^{*},\hfill \end{array}$$

thus finishing the proof. □

Remark. 

It follows from Theorems 6 and 7 that, if both conditions $1<{\mathcal{R}}_0^{\mathit{T}}<1+{\displaystyle \frac{\nu }{\widehat{y}}}$ and ${\mathcal{R}}_0^{\mathit{D}}>1+{\displaystyle \frac{C_A{x}^{*}}{B{\mu }_A}}$ hold, then both equilibrium points $E_W$ and $E^{*}$ would be locally asymptotically stable. However, those conditions cannot hold simultaneously. In fact, we have

Theorem 8. 

Let ${\mathcal{R}}_0^{\mathit{T}}>1$. Then, ${\mathcal{R}}_0^{\mathit{T}}<1+{\displaystyle \frac{\nu }{\widehat{y}}}\iff {\mathcal{R}}_0^{\mathit{D}}<1+{\displaystyle \frac{C_A{x}^{*}}{B{\mu }_A}}$.

Proof. 

It follows from (3.15), (3.20), and (3.22) that

$${\mathcal{R}}_0^{\mathit{T}}<1+{\displaystyle \frac{\nu }{\widehat{y}}}\iff y^{*}={\displaystyle \frac{\nu }{{\mathcal{R}}_0^{\mathit{T}}-1}}>\widehat{y}={\displaystyle \frac{\rho B{\mu }_A}{C_A(\gamma +{\mu }_N)}}\left({\mathcal{R}}_0^{\mathit{D}}-1\right).$$

(3.31)

Also, it follows from (3.17) and (3.22) that

$$\begin{array}{ccc}\hfill {\mathcal{R}}_0^{\mathit{D}}>1+{\displaystyle \frac{C_A{x}^{*}}{B{\mu }_A}}& \iff & {\left({\displaystyle \frac{2\gamma \rho }{{\mu }_A(\gamma +{\mu }_N)-1}}\right)}^2>1+{\displaystyle \frac{4\gamma \nu C_A}{B{\mu }_A^2({\mathcal{R}}_0^{\mathit{T}}-1)}}\hfill \\ & \iff & {\displaystyle \frac{{\rho }^2\gamma }{{(\gamma +{\mu }_N)}^2}}={\displaystyle \frac{\rho {\mu }_A{\mathcal{R}}_0^{\mathit{D}}}{\gamma +{\mu }_N}}>{\displaystyle \frac{\rho {\mu }_A}{\gamma +{\mu }_N}}+{\displaystyle \frac{\nu C_A}{B({\mathcal{R}}_0^{\mathit{T}}-1)}}\hfill \\ & \iff & {\displaystyle \frac{\rho {\mu }_A}{\gamma +{\mu }_N}}\left({\mathcal{R}}_0^{\mathit{D}}-1\right)>{\displaystyle \frac{\nu C_A}{B({\mathcal{R}}_0^{\mathit{T}}-1)}}.\hfill \end{array}$$

(3.32)

The theorem now follows from (3.31) and (3.32).

□

3.1.4. Global Stability of the Trivial Equilibrium

We can establish global stability conditions for the trivial equilibrium (3.4).

Theorem 9. 

If ${\mathcal{R}}_0^{\mathit{D}}<1$ and ${\mathcal{R}}_0^{\mathit{T}}<1$, the extinction equilibrium, $E_0=(0,0,0)$ is globally asymptotically stable in the positive octant, ${\mathbb{R}}_{+}^3$.

Proof. 

Consider the linear Lyapunov function

$$V(x,y,W)=a_{1}x+a_{2}y+a_{3}W,a_i>0,i=1,2,3.$$

Each term of $V(x,y,W)$ is nonnegative for $(x,y,W)\in {\mathbb{R}}_{+}^3\setminus \left\{E_0\right\}$ when $a_1,a_2,a_3\ge 0$. Differentiating V along the trajectories of system (3.4) yields

$$\dot{V}=a_1\left[\gamma y-{\mu }_{A}x-{\displaystyle \frac{C_A}{B}}{x}^2\right]+a_2\left[\rho x-({\mu }_N+\gamma +\xi W)y\right]+a_3\left[{\displaystyle \frac{\xi y}{\nu +y}}W-{\mu }_{W}W\right].$$

(3.33)

Regrouping the terms in (3.33),

$$\begin{array}{c}\hfill \dot{V}=\underset{\mathrm{coef}.\mathrm{of}x}{\underbrace{\left(-a_1{\mu }_A+a_2\rho \right)}}x+\underset{\mathrm{coef}.\mathrm{of}y}{\underbrace{\left(a_1\gamma -a_2({\mu }_N+\gamma +\xi W)\right)}}y-a_1{\displaystyle \frac{C_A}{B}}{x}^2\\ \hfill +a_3\left[{\displaystyle \frac{\xi y}{\nu +y}}-{\mu }_W\right]W.\end{array}$$

(3.34)

Let

$$a_2={\displaystyle \frac{a_1{\mu }_A}{\rho }}>0,$$

so that the linear x-term in (3.34) vanishes. Thus,

$$\dot{V}=a_1\gamma \left[1-{\displaystyle \frac{{\mu }_A({\mu }_N+\gamma +\xi W)}{\rho \gamma }}\right]y-a_1{\displaystyle \frac{C_A}{B}}{x}^2+a_3\left[{\displaystyle \frac{\xi y}{\nu +y}}-{\mu }_W\right]W.$$

(3.35)

Define

$$k_1:=a_1\gamma \left(\frac{{\mu }_A({\mu }_N+\gamma )}{\rho \gamma }-1\right),k_2:=a_1\frac{C_A}{B},k_3:=a_3({\mu }_W-\xi ).$$

Note that:

• Since $R_0^{\mathit{D}}<1$ is equivalent to ${\displaystyle \frac{{\mu }_A({\mu }_N+\gamma )}{\rho \gamma }}>1\Rightarrow a_1\gamma \left(\frac{{\mu }_A({\mu }_N+\gamma )}{\rho \gamma }-1\right)=k_1>0$
• $a_1,C_A,B>0\Rightarrow a_1\frac{C_A}{B}=k_2>0$
• Also, because $R_0^{\mathit{T}}={\displaystyle \frac{\xi }{{\mu }_W}}<1$ we have ${\mu }_W-\xi >0\Rightarrow a_3({\mu }_W-\xi )=k_3>0$.
  Use monotonicity bounds:
  ${\mu }_N+\gamma +\xi W\ge {\mu }_N+\gamma \Rightarrow 1-{\displaystyle \frac{{\mu }_A({\mu }_N+\gamma +\xi W)}{\rho \gamma }}\le 1-{\displaystyle \frac{{\mu }_A({\mu }_N+\gamma )}{\rho \gamma }}=-{\displaystyle \frac{k_1}{a_1\gamma }}$
  Also
  ${\displaystyle \frac{y}{\nu +y}}<1\Rightarrow \xi {\displaystyle \frac{y}{\nu +y}}-{\mu }_W\le \xi -{\mu }_W=-({\mu }_W-\xi )$

Therefore,

$$\dot{V}\le -k_{1}y-k_2{x}^2-k_{3}W\le 0.$$

The largest invariant set where $\dot{V}=0$ is exactly $E_0=(0,0,0)$. By LaSalle’s invariance principle, every trajectory with nonnegative initial conditions converges to $E_0$. Hence $E_0$ is globally asymptotically stable in ${\mathbb{R}}_{+}^3$.

□

3.2. Simulations

3.2.1. Reduced System

For our first simulations, corresponding to the reduced system (3.4), we use the parameter values displayed in Table 2, estimated from the literature or within a biologically relevant range —with the exception of some tree-related parameters that were set to zero to ensure a constant total number of shoots as required in (3.4). We then modify the values of some of the parameters so that they fall under the ranges of Theorems 5–7.

Table 2. Parameter values for simulating the reduced D. citri–T. radiata model

Parameter	Value	Parameter	Value
$\rho$	$6.52\times {10}^{-2}$	$C_A$	$2.07\times {10}^{-5}$
${\mu }_A$	1/32	${\mu }_N$	1/6
$\gamma$	1/6.11	${\mu }_W$	0.04
$\xi$	2	$\nu$	200

Table 2. Parameter values for simulating the reduced D. citri–T. radiata model

Parameter	Value	Parameter	Value
$\rho$	$6.52\times {10}^{-2}$	$C_A$	$2.07\times {10}^{-5}$
${\mu }_A$	1/32	${\mu }_N$	1/6
$\gamma$	1/6.11	${\mu }_W$	0.04
$\xi$	2	$\nu$	200

These result in the following basic reproduction numbers:

Table 3. Initial values for the reduced D. citri–T. radiata model.

Variable	Initial Value
$N\left(0\right)$	30
$A\left(0\right)$	0
$T_1\left(0\right)$	200
$T_2\left(0\right)$	100
$N_i\left(0\right)$	0.2
$A_i\left(0\right)$	150
${T_1}_i\left(0\right)$, ${T_2}_i\left(0\right)$	0
$W\left(0\right)$	1,000

Table 3. Initial values for the reduced D. citri–T. radiata model.

Variable	Initial Value
$N\left(0\right)$	30
$A\left(0\right)$	0
$T_1\left(0\right)$	200
$T_2\left(0\right)$	100
$N_i\left(0\right)$	0.2
$A_i\left(0\right)$	150
${T_1}_i\left(0\right)$, ${T_2}_i\left(0\right)$	0
$W\left(0\right)$	1,000

${\mathcal{R}}_0^{\mathit{T}}>1and{\mathcal{R}}_0^{\mathit{D}}>1$,

which correspond to Theorem (7). We ran a 60-year simulation using the following initial conditions:

The results of are displayed in Figure 1 and Figure 2, clearly corroborating Theorem 7.

Next, we modified the unit birth rate of T. radiata to $\xi =0.041$ so its basic reproduction number is ${\mathcal{R}}_0^{\mathit{T}}=1.025$ corresponding to Theorem (6), as corroborated in Figure 3 and Figure 4.

Finally, we restore $\xi =0.1$ and slightly modify the value of the per capita “oviposition” rate for D. citri to $\rho =0.0625$ to lower its basic reproduction number to ${\mathcal{R}}_0^{\mathit{D}}=0.991$ that corresponds to Theorem (5), as corroborated below in Figure 5, and Figure 6.

3.2.2. Full System

Because the model (2.1)–(2.2) was designed for a real agricultural context, we considered several different scenarios to simulate potential management practices. For all of them we assume:

• No transition between low-productivity and high-productivity trees (or vice versa)
• No tree eradication
• Effectiveness of agrochemical decays linearly to zero 30 days after application

For our second set of simulations, corresponding to the system (Section 2), we use the parameter values displayed in Table 4, estimated from the extant literature or chosen within a biologically relevant range. We ran a 1-year simulation using the initial conditions described in Table 5.

Table 4. Parameter values for simulating the extended D. citri-T. radiata model.

Parameter	Value	Source
$\rho$	0.995	Estimated value from [11]
$m_1$	42	Estimated value from [12]
${m_1}_i$	32	Assumed from [12]
$m_2$	21	Assumed from [12]
${m_2}_i$	21/2	Assumed from [12]
${\mu }_N$	0.02	Estimated value from [12]
${\mu }_A$	1/(48.25)	Estimated value from [13]
${\mu }_W$	0.04	Estimated value from [14]
${\mu }_{T_2}$	1/(10*365)	Estimated value from [15]
${\mu }_{{T_2}_i}$	1/(5*365)	Estimated value from [16]
$\gamma$	1/9	Estimated value from [13]
$\xi$	25	Estimated value from [17]
$C_A$	$2.07\times {10}^{-5}$	Assumed
${\theta }_T$	$7.00\times {10}^{-3}$	[18]
$\tau$	0	Assumed
$\lambda$, ${\lambda }_i$	0	Assumed
${\alpha }_1$, ${\alpha }_2$	0	Assumed
$\nu$	200	Assumed
$\eta$	${10}^{-6}$	Assumed

Table 4. Parameter values for simulating the extended D. citri-T. radiata model.

Parameter	Value	Source
$\rho$	0.995	Estimated value from [11]
$m_1$	42	Estimated value from [12]
${m_1}_i$	32	Assumed from [12]
$m_2$	21	Assumed from [12]
${m_2}_i$	21/2	Assumed from [12]
${\mu }_N$	0.02	Estimated value from [12]
${\mu }_A$	1/(48.25)	Estimated value from [13]
${\mu }_W$	0.04	Estimated value from [14]
${\mu }_{T_2}$	1/(10*365)	Estimated value from [15]
${\mu }_{{T_2}_i}$	1/(5*365)	Estimated value from [16]
$\gamma$	1/9	Estimated value from [13]
$\xi$	25	Estimated value from [17]
$C_A$	$2.07\times {10}^{-5}$	Assumed
${\theta }_T$	$7.00\times {10}^{-3}$	[18]
$\tau$	0	Assumed
$\lambda$, ${\lambda }_i$	0	Assumed
${\alpha }_1$, ${\alpha }_2$	0	Assumed
$\nu$	200	Assumed
$\eta$	${10}^{-6}$	Assumed

Table 5. Initial values for the extended D. citri–T. radiata model.

Variable	Initial Value
$N\left(0\right)$	192
$A\left(0\right)$	214
$T_1\left(0\right)$	100
$T_2\left(0\right)$	43
$N_i\left(0\right)$	96
$A_i\left(0\right)$	107
${T_1}_i\left(0\right)$	100
${T_2}_i\left(0\right)$	42
$W\left(0\right)$	0

Table 5. Initial values for the extended D. citri–T. radiata model.

Variable	Initial Value
$N\left(0\right)$	192
$A\left(0\right)$	214
$T_1\left(0\right)$	100
$T_2\left(0\right)$	43
$N_i\left(0\right)$	96
$A_i\left(0\right)$	107
${T_1}_i\left(0\right)$	100
${T_2}_i\left(0\right)$	42
$W\left(0\right)$	0

A total of 285 Tahiti acid lime trees (Citrus × latifolia Tanaka ex Q. Jiménez) are included, corresponding to standard planting density per hectare for Colombian cultivation [15]. The interventions we consider are one or both of releasing T. radiata and applying agrochemicals. Therefore, our simulation scenarios are the following four:

• Scenario 1, No control: (i.e., no management interventions)

  $${\delta }_N\left(t\right)={\delta }_A\left(t\right)={\delta }_W\left(t\right)=\Lambda \left(t\right)=0,\forall t\ge 0.$$
• Scenario 2,T. radiatarelease only: introduction of 400 females at the initial time, taking into account as a reference the prescribed release rates (individuals/ha, between male and female) [15], as shown in Figure .
• Scenario 3, Agrochemical application only: (${\delta }_N\left(t\right)$, ${\delta }_A\left(t\right)$, ${\delta }_W\left(t\right)$) shown in Figure .
• Scenario 4, Combined control: pairing initial agrochemical application (see Figure ) with T. radiata release of 100 females at day 30 (see Figure ).

Figure 7. Functions modeling the interventions

Figure 7. Functions modeling the interventions

The results of these simulation scenarios are shown in Figure 8, Figure 9, and Figure 10. In each of these figures, the right panels depict a zoom-in of the corresponding left panels near the origin to help understand the initial impact of the interventions, while the 1-year simulations displayed on the left panels are intended to show the long-term effect of the interventions. In Figure 8, we show the result of the simulations for the populations of uninfected D. citri nymphs and adults. In Figure 9, we show the result of the simulations for the populations of infected D. citri nymphs and adults. In Figure 10, we show the result of the simulations for the populations of T. radiata and for those of uninfected and infected trees. In all scenarios, the population of HLB-infected trees increases while the populations of uninfected and total trees decrease, because the model does not include any source for new healthy trees or recovery from HLB infection.

Figure 8. non-infectious D. citri populations

Figure 8. non-infectious D. citri populations

Figure 9. infected D. citri populations

Figure 9. infected D. citri populations

Figure 10. T. radiata and tree populations

Figure 10. T. radiata and tree populations

Scenario 1 serves as a baseline that underscores the devastating effect that HLB may have on the plantation (second panel of Figure 10). We see that after one year, the initial 142 infected trees would increase their numbers to 245 (72.5% increase), while the initial 143 uninfected trees would decrease to just 26 (81.8% decrease).

We see, in Figure 8 for the uninfected D. citri and in Figure 9 for the infectious D. citri, that just the release of 400 adult wasps of T. radiata when the D. citri nymph population is still small (here, one-half the number of released T. radiata, Scenario 2) results in a very quick extinction of the latter, in approximately 2 days in both figures (second panels of Scenario 2), which then depletes the D. citri adult population in about 8-9 months (first panels of Scenario 2), as well as the T. radiata population that cannot productively lay their eggs once the nymphs are depleted (first panel of Scenario 2 in Figure 10). After one year, this intervention is very successful in reducing the plantation’s losses, by reducing the proportional increase of infected trees from 72.5% when no control measures are implemented to 57.7%, as well as diminishing the 81.8% decrease in the number of uninfected trees to 66.4%.

On the other hand, we also see, in Figure 8 and Figure 9, that just the application of an agrochemical of high initial effectiveness at that same time when the D. citri nymph population is still small (numbering 200, Scenario 3) almost leads to the extinction of the latter in the 30 days that the agrochemical is still effective (second panels of Scenario 2 in both figures) and, consequently, we also see the adult D. citri population being almost led to extinction. However, after those 30 days, the D. citri nymph population takes off to exceedingly large numbers that peak after about 5 months (first panels of Scenario 3 in both figures) by which time the population of uninfected trees has fallen to approximately half of its original number (second panel of Scenario 3 in Figure 10), causing the decline of both nymph and adult populations of D. citri to decrease.

In Scenario 4, we considered the sequential use of both control measures rather than their combined used simultaneously, with the release of T. radiata 30 days after the application of agrochemicals, to avoid considering the negative impact that the use of the latter has on the population of the former (because we assume that the effectiveness of the agrochemical is negligible after 30 days). In this scenario the realease of T. radiata involves just 100 adult wasps, instead of the 400 used in Scenario 2. Yet, their timely release when the populations of D. citri had been decimated by the agrochemical, is enough to “finish the job". In fact, at day 30 in Scenario 4, the ratio of T. radiata to D. citri nymphs is considerably larger than the 2:1 ratio at the start of Scenario 2, thus making the wasps much more effective in extinguishing the D. citri population (and henceforth self-eliminating too!).

Comparing the second panels of Scenarios 3 and 4 in Figure 8 and Figure 9, we see that the plots are identical during the first 30 days because the two scenarios are actually the very same during those 30 days. However, the delay in the release of the wasps when compared to Scenario 2, as expected, still significantly improves the results after one year in Scenario 3, both in terms of the numbers of healthy and infected trees and in the D. citri populations but cannot match the outcomes of Scenario 2, either for the trees or for the psyllyd populations.

The seeming resurgence of D. citri adults after approximately 2-3 months (first panel of Figure 9, Scenarios 1 and 3) is a numerical artifact due to the positivity of solutions of our ODE system and the increased number of infected trees upon whose shoots the adult psyllid feeds and lays eggs.

4. Discussion

We proposed a 9-dimensional ODE model to describe the dynamics of interactions and infection of D. citri, T. radiata, and lemon trees with the bacterium Huanlongbing. The model incorporates control interventions by agrochemicals (pesticides to contain the D. citri population) and by the parasitoid of D. citri, T. radiata.

We reduced the system to a 3-dimensional system that omits the dynamics of the tree population and combines uninfected and infected populations of D. citri. That is, the reduced model described the dynamics of interaction of D. citri with its parasitoid T. radiata.

For the reduced system we carried out a fairly complete mathematical analysis of well-posedness and long term dynamics. We proved that, for non-negative initial conditions, the intial value problem has a unique, non-negative, global-in-time solution that converges asymptotically to one of three equilibrium points of the system. The asymptotic behavior is determined by the two basic reproduction numbers of the interacting populations in threshold form.

We found explicit expressions for the only three biologically relevant equilibrium points, one being the extinction equilibrium, a second one, the “no-Tamarixia-equilibrium" (a boundary equilibrium, with positive D. citri population but no T. radiata), and a coexistence equilibrium.

We provided explicit formulas for ${\mathcal{R}}_0^{\mathit{D}}$, the basic reproduction number for the D. citri population, and for ${\mathcal{R}}_0^{\mathit{T}}$, the basic reproduction number for the T. radiata population. We proved that the local asymptotic stability of equilibria operates exactly as one would expect: when both basic reproduction numbers are smaller than 1, the extinction equilibrium is the only one in the non-negative octant, and it is locally asymptotically stable. As ${\mathcal{R}}_0^{\mathit{D}}$ increases from below to above 1, the “no-Tamarixia-equilibrium", $E_W$, bifurcates from the extinction equilibrium into the non-negative octant and it is locally asymptotically stable as long as ${\mathcal{R}}_0^{\mathit{T}}$ is smaller than an explicitly computable critical value that is larger than 1. As ${\mathcal{R}}_0^{\mathit{T}}$ increases from below that critical value to above it, $E_W$ becomes unstable and the coexistence equilibrium (which bifurcated from $E_W$ when ${\mathcal{R}}_0^{\mathit{T}}$ crossed from below to above 1) becomes locally asymptotically stable.

We also proved global stability of the trivial equilibrium via Lyapunov functions and LaSalle’s invariance principle, establishing that the threshold conditions for local asymptotic stability actually result in global stability within the non-negative octant of ${\mathbb{R}}^3$.

We performed simulations of the full 9-dimensional model in four scenarios, to get a sense of the impact of the two types of control measures. The first scenario is without control and serves as a baseline. The second one introduces T. radiata only while the third one uses only one application of a high-efficiency agrochemical pesticide. The last scenario uses the agrochemical first, followed by the release of a reduced number of T. radiata after 30 days, which is the estimated duration of the pesticide’s effect.

Based on the simulations, it seems that the preferred control method is T. radiata, but a further cost anaylisis needs to be performed. The reason that T. radiata is so effective hinges upon its mode of reproduction. The fact that they lay one egg into a D. citri nymph that will not evolve to adulthood, makes this wasp species highly effective in controlling –and even eliminating– the D. citri population.