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.

Introduction

The forest has a significant role in human life, therefore, it is necessary to put in place safety measure in order to protect the trees from being infected with diseases. The trees do not only provide greenery to the environment but also provide pleasant atmosphere for human community. The pine wilt disease (PWD) is one of the major threat to the forest and ecosystem. It is one of the dramatic disease, that kills, pine trees within a very short time and one of the symptom of this diseases is the red dish brown foliage. The North American native pinewood nematode (Bursaphelenchus xylophilus) which is vectored by Monochamus alternatus (Japanese pine sawyer beetle, a species that transmits the nematode to healthy trees) is one of the causes of pine wilt disease[1-5]. This disease has been known since 1900s in Japan[6], and 1980s[7] in China and causes great economic and environmental loss to ecosystems worldwide.

This disease has three organisms: a pinewood nematode, a gymnosperm host, and an insect vector. During primary transmission, dauer juveniles (JIV stage) of Bursaphelenchus xylophilus are carried phoretically in the tracheae of their beetle host to young twigs of susceptible trees, where they enter through resin canals in wounds made during maturation feeding by the insect[6, 8]. The dead or dying conifer is a suitable breeding host for the next generation of Monochamus spp. vectors. Nematodes brought into the conifer during oviposition of the beetle (secondary transmission) will moult from the dauer juvenile stage and enter the propagative phase to grow and reproduce on secondary fungi that are present as the pine host dies. Then the cerambycid eggs hatch and develop through several larval instars while producing galleries, at first in the inner bark, cambium and outer sapwood and later in deeper woody tissue[8-10]. The pine wilt disease is caused by Bursaphelenchus xylophilus [11, 12] and the vector is the pine sawyer beetles (Monochamus alternatus). The nematode are scattered by the vector beetles over pine healthy trees from space, while in the vector population the direct transmission can happen while vectors are mating[13, 14]. During breeding, mature beetles use healthy tree twigs mainly for feeding purposes while they focus on the infected trees only for copulation and oviposition[15]. The transmission of Bursaphelenchus xylophilus horizontally between heterosexual vectors promotes multiple infections. The PWD’s first outbreak in 1905 occurred in Japan, the disease has spread throughout the country in the 1970s and only excluding the northern parts of Japan[16]. After a decade, PWD had spread to many parts of Asia, such as China, Taiwan, Hong Kong and South Korea and in 1999, the disease hit Europe (Portugal)[17]. Today, PWD has become one of the major threats all over the world to forests ecosystems.

Pine wilt disease is very disastrous and can within weeks kill affected trees. The death of the trees are actually caused by a microscopic pine wood nematodes (Bursaphelenchus xylophilus) and the management of pine wilt disease is limited to prevention primarily. There are no cures for the disease once a susceptible tree becomes infected with pinewood nematode and these parasitic nematodes do not associate with the plant roots, but to the aboveground parts of the pine trees. The spread of the nematode from trees to trees occur via the pine sawyer 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-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-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 (t), pine trees that are already exposed, E , and the pine trees that are infected I at any time t, with .

The total population of vector (beetles) is denoted by N , which is categorized further into subclasses, namely, the susceptible beetles, S , the exposed vector beetles, E (t) (which not carrying pinewood nematode), and infected vector beetles, I (t)(which have the ability to carry pinewood nematode) at any time t, with .

