Mathematical modeling and stability analysis of Pine Wilt Disease with optimal control

This paper presents and examine a mathematical system of equations which describes the dynamics of pine wilt disease (PWD). Firstly, we examine the model with constant controls. Here, we investigate the disease equilibria and calculate the basic reproduction number of the disease. Secondly, we incorporate time dependent controls into the model and then analyze the conditions that are necessary for the disease to be controlled optimally. Finally, the numerical results for the model are presented.

beetle (Monochamus spp.) as they feed on healthy trees branches, the nematodes would then emerge and thereafter enter through the feeding wounds caused by the beetles into the trees. Infected trees thereafter are killed by nematodes feeding on cells surrounding the resin ducts.
During oviposition (egg laying), the adult sawyers are by preference only attracted to the dead or dying trees. And the tree age also influences the risk of pine wilt. Trees that are 10 years and above are the most affected, beetles feed and infest one or more trees with nematodes and then move to other susceptible area. But very often the disease only affects a tree in a group. This is due to the nematodes movement, which is not based on physical contact or water or through grafts of roots 18,19 .
Some mathematical models have been presented on the dynamics of PWD. A mathematical model presented in ref. 20, examined the stability of pine wild disease and application of optimal control technique, where the population of pine trees were divided into two categories, that is, susceptible pine trees and infected pine trees, while the vector population (beetles) were divided into two classes; susceptible vector and infected vector 20 . Also, K. S. Lee and A. A. Lashari 21 introduced a mathematical model that incorporated the exposed class in the pine trees population with a detailed discussion made on the stability and optimal control of PWD. Also, M. Ozair presented a mathematical model on the dynamics of PWD by dividing the host pine trees and vector beetles into susceptible and infected classes with nonlinear incidence and horizontal transmission 22 . Recently, K. S. Lee and D. Kim introduced a mathematical model that describe the dynamics of PWD by presenting its global stability with nonlinear incidence rates 23 .
It is necessary to make the community and forest safe from infectious diseases, with different control measure used. Such as optimal control technique and mathematical models are widely used to understand the complicated nonlinear phenomena of infectious diseases with different purposes [24][25][26][27] which is helpful in analyzing biological models. In the literature, various articles on the dynamics of infectious disease with optimal control strategies are presented [28][29][30][31] and the reference there in.
The aim of this paper is to present a mathematical model on the dynamics of Pine trees and vector (beetles) population. We first presented the detail mathematical study of the model, which is the local and global stability, asymptotical stability and backward bifurcation phenomena. Then, we applied the optimal control technique to minimize the population of exposed pine trees, infected pine trees, susceptible vector, exposed vector and infected vector (beetles) as well as to maximize the population of susceptible pine trees. Different control strategies have been presented in order to reduce infection in the population of pine trees.
In the next, section Basic Model Formulation we give a detail analysis about the mathematical formulation of PWD. In section Stability analysis disease free a brief mathematical results are presented for disease free case and a backward bifurcation analysis. The stability analysis of the model at endemic equilibrium is presented in Section stability endemic equilibrium. Further, we apply the optimal control technique in section Optimal control problem while the numerical results and conclusion are presented in sections Numerical results and Conclusion respectively

