Dynamical features of pine wilt disease through stability, sensitivity and optimal control

This work investigates the dissemination mechanism of pine wilt disease. The basic reproduction number is computed explicitly, and an ultimate invariable level of contagious hosts and vectors, without and with disease, is discussed by using this number. Highly effective techniques, Lyapunov functional and graph theoretic, are utilised to obtain the ultimate constant level of the whole population. The idea of complete disease eradication and reduction of endemic level is explored through the utilisation of two efficient methods. Using sensitivity analysis approach, necessary control measures are suggested to overcome the disease. Using the literature data, the robustness of control strategies is shown graphically.


Introduction
Pine trees are evergreen conifers. These are some of the famous plants around the world because of trade of timber (Fig. 1). Manufacturing of turpentine, pulp rosin and paper explores the value of these trees. Some species have seeds (called pine nuts) that are used for cooking and baking.
Pine wilt disease (PWD), the disease of pines, is so nasty that after getting the disease the infected tree becomes dead within a few weeks to a few months. The usual symptom is needles' discolouration. Their colour first changes to greyish green and finally it turns out to be totally brown. Infected pine trees cease oleoresin after succumbing to the disease (Fig. 2). It was reported that PWD first occurred in Japan in the early twentieth century [1]. It was observed in the Japanese prefecture Hyogo in 1921. In 1930, this disease appeared in distant places of Kyushu. This malady ultimately spread to the shore of Seto Inland Sea [2].
Bursaphelenchus xylophilus, which has an international common name "pine wood nematode (PWN)", is the casual pathogen of this malady (Fig. 3). This nematode exists in the upper portions of pine trees. Disturbed water conditions, through small air pockets, are created by Bursaphelenchus xylophilus in the xylem of pine tree, which prevents water  Healthy and infected pine trees movement. As a result, the tree becomes dead soon [3]. After more than 60 years, this pathogen has been discovered. It was first reported as the casual agent in Japanese native pines in 1971 [4]. Numerous genetic varieties of PWN were studied, and it was concluded that it came to Japan from the United States. From this place, it spread to China, Portugal, Korea and Spain [5].
PWN is not capable of moving itself from one tree to another. Bark beetles (pine sawyer) serve as the transporting agent (Fig. 4). Monochamus alternatus, commonly known as the Japanese pine sawyer, is a species of beetles that belong to the family of Cerambycidae. These are the vectors for the transmission of PWN. When nematode approaches the dispersal stage, it attaches to the insects' respiratory system, and when adult beetles feed on healthy pine trees, this Bursaphelenchus xylophilus is transferred into the tree.
Transmission mechanism of PWD occurs in three ways: First, when infected pine sawyers fly and start maturation feeding. These infected pine sawyers feed on the twigs of a pine tree, and during this process the nematode is transferred into the tree. It is called the primary transmission. Second transmission occurs during the egg laying activities of the pine sawyer on a freshly cut, dying or dead pine tree [6]. Third transmission, which may be  called horizontal transmission, occurs during the mating process of the pine sawyer [7]. The PWD cycle can be explained through Fig. 5.
Modelling is the branch of science in which we use inspirational approaches of applied mathematics equipped with inspiration. In order to analyse the dynamics of an ecosystem and quantitative prediction, mathematical treatment is an indispensable tool. Mathematical modelling has become an essential tool which enables one to comprehend the causes, patterns and effects of the disease. Mathematicians, through investigations and prediction of an epidemic, provide an invaluable assistance to health care agencies who attempt and allocate resources to control an infectious disease or epidemic in the community.
Various mathematical models have been proposed for investigation of Bursaphelenchus xylophilus transmission. A simple deterministic model for wilt disease of pines has been presented by Yoshimura et al. [8]. The model demonstrates disease epidemic through three important factors, namely, initial pine density, beetles' density and the capability of elimination of the pine sawyer. Takasu et al. [9] presented a mathematical model based on integro-change equation to explain the range expansion of PWD. Mechanistic contact of beetles and pine trees is supposed in [10]. In this article, it was shown that the Alee effect appears from the reality that the beetles have contact with the pine tree at least twice to reproduce successfully. Shi and Song [11] presented the mathematical model of PWD and analysed it qualitatively by considering two types of transmissions: maturation feeding of infectious beetles on healthy pine trees and oviposition of the pine sawyer on dying or dead trees due to the pine wilt infection. However, another transmission that happens through copulation of pine sawyers, as suggested by Togashi and Arakawa, has not been discussed in their work [12]. Other mathematical models based on ordinary and fractional derivatives are comprehensively explained in [13][14][15][16][17][18][19][20][21][22][23][24][25].
This work is based on the reconsideration of model studied in [26] with the addition of the following distinguished aspects: • The exploitation rate of tainted pines is higher than all other compartments containing pine trees. • Transmission of nematode in vector population occurs during copulation. • Formulation of optimal control problem through sensitivity analysis of the threshold parameter, pine trees having nematodes and beetles having nematodes. This work has three fold intentions: 1. To reveal the stability of equilibria. 2. To recognise those factors, for basic reproduction number and infectious compartments, that have high impact. 3. To formulate the optimal control techniques by making use of sensitivity analysis.

