Stability Analysis and Optimal Control Strategy for Prevention of Pine Wilt Disease

and Applied Analysis 3 System (3) has always a disease-free equilibrium E 0 = (Λ h /μ 1 , 0, 0) and basic reproduction number R 0 = Λ h Λ Vγ(αφ + βψθ)/μ 2 1 μ 2 2 . If R 0 > 1, system (3) has a unique endemic equilibrium E ∗ = (S ∗ h , I ∗ h , I ∗ V ) in the interior of Γ, where

Pine trees infected by PWD usually die within a few months.The first observable symptom is the lack of resin exudation from bark wounds.The foliage then becomes pale green, then yellow, and finally reddish brown as the tree succumbs to the disease.It is well known that there are three transmission pathways of PWD.One occurs during maturation feeding, when the nematode is transferred from insect vectors to healthy pine trees via the insect feeding on wounds [4].The second occurs during oviposition of the mature female on dead, dying, or recently cut pine trees via the oviposition wounds [6].The third occurs during mating, that is, as mature males search for females in bark wounds, such as the oviposition wounds.This is referred to as horizontal transmission [14].
Vector-borne diseases are infectious diseases caused by viruses, bacteria, protozoa, or rickettsia.They are primarily transmitted by disease transmitting biological agents (anthropoids), called vectors, who carry the disease without getting it themselves.For example, the most prevalent vector-borne disease is malaria, whose vectors are mosquitoes.Numerous cases of vector-borne disease are known for plants.Among them, there are some important wilting diseases of trees, such as PWD and the red ring disease of palms.In recent years, many mathematical models have been used to investigate the transmission dynamics of the PWD [15][16][17][18].Mathematical models have become important tools in analyzing the spread and control of infectious diseases.Thus, we present 2 Abstract and Applied Analysis a mathematical model to describe the vector-host interaction between the pine trees and nematode-carrying pine sawyers.The objective of this study is to investigate global dynamics and control the spread of the PWD based on the proposed mathematical model.
In this paper, we propose a mathematical model of PWD using ordinary differential equations, with mass action incidence, and we consider a controlled PWD transmission model to prevent the spread of disease using optimal control strategies.Generally, PWD control has focused on the elimination of pine sawyer larvae inhabiting pine wilt trees (eradication of infected host pine trees), either by winter fumigation or by controlling the adult sawyers using insecticide over the summer.We apply two control parameters: the injection of nematicide to the host pine tree and the aerial spraying of pesticide to kill adult pine sawyers.
For this, we first introduce control variables representing the optimal treatment for infected hosts and vectors.We formulate an optimal system for the mathematical model of PWD caused by adult pine sawyers carrying pinewood nematodes.
Furthermore, we show the existence of an optimal control for this problem and analyze the system using optimal control techniques.
This paper is organized as follows.In Section 2, we derive the ordinary differential equations that describe the interactions between the host pine tree and vector forest beetles populations.In Section 3, the control problem is formulated.In Section 4, we derive the necessary conditions for an optimal control and the corresponding states using Pontryagin's Maximum Principle.Section 5 presents the results of numerical simulations, and our conclusions are given in Section 6.

