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

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 ℛ₀ and the other value ℛ_𝑐 which is larger than ℛ₀. If ℛ₀ and ℛ_𝑐 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.

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 m³/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 beetle Monochamus 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.(𝐻₁) The population of pine is divided into two classes: pine which is normal and susceptible and pine which has been infected by Bursaphelenchus xylophilus. 𝑆(𝑡), 𝐼(𝑡) express the number of pine which is normal and susceptible and pine which is infected by Bursaphelenchus xylophilus at time 𝑡, respectively.(𝐻₂) 𝑋(𝑡), 𝑌(𝑡) Are the number of Monochamus alternatus which does not carry Bursaphelenchus xylophilus Nickle and which carries Bursaphelenchus xylophilus Nickle at time 𝑡, respectively.(𝐻₃)  The pine is infected by the Monochamus alternatus which carries Bursaphelenchus xylophilus Nickle, and this contact process is assumed to follow the standard incident rate 𝛽(𝑌/(𝑋+𝑌)); 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 𝑘(𝐼/(𝑆+𝐼)).

The dynamic of the whole model is indeed determined by the four-dimensional system involving only the pine and Monochamus alternatus. And the model for the transmission of the virus in the pine Monochamus alternatus cycle reads 𝑑𝑆𝑌𝑑𝑡=𝑏−𝛽𝑋+𝑌𝑆−𝛼𝑆,𝑑𝐼𝑌𝑑𝑡=𝛽𝑋+𝑌𝑆−𝛿𝐼,𝑑𝑋𝐼𝑑𝑡=𝑎−𝑘𝑆+𝐼𝑋−𝜇𝑋,𝑑𝑌𝐼𝑑𝑡=𝑘𝑆+𝐼𝑋−𝜇𝑌, (1) where 𝑎 is the constant input rate of Monochamus alternatus; 𝑏 is the constant increase rate of pine; 𝜇 is the mortality of Monochamus alternatus; 𝑘 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) It follows from the system (1) that 𝐴(𝑡) satisfies the following differential equation: 𝑑𝐴𝑑𝑡=𝑎−𝜇𝐴. (3) This leads to 𝐴(𝑡)→(𝑎/𝜇) as 𝑡→∞. Thus, the system (1) is reduced to the following three-dimensional system: 𝑑𝑆𝜇𝑑𝑡=𝑏−𝑎𝛽𝑌𝑆−𝛼𝑆,𝑑𝐼=𝜇𝑑𝑡𝑎𝛽𝑌𝑆−𝛿𝐼,𝑑𝑌=𝑑𝑡𝑘𝐼𝑎𝑆+𝐼𝜇−𝑌−𝜇𝑌. (4)

The initial conditions of the system (4) are assumed as follows: 𝑆(0)≥0,𝐼(0)≥0,𝑌(0)≥0. (5)

This paper is organized as follows. In Section 2, 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 )∣𝑏𝐷=(𝑆,𝐼,𝑌𝛿𝑏≤𝑆+𝐼≤𝛼𝑎,0≤𝑌≤𝜇. (6)

Proposition 1. All solutions (𝑆,𝐼,𝑌) of the system (4) are nonnegative. Moreover, 𝐷 is a global attractor in ℝ³₊ 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 that 0≤𝑌≤𝑎/𝜇. It follows from the first two equations of the system (4) that ()≤𝑑()𝑏−𝛿𝑆+𝐼𝑆+𝐼()𝑑𝑡≤𝑏−𝛼𝑆+𝐼 (7) for all 𝑡≥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., 𝐼=𝑌=0), 𝐸₀(𝑆₀,0,0), where 𝑆₀=𝑏/𝛼. Let 𝑥=(𝐼,𝑌,𝑆)^⊤. Then system (4) can be written as 𝑑𝑥(𝑥)(𝑥),𝑑𝑡=ℱ−𝒱 (8) where ℱ(𝑥)=𝜇𝑎𝛽𝑌𝑆𝑘𝐼𝑎𝑆+𝐼𝜇0,𝒱(𝑥)=−𝜇−𝑌𝛿𝐼𝜇𝑌𝑏−𝑎.𝛽𝑌𝑆−𝛼𝑆 (9) We can get 0𝜇𝐅=𝑎𝛽𝑆−𝑘𝐼𝑎/𝜇−𝑌(𝑆+𝐼)²+𝑘𝑎/𝜇−𝑌−𝑆+𝐼𝑘𝐼,.𝑆+𝐼𝐕=𝛿00𝜇 (10) Submitting 𝐸₀ into 𝐅, then 0𝜇𝐅=𝑎𝛽𝑆₀𝑘𝑎𝛼0𝜇𝑏 (11) and giving 𝐕^−𝟏=1𝛿001𝜇,𝐅𝐕^−𝟏=0𝛽𝑏𝑎𝛼𝑘𝑎𝛼0𝑏𝛿𝜇 (12)(13) is the next generation matrix for the system (4). It then follows that the spectral radius of matrix 𝐅𝐕^−𝟏 is 𝜌(𝐅𝐕^−𝟏)=𝑘𝛽/𝛿𝜇. According to Theorem  2 in [8], the basic reproduction number of model (4) is ℛ₀=𝑘𝛽.𝛿𝜇 (14)

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 equilibrium 𝐸₀ is locally asymptotically stable for ℛ₀<1 and unstable for ℛ₀>1.

