Abstract
We propose a mathematical model of pine wilt disease (PWD) which is caused by pine sawyer beetles carrying the pinewood nematode (PWN). We calculate the basic reproduction number and investigate the stability of a disease-free and endemic equilibrium in a given mathematical model. We show that the stability of the equilibrium in the proposed model can be controlled through the basic reproduction number . We then discuss effective optimal control strategies for the proposed PWD mathematical model. We demonstrate the existence of a control problem, and then we apply both analytical and numerical techniques to demonstrate effective control methods to prevent the transmission of the PWD. In order to do this, we apply two control strategies: tree-injection of nematicide and the eradication of adult beetles through aerial pesticide spraying. Optimal prevention strategies can be determined by solving the corresponding optimality system. Numerical simulations of the optimal control problem using a set of reasonable parameter values suggest that reducing the number of pine sawyer beetles is more effective than the tree-injection strategy for controlling the spread of PWD.
1. Introduction
Pine wilt disease (PWD) is caused by the pinewood nematode Bursaphelenchus xylophilus. The vector for this parasite is the Japanese pine sawyer beetle Monochamus alternatus, which disperses the nematode to healthy trees [1–8]. PWD is lethal to healthy pine trees such as Pinus densiflora S. et Z., P. thunbergii Parl., and P. sylvestris L. [9]. The disease (PWD), caused by the pinewood nematode (PWN), Bursaphelenchus xylophilus, was first observed in 1905 in Japan [10]. By the 1980s, the PWD epidemic had spread to other Asian countries such as China, Taiwan, Hong Kong, and Korea. It subsequently reached Europe (Portugal) in 1999 [11–13]. Thus, PWD is the most serious global threat to pine forest ecosystems worldwide.
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–18]. Mathematical models have become important tools in analyzing the spread and control of infectious diseases. Thus, we present 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.
2. 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 , and infected vector beetles that carry the PWN . The total population sizes at time for the host pine trees and vector beetles are denoted by and , 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 . 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 . 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 . A constant emergence rate of the vector pine sawyer beetle is denoted by , 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 . Finally, , denote the natural death rate of the host population and the natural death rate of the vector population, respectively.
Since the model deals with host pine trees and vector beetle populations, it can easily be seen that every state variable will remain nonnegative for nonnegative initial conditions, that is, for all , , , , and . And we make some reasonable technical assumptions on the parameters of the model, namely, , , , , , , , , and . The above systems for the host population and the vector are also equipped with initial conditions as follows: , , , , and . The total host population dynamics is given by . The given initial conditions make sure that . The total dynamics of vector population is . We know that , .
Obviously, is positively invariant, and system (1) is dissipative and the global attractor is contained in . On , we have . Hence, we will study the following three-dimensional nonlinear system:
For system (3), the region is positively invariant.
System (3) has always a disease-free equilibrium and basic reproduction number .
If , system (3) has a unique endemic equilibrium in the interior of , where
Mathematical analysis of the system (3) is provided in Appendices A and B.
3. A Model for Optimal Control of PWD
In the host pine tree population, the force of infections is reduced by , where measures the level of successful prevention efforts. It follows that control variable 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 . In the vector adult beetles population, the control variable 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 . Further, we assumed that the mortality rate of the vector population increases at a rate proportional to , where 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 , , , and .
4. Analysis of Optimal Control
We define the objective functional 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 , . For the optimal control problem of the given system, we consider the control variables relative to the state variables , , , and where both control variables are bounded and measured with In the objective functional the quantities , represent the weight constants of the infected host population and total vector population, respectively. In the objective functional and are weight factors for reduction of vector-host contacts and vector control, respectively. The terms and 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 , such that 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.The set of controls and corresponding state variables is nonempty.The control set is convex and closed.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.There exist constants and 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 . We note that the solutions are bounded. Our control set satisfies condition . Since our state system is bilinear in , , the right-hand side the state system (5) satisfies condition , using the boundedness of the solutions. Note that the integrand of our objective functional is convex. Moreover we have the last condition needed where , , and . Thus we obtain the following.
Theorem 1. Given the objective functional , where subject to the system (5) with initial conditions, then there exists an optimal control such that .
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 , , and , to obtain In order to find the necessary conditions for this optimal control pair, we apply Pontryagin’s Maximum Principle [21] as follows.
If is an optimal solution of an optimal control problem, then there exists a nontrivial vector function 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 , and solutions , , , and of the corresponding state system (5), there exist adjoint variables , , , and satisfying with transversality conditions Furthermore, , are represented by
Proof. To determine the adjoint equations and the transversality conditions we use the Hamilton . We differentiate the Hamiltonian with respect to , , , . Thus we obtain with transversality conditions Using the optimality conditions we have and the property of the control space we can drive
5. 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 , , , and initial state variables , , , and 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 and );(ii)strategy 2: physical or chemical treatment of wilt pines to kill their larvae or eradication of vector beetles by aerial pesticide spraying ( control alone);(iii)strategy 3: precaution effort by tree-injection of a nematicide ( control 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).
5.1. Strategy 1: Optimal Use of Tree-Injection of a Nematicide () and Eradication of Vector Beetles ()
The computed results for optimal control under strategy 1 (, ) are shown in Figures 1, 2, 3, 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 , and infected vectors . With this strategy, an increase in the number of susceptible hosts is observed in Figure 1.
5.2. Strategy 2: Optimal Use of Eradication of Vector Beetles ()
The computed results for optimal control under strategy 2 (, ) 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.
5.3. Strategy 3: Optimal Use of Tree-Injection of a Nematicide ()
The computed results for optimal control under strategy 3 (, ) 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 , 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.
(a)
(b)
5.4. 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 , are shown in Figures 14 and 16 (, , , fixing ) and Figures 15 and 17 (fixing , , , and ). In Figure 16, the value of the weight constant 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 is varied. As the value of weight constants increases, control functions decrease, and increasing the cost of vector control, the population of infected vector is increased.
6. 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 , the disease-free equilibrium is globally asymptotically stable in , and thus the disease always dies out. However, if , 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 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.
Appendices
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 , then the disease-free equilibrium is locally asymptotically stable.
Proof. Linearize the system (3) around the disease-free equilibrium . The matrix of the linearization at is given by
The characteristic equation of this matrix is given by det, where is the unit matrix. Expanding the determinant into a characteristic equation we obtain the following equation:
So, and , where
Thus, according to the Routh-Hurwitz criteria [27], the quadratic equation has only roots with negative real parts. Hence, the disease-free equilibrium 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 , then the disease-free equilibrium is globally asymptotically stable in .
Proof. We define the following Lyapunov function on .
where
Clearly, along the solution of the system (3) and is zero if and only if both and are zeros. The derivative of along the solutions of (3) is
We know that
Thus
We see that if , the derivative if and only if , while in the case of the derivative if and only if or . Consequently, the largest compact invariant set in , when , is the singleton . Hence by Lasalle’s Invariance Principle [28], 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 real matrix. If , , and 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 , 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:
From the Jacobian matrix , the second additive compound matrix is given by
From the Jacobian matrix , we have
And we know that system (3) is given
Thus,
since we know that
So
Computing directly the determinant of and from (B.4) and (B.6) we can get
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 is greater than the unity. For this, we will prove the following result.
Theorem B.3. If , then system (3) is uniformly persistent; that is, there exists (independent of initial conditions), such that , , and .
Proof. Let be a semidynamical system (3) in , a locally compact metric space, and . The set is a compact subset of and is positively invariant set of system (3). Let be defined by and set , where is sufficiently small so that . Assume that there is a solution such that, for each , we have . Let us consider the following:
where is a sufficiently small constant so that . By a direct calculation, we have
Let
Thus, we have
The inequality (B.12) implies that as . However, is bounded on the set . According to Theorem 1 in [32], we complete the proof of Theorem 2.
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–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 , then the endemic equilibrium is globally asymptotically stable in .
Since system (3) is competitive and persistent, and is locally asymptotically stable if . 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 , 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 is asymptotically stable, where is the second additive compound matrix of the Jacobian . The Jacobian of system (3) is given by For the solution , (B.15) becomes In order to prove that (B.15) is asymptotically stable, we will use the following Lyapunov function which is similar to the one found in [31] for the SEIR model. Let where is the norm in defined by Since the given system is uniformly persistence, we obtain that the orbit of remains at a positive distance from the boundary of . Therefore, we have , with . Hence, the function is well defined along and Along the positive solution of the system, becomes Similarly as was done in [31], we found the following inequalities: We know that Therefore Thus we get where Since we have Hence we find Therefore, by Gronwall’s inequality, we obtain which implies that as . This implies that the given linear system is asymptotically stable and the periodic solution is asymptotically orbitally stable.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2013R1A6A3A01060805).