Model Formulation
In this paper, we consider a four-dimensional model, which consists of the populations of susceptible host pine trees  ℎ , infected host pine trees  ℎ , susceptible vector beetles that do not carry the PWN  V , and infected vector beetles that carry the PWN  V .The total population sizes at time  for the host pine trees and vector beetles are denoted by  ℎ and  V , respectively.In the host pine trees populations, Pine trees infected with the PWD never recovers from an infected states, and most infected pine trees die within a year of infection.Thus we do not consider the number of immune hosts  ℎ in the model.Thus  ℎ =  ℎ +  ℎ and  V =  V +  V .The model is given by the following system of differential equations: Susceptible host pine trees are recruited at a constant rate Λ ℎ .Parameter  represents the transmission rate per contact with infected vectors during maturation feeding.Parameter  is the average number of contacts per day with vector adult beetles during the maturation feeding period.We denote the incidence of new infections via this route by the mass actions term  ℎ  V .The transmission probability that infected beetles transmit a nematode by oviposition is denoted by , and the average number of contacts per day when adult beetles oviposit is denoted by .The parameter  is the probability that susceptible host pine trees cease oleoresin exudation without being infected by the nematode.In the model, the rate of transmission through oviposition is denoted by , and the incidence of new infections via this route is given by the mass action term  ℎ  V .A constant emergence rate of the vector pine sawyer beetle is denoted by Λ V , and parameter  denotes the rate at which adult beetles carry the PWN when the beetles escape from dead trees.The incidence terms for vector populations are given by  ℎ  V .Finally,  1 ,  2 denote the natural death rate of the host population and the natural death rate of the vector population, respectively.

A Model for Optimal Control of PWD
In the host pine tree population, the force of infections is reduced by 1 −  1 , where  1 ∈ [0, 1] measures the level of successful prevention efforts.It follows that control variable  1 represents the use of drugs or vaccine which are alternative preventive measures to minimize or eliminate PWN infection (e.g., tree-injection of a nematicide).In system (1), we modified the recruitment rate in susceptible vector populations by including density effects.Further, we replace the previous recruitment rates with Λ V → Λ V  V .In the vector adult beetles population, the control variable  2 ∈ [0, 1] represents the level of adulticide used for vector control such as aerial spraying of pesticide.It follows that the reproduction rate of the vector population is reduced by a factor of 1 −  2 .Further, we assumed that the mortality rate of the vector population increases at a rate proportional to  2 , where  0 > 0 is a rate constant.Considering these assumptions and extensions, it follows that the extended model is described by an initial value problem with a system of four differential equations: with initial conditions  ℎ (0) ≥ 0,  ℎ (0) ≥ 0,  V (0) ≥ 0, and  V (0) ≥ 0.

Analysis of Optimal Control
We define the objective functional 2 ) subject to the state system given by (5).The objective of our work is to minimize the infected host pine tree and the total number of vector population and the cost of implementing the control by using possible minimal control variables  1 ,  2 .For the optimal control problem of the given system, we consider the control variables  = ( 1 ,  2 ) ∈  relative to the state variables  ℎ ,  ℎ ,  V , and  V where both control variables are bounded and measured with In the objective functional the quantities  1 ,  2 represent the weight constants of the infected host population and total vector population, respectively.In the objective functional  1 and  2 are weight factors for reduction of vector-host contacts and vector control, respectively.The terms and ( 2 /2) 2 2 describe the costs associated with prevention of vector-host contacts and vector control, respectively.
The cost associated with the first control could come from the use of tree-injection of a nematicide.The cost associated with the second control could come from physical or chemical treatment of wilt pines to kill their larvae or eradication of vector beetles by aerial pesticide spraying.
Our aim is to find control functions subject to system (5), where the control set is defined as We note that the existence of an optimal control pair can be proved by using results from Fleming and Rishel [19].We know easily that the system of equations given by ( 5) is bounded from above by a linear system.According to the theorem in [19], the solution exists if the following hypotheses are met.
( 1 ) The set of controls and corresponding state variables is nonempty.( 2 ) The control set  is convex and closed.
( 3 ) Each right-hand side of the state system is bounded by a linear functional in the state and control and can be written as a linear function of  with coefficients depending on time and state.( 4 ) There exist constants  1 ,  2 > 0 and  > 1 such that the integrand (, , ) of the objective functional  is convex and satisfies In order to verify these conditions, we use a result by Lukes [20] to give the existence of solutions of the state system (5) with bounded coefficients, which gives ( 1 ).We note that the solutions are bounded.Our control set satisfies condition ( 2 ).Since our state system is bilinear in  1 ,  2 , the righthand side the state system (5) satisfies condition ( 3 ), using the boundedness of the solutions.Note that the integrand of our objective functional is convex.Moreover we have the last condition needed where  1 ,  2 > 0,  1 ,  2 ,  1 ,  2 > 0, and  > 1.Thus we obtain the following.
]} subject to the system (5) with initial conditions, then there exists an optimal control Lagrangian for a problem discusses how the techniques come and Hamiltonian helps in solving the adjoint variable, which we need to construct for the optimal control problem.In order to find an optimal solution pair, first we should construct the Lagrangian and Hamiltonian for the optimal control problem (5).In fact, the Lagrangian of the optimal problem is given by We seek for the minimal value of the Lagrangian.To do this, we define the Hamiltonian  for the control problem as taking  = ( ℎ ,  ℎ ,  V ,  V ),  = ( 1 ,  2 ), and  = ( 1 ,  2 ,  3 ,  4 ), to obtain In order to find the necessary conditions for this optimal control pair, we apply Pontryagin's Maximum Principle [21] as follows.
If ( * 1 ,  * 2 ) is an optimal solution of an optimal control problem, then there exists a nontrivial vector function () = ( 1 (),  2 (), . . .  ()) satisfying the following equalities.The state equation is the optimality condition and the adjoint equation Now we apply the necessary conditions to the Hamiltonian .
Theorem 2. Given optimal controls  * 1 ,  * 2 and solutions  * ℎ ,  * ℎ ,  * V , and  * V of the corresponding state system (5), there exist adjoint variables  1 ,  2 ,  3 , and  4 satisfying with transversality conditions Furthermore,  * 1 ,  * 2 are represented by Proof.To determine the adjoint equations and the transversality conditions we use the Hamilton .We differentiate the Hamiltonian  with respect to  ℎ ,  ℎ ,  V ,  V .Thus we obtain with transversality conditions Using the optimality conditions we have and the property of the control space  we can drive