Theorem 3. The disease-free equilibrium 𝐸₀ is globally asymptotically stable if ℛ_𝑐<1, where ℛ_𝑐=𝑘𝛽/𝛼𝜇.

Proof. Define a Lyapunov function 𝐿 of the system (4) as follows: 𝐿=𝛼𝑎𝑏𝛽𝐼+𝑌. (15) Its derivative along a solution of the system (4) is 𝑑𝐿=𝑑𝑡𝛼𝑎𝜇𝑏𝛽𝑎𝛽𝑌𝑆−𝛼𝑎𝛿𝑏𝛽𝐼+𝑘𝑎𝜇(𝑆+𝐼)𝐼−𝑘𝐼=𝑆+𝐼𝑌−𝜇𝑌𝛼𝜇𝑏𝑆−𝑘𝐼𝑆+𝐼−𝜇𝑌+𝑎𝑘𝜇()−𝑆+𝐼𝛼𝑎𝛿𝐼≤𝑏𝛽𝜇−𝑘𝐼𝑆+𝐼−𝜇𝑌+𝑎𝑘𝛿−𝜇𝑏𝛼𝑎𝛿𝐼≤𝑏𝛽𝑎𝛿𝑘𝛽−𝛼𝜇𝜇𝑏𝛽𝐼. (16) It is clear from (16) that for ℛ_𝑐≤1, 𝑑𝐿/𝑑𝑡≤0 as 𝑏/𝛿≤𝑆+𝐼≤𝑏/𝛼. Furthermore, if 𝑀 is the set of solutions of the system where 𝑑𝐿/𝑑𝑡=0, then the Lyapunov-Lasalle Theorem [9] implies that all paths in 𝐷 approach the largest positively invariant subset of the set 𝑀. Here, 𝑀 is the set where 𝐼=0. On the boundary of 𝐷 where 𝐼=0, we have 𝑌=0 and 𝑑𝑆/𝑑𝑡=𝑏−𝛼𝑆. So 𝑆→𝑏/𝛼 as 𝑡→∞. Hence, all solution paths in 𝐷 approach the disease-free equilibrium 𝐸₀.

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 ℛ₀ and ℛ_𝑐 to less than 1, then the pine wilt disease will disappear. Obviously, ℛ₀<ℛ_𝑐. The condition of global stability of disease-free equilibrium is stronger than that of local stability. In fact, since ℛ₀<1<ℛ_𝑐, it is possible to be no outbreak appearance, see Figure 2(a). ℛ_𝑐>1 indicates that pine wilt disease could occur, while ℛ₀<1 shows that the outbreak will stop since the infected pine decreases to the disease-free equilibrium. We know that ℛ₀ and ℛ_𝑐 depend on 𝜇, the mortality of Monochamus alternatus. Thus, if we are able to increase the value of 𝜇, ℛ₀ and ℛ_𝑐 will decrease.

(a)

