# Analysis of the Mathematical Model for the Spread of Pine Wilt Disease

Authors: Xiangyun Shi; Guohua Song
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

This paper formulates and analyzes a pine wilt disease model. Mathematical analyses of the model with regard to invariance of nonnegativity, boundedness of the solutions, existence of nonnegative equilibria, permanence, and global stability are presented. It is proved that the global dynamics are determined by the basic reproduction numberℛ0 and the other value ℛc which is larger than ℛ0. If ℛ0 and ℛc are both less than one, the disease-free equilibrium is asymptotically stable and the pine wilt disease always dies out. If one is between the two values, though the pine wilt disease could occur, the outbreak will stop. If the basic reproduction number is greater than one, a unique endemic equilibrium exists and is globally stable in the interior of the feasible region, and the disease persists at the endemic equilibrium state if it initially exists. Numerical simulations are carried out to illustrate the theoretical results, and some disease control measures are especially presented by these theoretical results.


---

## Body

## 1. Introduction

Pine wilt is a dramatic disease of pine caused by the pinewood nematode (Bursaphelenchus xylophilus), which constitutes a major threat to forest ecosystems worldwide, from both the economical point of view and the environmental (landscape) perspective [1]. Pine trees are affected by pine wilt disease, wilt, and usually die within a few months. Symptoms of pine wilt disease usually become evident in late spring or summer. The first observable symptom is the lack of resin exudation from barks wounds. The foliage then becomes pale green, then yellow, and finally reddish brown when the tree succumbs to the disease. The wood in affected trees is dry and totally lacks resin.In [2], Evans et al. reviewed the principle of the Bursaphelenchus xylophilus transmission and disease dissemination. The Bursaphelenchus xylophilus is transmitted from pine tree to pine tree by a bark beetle called the pine sawyer (Monochamus alternatus) either when the sawyer beetles are fed on the bark and phloem of twigs of susceptible live trees (primary transmission) or when the female beetles lay eggs (oviposition) in freshly cut timber or dying trees (secondary transmission). Nematodes introduced during primary transmission can reproduce rapidly in the sapwood, and a susceptible host can wilt and die within weeks of being infested if conditions are favorable to disease development. In the summer, adult pine sawyers emerge from pine trees and fly to new trees. If the beetle is carrying the pine wood nematode, it spreads to the new trees (see Figure 1).Figure 1
Schematic representation of the interrelationships between the pinewood nematode,Bursaphelenchus xylophilus, and its insect vector (Monochamus alternatus) (adapted from Evans et al. [2]).Pine wilt disease causes significant economic losses in natural coniferous forests in Eastern Asia (especially Japan, China, and South Korea) and Western Europe (especially Portugal). As such, the pine wilt disease is among the most important pests included in the quarantine lists of many countries around the world [3]. In the beginning of the 20th century, pine wilt disease was first reported in Japan. Up to the present day, pine wilt disease has become the major ecological catastrophe of pine forests in Japan. For example, it has losses reaching over 2 million m3/year in the 1980s. Since then, it has spread to other Asian countries and regions such as China, Taiwan, and Korea, causing serious losses and economic damage. The pinewood nematode was first detected in Portugal’s Setúbal region in 1999. Immediately, several governments of the European Union prompted actions to assess the extent of the nematode’s presence and to contain Bursaphelenchus xylophilus and Monochamus alternatus in an area with a 30 km radius in the Setbal Peninsula, 20 km south of Lisbon [4]. Despite the dedicated and concerted actions of government agencies, this disease continues to spread. In 2006, a new strategy which is the coordination of the national program for the control of the pinewood nematode for the control and ultimately the eradication of the nematode has been announced in Portugal [4].Experience from control actions in Japan included aerial spraying of insecticides to control the insect vector (the cerambycid beetleMonochamus alternatus), injection of nematicides to the trunk of infected trees, slashing and burning of large areas out of control, beetle traps, biological control, and tree breeding programs [1]. These actions allowed not only some positive results, but also unsuccessful cases due to the pinewood nematode spread and virulence. Other Asian countries also followed similar strategies, but the nematode is still spreading in many regions.In recent years, many attempts have been made to develop realistic mathematical models for investigating the transmission dynamics of the pine wilt disease, see [5–7]. Previous studies have worked on modeling of population dynamics of the vector beetle (Monochamus alternatus) and the pine tree to explore expansion of the disease using an integro-difference equation with a dispersal kernel that describes beetle mobility. In this paper, we revisit these previous models but retain individuality by building a differential system model to study the relation between pine tree, Bursaphelenchus xylophilus, and its insect vector Monochamus alternatus in the long run. The following assumptions are made in formulating the mathematical model.(H1)
The population of pine is divided into two classes: pine which is normal and susceptible and pine which has been infected by Bursaphelenchus xylophilus.S(t), I(t) express the number of pine which is normal and susceptible and pine which is infected by Bursaphelenchus xylophilus at time t, respectively.(H2)
X(t), Y(t) Are the number of Monochamus alternatus which does not carry Bursaphelenchus xylophilus Nickle and which carries Bursaphelenchus xylophilus Nickle at time t, respectively.(H3)
The pine is infected by theMonochamus alternatus which carries Bursaphelenchus xylophilus Nickle, and this contact process is assumed to follow the standard incident rate β(Y/(X+Y)); the Monochamus alternatus which does not carry the Bursaphelenchus xylophilus Nickle is infected by Bursaphelenchus xylophilus Nickle through the infected pine, and this contact process is also assumed to follow the standard incident rate k(I/(S+I)).The dynamic of the whole model is indeed determined by the four-dimensional system involving only the pine andMonochamus alternatus. And the model for the transmission of the virus in the pine Monochamus alternatus cycle reads
(1)dSdt=b-βYX+YS-αS,dIdt=βYX+YS-δI,dXdt=a-kIS+IX-μX,dYdt=kIS+IX-μY,
where a is the constant input rate of Monochamus alternatus; b is the constant increase rate of pine; μ is the mortality of Monochamus alternatus; k is the probability that Monochamus alternatus carries Bursaphelenchus xylophilus Nickle; β is the probability that pine is infected by Bursaphelenchus xylophilus; α is the exploitation percent of pine which is normal and susceptible; δ is the percent isolated and felled of pine which has infected Bursaphelenchus xylophilus. Here we suppose that δ>α. All parameters are assumed to be positive.Let(2)A=X+Y.
It follows from the system (1) that A(t) satisfies the following differential equation:
(3)dAdt=a-μA.
This leads to A(t)→(a/μ) as t→∞. Thus, the system (1) is reduced to the following three-dimensional system:
(4)dSdt=b-μaβYS-αS,dIdt=μaβYS-δI,dYdt=kIS+I(aμ-Y)-μY.The initial conditions of the system (4) are assumed as follows:
(5)S(0)≥0,I(0)≥0,Y(0)≥0.This paper is organized as follows. In Section2, we present some preliminaries such as the positivity and the boundedness of solutions. In Section 3, we firstly calculate the basic reproduction number of system (4). Then we obtain the local and global stability of the disease-free equilibrium. In Section 4, we present the local and global stability of the endemic equilibrium of the system (4). In Section 5, we conclude with some numerical simulations and discussions.


