Bio-inspired Analytical Heuristics to Study Pine Wilt Disease Model.

This paper portrays the dynamics of pine wilt disease. The specific formula for reproduction number is accomplished. Global behavior is completely demonstrated on the basis of the basic reproduction number [Formula: see text]. The disease-free equilibrium is globally asymptotically stable for [Formula: see text] and in such a case, the endemic equilibrium does not exist. If [Formula: see text] exceeds one, the disease persists and the unique endemic equilibrium is globally asymptotically stable. Global stability of disease-free equilibrium is proved using a Lyapunov function. A graph-theoretic approach is applied to show the global stability of the unique endemic equilibrium. Sensitivity analysis has been established and control strategies have been designed on the basis of sensitivity analysis.

Introduction

Forests play a vital role in human life, so it is essential to set up safety measures for the protection of infected trees. The trees do not just give greenery to nature, but also supply a charming climate for humans. Pine wilt disease (PWD) rapidly spreads in the tree, and within a few weeks, the tree starts exhibiting the symptoms and moves towards death quickly. Now, this disease has become a threat to the forest ecosystem and to pine forests. It is a very dangerous and wild disease among pine trees with huge economic losses which disturb the development of social sustainability. There is no way to cure the pine tree once infected. Consequently, prevention is the standard technique[1].

The early occurrence of PWD goes back to . During the the annual loss of pines increased from to , and during 1940s destructions were estimated at volumes of every year. The struggle to control the disease was abandoned during World War II, which resulted in an expeditious increase of infected pine trees and an immense loss of in 1948. Considerable efforts lessened the annual loss to during the 1950s. Elimination of dead trees through tumbling and blazing was used as the basic control method at that time. In the , annual losses exceeded million again. Prevalent growth in the infested region was identified as one of the most characteristic features of the epidemic at that time. Considerable timber loss, million , was recorded in [2].

Bursaphelenchus xylophilus is the nematode that is the main cause of this disease. It was first described in 1971 as the causal agent of PWD of native pines in Japan[3]. Monochamus alternatus (pine sawyer beetle) is the vector for this parasite. The initial notable symptom of PWD is the reduction in the emission of oleoresin from wounds of bark. A further indication is the conversion of needle color. It becomes light grayish green to yellowish-green then yellowish-brown. Lastly, fully brown when tree capitulates to the malady[4,5].

Three types of transmission of PWD are discovered:

First, when infected adult beetles fly to susceptible pine trees and start maturation feeding, transmit nematode in it, this is known as the primary transmission[6].

Secondary transmission develops when mature females lay eggs on dead or dying, newly cut pine trees[7].

The third is the horizontal transmission of the nematode. It happens through the mating of male and female bark beetles[8].

Numerous mathematical models have been designed to analyze the transmission dynamics of PWD. In 1999, Yoshimura et al.[9] studied that there is a minimum density of pine under which the invasion of the disease reduces and a minimum density of pine increases unjustifiably with increment in the termination rate. Takasu et al.[10] analyzed the reliance on parasite rate of range expansion on the extermination rate of the beetle, primary pine density, and the vector dispersal capacity. Mechanistic interactions are considered in[11], and it is demonstrated that the Allee effect can appear in the beetle dynamics due to the need for beetles to contact pine trees at least twice to imitate auspiciously. Mathematical models based on ordinary differential equations have been analyzed in[12-16]. Shi and Song[17] studied the mathematical model for the transmission of pine wilt disease by considering two routes. Firstly, healthy pine trees are infected through the maturation feeding of infectious bark beetles and secondly, when mature females lay eggs on infected dead or dying trees. However, they did not consider the third transmission of nematodes among bark beetles through mating that has been suggested by Togashi and Arakawa[18].

In this paper, we revisit the model given in[17] and included the following features: The aim of this paper is to study new features of the models:

Transmission occurs in the vector population through the mating.

The optimal control problem has been designed by performing a sensitivity analysis of basic reproduction number , infected host and infected beetles

To show the global behavior of equilibria.

Identification of sensitive parameters for reproduction number and endemic levels of infectious classes.

Design the optimal control strategy on the basis of sensitivity analysis.

Model Framework with Flow Diagram

