Stability analysis of pine wilt disease model by periodic use of insecticides

ABSTRACT This work is related to qualitative behaviour of an epidemic model of pine wilt disease. More precisely, we proved that the reproductive number has sharp threshold properties. It has been shown that how vector population can be reduced by the periodic use of insecticides. Numerical simulations show that epidemic level of infected vectors becomes independent of saturation level by including the transmission through mating.


Introduction
Pine wilt disease (PWD), caused by the pinewood nematode (PWN) Bursaphelenchus xylophilus is a dramatic disease because it usually kills affected trees within a few weeks to a few months. It causes major losses in coniferous forests in eastern Asia, in particular Japan, China, South Korea and Western Europe, e.g. Portugal. As such, the PWN is among the most important pests included in the quarantine lists of many countries around the world [19].
Although PWN supposedly originated in North America, the first occurrence of PWD was recorded in Japan in 1905. Since then, this disease spread to China, Korea, Mexico and Portugal, in the process becoming a major ecological catastrophe with serious economic losses to the pine forest industry. Pine wilt particularly kills Scots pine. However, some other pine species such as Austrian (Pinus nigra), jack (P. banksiana), mogo (P. Mugo), red (P. Resinosa) pines are occasionally killed by pine wilt.
In living pine trees, PWN rarely invade the xylem conduits but distribute in the cortical and xylem resin canals. They can then move through the resin canals to reach almost all parts of the stem and branches fairly rapidly. In infected pine trees, both virulent PWN and avirulent PWN cause some localized embolisms in the xylem. The virulent PWN then significantly multiply in the pine stem, and water potential of the leaf abruptly decreases inducing death [18].
The PWN is not able to migrate from tree-to-tree by itself. The main vector insect is a sawyer beetle (Monochamus alternatus). The sawyer beetle is infected with the PWN while it transforms from pupa to adult in the wood of the damaged pine tree. Once a beetle emerges and it flies to the healthy crown of pine tree for sufficient feeding. A beetle can holds about 270,000 worms of nematodes [7].
There is no treatment for pine wilt once a susceptible tree gets to be plagued with PWN. Consequently, the standard technique is one of prevention. The best system is to abstain from planting non-native pine species in areas in ranges where the normal summer temperature surpasses 20 o C (68 o F). In the areas of earlier existence of non-native pines consideration ought to be taken to water and fertilize them particularly during times of dry season. We can also take a set of techniques for specialized treatment of timber to increase the mortality of infected long horned beetle; removal of dying branches on affected trees or dead and dying trees and fumigation of nematode-infested lumber.
Mathematical modeling is very effective tool to explain the process of transmission of a disease. In this way, we can determine those factor that have significant impact on the spread of the disease and we can design different control strategies to restraint the spread of infection. Heylman et al. [5] developed a mathematical model that allows simulation of oxygen distribution in a bone defect as a tool to explore the likely effects of local changes in cell concentration, defect size or geometry, local oxygen delivery with oxygen-generating biomaterials (OGBs), and changes in the rate of oxygen consumption by cells within a defect. Pimentel et al. [17] discussed a model of normality based vital-sign data from patients who had a normal recovery was constructed using a kernel density estimate and tested with abnormal data from patients who deteriorated sufficiently to be re-admitted to the intensive care unit.
Many mathematical models have been proposed and analysed for vector -borne disease [14,15,20,21,24,26,28]. Numerous ventures have been made to develop realistic mathematical models for exploring the transmission dynamics of the PWD. Yoshimura et al. [27] demonstrated that there is a minimum pine density below which the disease dependably comes up short in incursion and the minimum pine density increases disproportionately with the increase in the abolition rate. Takasu et al. [23] investigated the explored the reliance of the parasites rate of extent extension on the destruction rate of the beetle, the initial pine density, and the beetle dispersal capacity. Mechanistic interactions are considered in [22], and it is proved that the Allee impact can emerge in the beetle dynamics because of the necessity for beetles to contact pine trees at least twice to replicate successfully. The incubation period after which a tree contacted by a first beetle becomes ready for beetle oviposition by later beetles is crucial for the emergence of this Allee effect. Lee and Lashari [8] discussed a mathematical model with nonlinear incidence rates to describe the host-vector interaction between pines and pine sawyers carrying nematode by considering the role of incubation period during disease transmission. However, they did not include direct mode of transmission. Recently, Ozair et al. [16] developed a PWD model considering a direct mode of transmission with standard incidence rate. They did not consider those pine trees that has been infected by the nematode but still sustains the ability for oleoresin exudation.
In this paper, we extended the mathematical model discussed in [13] by including the role of incubation period during disease transmission. The objective of this paper is twofold. The first is to carry out the stability analysis and second to observe that how the disease can be eradicated by periodic use of insecticides. The rest of the paper is organized as follows. In second section, the description of the extended mathematical model is presented. The stability of disease-free equilibrium and the stability of endemic equilibrium are investigated in third and fourth sections, respectively. In Section 5, the global stability of endemic equilibrium is proved using the geometric approach method. Discussions and simulations are done in the last section.