Model formulation and flow diagram
We construct the mathematical model by considering two populations, host and vector. The host population consists of susceptible, exposed, asymptomatic and infectious carriers represented by S h (t), E h (t), A h (t) and I h (t), respectively. There are two compartments of vector population, susceptible and infectious vectors. These are denoted by S v (t) and I v (t), respectively. Transmission of nematode among vectors and the pine sawyer arises during maturation feeding of the pine sawyer and oviposition of these vectors. However, the pine sawyer may also attain infection through copulation. Suppose that N h and N v denote the whole host and vector population at any time t, respectively. So, N h = S h + E h + A h + I h and Suppose that the planting rate of healthy pines is λ h and the recruitment rate of a mature pine sawyer is λ v . The exploitation rate of pine trees (except infectious pines, which is denoted by δ) is μ h and the fatality rate of vectors is μ v . The transferring rate of nematodes through an infectious pine sawyer during maturation feeding is α. The mean number of daily contacts of fully developed infectious vectors with the susceptible pines is φ. Suppose that the pine sawyer makes a mean number of ψ contacts with the tree for reproducing themselves and β is the transmission rate of nematode through oviposition.
The parameter θ shows the probability that a dead tree was not infected with the nematode and ceases oleoresin due to natural causes. Let γ denote the emergence rate of an infectious pine sawyer from an expired tree. Exposed trees move to either the infectious or the asymptomatic carrier class. These trees move with the per capita rate k, the fraction ρ (0 < ρ < 1) of them move to the infectious compartment.
We design a system of nonlinear ODEs, based on the above stated assumptions, as follows: Figure 6 shows the flow sheet of the model. All the above variables will be nonnegative, because we are dealing with the population, and the whole populace of both species has the following form of differential equations: The global attractor for (1) exists in the set , and the system is dissipative.

Equilibria and basic reproduction number
Two kinds of equilibria, disease free equilibrium (DFE) and endemic equilibrium (EE), exist for the above presented model. We denote DFE by E o and EE by E * , respectively. Disease free equilibrium can be easily computed by putting all the diseased classes and right-hand sides of each class of system (1) equal to zero. We get . The whole dynamics of the disease is based on the threshold parameter known as the basic reproduction number. It identifies the status of the disease, i.e. whether the disease is spreading in the community or not. It is termed "the mean number of secondary infections brought about by an infectious individual during its complete course of infection". We find the expression of basic reproduction number with the help of the next generation matrix method.
Let x = (y 1 , y 2 , . . . , y n ) T ∈ R n and y = (z 1 , z 2 , . . . , z m ) T ∈ R m denote the population in the infected and non-infected classes, respectively.F = (F 1 ,F 2 , . . . ,F n ) T andV = (V 1 ,V 2 , . . . ,V n ) T , whereF l shows a new infection rate in the lth infected compartment,V l symbolises the transition terms, for example, recovery and fatality in the lth infected class [27]. The dominant eigenvalue of the matrixFV -1 is called the basic reproduction number. Contaminated compartments are E h , A h , I h and I v . Therefore Thus, which gives As R • is the dominant eigenvalue ofFV -1 , hence is the endemic equilibrium for the reproduction number greater than one and and we can get the value I * v by solving the following second degree equation: The changes in signs of coefficients A, B and C tell us about the zeros of the above quadratic equation. We can see that A is always positive and C < 0 if and only if R • > 1. Paying no attention to the sign of B, we conclude that there is a unique positive value of I * v because there is only one sign change for R • > 1. Hence the unique EE

Global behaviour of equilibria
We explore the global behaviour of equilibria by using Lyapunov functional theory and graph theory. The comprehensive procedure of global stability using graph theory has been omitted here. We recommend the keen reader to go through [28]. Proof Suppose that we have the Lyapunov function It can be easily observed that L ≤ 0 when R • ≤ 1 and L = 0 for E h = I h = I v = 0. Thus, the maximal dense unvarying set, for a system of equations (1), contains a single component which is E • . With the help of Lasalle's invariance principle [29], we can say that E o is globally asymptotically stable in .

Theorem 2
The endemic equilibrium point (E * ) exhibits global asymptotic stability in for R • > 1.
Proof Suppose that Taking the rate of change of above D s w.r.t t and by making use of the result 1 + ln xx ≤ 0 for positive x, the derivatives can be expressed as follows: and similarly, Figure 7 shows a weighted digraph with ten arcs and six vertices. One can observe that Solving the above equations, we get is the Lyapunov function for system (1). Utilising LaSalle's invariance principle and the Lyapunov functionD [29], E * is GAS in the set ( ).