Numerical Results
In this section, we show the numerical simulations of the impacts of the optimal control strategy in the proposed the PWD model.The optimality system is solved using the Runge-Kutta fourth order scheme.The optimal strategy is obtained by solving the state and adjoint systems and the transversality conditions.First we start to solve the state ( 16) using the Runge-Kutta fourth order forward in time with a guess for the controls over the simulated time.Then, using the current iteration of the state equations, the adjoint equations in system ( 16) are solved by a backward method with the transversality conditions (17).Then the controls are updated by using a convex combination of the previous controls and the value from the characterizations (18).This process is repeated and iterations stopped if the values of the unknowns at the previous iterations are very close to the ones at the present iterations [21].We have chosen the same set of weight factors  1 = 0.01,  2 = 0.001,  1 = 50, and  2 = 200 initial state variables  ℎ (0) = 300,  ℎ (0) = 50,  V (0) = 65, and  V (0) = 20 to illustrate the effect of different optimal control strategies on the spread of pine wilt disease in a population.
In this section, we investigate numerically the effect of the following optimal control strategies on the spread of pine wilt disease.Three different control strategies are explored.Then we look at the following three alternatives: (i) strategy 1: precaution effort by tree-injection of a nematicide and eradication of vector beetles (controls  1 and  2 ); (ii) strategy 2: physical or chemical treatment of wilt pines to kill their larvae or eradication of vector beetles by aerial pesticide spraying ( control  2 alone); (iii) strategy 3: precaution effort by tree-injection of a nematicide ( control  1 alone).
Baseline values of model parameters used for simulations are presented in Table 1.Owing to the absence of data, some of the other parameters associated with given model are assumed (see Table 1).

Strategy 1: Optimal Use of Tree-Injection of a Nematicide (𝑢 1 ̸
= 0) and Eradication of Vector Beetles ( 2 ̸ = 0).The computed results for optimal control under strategy 1 and 4. In Figures 2, 3, and 4, we can see that the control strategy results in a decrease in the number of infected hosts  ℎ , susceptible vectors  V , and infected vectors  V .With this strategy, an increase in the number of susceptible hosts  ℎ is observed in Figure 1.