Model formulation
In this section, we formulate a continuous mathematical model for the transmission of vector-borne disease according to the basic rules of mathematical modelling in epidemiology. The common transmission of nematodes among pine trees and bark beetles occur during maturation feeding of infected vectors. The pine sawyers have PWN when it emerges from infected pine trees. However, the beetles may also get infected directly during mating. The following assumptions are made about the coexistence of the populations in the same environment.
Pine Trees: The trees that have a potential to be infected by the nematode and can exude oleoresin which acts as a physical barrier to beetle oviposition and beetles cannot oviposit on them.
Exposed E h (t): Those infected trees that still sustain the ability for oleoresin exudation.
The oleoresin exudation ability have been lost by these trees and also beetles can oviposit on them.
Bark Beetles: Susceptible S v (t): Adult beetles which do not carry PWN. Infected I v (t): The number of infective adult beetles which carry PWN.
The total pine population is N h (t) = S h (t) + E h (t) + I h (t) and total vector population is N v (t) = S v (t) + I v (t). We assume that pine and beetle populations have constant recruitment rates h and v , respectively and mortality rates given by μ h and μ v , respectively. We assume that the transmission rate per contact during maturation feeding is δ 1 and the rate at which infected vectors transmit the nematode via oviposition is δ 2 . We also assume that there is a probability for susceptible pines to die of natural causes and cease oleoresin exudation without being infected by the nematode which is denoted by θ . The rate at which adult beetles carry the PWN when they emerge from dead trees is β 2 . The transmission between susceptible pines and infected vectors occurs when infected beetles lay eggs on those dead pines that die of natural causes or through the maturation feeding of infected vectors; the incidence terms for these transmissions are δ 2 θ S h I v /(1 + α 1 I v ) and δ 2 θS h I v /(1 + α 1 I v ), respectively. The transmission between susceptible vectors and infected hosts occurs when adult beetles emerge from dead pine trees. This transmission is denoted by β 2 S v I h /(1 + α 2 I h ). The parameters α 1 and α 2 determine the level at which the infection is saturated. The beetles transmit nematodes directly through mating. Then incidence term for this transmission is β 3 S v I v , where β 3 is the transmission rate among beetles during mating. The parameter β 1 denotes the transfer rates between the exposed and the infectious.
According to assumptions above, the model is given by the following system of ODEs: The total population sizes of both the species satisfy the equations: We can easily see that the total population sizes for both host and vector populations are asymptotically constant such as In our model, we can assume without any loss of generality that Thus, we study the following system of differential equations: The invariant region for the system (4) is where R 3 + denotes the nonnegative cone of R 3 including its lower dimensional faces.

Disease-free equilibrium and its stability
Direct calculations show that the system (4) always has the disease-free equilibrium point given by E 0 (0, 0, 0). The dynamics of the disease are described by the quantity . R 0 is the critical threshold of model (4) that is called the basic reproduction number in the epidemic model. Using Theorem 2 in [3], at first, the following results are established.
Proof: We linearize the system (4) around the disease-free equilibrium E 0 . The matrix of the linearization at E 0 is given by The characteristic equation of this matrix is given by det(λI − J(E 0 )), where I is the 3 × 3 unit matrix. Expanding the determinant into a characteristic equation, we obtain the following equation: where These three eigenvalues have negative real part if they satisfy the Routh-Hurwitz criteria [1], such that a i > 0 for i = 1,2,3, with a 1 a 2 > a 3 . If R 0 < 1, then According to the Routh-Hurwitz criteria, the disease-free equilibrium E 0 of the model (4) is locally asymptotically stable. Now, we study the global behaviour of the disease-free equilibrium for system (4).
Proof: We construct the following Lyapunov function: Its derivative along the solutions to the system (4) is  [6], then it implies that E 0 is globally asymptotically stable in .