## 2. Positivity and Boundedness of Solutions

It is important to prove that the solutions of the system (4) are positive and bounded with the positive initial conditions (5) for the biology meaning.Let(6)D={(S,I,Y)∣bδ≤S+I≤bα,0≤Y≤aμ}.Proposition 1.
All solutions(S,I,Y) of the system (4) are nonnegative. Moreover, D is a global attractor in ℝ+3 and positively invariant for the system (4).Proof.
The first statement is trivial. It easily follows from the argument for reduction in the last equation that0≤Y≤a/μ. It follows from the first two equations of the system (4) that
(7)b-δ(S+I)≤d(S+I)dt≤b-α(S+I)
for all t≥0, then the second statement follows immediately.


## 3. Stability of the Disease-Free Equilibrium

For all infectious diseases, the basic reproduction number, sometimes called the basic reproductive rate or the basic reproductive ratio, is one of the most useful threshold parameters which characterize mathematical problems concerning infectious diseases. This metric is useful because it helps determine whether or not an infectious disease will spread through a population. In this section, we will calculate the basic reproduction number of the system (4). Moreover, we will obtain the local and global stability of the disease-free equilibrium.It is easy to see that the system (4) always has a disease-free equilibrium (the absence of infection, i.e., I=Y=0), E0(S0,0,0), where S0=b/α.  Let x=(I,Y,S)⊤. Then system (4) can be written as
(8)dxdt=ℱ(x)-𝒱(x),
where
(9)ℱ(x)=(μaβYSkIS+I(aμ-Y)0),𝒱(x)=(δIμY-(b-μaβYS-αS)).
We can get
(10)F=(0μaβS-kI((a/μ)-Y)(S+I)2+k((a/μ)-Y)S+I-kIS+I),V=(δ00μ).
Submitting E0 into F, then
(11)F=(0μaβS0kaαμb0)
and giving
(12)V-1=(1δ001μ),(13)FV-1=(0βbaαkaαbδμ0)
is the next generation matrix for the system (4). It then follows that the spectral radius of matrix FV-1 is ρ(FV-1)=kβ/δμ. According to Theorem  2 in [8], the basic reproduction number of model (4) is
(14)ℛ0=kβδμ.In the following, we will discuss the local and global stability of the disease-free equilibrium. From above and [8], we can obtain the following theorems.Theorem 2.
The disease-free equilibriumE0 is locally asymptotically stable for ℛ0<1 and unstable for ℛ0>1.Theorem 3.
The disease-free equilibriumE0 is globally asymptotically stable if ℛc<1, where ℛc=kβ/αμ.Proof.
Define a Lyapunov functionL of the system (4) as follows:
(15)L=αabβI+Y.
Its derivative along a solution of the system (4) is
(16)dLdt=αabβμaβYS-αaδbβI+kaμ(S+I)I-kIS+IY-μY=(αμbS-kIS+I-μ)Y+(akμ(S+I)-αaδbβ)I≤(μ-kIS+I-μ)Y+(akδμb-αaδbβ)I≤aδ(kβ-αμ)μbβI.
It is clear from (16) that for ℛc≤1, dL/dt≤0 as b/δ≤S+I≤b/α. Furthermore, if M is the set of solutions of the system where dL/dt=0, then the Lyapunov-Lasalle Theorem [9] implies that all paths in D approach the largest positively invariant subset of the set M. Here, M is the set where I=0. On the boundary of D where I=0, we have Y=0 and dS/dt=b-αS. So S→b/α as t→∞. Hence, all solution paths in D approach the disease-free equilibrium E0.Remark 4.
This above result is of outmost importance because it shows that if at any time, through appropriate interventions, we are able to lowerℛ0 and ℛc to less than 1, then the pine wilt disease will disappear. Obviously, ℛ0<ℛc. The condition of global stability of disease-free equilibrium is stronger than that of local stability. In fact, since ℛ0<1<ℛc, it is possible to be no outbreak appearance, see Figure 2(a). ℛc>1 indicates that pine wilt disease could occur, while ℛ0<1 shows that the outbreak will stop since the infected pine decreases to the disease-free equilibrium. We know that ℛ0 and ℛc depend on μ, the mortality of Monochamus alternatus. Thus, if we are able to increase the value of μ, ℛ0 and ℛc will decrease.(a) Time series diagram of the system (4), where the blue and green continuous curves express the number of pines which is normal and susceptible S(t) and pine which is infected Bursaphelenchus xylophilus I(t), respectively. And the red continuous curves expresses the number of Monochamus alternatus which carries Bursaphelenchus xylophilus Nickle. The pine wilt disease will not break out, where b=0.25,μ=0.05,β=0.006,a=0.003,α=0.002,δ=0.01 and k=0.05. In this case, ℛ0=0.7746 and ℛc=1.7321. (b) Time series diagram of the system (4), where the blue and green continuous curve express the number of pines which is normal and susceptible S(t) and pine which is infected Bursaphelenchus xylophilus I(t), respectively. And the red continuous curve express the number of Monochamus alternatus carry Bursaphelenchus xylophilus Nickle. the endemic equilibrium E*(S*,I*,Y*) is globally asymptotically stable, where b=0.25,μ=0.0005,β=0.006,a=0.003,α=0.002,δ=0.008 and k=0.05. In this case, ℛ0=8.6603 and ℛc=17.3205. (c) Time series diagram of the system (4), where the blue and green continuous curves express the number of pines which is normal and susceptible S(t) and pine which is infected Bursaphelenchus xylophilus I(t), respectively. And the red continuous curve expresses the number of Monochamus alternatus which carries Bursaphelenchus xylophilus Nickle. The disease-free equilibrium E0(S0,0,0) of the system (4) is globally asymptotically stable, where b=0.25,μ=0.5,β=0.006,a=0.003,α=0.002,δ=0.008 and k=0.05. In this case, ℛ0=0.2739 and ℛc=0.5477.
(a)
(b)
(c)