We formulate the mathematical model comprises of susceptible host pine trees , infected host pine trees , susceptible vector beetles and infected beetles . The usual mode of nematode transmission between pine trees and bark beetles develops as a result of maturation feeding. The pine sawyers have pinewood nematodes when they emerge from infected pines. However, the vectors may also be directly infected during mating. The population of host pine trees at time is represented by , and whole vector adult beetles are indicated by . Thus, and

Let be the constant recruitment rate of Monochamus alternatus, be the increase rate of pines, be the fatality rate of Monochamus alternatus. Pines are infected by the Monochamus alternatus having Bursaphelenchus xylophilus Nickle. This contact process is assumed to follow the bilinear incident rate . Bark beetles, which do not carry nematode is infected by Bursaphelenchus xylophilus Nickle through infectious pines and this contact process is represented by the incidence rate . Monochamus alternatus is infected through mating and the contact process is denoted by Let be the felling rate of susceptible pines, be the exploitation rate of pines which have Bursaphelenchus xylophilus. Here we assume that .

With these assumptions, the model can be demonstrated by the following system of non-linear differential equations: The flow diagram of the model (1) is shown in Fig. 1.

Figure 1

Flow chart of the model (1).

As the model deals with beetle and tree populations, it is clear that all the state variables will be non-negative for non-negative initial conditions, that is, . The total population of both the species satisfies the equations: Clearly, = is positively invariant, the global attractor is contained in and the model (1) is dissipative.

Equilibria and Basic Reproduction Number

Direct calculations show that the model (1) always attains a disease-free equilibrium (DFE) given by . The dynamics of the malady is characterized by the basic reproduction number that is defined to be the average number of secondary infections, which are accomplished by an infected individual in a totally susceptible population. This enables us to find whether an infectious disease will prevail through population or not. We use the next-generation method to calculate the basic reproduction number.

Let and be the population in disease and non-disease compartments, respectively, and where shows the rate of new infections in ith disease class and shows the transition terms, for example recovery and death in the ith disease class[19]. is acquired by choosing the dominant eigenvalue (spectral radius) of , where and The infected compartments are and Therefore, Thus, which gives Since is the spectral radius of , therefore, For , the system (1) has a unique endemic equilibrium (EE)  = , where and is uniquely calculated from the equation where Noting that while the sign of coincides with that of so that if and only if . Thus we conclude that if , the model (1) has unique EE given as .

Stability of Equilibria

We investigate the global behavior of the system (1) around DFE. For this, we construct the Lyapunov function.

Theorem 1.

For , the DFE of the system (1) is globally asymptotically stable in  .

Proof. Consider the Lyapunov function Taking time derivative of , we have if and only if . Furthermore, if and only if . Thus, it is easy to check that the largest compact invariant set for the system (1) is singleton . By Lasalle’s Invariance Principle[21], is globally asymptotically stable (GAS) in .□

We now will show that the unique EE is GAS using the graph-theoretic approach. For further terminologies, the readers are referred to[20].

Theorem 2.

The unique positive disease-present equilibrium point has global asymptotic stability in the interior of .Proof. Let Differentiating and using the inequality for , we have Similarly, A weighted digraph is constructed with four vertices and six arcs, which is shown in Fig. 2:

Figure 2

Weighted diagraph.

It is easy to see that there are three cycles and along each cycle , . Thus, by Theorem (3.5) from[20], there exist , , , , such that . From Theorem (3.3) given in[20], it follows that Thus, is the Lyapunov function of the system (1). Using this Lyapunov function and LaSalle’s invariance principle[21], it follows that is GAS in int ().□

Sensitivity Analysis

The main thing for an infectious disease is to study its capability to enter in a population. To check that which factors are responsible for the expanse and the existence of the disease, we shall carry out the sensitivity analysis. In this analysis, we explain the following features:

How the endemic level of host and vector population is affected by changing the values of parameters?

The measurement of the corresponding change in the state variable when a parameter changes.

Effect of parameters on the endemic level of infectious hosts

Figure 3 shows the endemic level of infectious hosts for different values of the parameters. Graphs explore that parameters are directly related to and inversely related to . It can also be observed that considerable change occurs in the value of by an increment or reduction of and . However, when we calculate the percentage difference of the values of corresponding to the percentage difference of the parameter values, the most sensitive parameter is , which is the transmission rate of the Bursaphelenchus xylophilus Nickle through the maturation feeding of infectious Monochamus alternatus. The value of this percentage difference of the endemic level is . It means that we should focus on the reduction of the transmission rate of nematode while designing the control policy. Percent increment or reduction in the value of corresponding to all parameters is given in Table 1.