Basic Model Formulation
The total population of pine wood trees is denoted by N(t) and we subdivide into four subclasses; the pine trees that are susceptible, S H (t), pine trees that are already exposed, E H , and the pine trees that are infected I H at any H . The total population of vector (beetles) is denoted by N V , which is categorized further into subclasses, namely, the susceptible beetles, S V , the exposed vector beetles, E V (t) (which not carrying pinewood nematode), and infected vector beetles, I V (t)(which have the ability to carry pinewood nematode) at any time t, We are not interested to include the class R H (t) for pine wood trees population, this is because the infected pine wood tree dies in a year or may in a subsequent next years.
With the above assumptions, we present the model and their schematic diagram in Fig. 1 is follows: The parameter Λ H denotes the recruitment rate of pine trees that are susceptible and the parameter κ 1 stands for the contact rate during maturation. The average number of contacts during maturation per day with vector adult beetles is denoted by ψ. The mass action term κ 1 ψS H I V represents the incidence rate. The parameter κ 2 , denote the probability that a nematode is being transmitted through oviposition by an infected beetle, and the average number of contacts per day when adult beetles oviposit is denoted by φ. The probability that pine trees that are susceptible cease oleoresin exudation without infected by the nematode is represented by α. We show the transmission through oviposition by κ 2 φα and hence, κ 2 φαS H I V denotes the number of new infections. The rate of progression from pine trees that are exposed to trees that are infected and the natural death rate of pine trees respectively denoted by δ and d 1 .
The vector pine beetle emergence rate is denoted by Λ V , while η measure the rate at which adult beetles that are escaping from dead trees carry the PWN with them and so the transmission via this route is denoted by ηS V I H . The transfer rate from E v to I v , the natural death rate and disease induced death rate for the vector beetles population are respectively denoted by μ, d 2 and γ. The model (1) presents the pine wood trees and vector beetles populations, and it is understood that all the variables involved with parameters are nonnegative for nonnegative initial conditions. So, the initials conditions for the model (1) is follows as: When t → ∞, the total dynamics of pine wood trees and beetles approaches Thus, the biological feasible region for model (1) is which is positively invariant and the global is attracted in Ψ.

Basic reproduction number.
Here, we will compute the basic reproduction number 0  of system (1). The concept of next generation matrix and basic reproduction number in refs 32, 33 will be used to obtain 0  for the proposed model (1 The inverse of V equals Thus, the next generation matrix of system (2) is So the basic reproduction number is represent respectively, the spectral radius, disease free equilibrium of trees and disease free equilibrium of vector. In the following, we show that 0  is the key threshold parameters whose values completely characterize the global dynamics of system (1).

Stability Analysis of Disease free Equilibrium
Calculated the eigenvalues of J 4 , , the disease-free equilibrium  0 of system (1) is locally asymptotically stable.

Analysis of Backward Bifurcation
In this subsection, we analyze the existence of Backward bifurcation for the model (1). To analyze this, we use the center manifold theory as described in Castillo-Chavez and Song (2004) (Theorem 4.1) 35 , which is reproduced here below for convenience. (2004)) 35 Consider the following general system of ordinary differential equations with a parameter φ.

Theorem 0.2. (Theorem 4.1 of Castillo-Chavez and Song
) and assume The local dynamics of the system (3) around 0 is totally determined by the signs of a and b.
0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when φ <  0 1 , 0 is unstable and there exists a negative and locally asymptotically stable equilibrium; 1 , 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; (iii) a > 0, b < 0. When φ < 0 with φ  1 is unstable, and there exists a locally asymptotically stable negative equilibrium; when φ <  0 1 , 0 is stable, and a positive unstable equilibrium appears; (iv) a < 0, b > 0. When φ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.
Particularly, if a > 0 and b > 0, then a backward bifurcation occurs at φ = 0. If we choose κ 1 as a bifurcation parameter, then at Then, we make the following change of variables 5 and I V = x 6 . In addition, 3 4 5 6 , the PWD model can then be written in the form To follow the above method, we find the Jacobian matrix evaluated at disease free equilibrium 0  is given by The Jacobian matrix  J( ) 0 has a simple zero eigenvalue evaluated at κ ⁎ 1 . The Jacobian matrix  J( ) 0 evaluated at κ ⁎

Computation of a:
To compute a, we find the following partial derivatives we obtain

Computation of b:
To compute b, we find the following partial derivatives Obviously, the coefficient b is positive always so that, according to Theorem (0.2), it is the sign of the coefficient a which decides the local dynamics around the disease-free equilibrium for κ κ , then the disease-free equilibrium 0  is globally asymptotically stable on Ψ.
Proof: Consider the following lyapunov function: The derivative of L along the solution of model (1) is where w i , for = … i 1, 2, 6 are positive constants to be chosen later.
Now cho osing t he const ants such as and simplifying, we obtain  (1) is globally asymptotically stable.