## 4. Global Stability of the Endemic Equilibrium

In this section, we will discuss the local and global stability of the endemic equilibrium. The endemic equilibriumE*(S*,I*,Y*) of the system (4) can be deduced by the following equations:
(17)b-μaβY*S*-αS*=0,μaβY*S*-δI*=0,kI*S*+I*(aμ-Y*)-μY*=0.Clearly, the system (4) has a unique endemic equilibrium E*(S*,I*,Y*) when ℛ0>1 and no endemic equilibrium when ℛ0≤1, where S*=b(k+μ)/(kβ+kα-μδ+μα), I*=(kβ-μδ)b/(kβ+kα-μδ+μα)δ, and Y*=(kβ-μδ)a/(k+μ)μβ.In the following, we will consider the locally asymptotical stability of the positive equilibrium whenℛ0>1.Theorem 5.
Ifℛ0>1, the endemic equilibrium E* is locally asymptotically stable.Proof.
Jacobian matrix of the system (4) is evaluated at the endemic equilibrium E*:
(18)J(E*)=(-μaβY*-α0-μaβS*μaβY*-δμaβS*-μY*S*+I*μY*I*-μY*S*+I*-kI*S*+I*-μ).
The eigenvalue problem for the Jacobian matrix (18) provides the characteristic equation
(19)λ3+a1λ2+a2λ+a3=0,
where the coefficients ai (i=1,2,3) are
(20)a1=(kβ-μδ)×[ab2k(αδk+βμ2+2kβμ+αkβ+2δkβ+kβ2+αδμ+δ2k+αβμ+k2β)]×((k+μ)(kβ+kα-μδ+μα)2δ2a(S*+I*)I*)-1,a2=(kβ-μδ){ab2k[(δk+δ2+μβ+δβ+kβ+μδ)(kβ-μδ)×(kβ-μδ)+αδ2k+αk2β+αδβμ+αδ2μ+μ2αβ(kβ-μδ)+αδkβ+2αkβμ{ab2}]}×((k+μ)(kβ+kα-μδ+μα)2δ2a(S*+I*)I*)-1,a3=(kβ-μδ)2ab2δ(kβ+kα-μδ+μα)a(S*+I*)I*.
Note thata1>0, a2>0, and a3>0 if ℛ0>1.
Expressinga1a2-a3 in terms of α, we have
(21)a1a2-a3=d1α2+d2α+d3α,
where
(22)d1=k(δ+β)(k+μ)2(βμ+kβ+δβ+δ2),d2=(k+μ)(kβ2μ3+μ3δ2β+μ3δ3+3k2β2μ2+μ2δ3k+μδ2k2β+3μδk2β2+3k3β2μ+k4β2+4k3β2δ+2k3β3+2k3βδ2+δ4k2+(4kβ-δμ)kδ3+(5kβ-2δμ)kδ2β+(2kβ-δμ)kβ2(δ+μ)),d3=(δ+β)(kβ-μδ)×(β2δk2+β2k2μ+β2k3+βk4+3βk3μ+2βδk3+2βk2δμ+2δ2k2β+3k2βμ2+kβμ3+δ2k3+2δ2k2μ+k2δ3+2δ2kμ2+δ2μ3).
Clearly, d1,d2,d3>0 since ℛ0>1. Then it becomes obvious from the expression for d1,d2, and d3 that a1,a2,a3>0. Hence, thanks to the Routh-Hurwitz criterion all eigenvalues of J(E*) have negative real part, and consequently E* is locally asymptotically stable.
We have shown thatℛ0>1 implies the existence and uniqueness of the endemic equilibrium E*. In the following, we provide sufficient conditions leading to a globally asymptotically stable infected steady state when ℛ0>1. The stability analysis of E* will be here performed through the geometric approach to global stability due to Li and Muldowney [10–12]. Firstly, we will summarize the main facts related to our research.
Let us consider the system of differential equations(23)dXdt=F(X),X∈D,
where D is an open subset on ℝ3 and F is twice continuously differentiable in D. The noncontinuable solution of (23) satisfying X(0)=X0 is denoted by X(t,X0), the positive (negative) semiorbit through X0 is denoted by ϕ+(X0)(ϕ-(X0)), and the orbit through X0 is denoted by ϕ(0)=ϕ-(X0)∪ϕ+(X0). We use the notation ω(X0)(α(X0)) for the positive (negative) limit set of ϕ+(X0)(ϕ-(X0)), provided the latter semiorbit has compact closure in D.
The system (23) is competitive in D [13–16] if, for some diagonal matrix H=diag(ϵ1,ϵ2,ϵ3), where ϵi is either 1 or -1, H(DF(X))H has nonpositive off-diagonal elements for X∈D, where DF(X) is the Jacobian of (23). It is shown in [16] that if D is convex the flow of such a system preserves for t<0 the partial order in ℝ3 defined by the orthant
(24)K1={(X1,X2,X3)∈ℝ3∣ϵiXi≥0}.
Hirsch [13] and Smith [15, 16] proved that three-dimensional competitive systems that live in convex sets have the Poincare-Bendixson property [17]; that is, any nonempty compact omega limit set that contains no equilibria must be a closed orbit.
We recall additional definitions that we will use later. We first recall the basic definitions in [18]. Suppose that (23) has a periodic solution x=p(t) with minimal period ω>0 and orbit γ={p(t):0≤t≤ω}. This orbit is orbitally stable if and only if, for each ε>0, there exists a δ1>0 such that any solution x(t), for which the distance of x(0) from γ is less δ, remains at a distance less than ε from γ for all t≥0. It is asymptotically orbitally stable, if the distance of x(t) from γ also tends to 0 as t goes to ∞. This orbit γ is asymptotically orbitally stable with asymptotic phase if it is asymptotically orbitally stable and there exists a b>0, such that, any solution x(t), for which the distance of x(0) from γ is less than b, satisfying |x(t)-p(t-τ)|→0 as t→+∞ for some τ which may depend on x(0).Definition 6.
We say that the system (23) has the property of stability of periodic orbits if and only if the orbit of any periodic solution γ(t), if it exists, is asymptotically orbitally stable.The following lemma is the main tool to prove the global stability of the endemic equilibrium with disease.Lemma 7 (see [19]).
Assume thatn=3 and D is convex and bounded. Suppose that (4) is competitive and persistent and has the property of stability of the periodic orbits. If x0 is the only equilibrium in int (D) and if it is locally asymptotically stable, then it is globally asymptotically stable in int (D).In order to apply this lemma to prove the globally asymptotically stability of the endemic equilibrium, we will prove the persistence of the system (4).Theorem 8.
On the boundary ofD, the system (4) has only one ω-limit point which is the equilibrium E0. Moreover, for ℛ0>1, E0 cannot be the ω-limit of any orbit in int(D).Proof.
The vector field is transversal to the boundary ofD except in the S-axis, which is invariant with respect to (4). On the S-axis we have
(25)dSdt=b-αS,
which implies that S→b/α as t→∞. Therefore, E0 is the only ω-limit point on the boundary of D.
To prove the second part of the position, we consider the function(26)V1=Y+αabβ1+ℛ022I,
the derivative of which along solutions is given by
(27)dV1dt=dYdt+αabβ1+ℛ022dIdt=kS+IaμI-kYIS+I-μY+αμb1+ℛ022βYS-αaδbβ1+ℛ022I≥(ka(b/α)μ-αaδbβ1+ℛ022)I+(αμb1+ℛ022S-kb/δI-μ)Y≥(1-1ℛ021+ℛ022)bμkaαI+(S-kδαμ21+ℛ02I-21+ℛ02bα)bαb21+ℛ02Y.
Since ℛ0>1, then ((1/ℛ02)+1)/2<1 and 2/(1+ℛ02)<1. Therefore, there exists a neighborhood U of E0 such that for (S,I,Y)∈U∪int(D) the expression inside the brackets is positive. In this neighborhood, we have V1>0 unless I=Y=0. Moreover, the level sets of V1 are the planes
(28)Y+αabβ1+ℛ022I=C,
which go away from the S-axis as C increase. Since V1 increases along the orbits starting in U∪int(D), we conclude that they go away from E0. This proves the proposition and therefore the persistence of system (4) when ℛ0>1.
By looking at the Jacobian matrix of the system (4) and choosing the matrix H as
(29)H=(-10001000-1),
we can see that the system (4) is competitive in D, with respect to the partial order defined by the orthant K1={(S,I,Y)∈ℝ3:S≤0, I≤0, Y≤0}. Our main results will follow from this observation and the above theorems.Theorem 9.
Ifℛ0>1, then the positive equilibrium E* of the system (4) is globally asymptotically stable.Proof.
Since the system (4) is competitive permanent if ℛ0>1 holds and the only equilibrium point E* of the system (4) is locally asymptotically stable. Furthermore, in accordance with Lemma 7, Theorem 9 would be established if we show that the system (4) has the property of stability of periodic orbits. In the following, we prove it.
Letp(t)=(S(t),I(t),Y(t)) be a periodic solution whose orbit γ is contained in int(R+3) and suppose that its minimal period is ω>0. The second compound equation is the following periodic linear system:
(30)dZ(t)dt=∂f[2]∂x(p(t))Z(t),
where Z=(Z1,Z2,Z3)⊤ and ∂f/∂x is derived from the Jacobian matrix of the system (4) and defined as follows: (31)∂f[2]∂x=(-μaβY-α-δμaβSμaβS0-k((a/μ)-Y)I(S+I)2+k((a/μ)-Y)S+I-μaβY-α-kIS+I-μ0k((a/μ)-Y)I(S+I)2μaβY-δ-kIS+I-μ).For the solution p(t), (30) becomes
(32)Z˙1(t)=-(μaβY+α+δ)Z1+μaβS(Z2+Z3),Z˙2(t)=(-k((a/μ)-Y)I(S+I)2+k((a/μ)-Y)S+I)Z1-(μaβY+α+kIS+I+μ)Z2,Z˙3(t)=k((a/μ)-Y)I(S+I)2Z1+μaβYZ2-(δ+kIS+I+μ)Z3.
To prove that (32) is globally asymptotically stable, we shall use the following Lyapunov function:
(33)L(Z1,Z2,Z3;S(t),I(t),Y(t))=sup{|Z1|,IY(|Z2|+|Z3|)}.
Function (33) is positive, but not differentiable everywhere. Fortunately, this lack of differentiability can be remedied by using the right derivative of L(t), denoted as D+L(t). Then we have the following equalities:
(34)|Z˙1(t)|≤-(μaβY+α+δ)|Z1|+μaβS(|Z2|+|Z3|)≤-(μaβY+α+δ)|Z1|+μaβSYI(IY|Z2|+|Z3|),|Z˙2(t)|≤(-k((a/μ)-Y)I(S+I)2+k((a/μ)-Y)S+I)|Z1|-(μaβY+α+kIS+I+μ)|Z2|,|Z˙3(t)|≤k((a/μ)-Y)I(S+I)2|Z1|+μaβY|Z2|-(δ+kIS+I+μ)|Z3|.
Therefore,(35)D+(IY(|Z2(t)|+|Z3(t)|))=(I˙I-Y˙Y)IY(|Z2|+|Z3|)+IYD+(|Z2|+|Z3|)≤(I˙I-V˙V)IY(|Z2|+|Z3|)+IY[k((a/μ)-Y)S+I|Z1|-(α+kIS+I+μ)(|Z2|+|Z3|)]≤(I˙I-Y˙Y-α-kIS+I-μ)(|Z2|+|Z3|)+IYk((a/μ)-Y)S+I|Z1|.Define(36)g1(t)=-(μaβY+α+δ)+μaβSYI,g2(t)=IYk((a/μ)-Y)S+I+(I˙I-Y˙Y-α-kIS+I-μ).
Rewriting the last two equations of (4) as
(37)I˙I+δ=μaβSYI,Y˙Y+kIS+I+μ=IYka(S+I)μ,
then
(38)g1(t)=-(μaβY+α+δ)+μaβSYI=I˙I-μaβY-α,g2(t)=IYk((a/μ)-Y)S+I+(I˙I-Y˙Y-α-kIS+I-μ)=I˙I-kIS+I-α.Thus, we obtain(39)D+L(t)≤sup{g1(t),g2(t)}L(t).
It follows from (38) that g1(t)≤(I˙/I)-α, and thus, from (39) and Gronwall’s inequality, we obtain
(40)L(t)≤L(0)Ie-αt.
Since b/δ≤I≤b/α in int(D), which implies that L(t)→0 as t→∞, therefore, (Z1(t),Z2(t), Z3(t))→(0,0,0) as t→∞. This implies that the linear system (32) is asymptotically stable, and therefore, the periodic solution is asymptotically orbitally stable. This proves the theorem.