Figure 3

Variation in the endemic level of I∗ with respect to model parameters.

Table 1

Sensitivity of for all model parameters.

Parameters	Initial	Final	Difference	Percentage	Initial	Final	Difference	Percentage	C_5∕91.11111111	C_9*C_10
		Value		Difference	Value of I^*	Value of I^*		Difference
a	10	27	17	170	220	620	400	181.8181818	1.865853659	339.2461197
β	0.001	0.08	0.079	7900	5	245	240	4800	86.70731707	416195.122
α	0.9	0.01	−0.89	−98.88888889	10	240	230	2300	−1.085365854	−2496.341463
δ	0.5	0.03	-0.47	−94	7	304	297	4242.857143	−1.031707317	−4377.38676
β_1	0.0001	00.05	0.0499	49900	206	210	4	1.941747573	547.6829268 k8	1063.461994
b	9	25	16	177.7777778	214	236	22	10.28037383	1.95122	20.05927
k	0.005	0.025	0.02	400	224	256	32	14.28571429	4.390243902	62.71777003
μ	0.9	0.08	−0.82	-91.11111111	224	246	22	9.821428571	−1	−9.821428571

Effect of parameters on the endemic level of infectious bark beetles

It can be observed from the (Fig. 4) that the endemic level of infectious bark beetles is greatly influenced by the parameters and . But the calculation of percentage difference of the values of corresponding to percentage difference of the parametric values shows that the most sensitive parameter is again . The value of this difference is . Percent increment or reduction in the value of corresponding to all parameters is given in Table 2.

Figure 4

Variation in the endemic level of with respect to model parameters.

Table 2

Sensitivity of for all model parameters:.

Parameters	Initial	Final	Difference	Percentage	Initial	Final	Difference	Percentage	C_5/91.11111111	C_9∗C_10
		Value		Difference	Value of Y^*	Value of Y^*		Difference
a	10	28	18	180	10	13	3	30	1.97561	59.26829
β	0.001	0.08	0.079	7900	1	10	9	900	86.7073	78036.59
α	0.9	0.01	-0.89	-98.88888889	1	9.5	8.5	850	-1.085366	-922.561
δ	0.4	0.02	-0.38	-95	2	11	9	450	-1.042683	-469.207
β_1	0.05	0.001	-0.049	-98	15	17	2	13.33333333	-1.075609756	-14.34146341
b	10	26	16	160	7	19	12	171.4285714	1.756097561	301.0452962
k	0.001	0.01	0.009	900	3	12	9	300	9.87804878	2963.414634
μ	0.9	0.08	-0.82	-91.11111111	10	170	160	1600	-1	-1600

Sensitivity indices

The sensitivity index helps to measure the relative change in the state variable when a parameter changes. The normalized forward sensitivity of a variable to the parameter is defined as the ratio of relative change in the variable to the relative change in the parameter[22]. When a variable is a differentiable function of the parameter, the sensitivity index may be defined in terms of partial derivatives.

Definition 1.

The normalized forward sensitivity index of a variable , that depends differentiability on a parameter , is defined as

Sensitivity indices of

The analytical expressions for the sensitivity indices of can be obtained by using the above definition. For example, the sensitivity index of with respect to is The above sensitivity index has no obvious mathematical structure so that we could decide whether is an increasing function of or decreasing function. However, when we use the numerical values of the parameters given in Table 2, we see that the sensitivity index of of the parameter is . It means that when we increase the value of by , the value of increases by almost . Looking at the sensitivity indices of with respect to all the model parameters given in Table 3, we observe that the most sensitive parameter is (the mortality rate of bark beetles), which has the highest value  . Negative sign shows that is a decreasing function of . An increment in the value of by , the reduction in occurs by .

Table 3

Description of the model parameters and nominal values.

Parameters	Explanation	Value
a	Recruitment rate of Monochamus alternatus	10
b	Increase rate of Pines	14
β	Transmission rate of nematode through the maturation feeding of infected bark beetles	0.03
α	Exploitation rate of Susceptible pines	0.03
δ	Exploitation rate of infected pines	0.04
k	Trnasmission rate of Bursaphelenchus xylophilus through the contact of infectious pines and susceptible bark beetles	0.0058
β_1	Trnasmission rate of nematode as a result of mating	0.01
μ	Mortality rate of bark beetles	0.9