Strategy 2:
Optimal Use of Eradication of Vector Beetles ( 2 ̸ = 0).The computed results for optimal control under strategy 2 ( 1 = 0,  2 ̸ = 0) are shown in Figures 5, 6, 7, and 8. Figures 6, 7, and 8 show a significant difference in the number of infected hosts, susceptible vectors, and infected vectors, respectively, compared to the situation where there is no control.Our numerical results suggest that tree-injection of a nematicide is not effective if we want to reduce the spread of PWD.

Strategy 3: Optimal Use of Tree-Injection of a Nematicide
( 1 ̸ = 0).The computed results for optimal control under strategy 3 ( 1 ̸ = 0,  2 = 0) are shown in Figures 9, 10, 11, and 12. Figures 10 and 12 show no significant difference in the number of infected hosts  ℎ and infected vectors  V , respectively, compared to the situation where there is no control.This shows that the eradication of vector beetles is more effective than injecting trees with nematicide.Figure 13 represents the optimal controls.

Effects of Weight Constants.
A sensitivity analysis is carried out by studying the adequacy of our simulations in relation to the weight constants.The role of weight constants is explored by assessing their quantitative impact in terms of the number of infected cases.Comparative results with the implementation of strategy 1 under different weight constants on the controls  1 ,  2 are shown in Figures 14 and 16 ( 1 = 50,  1 = 5,  1 = 2, fixing  2 = 200) and Figures 15 and 17 (fixing  1 = 50,  2 = 200,  2 = 100, and  2 = 50).In Figure 16, the value of the weight constant  1 is varied.As the value of weight constants increases, control functions decrease, and increasing the cost of successful prevention efforts, the population of infected host  ℎ is increased.In Figure 17, the value of the weight constant  2 is varied.As the value of weight constants increases, control functions decrease, and increasing the cost of vector control, the population of infected vector  V is increased.

Conclusions
In this paper, we proposed a host vector transmission model to study the dynamics of PWD.A mathematical analysis was carried out for a model in which PWD is caused by the PWN, a parasite hosted by the pine sawyer beetle.The global dynamics of the model showed the existence and stability of disease-free and endemic equilibria.We proved that if  0 ≤ 1, the disease-free equilibrium  0 is globally asymptotically stable in Γ, and thus the disease always dies out.However, if  0 > 1, the unique endemic equilibrium  * exists and is globally asymptotically stable, indicating that PWD persists at the endemic equilibrium if it is initially present.We have

𝜓
The number of contacts during the oviposition period 0.41 day −1 Assumed

𝜃
The probability in which susceptible host pine is not infectious by nematode and ceases oleoresin exudation naturally 0.0000301 day considered the optimal control strategy for preventing the spread of PWD.The conditions for optimal control of the disease were derived and analyzed, and the injection of nematicide to infected hosts was compared with the aerial spraying of insecticides to control infected vectors.From the control plots, we showed that the number of infected hosts and the total vector population decreased in the optimal system.Control programs that follow these strategies can effectively reduce the spread of PWD in forest ecosystems.Numerical simulations illustrated the effectiveness and efficiency of the proposed control problem more effective than  the tree-injection strategy for controlling the spread of PWD in the forest ecosystem.