## 5. Numerical Simulations and Disease Control

In this paper, we investigate the dynamical behavior of a pine wilt disease model that incorporates a standard incidence rate. Qualitative analysis of the model is presented. The model has two equilibria, the disease-free equilibrium and endemic equilibrium. The behavior of the system (4) near each equilibrium has been studied. Threshold value of the relative basic reproductive number ℛ0 which determines the spread of infection has been worked out.We obtain the basic reproduction number of the system (4), ℛ0=kβ/δμ. Furthermore, we obtain that the disease-free equilibrium E0(S0,0,0) of the system (4) is globally asymptotically stable when ℛ0<ℛc<1 (see Figure 2(c)). And if ℛ0>1, the endemic equilibrium E*(S*,I*,Y*) is globally asymptotically stable (see Figure 2(b)). That is to say, if ℛ0 and ℛc are lower than the threshold 1, the system of forest insect pest appears in disease-free state, that is, Bursaphelenchus xylophilus will be eradicated ultimately, which is what we hope; if ℛc>1, while ℛ0<1, though the pine wilt disease could occur, the outbreak will stop since the infected pine decreases to the disease-free equilibrium. If ℛ0 is larger than 1, the system of forest insect pest have insect pests equilibrium, that is Bursaphelenchus xylophilus will not disappear, and with the time to develop, susceptible pine, infected pine, Bursaphelenchus xylophilus, and its vector will be extend to a stable state, that is, pine wilt disease will in an endemic steady state.By simple calculation, we can find∂ℛc/∂μ<0. Hence, we can make ℛc<1 by increasing μ. Also from Figure 2, we can find that ℛc=17.3205>1 in Figure 2(b) and the pine wilt disease is not disappearing. As μ increases from 0.005 to 0.5, we calculate the ℛc=0.5477<1 in Figure 2(c), keeping all the other parameters the same as in Figure 2(b). Thus, the disease-free equilibrium E0(S0,0,0) of system (4) is globally asymptotically stable, and Bursaphelenchus xylophilus will be eradicated ultimately. From the Remark 4, we know that the pine wilt disease will not break out when ℛ0<1. So we only need to lower ℛ0 to less than 1 to control the pine wilt disease. We can also find ∂ℛ0/∂δ<0 and ∂ℛ0/∂μ<0 by a simple calculation. Hence, we can make ℛ0<1 by increasing μ or (and) δ. Since μ and δ are the mortality of Monochamus alternatus and the percent isolated and felled of pine which has been infected by Bursaphelenchus xylophilus, respectively, we can increase μ or (and) δ by various control efforts in reality. On the one hand, we can take a set of measures to increase the mortality of Monochamus alternatus, for example, setting out beetle traps, setting vertical wood traps, releasing the natural enemies of Bursaphelenchus xylophagous (such as releasing the Scleroderma guanior), using chemicals to kill sawyer beetles to reduce the number of longhomed beetle. On the other hand, we can also take a set of techniques for specialized treatment of timber to increase the mortality of infected longhorned beetle; removal of dying branches on affected trees or dead and dying trees and fumigation of nematode-infested lumber.But it is difficult to makeℛ0<1 by increasing μ or (and) δ dramatically in practice. In fact, although these above measures are very effective to control pine wilt disease, they have not yet been practiced to eradicate pine wilt disease ultimately because these measures require more cost and labor and even entail danger of forest fires that most owners of forests hesitate to use these measures.From Theorem8, we can know that the endemic equilibrium E* is globally asymptotically stable when ℛ0>1. That is to say, if we do not control the pine wilt disease, then the disease will be prevalent, and it will achieve a balance in the ecological environment finally. But it will suffer large economic losses if we let the disease spread at will. If we do not control the parameters well, the pine which has been infected by Bursaphelenchus xylophilus at the endemic equilibrium will reach a significant number. As the infected wood can not be used as wood products, the loss of afforestation will be large. We have to think of any other ways to control pine wilt disease. We can control the number of infected pine trees in the infected equilibrium point by controlling some parameters.We can also get∂I*/∂δ=-b[(μδ-βk)2+βk2α+βkμα]/(βk+kα-μδ+μα)2δ2<0 and ∂I*/∂μ=-bkα(δ+β)/(βk+kα-μδ+μα)2δ<0. I* is a decreasing function of parameter δ (or μ). Hence, we can decrease I* by increasing δ (or μ) properly. In this case, we need not increase δ (or μ) large enough to make ℛ0<1 to eradicate the pine wilt disease, and we only need to control the infected pine in a lower level stable state (see Figure 3). In Figure 3(a), we illustrate the effect of varying μ on I*. Keeping all the other parameters the same as in Figure 2(b), we plot I(t) for different μ (μ=0.0005,0.001,0.0015,0.002,0.0025). In Figure 3(b), we illustrate the effect of varying δ on I*. Keeping all the other parameters the same as in Figure 2(b), we plot I(t) for different δ (δ=0.008,0.016,0.024,0.032,0.04). The numerical simulations show that the infected pine I(t) approaches to the steady state I*, and I* decreases as μ or δ increases. Though the Bursaphelenchus xylophilus exists, the pine wilt disease spreads only in a low level, which can lower our economic losses to a large extent, and the pine wilt disease will not break out in the long time.(a) Simulation results showing the effect of decreasing infected pine trees to parameterμ for μ=0.0005,0.001,0.0015,0.002,0.0025, and other parameter values are as given in Figure 2(b). (b) Simulation results showing the effect of decreasing infected pine trees to parameter δ for δ=0.008,0.016,0.024,0.032,0.04, and other parameter values are as given in Figure 2(b).
(a)
(b)



---

*Source: 184054-2013-03-10.xml*