Sensitivity indices of and

The endemic level of infectious pines given in (3) is represented in the form of endemic level of infectious bark beetles and it can be obtained by the quadratic equation (4). Therefore, we do not have an explicit expression for the endemic equilibrium. However, by using the values of parameters given in table 3, we can calculate the sensitivity indices of and .

The calculated values of sensitivity indices of given in Table 4, express that the most sensitive parameter for the endemic level of infectious pines is the parameter . It is observed that the increment in the value of by , would decrease almost . The endemic level of infectious hosts is also an increasing function of . By decreasing the value of by , would decrease almost .

Table 4

Values of sensitivity indices of and .

Parameters	Sensitivity index of R_∘	Sensitivity index of I^*	Sensitivity index of Y^*
a	0.51	0.0056	0.0000088
b	0.49	−0.31	1.0
β	0.49	1.45	−0.00070
α	−0.49	−0.0018	−0.0016
δ	−0.49	−1.14	−1.002
k	0.49	−0.31	−0.00049
β_1	0.012	0.45	0.00071
μ	−1.012	−1.60	−0.0025

Sensitivity indices of the endemic level of infectious bark beetles , given in Table 4, show that the parameters and play a vital role in the enhancement or reduction in the values of . is an increasing function of and a decreasing function of . Increment in the value of by , reduction in the value of occurs by almost . Contrary to this, by decreasing the value of   , decreases by almost . To decrease the endemic level of infectious bark beetles, we cannot decrease the input rate of pines. However, we should focus on the felling rate of infectious pines in order to achieve the decreased endemic level of infectious bark beetles.

Optimal Control

Sensitivity analysis urges to modify the model (1) to include some control measures for the specific counteractive action. The model including three controls and is given as: The control represents the utilization of nematocide infused into the stem of healthy trees, shows the increase in the slashing rate of infected pine trees with the goal that bark beetle could not be able to oviposit on them and use of adulticide applied for vector control, for example, aerial spraying of pesticide, is expressed by the control function . To explore the ideal level of efforts to rein the ailment, we construct the objective functional . Its purpose is to reduce the infectious trees and the cost of applied k8controls. One has where represent the positive weights. We will likely limit the number of infected trees figure while cutting down the expenditure of controls and with the aid of the above-mentioned objective functional. We shall find an optimal control ,  and such that where is the control set. The solution of this optimal control (OC) problem and derivation of necessary conditions are obtained by using the Pontryagin’s Maximum Principle[23].

Optimal Control Existence

Optimal Control existence can be demonstrated through a well known classical result: according to[24], we must inspect that the following axioms are satisfied: where and . Thus we have the following theorem:

(H1) Controls set and set of state variables are nonempty.

(H2) The admissible control set is convex and closed.

(H3) R.H.S of the state system is bounded by a linear function of the state variables and controls.

(H4) The objective functional has convex integrand on and is bounded below by , where and .

(H5) The existence of solutions for the system (5) is established by using the result given by Lukes ([25], Th 9.2.1, p 182). In this way, we confirm the above hypothesis. is accomplished because the coefficients are bounded. The boundedness of solutions shows that the set of controls fulfills . Since the system of equations is bilinear in , and solutions are bounded. Hence, R.H.S of (5) meets the criteria . The last condition is also satisfied because the integrand of the objective function is convex.

Theorem 3.

For the objective functional where subject Eq.(5) with initial conditions, there exists an optimal control so that .

The optimal solution can be attained by finding the Lagrangian as well as Hamiltonian for the system (5). Its Lagrangian is We have to set up the minimal value of the Lagrangian. For this purpose, we construct the Hamiltonian for OC problem as follows:

Let us take (), () and (), then we have

The optimality system

We apply the Pontryagin’s Maximum Principle[23] for finding the essential conditions for this optimal control. It is done as follows:

For the optimal solutionof OC problem (5), a non-zero vector function exists that satisfies the subsequent conditions. The state equation, optimality condition and adjoint equation are given as respectively The essential conditions applied to the Hamiltonian give the following result:

Theorem 4.