Stability Endemic Equilibrium
In this subsection, we investigate the stability results for the endemic case. The endemic equilibrium of the model (1) denoted by = ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ S E I S E I ( , , , , , ) and is given by     and δ γ µ 3 and κ ψ κ φα = + P ( ) 1 2 . In equation (7) one of the root −d 2 is clearly negative, for the remaining the roots we use the following Routh-Hurtwiz conditions. c i > 0 for Here, all c i > 0 for = … i 1, 2, 5 and the parameters are positive. Thus, the given conditions in the theorem above ensure the local endemic stability of the endemic equilibrium of the model (1).

Global stability of Endemic Equilibrium
Here, we present the global stability of the system (1) at  1 by using the approach presented in refs 36-38. At the steady state the model (1) at  1 is The above equations will be used in the following equations (10-15).
Proof: Consider the lyapunove function: The derivative of L along the solutions of system (1) is By direct calculations, we have that: It follows from (10-15) In equation (16), One can see that the largest invariant subset, =  L 0 is  1 . So, by LaSalle's invariance Principle 39 , 1  is globally asymptotically stable whenever > 1 0  .

Optimal Control Problem
Here, we incorporate control measures into the system (1) by including the density effects respectively to modifying the recruitment rate of the susceptible classes, that is, Λ → Λ + cN H H H and Λ → Λ N V V V , where the constant c represent the density impact on recruitment rate 40 . Our main goal is to investigate the best control strategies with minimum cost of implementation. The force of infections in pine trees population is reduce by (1 − u 1 ), where u 1 measures the effort due to the precautionary measures, such as nematode tree-injection and vaccination. In order to prevent infection of the population of pine trees, a nematicide-injection preventative control is used. The deforestation of infected trees effort is denoted by the control variable u 2 . The destruction and removal of infected pine trees can drastically reduce further infection. This will further ensure that eggs, larvae and pupa that are inhabiting the pines are destroyed. While the control variable u 3 is the effort due to eradication through aerial insecticide spraying. Here, we consider that the death rate of the population of adult vector (beetle) increases as u 3 increases. The factor (1 − u 3 ) measure the reduction in the population of adult beetles. We present the control model based on the assumptions and extensions made above,  subject to non-negative initial conditions. In system (17), we use three control variables = ∈ u t u u u ( ) ( , , ) 1 2 3  , which is define briefly above. The control variables  = ∈ u t u u u ( ) ( , , ) Here, for the control problem we defined the objective functional as 1  2  3  0  1  2  1 1  2  2 2  2  3 3  2 subject to the control system (17). The constants in (19), C 1 , C 2 , B 1 , B 2 and B 3 are the weight and balancing constants. These constants C 1 , C 2 and B 1 , B 2 , B 3 measure respectively the relative cost of the interventions over the interval [0, T]. In order to find an optimal control, ⁎ ⁎ ⁎ u u u , , where  is defined in (18) and subject to control system (17) with non-negative initial conditions. Next, we use the well known Pontryagin's Maximum Principle to find the solution to the control problem and derive its necessary conditions. Here, we show the existence for the control system (control problem). Let the variables S H (t), E H (t), I H (t), S V (t), E V (t) and I V (t) be represent the state variables with control associated control variables u 1 (t), u 2 (t) and u 3 (t). For existence, we rewrite the control problem (17) in the following form: and ′  represents the derivative with respect to time t. The system (21) is a nonlinear system with bounded coefficient. By setting The term F  ( ) on right hand side of (22) satisfies where Z is the positive constant, Z = max(Z 1 , Z 2 , Z 3 , Z 4 , Z 5 , Z 6 ) is independent of the state variables. Also we have Thus, it follows that the function G( )  is uniformly Lipschitz continuous. From the definition of control variables we can see that a solution of the system (17) exists 41 .
Existence of the control problem. Following the result in ref. 42 we prove the existence of an optimal control problem. In control system (17), the equations are obviously bounded above and so the result in ref. 42 can be applied to the model (17), since the state variables and the set of control variables are nonempty, the set  of the control variables is closed and convex. The right hand side of each equation in control problem (17) is continuous, bounded above by a sum of the bounded control and state, and can be written as a linear function of u with coefficients depending on time and state and there exists constants l 1 , l 2 > 0 and m > 1 such that the integrand L(y, u, t) of the objective functional J is convex and satisfies To satisfy the conditions mentioned above, we use the result given in ref. 43 to establish the existence of (17). The state variables and the set of control is obviously bounded and nonempty. The solutions are bounded, and convex.
The system is bilinear in control variables (since the solutions are bounded). The last one can be verified as and λ = (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 ), we have: Solution to the Optimal control problem. To obtain the optimal solution of the control system (17), we use the well-known Pontryagin's Maximum Principle 44 : Letting ⁎ u 1 , ⁎ u 2 and ⁎ u 3 represent the solutions to the control problem (17), then, there exists the adjoint variables λ i for i = 1, 2, 3, 4, 5, 6 satisfying the following conditions given below.  Table 1. Values of parameter used in numerical simulation of the optimal control system.

