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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{R}}}_{{\boldsymbol{o}}}$$\end{document}Ro. The disease-free equilibrium is globally asymptotically stable for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{R}}}_{{\boldsymbol{o}}}{\boldsymbol{ < }}{\bf{1}}$$\end{document}Ro<1 and in such a case, the endemic equilibrium does not exist. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{R}}}_{{\boldsymbol{o}}}$$\end{document}Ro 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.


Bio-inspired Analytical Heuristics to Study pine Wilt Disease Model
Muhammad ozair 1 , takasar Hussain 1* , Aziz Ullah Awan 2 , Adnan Aslam 3 , Riaz Ahmad Khan 4 , farhad Ali 5 & fatima tasneem 1 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 R o . The disease-free equilibrium is globally asymptotically stable for R < 1 o and in such a case, the endemic equilibrium does not exist. If R o 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.
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 1905. During the s 1930 the annual loss of pines increased from m 30,000 3 to m 200,000 3 , and during 1940s destructions were estimated at volumes of m 400,000 3 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 m 1,230,000 3 in 1948. Considerable efforts lessened the annual loss to m 400,000 3 during the 1950s. Elimination of dead trees through tumbling and blazing was used as the basic control method at that time. In the s 1970 , annual losses exceeded 1 million m 3 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, .
2 4 million m 3 , was recorded in 1979 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][13][14][15][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: • 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 R o , infected host I * and infected beetles Y . * • The aim of this paper is to study new features of the models: 1. To show the global behavior of equilibria. 2. Identification of sensitive parameters for reproduction number and endemic levels of infectious classes. 3. 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 S t ( ), infected host pine trees I t ( ), susceptible vector beetles X t ( ) and infected beetles Y t ( ). 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 t is represented by , and whole vector adult beetles are indicated by . Thus, = + S I  and  = + . X Y Let a be the constant recruitment rate of Monochamus alternatus, b 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 YS β . 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 kIX. Monochamus alternatus is infected through mating and the contact process is denoted by XY 1 β . 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.
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, ≥ t 0. The total population of both the species satisfies the equations: 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.
www.nature.com/scientificreports www.nature.com/scientificreports/ where  F i shows the rate of new infections in ith disease class and  V i shows the transition terms, for example recovery and death in the ith disease class 19 .  R is acquired by choosing the dominant eigenvalue (spectral radius) of FV The infected compartments are I and . Y Therefore, Thus,  www.nature.com/scientificreports www.nature.com/scientificreports/ For  > R 1, the system (1) has a unique endemic equilibrium (EE) E * = S I X Y ( * , * , * , * ), where 1 and Y * is uniquely calculated from the equation Noting that > A 0 while the sign of C coincides with that of (1) has unique EE given as E S I X Y * ( * , * , * , * ) = .

Stability of Equilibria
We investigate the global behavior of the system (1) around DFE. For this, we construct the Lyapunov function.
Taking time derivative of L, we have □ We now will show that the unique EE is GAS using the graph-theoretic approach. For further terminologies, the readers are referred to 20 .  Differentiating and using the inequality * * ln * ln * ln * * * * * * ln * * ln * * * * : A weighted digraph is constructed with four vertices and six arcs, which is shown in Fig. 2: It is easy to see that there are three cycles and along each cycle G G G 0 31 43 14 Thus, www.nature.com/scientificreports www.nature.com/scientificreports/ is the Lyapunov function of the system (1). Using this Lyapunov function and LaSalle's invariance principle 21 , it follows that E * 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 a b k , , , , 1 β β are directly related to I * and inversely related to α δ µ , , . It can also be observed that considerable change occurs in the value of I * by an increment or reduction of β α a, , and δ. However, when we calculate the percentage difference of the values of I * 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 416195. 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 I * corresponding to all parameters is given in Table 1.

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 β 1 and µ. But the calculation of percentage difference of the values of Y * corresponding to percentage difference of the parametric values shows that the most sensitive parameter is again β. The value of this difference is 8890. Percent increment or reduction in the value of Y * corresponding to all parameters is given in Table 2.
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.
The above sensitivity index has no obvious mathematical structure so that we could decide whether R  is an increasing function of a or decreasing function. However, when we use the numerical values of the parameters given in Table 2, we see that the sensitivity index of  R of the parameter a is . 0 51. It means that when we increase the value of a by 10%, the value of  R increases by almost 5%. Looking at the sensitivity indices of  R 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 1 012 − . . Negative sign shows that  R is a decreasing function of µ. An increment in the value of µ by 10%, the reduction in  R occurs by 10%.
Sensitivity indices of I * and Y * . The endemic level of infectious pines given in (3) is represented in the form of endemic level of infectious bark beetles Y * 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 I * and Y * .
The calculated values of sensitivity indices of I * 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 10%, I * would decrease almost 16%. The endemic level of infectious hosts is also an increasing function of β. By decreasing the value of β by 10%, I * would decrease almost 14%. Sensitivity indices of the endemic level of infectious bark beetles Y * , given in Table 4, show that the parameters b and δ play a vital role in the enhancement or reduction in the values of Y * . Y * is an increasing function of b and a decreasing function of δ. Increment in the value of δ by 10%, reduction in the value of Y * occurs by almost 10%. Contrary to this, by decreasing the value of b 10%, Y * decreases by almost 10%. 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 u u , 1 2 and u 3 is given as: The control u t ( ) 1 represents the utilization of nematocide infused into the stem of healthy trees, u t ( ) 2 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 u t ( ) 3 . To explore the ideal level of efforts to rein the ailment, we construct the objective functional J. Its purpose is to reduce the infectious trees and the cost of applied k8controls. One has 2 represent the positive weights. We will likely limit the number of infected trees figure while cutting down the expenditure of controls u t u t ( ), ( ) 1 2 and u t ( ) 3 with the aid of the above-mentioned objective functional. We shall find an optimal control u * 1 , u * 2 and u * 3 such that J u u u min J u u u u u u U ( * , * , * ) { ( , , ), ( , , ) }, 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:      . The optimal solution can be attained by finding the Lagrangian as well as Hamiltonian for the system (5). Its Lagrangian is    www.nature.com/scientificreports www.nature.com/scientificreports/ Let us take  = (S I X Y , , , ), λ = ( , , , λ λ λ λ ) and U = (u u u , , 3 ), then we have 1  2  1   3  2  1  1   2  1  0 2   3  3  1  1  1 3   4  1  1  1 3 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 solution u u u ( * , * , * ) 1 2 3 of OC problem (5), a non-zero vector functionλ t ( ) exists that satisfies the subsequent conditions. The state equation, optimality condition and adjoint equation are given as respectively dx dt H t u u u t ( ( , * , * , * , ( ))), 1 2 3 ( ( , * , * , * , ( ))), The essential conditions applied to the Hamiltonian H give the following result: Theorem 4. For optimal controls u u u * , * , * 1 2 3 and solutions S I X Y * , * , * , * of the system (5), variables λ λ … , , 1 4 occurs having  . The conditions of optimality and the property of control set U give the result (2020) 10:3534 | https://doi.org/10.1038/s41598-020-60088-1 www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ Experimental results of optimal control. We now solve the OC problem numerically and show some results. Since the initial values of the state variables S I X , , and Y 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 ≠ u ( 0 ) 1 and slashing of infected pines ( ≠ u 0 2 ): In this strategy, only two controls, nematocide and slashing of infected pine trees, are used to optimize the objective functional J, 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.
Use of nematocide ( ≠ u 0 1 ) and aerial spraying ( ≠ u 0 3 ): 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.
Slashing of infected pines (u 0 2 ≠ ) and aerial spraying (u 0 3 ≠ ): In this policy, slashing of infected pines and aerial spraying are used to optimize the objective functional J. 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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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.
Application of three controls. In this strategy, all control measures are used to optimize the objective functional J. 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. 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 I * and Y * corresponding to the percentage difference of the parameter values, it is observed that the most sensitive parameter is β for I * as well as for Y * . The calculation of sensitivity indices of R 0 has shown that the most sensitive parameter for R 0 is the mortality rate of bark beetles. The extreme values of the sensitivity indices of I * and Y * are 1 60 − . and − . 1 002 corresponding to the parameter values µ and δ, respectively. These sensitivity indices show that by increasing the mortality rate of bark beetles by 10%, I * would decrease almost 16%. Similarly, index − .
1 002 shows that the increment of the slashing of infected pines by 10% results reduction in the endemic level of Y * by 10%. 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.