Application of the model
Now the given model (1) is applied on real data containing infectious pines in order to check the robustness of the model. Numerical values obtained through the calibration of model outcomes with the real number of spoiled pines will be used for the identification of crucial factors for the transmission of the disease as well as to design the control strategies to reduce the malady from tree community. We used the data of the real number of spoiled pines that were infected in Korea from 2010 to 2019 and were given in [30]. Explanation of the parameters used in model (1) and their estimated values are shown in Table 1. It can be seen from Fig. 8 that the model outcomes are in great ascension with the real figure of pines that become infected in Korea during the period of 10 years.

Sensitivity analysis
The primary concern about PWD is to ponder the potential of this malady in the community. Sensitivity analysis will be performed to comprehend the main factors which play the key role in dissemination or persistence of the disease in the community. The complete eradication of PWD from tree community is an ideal situation. It means that we cannot make R • less than unity absolutely. However, attention can be paid to those parameters that cause high infection level. They can be identified by observing the variation in a constant level of the infectious compartments by swapping parameter values. We can focus specially on those parameters due to which huge fluctuation in the endemic state of contagious pine trees and pine sawyers occurs.

Sensitivity indices of R 0
For brief introduction about the calculation of sensitivity indices, one can consult [31]. We use its definition here in the form of partial derivative. It is given as follows: Using the above definition, mathematical expressions to calculate the sensitivity indices of R 0 are obtained. For instance, R 0 has the sensitivity index with respect to δ and μ v as follows: . Mathematical expressions given in the above two equations do not provide any information on how much difference happens in the value of basic reproduction number when values of δ and μ v are changed. However, after utilisation of fitted values of all the parameters mentioned in Table 1, we examine that R 0 has sensitivity index -0.276506 with respect to δ. It means that deceleration in the basic reproduction number by almost 3% occurs as a result of enhancement of falling rate of infectious pines by 10%. If we have a look at all the sensitivity indices of the threshold parameter (R 0 ) given in Table 2, we inspect that the most responsive parameter for R 0 is the fatality rate of pine sawyer. It has value -1.44699. This value indicates that the reproduction number is going to decelerate by almost 15% if we amplify the fatality rate of pine sawyer by 10%.

Parameters' influence on infectious pines and pine sawyer
The most important parameters by which R • can be made less than unity have been calculated in the above subsection. It seems an impractical task to get rid of PWD completely from pine tree community. However, the factors that help to reduce the infection may be identified. This objective can be obtained by the calculation of a relative change in the number of infectious hosts and vectors with the amplification of difference in parameter values.

Reduction of infectious pines
We can see from Figures 9-22 that a constant level of infectious pine trees is extremely influenced by the parameters α, β, δ, λ h and μ h . However, the computation of percentage difference of I * h shows that the most sensitive parameter is the probability of susceptible pines not to be infectious (θ ). Its value is 8765.984. It means that we should cut those trees that stop oleoresin as soon as possible. Percent enhancement in the number of infectious pines with respect to all parameters given in the model has been calculated in Table 3. Figures 23-36 show that if we increase the fatality rate of pine sawyers, exploitation rate of infectious pines and decrease the transmission probability through maturation feeding, the number of infectious pine sawyers can be reduced significantly. We also observe from Table 4 that the removal of adult pine sawyers is essential to decreasing the infection. This observation tells us that we should focus on the minimum possible emergence of pine sawyers. It can be made achievable through the use of pesticide spray into newly cut timber or bark of dead trees.

Optimal control
The observations of former section indicates that we should target the appearance of mature pine sawyers, ousting the spoiled pines rapidly and demise of pine sawyers. Now, we have to redesign model (1) to get the results of some control strategies to overcome this disease. Model (1) gets the following shape after the inclusion of these controls: The control u 1 (t) expresses the defensive measures like the application of nematocide inserted into the cane of well fit pines and vaccination for the reduction of the strength of infection. Control u 2 (t) means ousting the spoiled pines quickly so that pine sawyers have no option to lay eggs on them. To keep the number of vectors in control, we are using u 3 as control which is the application of adulticide. Objective functional J has been formulated in the way to cut down this malady. The main objective of its construction is to minimise the spoiled pines and also to reduce the expenses of implemented controls u 1 , u 2 and u 3 . One has where M 1 , M 2 are the weights having positive values. The main objective of applying the controls is to make the effort to retrench the number of contaminated pines and vectors keeping in mind that the cost of applied controls should be minimum. It is done with the help of the above objective functional. Controls u * 1 , u * 2 and u * 3 can be found in such a way thatÕ The set U is called control set. It is defined as We have a very useful result in the form of a principle called Pontryagin's maximum principle and given in [32]. We will use it to get the mandatory states and the solution of (3).