Endemic equilibrium and its stability
Here, we study the existence and stability of the endemic equilibrium points. By straightforward computation, if R 0 > 1 then the host-vector model system (4) has a unique endemic equilibrium given by and I * v is the root of the following equation: with From Equation (11), we see that R 0 > 1 if and only if C < 0. Since A > 0, Equation (10) has a unique positive root in feasible region. If R 0 < 1, then C > 0. Also, it can be easily seen that B > 0 for R 0 < 1. Thus, by considering the shape of the graph of Equation (10) and noting that C > 0, we have that there will be zero positive endemic equilibrium in this case. Therefore, we can conclude that if R 0 < 1, Equation (10) has no positive root in the feasible region. If R 0 > 1, (9) has a unique positive root in the feasible region. This result is summarized below.  (9) and (10).
In order to investigate the stability of the endemic equilibrium, the additive compound matrices approach as in [9,12] is used. We will linearize system (4) about an endemic equilibrium E * and get the following Jacobian matrix: where, From the Jacobian matrix J(E * ), the second additive compound matrix is given by The following lemma stated and proved in McCluskey and van den Driessche [11] is used to demonstrate the local stability of endemic equilibrium point E * . Using the above lemma, we will study the stability of the endemic equilibrium.
Hence, by Lemma, the endemic equilibrium (E * ) of the model (4) is locally asymptotically stable in .

Global stability of endemic equilibrium
We now prove the global stability of the endemic equilibrium E * , when the reproduction number R 0 is greater than the unity. For this, first we will prove the following result.
Theorem 5.1: If R 0 > 1, then system (4) is uniformly persistent; that is, there exists c > 0 (independent of initial conditions), such that A uniform persistence result given in [4] requires the following hypothesis (H) to be satisfied.
Hypothesis Here, E is a closed positively invariant subset of X on which a continuous flow F is defined, N is the maximal invariant set of ∂F on ∂E, and {N} are pairwise disjoint closed invariant subsets of ∂E which cover N. Now, we prove the following result. Proof: Suppose R 0 > 1. We show that system (4) satisfies all the conditions of Theorem 4.3 in [4]. Choose X = R 3 and E = . The maximal invariant set N of ∂ is the singleton {E 0 } and it is isolated. Therefore, Hypothesis (H) holds for system (4). The flow induced by f (x) is point dissipative ∂ by the positive invariance of . Define W + (N) = {y ∈ X : + (y) ⊂ N}, where + (y) is the omega limit set of y. When R 0 ≤ 1, we have that E is contained in the set W + (N) and for R 0 > 1, W + (N) = φ. Therefore, the uniform persistence of system (4) is equivalent to E 0 being unstable, and the theorem is proved.
Here, we use the geometrical approach of Li and Muldowney [10] to investigate the global stability of the endemic equilibrium E * in the feasible region . We have omitted the detailed introduction of this approach and we refer the interested readers to see [10]. For the applications of the Li and Muldowney approach to host-vector models [2,25]. We summarize this approach as follows.
Consider a C 1 map f : x → f (x) from an open set D ⊂ R n to R n such that each solution x(t, x 0 ) to the differential equation is uniquely determined by the initial value x(0, x 0 ). We have following assumptions: there exists a compact absorbing set K ⊂ D; (H 3 ) Equation (14) has unique equilibriumx in D.
Let P : x → P(x) be a nonsingfular n 2 × n 2 matrix-valued function which is C 1 in D and a vector norm |.| on R N , where N = n 2 . Let μ be the Lozinskiȋ measure with respect to the |.|. Define a quantityq 2 as where B = P f P −1 + PJ [2] P −1 , the matrix P f is obtained by replacing each entry p of P by its derivative in the direction of f, (p ij ) f , and J [2] is the second additive compound matrix of the Jacobian matrix J of Equation (4). The following result has been established in [9].