We are not interested to include the class R (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:

Figure 1

Flow chart for the transmission for the Pine wilt disease PWD.

The parameter Λ 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 I 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 I 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 Λ, 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 I . The transfer rate from E to I , 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: , , , , and .

The total population of pinewood trees is . The total population of vector beetles is .

When t → ∞, the total dynamics of pine wood trees and beetles approaches . Thus, the biological feasible region for model (1) iswhich is positively invariant and the global is attracted in .

The disease free equilibrium of the system (1) denoted by and is given by .

Basic reproduction number

Here, we will compute the basic reproduction number of system (1). The concept of next generation matrix and basic reproduction number in refs 32, 33 will be used to obtain for the proposed model (1). Hence, we define the new vector contains only the infected variables. Consider the following system:In the computation of [32], the necessary matrices involved are obtained as follows:The inverse of V equalsThus, the next generation matrix of system (2) isSo the basic reproduction number iswhere ρ, , represent respectively, the spectral radius, disease free equilibrium of trees and disease free equilibrium of vector. In the following, we show that is the key threshold parameters whose values completely characterize the global dynamics of system (1).

Stability Analysis of Disease free Equilibrium

LetDefine , so s(M) is a simple eigenvalue of M with a positive eigenvector[34]. It follows from ref. 32, that two equivalences hold: , .

Theorem 0.1. If , then the disease-free equilibrium is locally asymptotically stable on .

Proof: To show this results, we check the hypothesis present in ref. 32, namely (A 1)–(A 5). The hypothesis (A 1)–(A 4) can be verified easily, while A 5 could be satisfied if all the eigenvalues of the 6 × 6 matrixhave negative real parts, where J 3 = −F,Calculated the eigenvalues of J 4,If , then s(M) < 0 and , the disease-free equilibrium 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.

Theorem 0.2. (Theorem 4.1 of Castillo-Chavez and Song (2004))[35] Consider the following general system of ordinary differential equations with a parameter ϕ. where 0 is an equilibrium point of the system (that is, ) and assume

is the linearization matrix of (3) around the equilibrium point 0 with ϕ evaluated at 0. Zero is a simple eigenvalue of A and other eigenvalues of A have negative real parts;

Matrix A has a right eigenvector w and a left vector v (each corresponding to the zero eigenvalue).

Let f be the k-th component of f and The local dynamics of the system (3) around 0 is totally determined by the signs of a and b.

a > 0, b > 0. When ϕ < 0 with , 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when , 0 is unstable and there exists a negative and locally asymptotically stable equilibrium;

a < 0, b < 0. When ϕ < 0 with , 0 is unstable; when , 0 is locally asymptotically stable, and there exists a positive unstable equilibrium;

a > 0, b < 0. When ϕ < 0 with is unstable, and there exists a locally asymptotically stable negative equilibrium; when , 0 is stable, and a positive unstable equilibrium appears;

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 , we haveThen, we make the following change of variables S  = x 1, E  = x 2, I  = x 3, S  = x 4, E  = x 5 and I  = x 6. In addition, using vector notation , the PWD model can then be written in the form , with , as shown below:To follow the above method, we find the Jacobian matrix evaluated at disease free equilibrium is given byThe Jacobian matrix has a simple zero eigenvalue evaluated at . The Jacobian matrix evaluated at has right and left eigenvector denoted by and respectively, and obtainand

Computation of a:

To compute a, we find the following partial derivativeswe obtain

Computation of b:

To compute b, we find the following partial derivativeswe obtainObviously, 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 .

Theorem 0.3. If , then the disease-free equilibrium is globally asymptotically stable on .

Proof: Consider the following lyapunov function:The derivative of L along the solution of model (1) iswhere w , for are positive constants to be chosen later.Using and in equation (6), we getNow choosing the constants such as w 1 = w 2 = δ, w 3 = (d 1 + δ), , and simplifying, we obtainThus, L′(t) is negative for and zero if and only if , , E  = I  = 0 and E  = I  = 0. Therefore the largest compact invariant set in is the singleton set . So, the model (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 and is given by

Lemma 0.1. A unique endemic equilibrium exists for model (1) whenever .

Theorem 0.4. The endemic equilibrium of the model (1) is locally asymptotically stable when and the conditions

.

are satisfied.

Proof: The jacobian matrix evaluated at isThe associated characteristics equation of iswhereand and . 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  > 0 for whereHere, all c  > 0 for 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 by using the approach presented in refs 36–38. At the steady state the model (1) at isThe above equations will be used in the following equations (10–15).

Theorem 0.5. If , then the endemic equilibrium is globally asymptotically stable.

Proof: Consider the lyapunove function:The derivative of L along the solutions of system (1) isBy direct calculations, we have that: andIt follows from (10–15)In equation (16),One can see that the largest invariant subset, is . So, by LaSalle’s invariance Principle[39], is globally asymptotically stable whenever .

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, and , 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 , which is define briefly above. The control variables associate to the variables S , E , I , S , E and I are bounded and measure withHere, for the control problem we defined the objective functional assubject 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, , such thatwhere 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 (t), E (t), I (t), S (t), E (t) and I (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:whereand and represents the derivative with respect to time t. The system (21) is a nonlinear system with bounded coefficient. By settingThe term F on right hand side of (22) satisfieswhere 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 havewhere Z = Z 1 + Z 2 + Z 3 + Z 4 + Z 5 + Z 6 + ||K|| < ∞. 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 satisfiesTo 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 aswhere C 1, C 2, B 1, B 2, B 3, l 1, l 2 > 0 and m > 1. So, we get

Theorem 0.6. The objective functional (19) and the control set (18) subject to control system (17) there exists an optimal control such that .

In order to obtain the solution to the control problem, it is necessary to get the Lagrangian and Hamiltonian (17). Hence, the Lagrangian L for the control problem is expressed asand for the Hamiltonian H, by choosing , 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 , and represent the solutions to the control problem (17), then, there exists the adjoint variables λ for i = 1, 2, 3, 4, 5, 6 satisfying the following conditions given below.Here, by applying the necessary conditions to the Hamiltonian H, the following is obtained:

Theorem 0.7. The given optimal controls and the solutions of the state system, there exists the adjoint variables λ for i = 1, 2, 3, 4, 5, 6 with transversality conditions Further, the controls, are given by

Proof: To find the adjoint system (28) and the transversality conditions (29), we use H (26) by setting , , , , and and taking the time derivative of H with respect to S , E , I , S , E , I , we obtain (28). To find the optimal control characterization (30), we use , for i = 1, 2, 3.

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, B 2 = 0.3, B 3 = 3. 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 .

Table 1

Values of parameter used in numerical simulation of the optimal control system.

Notation	Value	References
ΛH	0.002021/day	Assumed
c	0.001241/day	Assumed
κ 1	0.00166/day	48
κ 2	0.0004/day	48
ψ	0.20/day	49
γ	0.00220/day	Assumed
d 1	0.0000301/day	50
δ	0.0133/day	Assumed
α	0.0032/day	assumed
η	0.00305/day	51
μ	0.01/day	Assumed
ΛV	0.0132652/day	Assumed
α o	0.21/day	Assumed
ϕ	0.0023/day	Assumed
d 2	0.011764/day	52

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

Figure 2

Comparison of both the systems: without and with control.

Figure 3

Comparison of both the systems: without and with control.

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.

Figure 4

Comparison of both the systems: without and with control.

Figure 5

Comparison of both the systems: without and with control.

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

Figure 6

Comparison of both the systems: without and with control.

Figure 7

Comparison of both the systems: without and with control.

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.

Figure 8

Comparison of both the systems: without and with control.

Figure 9

Comparison of both the systems: without and with control.

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 . Whenever the basic reproduction number , the disease free equilibrium is stable both locally and globally. If the basic reproduction number , 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 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].