For optimal controls and solutions of the system (5), variables occurs having

with transversality conditions . Additionally, are expressed as

Proof. The adjoint equations and the conditions of transversality are obtained by the Hamilton . Set and differentiate the Hamiltonian with respect to . Thus we have with transversality conditions . The conditions of optimality and the property of control set give the result

Experimental results of optimal control

We now solve the OC problem numerically and show some results. Since the initial values of the state variables and are given, therefore, forward RK-4 method is applied to solve the system of equations (5). The backward RK-4 method is applied to solve the adjoint system (6) because final values are given. Regarding the control variables, an initial guess is given. Time-varying control variables have been found by putting the values of the state and adjoint variables in the system (7). We consider that the optimal campaign continues for 120 days. We present the use of different controls.

Applications of two controls

Application of nematocide and slashing of infected pines (): In this strategy, only two controls, nematocide and slashing of infected pine trees, are used to optimize the objective functional , whereas third control (aerial spraying) is set to zero. Figure 5a shows that a considerable decrease occurs in the population of infectious pines and according to Fig. 5b small decrease appears in the population of vectors. The control profile is given in Fig. 5c. From this control profile, we observe that nematocide has to apply continuously for 26 days with great animation and then we can decrease this effort. From Fig. 5c, it can also be seen that we should make a maximum effort to slash the infected pines for 18 days and then slow down this control up to 120 days.

Figure 5

Use of controls and .

Use of nematocide () and aerial spraying (): This control policy is based on the use of two controls, nematocide and aerial spraying. Figure 6a shows that infectious pines slightly increases first, and then began to decrease for 40 days. After this period the number of infectious pines increases up to 120 days. Figure 6b shows that infectious vectors decrease to zero after 7 days. Control profile is given in Fig. 6c shows that we have to use nematocide and aerial spraying for 120 days continuously. However, Fig. 6a,b show the non-suitable strategy.

Figure 6

Use of controls and .

Slashing of infected pines () and aerial spraying (): In this policy, slashing of infected pines and aerial spraying are used to optimize the objective functional . Fig. 7a shows that the application of this strategy minimizes infectious pines first but after 80 days they again start to increase. Infectious vectors decrease to zero as shown in Fig. 7b. The control profile showing in Fig. 7c informs that we have to apply this control strategy for 118 days with great animation. However, Fig. 7a shows that this is not an effective control strategy because the infectious pines increase ultimately.

Figure 7

Use of controls and .

Application of three controls

In this strategy, all control measures are used to optimize the objective functional . Figure 8a shows that infectious pines decrease to zero after 18 days. The number of infected bark beetles also decreases to zero after 7 days as shown in Fig. 8b. This is the best control strategy because from Fig. 8c, we see that aerial spraying should apply for 5 days and nematocide and slashing of infected pines should be applied for 15 days with full effort. Since the period of application of these controls with great animation is short, therefore, the cost of these applied controls is reasonable.

Figure 8

Use of three controls.

Conclusions

The deterministic model of pine wilt disease has been rigorously analyzed. In this article, the global behavior of equilibria has been discussed on the basis of basic reproduction. Global stability of disease-free equilibrium has been proved by constructing a suitable Lyapunov function whereas the global behavior of endemic equilibrium has been investigated by using the graph-theoretic approach. Two features have been discussed in the sensitivity analysis. Firstly, it has been observed that how the endemic level of infectious hosts and vectors is effected by varying the values of parameters. Secondly, the corresponding change in the state variable has been measured when a parameter changes. When we calculate the percentage difference of the values of and corresponding to the percentage difference of the parameter values, it is observed that the most sensitive parameter is for as well as for . The calculation of sensitivity indices of has shown that the most sensitive parameter for is the mortality rate of bark beetles. The extreme values of the sensitivity indices of and are and corresponding to the parameter values and , respectively. These sensitivity indices show that by increasing the mortality rate of bark beetles by , would decrease almost . Similarly, index shows that the increment of the slashing of infected pines by results reduction in the endemic level of by . The sensitivity analysis has shown us to design optimal control problem having three controls, namely, the use of nematocide implanted into the trunk of healthy trees, the increase in felling rate of affected pine trees so that bark beetle could not succeed to oviposit on them and the level of adulticide proposed for vector control such as aerial spraying of pesticide.