Theorem 5.3:
Suppose that H 1 , H 2 and H 3 hold, the unique endemic equilibrium E * is globally stable in ifq 2 < 0.
We choose a suitable vector norm |.| in R 3 and a 3 × 3 matrix-valued function Obviously, P is C 1 and nonsingular in the interior of . Linearizing system (4) about an endemic equilibrium E * gives the following Jacobian matrix: The second additive compound matrix of J(E * ) is given by The matrix B = P f P −1 + PJ [2] P −1 can be written in block form as with, Consider the norm in R 3 as: |(u, v, w)| = max(|u|, |v| + |w|), where (u, v, w) denotes the vector in R 3 . The Lozinskiȋ measure with respect to this norm is defined as μ(B) ≤ sup(g 1 , g 2 ), where From system (4), we can write Since B 11 is a scalar, its Lozinskiȋ measure with respect to any vector norm in R 1 will be equal to B 11 . Thus and g 1 will become , |B 12 | and |B 21 | are the operator norms of B 12 and B 21 which are mapping from R 2 to R and from R to R 2 , respectively, and R 2 is endowed with the l 1 norm. μ 1 (B 22 ) is the Lozinskiȋ measure of 2 × 2 matrix B 22 with respect to l 1 norm in R 2 . Hence, Thus, Since Equation (4) is uniformly persistent when R 0 > 1, so for T > 0 such that t > T implies E h (t) ≥ c, I h (t) ≥ c and (1/t) log E h (t) < μ h /4 for all (E h (0), I h (0), I v (0)) ∈ K.

Discussion and simulations
This paper presents a host-vector model for PWD with nonlinear incidence rate. The mathematical analysis is carried out for a model for forest insect pests with PWD. The global dynamics of the model are shown to be determined by the basic reproduction number R 0 . More specifically, by constructing suitable Lyapunov function, we proved that if R 0 ≤ 1, then disease-free equilibrium E 0 is globally asymptotically stable in , and thus the disease always dies out. If R 0 > 1, the unique endemic equilibrium E * exists and is globally asymptotically stable, so that the disease persists at the endemic equilibrium if it is initially present. We know that the basic reproduction number R 0 of the model is proportional to the total number of the host tree population available as oviposition sites for the vector beetles and the number of vector population and host infectious rates and vector infectious rate, respectively. The basic reproduction number R 0 is independent of α 1 and α 2 . However, numerical simulations indicate that when the disease is endemic, the steady-state value of infectious hosts I * h decreases as α 1 increases(see Figure 1). The steady-state value of infectious vectors I * v decreases as α 2 increases when β 3 = 0 whereas when we take β 3 = 0 small change occurs in the steady-state value of infectious vectors I * v as α 2 increases. It means that saturation level of infectious pines do not effect the   steady-state value of infectious vectors I * v when vectors are getting infected through mating (see Figures 2 and 3).
Approach to diminish the vector population and this control means boosting vector mortality. By simple calculation, we can also see that R 0 is a decreasing function of μ v . The reproductive number would decrease by a factor 1/κ if we increase the vector mortality by a factor κ > 1. For example, the average lifespan of bark beetles is six months. If we apply insecticides every three months and assume that all beetles disappear then the average life span decreases from six months to three months. This would imply that μ v increases twice, R 0 would decrease by 1 4 . Hence, the disease could be controlled by the application of insecticides every three months if R 0 /4 < 1. We assume that h = 0.0009041, v = 0.0002691, mu h = 0.0000301, β 1 = 0.0057142, β 2 = 0.00305, β 3 = 0.00385, δ 1 = 0.000332, δ 2 = 0.00166, θ = 0.02 and average life span of bark beetles is six months, we obtain R 0 = 11.9431. 1/κ 2 × 11.9431 < 1, or κ > 3. 4, which means that we should increase vector mortality almost three-and-half times. Again assuming the extreme case, i.e. almost all beetles disappear after each application, R 0 would be less than one if pesticide is sprayed almost after 50 days (Figure 4).

Disclosure statement
No potential conflict of interest was reported by the authors.