Numerical Simulations and Results
Here, we are investigating the numerical solutions of the system (17) and that of the model without control (1). For the solution of both systems (that is the control and without control) many method are available in the literature 45,46 . In this simulation, the cases without control population are labeled with bold line and the control cases by a dashed line. The constants (weight) values in the objective functional are C 1 = 0.01, C 2 = 0.0036, B 1 = 0.02, The parameters values used in the optimal control solution is chosen in such a way, that the number of infected trees, exposed trees, susceptible vector, exposed vector and infected vector decreased while the population of susceptible trees increased. These are given in Table 1. We adopted different control strategies to minimize the infection in pine trees population by considering , , , and .
In this strategy, we set the control variable u 1 = 0 zero and the rest of the control variables are non-zero, this effect can bee seen in Fig. 2, (sub-figures (a-c)) and Fig. 3 (sub-figures, (a-d)). The population of susceptible trees increased while the exposed, infected trees and vector population decreased.

Strategy 2:
In this strategy, we set the control variable u 2 = 0 and the rest of the control variables are non-zero, this impact can bee seen in Fig. 4, (sub-figures (a-c)) and Fig. 5 (sub-figures, (a-d)). The population of susceptible trees increased while the exposed, infected trees and vector population decreased. One can observe in Figs 2(c) and 4(c), that the infected trees population did not decrease compare to strategy 1.

Strategy 3:
In this strategy, we set the control variable u 3 = 0 and the rest of the control variables are non-zero, this effect can bee seen in Fig. 6, (sub-figures (a-c)) and Fig. 7 (sub-figures, (a-d)). The  population of susceptible trees increased sharply while the exposed, infected trees and vector population decreased. One can see that infected population significantly reduced, see Fig. 6(c).

Strategy 4:
In this strategy, we activate all the control. The figures for this strategy can be seen in Fig. 8, (sub-figures (a-c)) and Fig. 9 (sub-figures, (a-d)). The population of susceptible trees increased, while the exposed, infected trees and vector population decreased. Thus, by comparing all the control stratigies from 1-4, one can observe that the strategy 4 is the best strategy to control infections in the pine trees population.

Conclusion
We have considered a mathematical system of equation which describes the pine wilt disease. The analysis of the system is well established. The stability analysis of the disease free and endemic equilibria is presented on the basis of basic reproduction number 0  . Whenever the basic reproduction number  < 1 0 , the disease free equilibrium is stable both locally and globally. If the basic reproduction number  > 1 0 , then the endemic equilibrium equilibrium is stable both locally and globally. A bifurcation analysis of the model is presented and a control problem is formulated by using three control variables, two control for pine trees population and one control for vector population. The mathematical results for the control problem with different control strategies are presented and concluded that the strategy 4 is the best control strategy for cost minimization in the population of pine trees diseases. This work generalized the work presented in ref. 21 by incorporating the exploded class E V in the vector population. In ref. 21, the optimal control solution is investigated without defining the control strategies. In this new work, we suggests different control strategies for the eradication of infection in Pine trees. So, this work compared to the previous study is more helpful for the elimination of pine trees infections 47 .