Persistence of optimal control
A very powerful procedure is applied for existence and analysis of the optimal control. According to [33], the necessary hypotheses are as follows: (H 1 ) A set of control and corresponding described unknowns has at least one element.
(H 2 ) U must fulfil the two properties, convexity and closeness.
(H 3 ) The right-hand side of the described system is bounded by a linear function with respect to state and control.
(H 4 ) The integrand of the objective functional is convex on U and has lower where c 1 , c 2 > 0 and τ > 1. The existence of solutions of system (3) can be proved by utilising the result given in ( [34], Th. 9.2.1, p. 182). In this manner, we validate the above axioms. (H 1 ) is fulfilled because of finiteness of coefficients. The finiteness of solutions reflects that U attains the 2nd hypothesis. Note that system (3) is bilinear in u 1 , u 2 , u 3 . Because of bilinearity of controls and boundedness of solutions, R. H. S of (3) satisfy the 3rd hypothesis H 3 . As the integrand of the integral J is convex, we have the assurance of H 4 .
where D 1 , D 2 , D 3 , M 1 , M 2 , c 1 , c 2 > 0 and τ is greater than 1, hence we come up with the conclusion given below.

Theorem 3 For the objective functionalÕ
The optimal solution of system (3) can be obtained through the Lagrangian and Hamiltonian. We consider the Lagrangian for (3) as follows: We have to find the least possible value of the Lagrangian L given above. To do this, we introduce the Hamiltonian H as follows: Let us take X = (S h , E h , A h , I h , S v , I v ), μ = (μ 1 , μ 2 , μ 3 , μ 4 , μ 5 , μ 6 ) and U = (u 1 , u 2 , u 3 ), then we have

The optimality system
We have made use of Pontryagin's maximum principle, contained in [35], for finding necessary states for the control. This has been done in the following way: with boundary conditions μ 1 (T) = μ 2 (T) = · · · = μ 6 (T) = 0. The property of the control set U and optimality conditions give the result u * 1 = max min 1, One can also draw the conclusions, when control strategies have been used to limit the spread of PWD, numerically. It can be seen from Fig. 37 that the strength of healthy trees goes on increasing due to these controls. It seems true because after the application of controls all infected pine trees will be removed, as a result we have more tress in the susceptible class. Figure 38 shows the substantial depletion in exposed pines when we apply these controls. Figure 39 depicts the behaviour of asymptotic carrier trees when the control is applied, one can see a little decrease here. From Fig. 40, we can clearly see that there is a substantial decrease in spoiled pine trees due to these controls.
To get the information about the number of susceptible pine sawyers by the application of controls, we see Fig. 41; it shows no significant depletion. On the other hand, Fig. 42 shows the outstanding result of the use of these controls on a class of infected vectors: first they are approaching zero for ten days, then they increase moderately. This behaviour of

Conclusions
PWD has been analysed in this work through a deterministic model. First, we studied the model qualitatively and calculated the explicit formula for the basic reproduction number. Constant levels (infection free and disease present) of all the compartments have been calculated. We have also shown that if at any stage we are able to reduce the basic reproduction number below unity, then the disease will be completely eradicated from the community. If the policy makers of any country fail to get the basic reproduction number less than one, then we have identified the constant level of vector and host population. The actual data of infectious Korean pines from 2010 to 2010 have been used, and the model has been simulated. The values of parameters for which the model outcomes had the best match with the actual data have been identified. We performed sensitivity analysis so that an effective control strategy could be designed to remove the disease from pine tree community completely or to obtain the least figure of infectious pines. We computed the sensitivity indices of the threshold parameter by using the fitted values of the parameters. The parameters that play a significant role in reducing the threshold parameter to less than one have been identified. The effects of parameters on the constant level of contaminated trees and vectors have been observed if we fail to get the value of the basic reproduction number below unity. We varied the values of the parameters and noticed the correspond- To get the most effective parameters for the endemic position of contaminated trees and pine sawyers, we used the least value to normalise all the percentage changes.
For successful strategies to overcome the disease, we have done sensitivity analysis. We modified the model by inserting the control functions and calculated the optimal values of these functions. We simulated the model by using the fitted values of the parameters used in the model, and the efficacy of implemented strategies has been examined. It can be seen that control strategies are very effective in reducing the number of infectious pines and vectors. The same idea can be applied to any community whose real data is available. The model can also be modified by taking the fractional derivative instead of ordinary deriva-tives to gain a better insight about PWD. It will be included in the future research, and comparison of outcomes through ordinary and fractional derivatives will be performed.