(b)

(c)

Figure 2 

(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 𝑆(𝑡) and pine which is infected Bursaphelenchus xylophilus 𝐼(𝑡), 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 𝑏=0.25,𝜇=0.05,𝛽=0.006,𝑎=0.003,𝛼=0.002,𝛿=0.01 and 𝑘=0.05. In this case, ℛ0=0.7746 and ℛ𝑐=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 𝑆(𝑡) and pine which is infected Bursaphelenchus xylophilus 𝐼(𝑡), respectively. And the red continuous curve express the number of Monochamus alternatus carry Bursaphelenchus xylophilus Nickle. the endemic equilibrium 𝐸∗(𝑆∗,𝐼∗,𝑌∗) is globally asymptotically stable, where 𝑏=0.25,𝜇=0.0005,𝛽=0.006,𝑎=0.003,𝛼=0.002,𝛿=0.008 and 𝑘=0.05. In this case, ℛ0=8.6603 and ℛ𝑐=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 𝑆(𝑡) and pine which is infected Bursaphelenchus xylophilus 𝐼(𝑡), respectively. And the red continuous curve expresses the number of Monochamus alternatus which carries Bursaphelenchus xylophilus Nickle. The disease-free equilibrium 𝐸0(𝑆0,0,0) of the system (4) is globally asymptotically stable, where 𝑏=0.25,𝜇=0.5,𝛽=0.006,𝑎=0.003,𝛼=0.002,𝛿=0.008 and 𝑘=0.05. In this case, ℛ0=0.2739 and ℛ𝑐=0.5477.

4. Global Stability of the Endemic Equilibrium

In this section, we will discuss the local and global stability of the endemic equilibrium. The endemic equilibrium 𝐸^∗(𝑆^∗,𝐼^∗,𝑌^∗) of the system (4) can be deduced by the following equations: 𝜇𝑏−𝑎𝛽𝑌^∗𝑆^∗−𝛼𝑆^∗𝜇=0,𝑎𝛽𝑌^∗𝑆^∗−𝛿𝐼^∗=0,𝑘𝐼^∗𝑆^∗+𝐼^∗𝑎𝜇−𝑌^∗−𝜇𝑌^∗=0. (17)

Clearly, the system (4) has a unique endemic equilibrium 𝐸^∗(𝑆^∗,𝐼^∗,𝑌^∗) when ℛ₀>1 and no endemic equilibrium when ℛ₀≤1, where 𝑆^∗=𝑏(𝑘+𝜇)/(𝑘𝛽+𝑘𝛼−𝜇𝛿+𝜇𝛼), 𝐼^∗=(𝑘𝛽−𝜇𝛿)𝑏/(𝑘𝛽+𝑘𝛼−𝜇𝛿+𝜇𝛼)𝛿, and 𝑌^∗=(𝑘𝛽−𝜇𝛿)𝑎/(𝑘+𝜇)𝜇𝛽.

In the following, we will consider the locally asymptotical stability of the positive equilibrium when ℛ₀>1.

Theorem 5. If ℛ₀>1, the endemic equilibrium 𝐸^∗ is locally asymptotically stable.

Proof. Jacobian matrix of the system (4) is evaluated at the endemic equilibrium 𝐸^∗: 𝐽𝐸^∗=−𝜇𝑎𝛽𝑌^∗𝜇−𝛼0−𝑎𝛽𝑆^∗𝜇𝑎𝛽𝑌^∗𝜇−𝛿𝑎𝛽𝑆^∗−𝜇𝑌^∗𝑆^∗+𝐼^∗𝜇𝑌^∗𝐼^∗−𝜇𝑌^∗𝑆^∗+𝐼^∗−𝑘𝐼^∗𝑆^∗+𝐼^∗−𝜇. (18)
The eigenvalue problem for the Jacobian matrix (18) provides the characteristic equation 𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃=0, (19) where the coefficients 𝑎_𝑖 (𝑖=1,2,3) are 𝑎₁=×𝑘𝛽−𝜇𝛿𝑎𝑏²𝑘𝛼𝛿𝑘+𝛽𝜇²+2𝑘𝛽𝜇+𝛼𝑘𝛽+2𝛿𝑘𝛽+𝑘𝛽²+𝛼𝛿𝜇+𝛿²𝑘+𝛼𝛽𝜇+𝑘²𝛽×𝑘+𝜇𝑘𝛽+𝑘𝛼−𝜇𝛿+𝜇𝛼²𝛿²𝑎𝑆^∗+𝐼^∗𝐼^∗−1,𝑎₂=𝑘𝛽−𝜇𝛿𝑎𝑏²𝑘𝛿𝑘+𝛿²×+𝜇𝛽+𝛿𝛽+𝑘𝛽+𝜇𝛿𝑘𝛽−𝜇𝛿+𝛼𝛿²𝑘+𝛼𝑘²𝛽+𝛼𝛿𝛽𝜇+𝛼𝛿²𝜇+𝜇²×𝛼𝛽+𝛼𝛿𝑘𝛽+2𝛼𝑘𝛽𝜇𝑘+𝜇𝑘𝛽+𝑘𝛼−𝜇𝛿+𝜇𝛼²𝛿²𝑎𝑆^∗+𝐼^∗𝐼^∗−1,𝑎₃=𝑘𝛽−𝜇𝛿²𝑎𝑏²𝛿𝑘𝛽+𝑘𝛼−𝜇𝛿+𝜇𝛼𝑎(𝑆^∗+𝐼^∗)𝐼^∗. (20)
Note that 𝑎₁>0, 𝑎₂>0, and 𝑎₃>0 if ℛ₀>1.
Expressing 𝑎₁𝑎₂−𝑎₃ in terms of 𝛼, we have 𝑎₁𝑎₂−𝑎₃=𝑑₁𝛼²+𝑑₂𝛼+𝑑₃𝛼, (21) where 𝑑₁=𝑘𝛿+𝛽𝑘+𝜇²𝛽𝜇+𝑘𝛽+𝛿𝛽+𝛿²,𝑑₂=𝑘+𝜇𝑘𝛽²𝜇³+𝜇³𝛿²𝛽+𝜇³𝛿³+3𝑘²𝛽²𝜇²+𝜇²𝛿³𝑘+𝜇𝛿²𝑘²𝛽+3𝜇𝛿𝑘²𝛽²+3𝑘³𝛽²𝜇+𝑘⁴𝛽²+4𝑘³𝛽²𝛿+2𝑘³𝛽³+2𝑘³𝛽𝛿²+𝛿⁴𝑘²+4𝑘𝛽−𝛿𝜇𝑘𝛿³+5𝑘𝛽−2𝛿𝜇𝑘𝛿²𝛽+2𝑘𝛽−𝛿𝜇𝑘𝛽²,𝑑𝛿+𝜇₃=×𝛽𝛿+𝛽𝑘𝛽−𝜇𝛿₂𝛿𝑘²+𝛽²𝑘²𝜇+𝛽²𝑘³+𝛽𝑘⁴+3𝛽𝑘³𝜇+2𝛽𝛿𝑘³+2𝛽𝑘²𝛿𝜇+2𝛿²𝑘²𝛽+3𝑘²𝛽𝜇²+𝑘𝛽𝜇³+𝛿²𝑘³+2𝛿²𝑘²𝜇+𝑘²𝛿³+2𝛿²𝑘𝜇²+𝛿²𝜇³. (22) Clearly, 𝑑₁,𝑑₂,𝑑₃>0 since ℛ₀>1. Then it becomes obvious from the expression for 𝑑₁,𝑑₂, and 𝑑₃ that 𝑎₁,𝑎₂,𝑎₃>0. Hence, thanks to the Routh-Hurwitz criterion all eigenvalues of 𝐽(𝐸^∗) have negative real part, and consequently 𝐸^∗ is locally asymptotically stable.
We have shown that ℛ₀>1 implies the existence and uniqueness of the endemic equilibrium 𝐸^∗. In the following, we provide sufficient conditions leading to a globally asymptotically stable infected steady state when ℛ₀>1. The stability analysis of 𝐸^∗ 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) where 𝐷 is an open subset on ℝ³ and 𝐹 is twice continuously differentiable in 𝐷. The noncontinuable solution of (23) satisfying 𝑋(0)=𝑋₀ is denoted by 𝑋(𝑡,𝑋₀), the positive (negative) semiorbit through 𝑋₀ is denoted by 𝜙⁺(𝑋₀)(𝜙⁻(𝑋₀)), and the orbit through 𝑋₀ is denoted by 𝜙(0)=𝜙⁻(𝑋₀)∪𝜙⁺(𝑋₀). We use the notation 𝜔(𝑋₀)(𝛼(𝑋₀)) for the positive (negative) limit set of 𝜙⁺(𝑋₀)(𝜙⁻(𝑋₀)), provided the latter semiorbit has compact closure in 𝐷.
The system (23) is competitive in 𝐷 [13–16] if, for some diagonal matrix 𝐻=diag(𝜖₁,𝜖₂,𝜖₃), where 𝜖_𝑖 is either 1 or −1, 𝐻(𝐷𝐹(𝑋))𝐻 has nonpositive off-diagonal elements for 𝑋∈𝐷, where 𝐷𝐹(𝑋) is the Jacobian of (23). It is shown in [16] that if 𝐷 is convex the flow of such a system preserves for 𝑡<0 the partial order in ℝ³ defined by the orthant 𝐾₁=𝑋₁,𝑋₂,𝑋₃∈ℝ³∣𝜖_𝑖𝑋_𝑖≥0. (24)
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 𝑥=𝑝(𝑡) with minimal period 𝜔>0 and orbit 𝛾={𝑝(𝑡)∶0≤𝑡≤𝜔}. This orbit is orbitally stable if and only if, for each 𝜀>0, there exists a 𝛿₁>0 such that any solution 𝑥(𝑡), for which the distance of 𝑥(0) from 𝛾 is less 𝛿, remains at a distance less than 𝜀 from 𝛾 for all 𝑡≥0. It is asymptotically orbitally stable, if the distance of 𝑥(𝑡) from 𝛾 also tends to 0 as 𝑡 goes to ∞. This orbit 𝛾 is asymptotically orbitally stable with asymptotic phase if it is asymptotically orbitally stable and there exists a 𝑏>0, such that, any solution 𝑥(𝑡), for which the distance of 𝑥(0) from 𝛾 is less than 𝑏, satisfying |𝑥(𝑡)−𝑝(𝑡−𝜏)|→0 as 𝑡→+∞ for some 𝜏 which may depend on 𝑥(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 𝛾(𝑡), 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 that 𝑛=3 and 𝐷 is convex and bounded. Suppose that (4) is competitive and persistent and has the property of stability of the periodic orbits. If 𝑥₀ is the only equilibrium in int (𝐷) and if it is locally asymptotically stable, then it is globally asymptotically stable in int (𝐷).

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 of 𝐷, the system (4) has only one 𝜔-limit point which is the equilibrium 𝐸₀. Moreover, for ℛ₀>1, 𝐸₀ cannot be the 𝜔-limit of any orbit in int(𝐷).

Proof. The vector field is transversal to the boundary of 𝐷 except in the 𝑆-axis, which is invariant with respect to (4). On the 𝑆-axis we have 𝑑𝑆𝑑𝑡=𝑏−𝛼𝑆, (25) which implies that 𝑆→𝑏/𝛼 as 𝑡→∞. Therefore, 𝐸₀ is the only 𝜔-limit point on the boundary of 𝐷.
To prove the second part of the position, we consider the function 𝑉₁=𝑌+𝛼𝑎𝑏𝛽1+ℛ²₀2𝐼, (26) the derivative of which along solutions is given by 𝑑𝑉₁=𝑑𝑡𝑑𝑌+𝑑𝑡𝛼𝑎𝑏𝛽1+ℛ²₀2𝑑𝐼=𝑘𝑑𝑡𝑎𝑆+𝐼𝜇𝐼−𝑘𝑌𝐼𝑆+𝐼−𝜇𝑌+𝛼𝜇𝑏1+ℛ²₀2−𝛽𝑌𝑆𝛼𝑎𝛿𝑏𝛽1+ℛ²₀2𝐼≥𝑘𝑎()𝜇−𝑏/𝛼𝛼𝑎𝛿𝑏𝛽1+ℛ²₀2𝐼+𝛼𝜇𝑏1+ℛ²₀2𝑘𝑆−𝑌≥1𝑏/𝛿𝐼−𝜇1−ℛ²₀1+ℛ²₀2𝑏𝜇𝐼+𝑘𝑎𝛼𝑆−𝑘𝛿2𝛼𝜇1+ℛ²₀2𝐼−1+ℛ²₀𝑏𝛼𝑏2𝛼𝑏1+ℛ²₀𝑌. (27) Since ℛ₀>1, then ((1/ℛ²₀)+1)/2<1 and 2/(1+ℛ²₀)<1. Therefore, there exists a neighborhood 𝑈 of 𝐸₀ such that for (𝑆,𝐼,𝑌)∈𝑈∪int(𝐷) the expression inside the brackets is positive. In this neighborhood, we have 𝑉₁>0 unless 𝐼=𝑌=0. Moreover, the level sets of 𝑉₁ are the planes 𝑌+𝛼𝑎𝑏𝛽1+ℛ²₀2𝐼=𝐶, (28) which go away from the 𝑆-axis as 𝐶 increase. Since 𝑉₁ increases along the orbits starting in 𝑈∪int(𝐷), we conclude that they go away from 𝐸₀. This proves the proposition and therefore the persistence of system (4) when ℛ₀>1.
By looking at the Jacobian matrix of the system (4) and choosing the matrix 𝐻 as 𝐻=−10001000−1, (29) we can see that the system (4) is competitive in 𝐷, with respect to the partial order defined by the orthant 𝐾₁={(𝑆,𝐼,𝑌)∈ℝ³∶𝑆≤0, 𝐼≤0, 𝑌≤0}. Our main results will follow from this observation and the above theorems.

Theorem 9. If ℛ₀>1, then the positive equilibrium 𝐸^∗ of the system (4) is globally asymptotically stable.

Proof. Since the system (4) is competitive permanent if ℛ₀>1 holds and the only equilibrium point 𝐸^∗ 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.
Let 𝑝(𝑡)=(𝑆(𝑡),𝐼(𝑡),𝑌(𝑡)) be a periodic solution whose orbit 𝛾 is contained in int(𝑅³₊) and suppose that its minimal period is 𝜔>0. The second compound equation is the following periodic linear system: 𝑑𝑍(𝑡)=𝑑𝑡𝜕𝑓^[2]𝑝(𝑡)𝑍(𝑡),𝜕𝑥 (30) where 𝑍=(𝑍₁,𝑍₂,𝑍₃)^⊤ and 𝜕𝑓/𝜕𝑥 is derived from the Jacobian matrix of the system (4) and defined as follows: 𝜕𝑓^[2]=−𝜇𝜕𝑥𝑎𝜇𝛽𝑌−𝛼−𝛿𝑎𝜇𝛽𝑆𝑎−𝑘𝐼𝛽𝑆0𝑎/𝜇−𝑌(𝑆+𝐼)²+𝑘𝑎/𝜇−𝑌−𝜇𝑆+𝐼𝑎𝛽𝑌−𝛼−𝑘𝐼𝑘𝐼𝑆+𝐼−𝜇0𝑎/𝜇−𝑌(𝑆+𝐼)²𝜇𝑎𝛽𝑌−𝛿−𝑘𝐼𝑆+𝐼−𝜇. (31)For the solution 𝑝(𝑡), (30) becomes 𝑍̇₁(𝑡)=−𝜇𝑎𝛽𝑌+𝛼+𝛿𝑍₁+𝜇𝑎𝛽𝑆𝑍₂+𝑍₃,𝑍̇₂(𝑡)=−𝑘𝑎/𝜇−𝑌𝐼(𝑆+𝐼)²+𝑘𝑎/𝜇−𝑌𝑍𝑆+𝐼¹−𝜇𝑎𝛽𝑌+𝛼+𝑘𝐼𝑍𝑆+𝐼+𝜇₂,𝑍̇₃(𝑡)=𝑘𝑎/𝜇−𝑌𝐼(𝑆+𝐼)²𝑍₁+𝜇𝑎𝛽𝑌𝑍₂−𝛿+𝑘𝐼𝑍𝑆+𝐼+𝜇₃. (32)
To prove that (32) is globally asymptotically stable, we shall use the following Lyapunov function: 𝐿𝑍₁,𝑍₂,𝑍₃(𝑡)(𝑡)(𝑡)𝑍;𝑆,𝐼,𝑌=sup₁,𝐼𝑌𝑍₂+𝑍₃. (33)
Function (33) is positive, but not differentiable everywhere. Fortunately, this lack of differentiability can be remedied by using the right derivative of 𝐿(𝑡), denoted as 𝐷₊𝐿(𝑡). Then we have the following equalities: 𝑍̇₁(𝑡)≤−𝜇𝑎𝛽𝑌+𝛼+𝛿𝑍₁+𝜇𝑎𝛽𝑆𝑍₂+𝑍₃𝜇≤−𝑎𝛽𝑌+𝛼+𝛿𝑍₁+𝜇𝑎𝑌𝛽𝑆𝐼𝐼𝑌𝑍₂+𝑍₃,𝑍̇₂(𝑡)≤−𝑘𝑎/𝜇−𝑌𝐼(𝑆+𝐼)²+𝑘𝑎/𝜇−𝑌𝑍𝑆+𝐼¹−𝜇𝑎𝛽𝑌+𝛼+𝑘𝐼𝑍𝑆+𝐼+𝜇₂,𝑍̇₃(𝑡)≤𝑘𝑎/𝜇−𝑌𝐼(𝑆+𝐼)²𝑍₁+𝜇𝑎𝛽𝑌𝑍₂−𝛿+𝑘𝐼𝑍𝑆+𝐼+𝜇₃. (34)
Therefore, 𝐷₊𝐼𝑌𝑍₂(𝑡)+𝑍₃(𝑡=𝐼̇)𝐼−𝑌̇𝑌𝐼𝑌𝑍₂+𝑍₃+𝐼𝑌𝐷₊𝑍₂+𝑍₃≤𝐼̇𝐼−𝑉̇𝑉𝐼𝑌𝑍₂+𝑍₃+𝐼𝑌𝑘𝑎/𝜇−𝑌𝑍𝑆+𝐼¹−𝛼+𝑘𝐼𝑍𝑆+𝐼+𝜇₂+𝑍₃≤𝐼̇𝐼−𝑌̇𝑌−𝛼−𝑘𝐼𝑍𝑆+𝐼−𝜇₂+𝑍₃+𝐼𝑌𝑘𝑎/𝜇−𝑌𝑍𝑆+𝐼¹. (35)

Define 𝑔₁(𝑡)𝜇=−𝑎+𝜇𝛽𝑌+𝛼+𝛿𝑎𝑌𝛽𝑆𝐼,𝑔₂(𝑡)=𝐼𝑌𝑘𝑎/𝜇−𝑌+𝐼̇𝑆+𝐼𝐼−𝑌̇𝑌−𝛼−𝑘𝐼.𝑆+𝐼−𝜇 (36) Rewriting the last two equations of (4) as 𝐼̇𝐼𝜇+𝛿=𝑎𝑌𝛽𝑆𝐼,𝑌̇𝑌+𝑘𝐼𝐼𝑆+𝐼+𝜇=𝑌𝑘𝑎()𝜇,𝑆+𝐼 (37) then 𝑔₁(𝑡)𝜇=−𝑎+𝜇𝛽𝑌+𝛼+𝛿𝑎𝑌𝛽𝑆𝐼=𝐼̇𝐼−𝜇𝑎𝑔𝛽𝑌−𝛼,₂(𝑡)=𝐼𝑌𝑘𝑎/𝜇−𝑌+𝐼̇𝑆+𝐼𝐼−𝑌̇𝑌−𝛼−𝑘𝐼=𝐼̇𝑆+𝐼−𝜇𝐼−𝑘𝐼𝑆+𝐼−𝛼. (38)

Thus, we obtain 𝐷₊𝐿(𝑡)≤sup𝑔₁(𝑡),𝑔₂(𝑡)𝐿(𝑡). (39) It follows from (38) that 𝑔₁(𝑡)≤(𝐼/̇𝐼)−𝛼, and thus, from (39) and Gronwall’s inequality, we obtain 𝐿(𝑡)≤𝐿(0)𝐼𝑒^−𝛼𝑡. (40) Since 𝑏/𝛿≤𝐼≤𝑏/𝛼 in int(𝐷), which implies that 𝐿(𝑡)→0 as 𝑡→∞, therefore, (𝑍₁(𝑡),𝑍₂(𝑡), 𝑍₃(𝑡))→(0,0,0) as 𝑡→∞. 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 ℛ₀ which determines the spread of infection has been worked out.

We obtain the basic reproduction number of the system (4), ℛ₀=𝑘𝛽/𝛿𝜇. Furthermore, we obtain that the disease-free equilibrium 𝐸₀(𝑆₀,0,0) of the system (4) is globally asymptotically stable when ℛ₀<ℛ_𝑐<1 (see Figure 2(c)). And if ℛ₀>1, the endemic equilibrium 𝐸^∗(𝑆^∗,𝐼^∗,𝑌^∗) is globally asymptotically stable (see Figure 2(b)). That is to say, if ℛ₀ and ℛ_𝑐 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 ℛ_𝑐>1, while ℛ₀<1, though the pine wilt disease could occur, the outbreak will stop since the infected pine decreases to the disease-free equilibrium. If ℛ₀ 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 𝜕ℛ_𝑐/𝜕𝜇<0. Hence, we can make ℛ_𝑐<1 by increasing 𝜇. Also from Figure 2, we can find that ℛ_𝑐=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 ℛ_𝑐=0.5477<1 in Figure 2(c), keeping all the other parameters the same as in Figure 2(b). Thus, the disease-free equilibrium 𝐸₀(𝑆₀,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 ℛ₀<1. So we only need to lower ℛ₀ to less than 1 to control the pine wilt disease. We can also find 𝜕ℛ₀/𝜕𝛿<0 and 𝜕ℛ₀/𝜕𝜇<0 by a simple calculation. Hence, we can make ℛ₀<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 ℛ₀<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 Theorem 8, we can know that the endemic equilibrium 𝐸^∗ is globally asymptotically stable when ℛ₀>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 𝜕𝐼^∗/𝜕𝛿=−𝑏[(𝜇𝛿−𝛽𝑘)²+𝛽𝑘²𝛼+𝛽𝑘𝜇𝛼]/(𝛽𝑘+𝑘𝛼−𝜇𝛿+𝜇𝛼)²𝛿²<0 and 𝜕𝐼^∗/𝜕𝜇=−𝑏𝑘𝛼(𝛿+𝛽)/(𝛽𝑘+𝑘𝛼−𝜇𝛿+𝜇𝛼)²𝛿<0. 𝐼^∗ is a decreasing function of parameter 𝛿 (or 𝜇). Hence, we can decrease 𝐼^∗ by increasing 𝛿 (or 𝜇) properly. In this case, we need not increase 𝛿 (or 𝜇) large enough to make ℛ₀<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 𝐼^∗. Keeping all the other parameters the same as in Figure 2(b), we plot 𝐼(𝑡) for different 𝜇 (𝜇=0.0005,0.001,0.0015,0.002,0.0025). In Figure 3(b), we illustrate the effect of varying 𝛿 on 𝐼^∗. Keeping all the other parameters the same as in Figure 2(b), we plot 𝐼(𝑡) for different 𝛿 (𝛿=0.008,0.016,0.024,0.032,0.04). The numerical simulations show that the infected pine 𝐼(𝑡) approaches to the steady state 𝐼^∗, and 𝐼^∗ 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)

(b)

Figure 3 

(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).

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11171284 and 11071011), program for Innovative Research Team (in Science and Technology) in University of Henan Province (2010IRTSTHN006), Innovation Scientists and Technicians Troop Construction Projects of Henan Province (104200510011), Natural Science Foundation of the Education Department of Henan Province (2011B110028), and Key Project for the Education Department of Henan Province (13A110771).