A. The Disease-Free Equilibrium and Its Stability
In this section, we will investigate the disease-free equilibrium and its stability of the system (3).Theorem A.1.If  0 < 1, then the disease-free equilibrium  0 is locally asymptotically stable.
Proof.Linearize the system (3) around the disease-free equilibrium  0 .The matrix of the linearization at  0 is given by The characteristic equation of this matrix is given by det( − ( 0 )) = 0, where  is the 3 × 3 unit matrix.Expanding the determinant into a characteristic equation we obtain the following equation: So,  = − 1 and  2 +  1  +  2 = 0, where  Thus, according to the Routh-Hurwitz criteria [27], the quadratic equation has only roots with negative real parts.Hence, the disease-free equilibrium  0 of the system (3) is locally asymptotically stable.Now, we study the global behavior of the disease-free equilibrium for system (3).Theorem A.2.If  0 ≤ 1, then the disease-free equilibrium  0 is globally asymptotically stable in Γ.
Proof.We define the following Lyapunov function on Γ.where Clearly,  ≥ 0 along the solution of the system (3) and is zero if and only if both  ℎ and  V are zeros.The derivative of  along the solutions of (3) is We know that   We see that if  0 < 1, the derivative   = 0 if and only if  ℎ = 0, while in the case of  0 = 1 the derivative   = 0 if and only if  ℎ = 0 or  V = 0. Consequently, the largest compact invariant set in {( ℎ ,  ℎ ,  V ) ∈ Γ |   = 0}, when  0 ≤ 1, is the singleton { 0 }.Hence by Lasalle's Invariance Principle [28],  0 is globally asymptotically 10 stable in Γ.This completes the proof.

B. The Endemic Equilibrium and Its Stability
In this section, we will investigate the endemic equilibrium and its stability of the system (3).We first give the following results.Lemma B.1.Let  be a 3 × 3 real matrix.If tr(), det(), and det( [2] ) are all negative, then all eigenvalues of  have negative real part.
Proof.The proof of this lemma is given in [29].
Theorem B.2.If  0 ≥ 1, then the endemic equilibrium  * is locally asymptotically stable.Proof.In order to investigate the local stability of the endemic equilibrium, the additive compound matrices approach as in [30,31] is used.We will linearize system (3) about an endemic equilibrium  * and get the following Jacobian matrix: ) .
(B.1)  From the Jacobian matrix ( * ), the second additive compound matrix  [2] ( * ) is given by  [2]  And we know that system (3) is given since we know that Computing directly the determinant of  [2] ( * ) and from (B.4) and (B.6) we can get det ( [2] Hence by Lemma B.1, the endemic equilibrium  * of the model ( 3) is locally asymptotically stable in Γ.
Now we prove the global stability of the endemic equilibrium  * , when the reproduction number  0 is greater than the unity.For this, we will prove the following result.
Finally, we will investigate the global stability of the endemic equilibrium  * in the feasible region Γ.To do this we will use the results about the three dimensionless competitive systems that live in convex sets [33][34][35] and powerful theory of additive compound matrix to prove asymptotic orbital stability of periodic solutions [30].We easily can get that the system (3) is competitive in Γ, by looking at its Jacobian matrix and choosing the matrix  as with respect to the partial order defined by the orthant Theorem B.4.If  0 > 1, then the endemic equilibrium  * is globally asymptotically stable in int Γ.
Since system (3) is competitive and persistent, and  * is locally asymptotically stable if  0 > 1.Furthermore, in accordance with Theorem 2.1 and Theorem 2.2 [35], Theorem 600 would be established if we show that system (3) has the property of stability of periodic orbits.Hence we will prove the following theorem.
Theorem B.5.If  0 > 1, then system (3) has the property of stability of periodic orbits.
Proof.In accordance with the criterion given by Li and Muldowney [31], for the asymptotic orbital stability of a periodic orbit of a general autonomous system, it is sufficient to prove that the linear nonautonomous system   () =  ( [2] ( ()))  () (B.15) is asymptotically stable, where ( [2] ) is the second additive compound matrix of the Jacobian ().The Jacobian of system (3) is given by  () = ( − ( + )  V −  where ‖ ⋅ ‖ is the norm in  Similarly as was done in [31], we found the following inequalities: which implies that () → 0 as  → ∞.This implies that the given linear system is asymptotically stable and the periodic solution is asymptotically orbitally stable.
Control in the infected host individuals

Table 1 :
Parameters used for numerical simulation.