Analysis of Pine Wilt Disease Model with Nonlinear Incidence and Horizontal Transmission

Abstract

The deterministic pine wilt model with vital dynamics to determine the equilibria and their stability by considering nonlinear incidence rates with horizontal transmission is analyzed. The complete global analysis for the equilibria of the model is discussed. The explicit formula for the reproductive number is obtained and it is shown that the “disease-free” equilibrium always exists and is globally asymptotically stable whenever 𝑅₀≤1. Furthermore, the disease persists at an “endemic” level when the reproductive number exceeds unity.

1. Introduction

Pine wilt, a fatal disease of commonly planted pines brought on by the pinewood nematode (Bursaphelenchus xylophilus), causes changes to ecosystem and destructs the variety of ecosystem. Pine trees affected by pine wilt disease usually die within few months. Symptoms of pine wilt disease normally appear in late spring or summer. The most prominent symptom is the lack of resin exudation from barks wounds. The foliage becomes light grayish green, then becomes yellow, and finally it becomes reddish brown. The tree succumbs to the disease at this stage. The affected trees totally lack resin and their wood becomes dry.

The long-horned pine sawyer beetles (Monochamus alternatus) are the main culprits for the spread of pinewood nematodes from infected pines to healthy or stressed pines. When new adult beetles emerge in spring, they locate a living host tree to feed on the bark of the young branches and transfer nematodes to the healthy trees through the feeding wounds produced by these sawyers. This transmission is referred to as primary transmission. The transmission of the nematodes during egg-laying activities in freshly cut timber or dying trees is referred to as secondary transmission. Nematodes, introduced during primary transmission, migrate to the resin canals of their hosts and kill these cells rendering them ineffective due to which a susceptible host can wilt and die within weeks of being infested upon the availability of favorable conditions to disease development. The principle of the Bursaphelenchus xylophilus transmission and disease dissemination is reviewed by Evans et al. [1]. Pine wilt particularly kills Scots pine within few weeks to few months. Some other pine species as Austrian (Pinus nigra), jack (P. banksiana), mogo (P. mugo), and red (P. resinosa) pines are occasionally killed by pine wilt.

Mathematical modeling became a considerably important tool in the study of epidemiology because it helps us to understand the observed epidemiological patterns and disease control and provides understanding of the underlying mechanisms which influence the spread of disease and may suggest control strategies. The incidence of a disease is defined as new cases occurring per unit time. It plays a vital role in mathematical epidemiology. The classical epidemiological models are developed by the assumption of bilinear incidence rate 𝛽𝑆𝐼 and the standard incidence rate 𝛽𝑆𝐼/𝑁, where 𝛽 is the transmission probability per contact and 𝑆,𝐼 are susceptible and infected individuals, respectively. However, there are several reasons that require modification in these incidence rates. For example, the assumption of homogeneous mixing may be invalid and in this case a model having a particular form of nonlinear transmission may be incorporated for heterogeneous mixing and the necessary population structure.

The saturation effects may also require a nonlinear incidence rate because if the proportion of the infected population is high enough so that exposure to the disease agent is virtually certain, then transmission rate may react more slowly than linear in order to increase the number of infected individuals. Capasso and Serio [2] who studied the cholera epidemic spread in Bari in 1973 and introduced the saturated incidence rate 𝑔(𝐼)𝑆 in the epidemic model observed this effect. A variety of nonlinear incidence rates has been utilized in epidemic models [3–7]. An epidemic model with nonlinear incidence is proposed in [8], in which the authors described the dynamics of diseases spread by vector mosquitoes such as malaria, dengue, and yellow fever.

The incidence of pine wilt disease depends on beetles’ density because pine sawyer beetles are the source of transmission of pinewood nematode. This incidence may approach its saturation level at very high beetle densities. The adult female pine sawyer attempts to avert from erstwhile oviposition scants. It approaches another tree before the saturation point of oviposition is reached. Thus the isolation of infected individuals results in the decrease in the number of contacts between the susceptible and infected individuals at high infective levels. These observations inspire to consider nonlinearities in the incidence rates.

In this paper, based on the ideas posed in [9–11], a pine wilt disease model considering a direct mode of transmission as well as nonlinear incidence rate is formulated. The aim of this paper is to establish stability properties of equilibria and the threshold parameter 𝑅₀ that completely determines the existence of endemic or disease-free equilibrium. If 𝑅₀≤1, the disease-free equilibrium is globally asymptotically stable. If 𝑅₀>1, a unique endemic equilibrium exists and is globally asymptotically stable. The rest of the paper is organized as follows. In Section 2, the description of the extended mathematical model is presented. Section 3 is devoted to the existence of equilibria. In Section 4, the Lyapunov function theory is used to show global stability of disease-free equilibrium and geometric approach is used to prove global stability of endemic equilibrium in Section 5. Discussions and simulations are done in Section 6.

2. Model Description

The pine population, with total population size denoted by 𝑁_ℎ(𝑡), is subdivided into two mutually exclusive compartments: susceptible pine trees 𝑆_ℎ(𝑡) and infectious pine trees 𝐼_ℎ(𝑡). Thus, 𝑁_ℎ(𝑡)=𝑆_ℎ(𝑡)+𝐼_ℎ(𝑡). The emission of oleoresin from susceptible host pines behaves like a physical barrier for beetle oviposition. Beetles can oviposit on the infected pine trees because these trees cease oleoresin. Since there are no cures for pine wilt once a susceptible tree becomes infested with pinewood nematodes, the recovered class 𝑅_ℎ(𝑡) has not been considered.

The total vector population at any time 𝑡 is denoted by 𝑁_𝑣(𝑡)=𝑆_𝑣(𝑡)+𝐼_𝑣(𝑡), where 𝑆_𝑣(𝑡) denotes the susceptible adult beetles that do not have any pinewood nematode at time 𝑡 and 𝐼_𝑣(𝑡) denotes the infected adult beetles carrying pinewood nematode at time 𝑡 when they emerge from dead pine trees. After emergence from the dead tree, beetles choose a healthy tree for sufficient feeding and transmit nematodes into the tree. These nematodes move through the feeding wounds and approach the xylem of the tree. When beetles are in oviposition they choose dying or dead tree and transmit nematode when they lay eggs in slits in bark. Nematodes enter these slits, feed on wood cells or fungi, and reproduce themselves. Before beetle’s emergence from dead tree the nematodes attach with the tracheae of its respiratory system. The following assumptions are made in formulating the mathematical model.(i)The exploitation rate of pine trees infected with Bursaphelenchus xylophilus is greater than the normal and susceptible pine trees.(ii)The susceptible beetles receive nematodes directly from infectious ones through mating.(iii)Adult beetles emerging from infected trees have pinewood nematode.(iv)The infected vectors transmit the nematode during maturation feeding as well as via oviposition. Under these assumptions, the vector-host model with nonlinear incidence can be described by the following system of differential equations: 𝑑𝑆_ℎ𝑑𝑡=Π_ℎ−𝛿₁𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝛿₂𝜃𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝜇_ℎ𝑆_ℎ,𝑑𝐼_ℎ=𝛿𝑑𝑡¹𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣+𝛿₂𝜃𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝜔𝐼_ℎ,𝑑𝑆_𝑣𝑑𝑡=Π_𝑣−𝛽₁𝑆_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ−𝛽₂𝑆_𝑣𝐼_𝑣−𝜇_𝑣𝑆_𝑣,𝑑𝐼_𝑣=𝛽𝑑𝑡¹𝑆_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ+𝛽₂𝑆_𝑣𝐼_𝑣−𝜇_𝑣𝐼_𝑣, (1) where Π_ℎ is the constant increase rate of pines, Π_𝑣 is the constant input rate of vectors, and 𝜇_𝑣 is the mortality rate of vectors. The exploitation rate of susceptible pines is 𝜇_ℎ whereas the isolation and felling rate of infected pines is 𝜔. 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 𝛿₂𝜃𝑆_ℎ𝐼_𝑣/(1+𝛼₁𝐼_𝑣) and 𝛿₁𝑆_ℎ𝐼_𝑣/(1+𝛼₁𝐼_𝑣), respectively. The parameter 𝜃 is the probability by which susceptible pines die of natural causes and cease oleoresin exudation without being infected by the nematode, and 𝛿₂ indicates the rate at which infected vectors transmit the nematode via oviposition whereas 𝛿₁ denotes transmission rate per contact during maturation feeding. The transmission between susceptible vectors and infected hosts occurs when adult beetles emerge from dead pine trees. This transmission is denoted by 𝛽₁𝑆_𝑣𝐼_ℎ/(1+𝛼₂𝐼_ℎ), where 𝛽₁ is the rate at which adult beetles carry the pinewood nematode when they emerge from dead trees. The parameters 𝛼₁ and 𝛼₂ determine the level at which the infection is saturated. The beetles transmit nematodes directly through mating. The incidence term for this transmission is 𝛽₂𝑆_𝑣𝐼_𝑣, where 𝛽₂ is the transmission rate among beetles during mating. All parameters are assumed to be positive.

The total dynamics of vector population satisfy the following equation: 𝑑𝑁_𝑣𝑑𝑡=Π_𝑣−𝜇_𝑣𝑁_𝑣. (2) This leads to 𝑁_𝑣→Π_𝑣/𝜇_𝑣 as 𝑡→∞. Thus, the system (1) is reduced to the following system of differential equations: 𝑑𝑆_ℎ𝑑𝑡=Π_ℎ−𝛿₁𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝛿₂𝜃𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝜇_ℎ𝑆_ℎ,𝑑𝐼_ℎ=𝛿𝑑𝑡¹𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣+𝛿₂𝜃𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝜔𝐼_ℎ,𝑑𝐼_𝑣𝑑𝑡=𝛽₁Π_𝑣𝜇_𝑣−𝐼_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ+𝛽₂Π_𝑣𝜇_𝑣−𝐼_𝑣𝐼_𝑣−𝜇_𝑣𝐼_𝑣. (3) Considering ecological significance, we study system (3) in the closed set Ω={(𝑆_ℎ,𝐼_ℎ,𝐼_𝑣)∶Π_ℎ/𝜔≤𝑆_ℎ+𝐼_ℎ≤Π_ℎ/𝜇_ℎ,0≤𝐼_𝑣≤Π_𝑣/𝜇_𝑣}. It can be easily verified that Ω is positively invariant with respect to (3).

3. Existence of Equilibria

The dynamics of the disease are described by the threshold quantity 𝑅₀ which is called the reproduction number defined as “the average number of secondary infections produced by an infected individual in a completely susceptible population.” It is one of the most useful threshold parameters that characterizes mathematical problems related to infectious diseases. This metric helps to determine whether or not an infectious disease will spread through a population. The basic reproduction number of model (3) is given by 𝑅₀=𝛽₂Π_𝑣𝜇²_𝑣+𝛽₁Π_𝑣𝜇²_𝑣Π_ℎ𝜇_ℎ𝜔𝛿₁+𝜃𝛿₂. (4) Direct calculation shows that for 𝑅₀≤1, there is only disease-free equilibrium 𝐸₀(Π_ℎ/𝜇_ℎ,0,0) and for 𝑅₀>1, there is an additional equilibrium 𝐸^∗(𝑆^∗_ℎ,𝐼^∗_ℎ,𝐼^∗_𝑣) which is called endemic equilibrium, with 𝑆^∗_ℎ=Π_ℎ−𝜔𝐼^∗_ℎ𝜇_ℎ,𝐼^∗_ℎ=Π_ℎ𝐼^∗_𝑣𝛿₁+𝛿₂𝜃+𝛿₁𝛼₁+𝛿₂𝜃𝛼₁𝐼^∗_𝑣×𝛼₁𝜇_ℎ+𝛿₁+𝛿₂𝜃𝛼₁𝜔𝐼_𝑣^∗2+𝛿₁+𝛿₂𝜃+2𝛼₁𝜇_ℎ𝜔𝐼^∗_𝑣+𝜔𝜇_ℎ⁻¹, (5) and 𝐼^∗_𝑣 is the root of the following equation: 𝐴𝐼_𝑣^∗3+𝐵𝐼_𝑣^∗2+𝐶𝐼^∗_𝑣+𝐷=0, (6) where, 𝐴=Π_ℎ𝜇_𝑣𝛼₂𝛽₂𝛼₁𝛿₁+𝜃𝛼₁𝛿₂,𝐵=𝜃𝜔𝛼₁𝛽₂𝛿₂𝜇_𝑣+𝜃𝛼₁𝛽₁𝛿₂Π_ℎ𝜇_𝑣+𝜃𝛼₂𝛽₂𝛿₂Π_ℎ𝜇_𝑣+𝜔𝛼₁𝛼₁𝛽₂𝜇_ℎ𝜇_𝑣,𝐶=𝜔𝛽₂𝜇_ℎ𝜇_𝑣+𝛼₁𝛽₁𝛿₁Π_ℎΠ_𝑣+𝜃𝛼₁𝛽₁𝛿₂Π_ℎΠ_𝑣,𝐷=𝜔𝜇_ℎ𝜇²_𝑣1−𝑅₀. (7) From (7), we see that 𝑅₀>1 if and only if 𝐷<0. Since 𝐴, 𝐵, and 𝐶 are always positive, there will be zero or unique positive endemic equilibrium accordingly as 𝑅₀≤1 or 𝑅₀>1. Thus we have the following theorem.

Theorem 1. System (3) always has the infection-free equilibrium 𝐸₀. If 𝑅₀>1, system (3) has a unique endemic equilibrium 𝐸^∗(𝑆^∗_ℎ,𝐼^∗_ℎ,𝐼^∗_𝑣) defined by (5) and (6).

4. Stability of Disease-Free Equilibrium

Here, we analyze stability of disease-free equilibrium 𝐸₀(Π_ℎ/𝜇_ℎ,0,0) for system (3). The linearization of the system (3) at 𝐸₀ results in the following characteristic equation: −𝜇_ℎ−𝜆𝜆²+𝜆𝜔+𝜇_𝑣−𝛽₂Π_𝑣𝜇_𝑣+𝜔𝜇_𝑣1−𝑅₀=0. (8) The characteristic equation (8) has one eigenvalue −𝜇_ℎ. The other eigenvalues can be found by the equation 𝜆²+𝑎𝜆+𝑏=0, (9) where 𝑎=𝜔+𝜇_𝑣−(𝛽₂Π_𝑣/𝜇_𝑣) and 𝑏=𝜔𝜇_𝑣(1−𝑅₀).

We observe that the roots of the quadratic equation (9) have negative real parts if 𝑅₀<1. If 𝑅₀=1, one root of (9) is 0. This fact does not guarantee that all eigenvalues have negative real parts. It will only be possible in case of real roots. If 𝑅₀>1, one of the roots of (9) has positive real part. The above discussion leads to the following theorem.

Theorem 2. The disease-free equilibrium of system (3) is locally asymptotically stable in Ω if 𝑅₀<1 and it is unstable if 𝑅₀>1.

Now, we analyze the global behavior of the disease-free equilibrium 𝐸₀. The following theorem provides the global property of the system.

Theorem 3. If 𝑅₀≤1, then the infection-free equilibrium 𝐸₀ is globally asymptotically stable in the interior of Ω.

Proof. The following Lyapunov function is proposed to establish the global stability of disease-free equilibrium: 𝐿=𝛽₁Π_𝑣𝜔𝐼_ℎ+𝜇_𝑣𝐼_𝑣. (10) Taking the time derivative of 𝐿 along the solutions of (3), we have 𝐿^′=𝛽₁Π_𝑣𝜔𝐼^′_ℎ+𝜇_𝑣𝐼^′_𝑣=𝛽₁Π_𝑣𝜔𝛿₁𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣+𝛿₂𝜃𝑆_ℎ𝐼_𝑣1+𝛼₁𝐼_𝑣−𝜔𝐼_ℎ+𝜇_𝑣𝛽₁Π_𝑣𝜇_𝑣−𝐼_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ+𝛽₂Π_𝑣𝜇_𝑣−𝐼_𝑣𝐼_𝑣−𝜇_𝑣𝐼_𝑣≤𝛽₁Π_𝑣𝜔𝛿₁𝑆_ℎ𝐼_𝑣+𝛿₂𝜃𝑆_ℎ𝐼_𝑣−𝜔𝐼_ℎ+𝜇_𝑣𝛽₁Π_𝑣𝜇_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ−𝛽₁𝐼_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ+𝛽₂Π_𝑣𝜇_𝑣𝐼_𝑣−𝛽₂𝐼_𝑣𝐼_𝑣−𝜇²_𝑣𝐼_𝑣<𝛿₁+𝛿₂𝜃𝛽₁Π_𝑣𝜔Π_ℎ𝜇_ℎ𝐼_𝑣+𝜇_𝑣𝛽₂Π_𝑣𝜇_𝑣𝐼_𝑣−𝜇²_𝑣𝐼_𝑣−𝛽₁Π_𝑣𝐼_ℎ+𝛽₁Π_𝑣𝐼_ℎ−𝜇_𝑣𝛽₁𝐼_𝑣𝐼_ℎ1+𝛼₂𝐼_ℎ−𝜇_𝑣𝛽₂𝐼²_𝑣=𝐼_𝑣𝜇²_𝑣𝑅₀−1−𝜇_𝑣𝛽₁𝐼_ℎ1+𝛼₂𝐼_ℎ−𝜇_𝑣𝛽₂𝐼_𝑣≤0. (11) Thus 𝐿^′(𝑡) is negative if 𝑅₀≤1. When 𝑅₀<1, the derivative 𝐿^′=0 if and only if 𝐼_𝑣=0, while in the case 𝑅₀=1, the derivative 𝐿^′=0 if and only if 𝐼_ℎ=0 or 𝐼_𝑣=0. Consequently, the largest compact invariant set in {(𝑆_ℎ,𝐼_ℎ,𝐼_𝑣∈Ω),𝐿^′=0}, when 𝑅₀≤1, is the singleton 𝐸₀. Hence, by LaSalle’s invariance principle [12], 𝐸₀ is globally asymptotically stable in Ω. This completes the proof.

5. Stability of Endemic Equilibrium

In this section, we will discuss global stability of endemic equilibrium 𝐸^∗ in the feasible region Ω. This is done through the geometrical approach applied by Li and Muldowney [13]. We summarize this approach below.

Consider a 𝐶¹ map 𝑓∶𝑥↦𝑓(𝑥) from an open set 𝐷⊂𝑅^𝑛 to 𝑅^𝑛 such that each solution 𝑥(𝑡,𝑥₀) to the differential equation 𝑥^′=𝑓(𝑥) (12) is uniquely determined by the initial value 𝑥(0,𝑥₀). We have the following assumptions:(𝐻₁)𝐷 is simply connected;(𝐻₂)there exists a compact absorbing set 𝐾⊂𝐷;(𝐻₃)Equation (12) has unique equilibrium 𝑥 in 𝐷.Let 𝑃∶𝑥↦𝑃(𝑥) be a nonsingular (^𝑛₂) ×(^𝑛₂) matrix-valued function which is 𝐶¹ in 𝐷 and a vector norm |⋅| on 𝑅^𝑁, where 𝑁=(^𝑛₂). Let 𝜇 be the Lozinskiĭ measure with respect to the |⋅|. Define a quantity 𝑞₂ as 𝑞₂=limsup_𝑡→∞sup_𝑥0∈𝐾1𝑡^𝑡₀𝜇𝐵𝑥𝑠,𝑥₀𝑑𝑠, (13) where 𝐵=𝑃_𝑓𝑃⁻¹+𝑃𝐽^[2]𝑃⁻¹, the matrix 𝑃_𝑓 is obtained by replacing each entry 𝑝 of 𝑃 by its derivative in the direction of 𝑓, (𝑝_𝑖𝑗)_𝑓, and 𝐽^[2] is the second additive compound matrix of the Jacobian matrix 𝐽 of (12). The following result has been established by Li and Muldowney [13].

Theorem 4. Suppose that 𝐻₁, 𝐻₂, and 𝐻₃ hold; the unique endemic equilibrium 𝐸^∗ is globally stable in Ω if 𝑞₂<0.

Obviously Ω is simply connected and 𝐸^∗ is a unique endemic equilibrium for 𝑅₀>1 in Ω. To apply the result of the above theorem for global stability of endemic equilibrium 𝐸^∗, we first prove the uniform persistence of (3) when the threshold parameter 𝑅₀>1, by applying the acyclicity theorem (see [14]).

Definition 5 (see [15]). The system (3) is uniformly persistent; that is, there exists 𝑐>0 (independent of initial conditions), such that liminf_𝑡→∞𝑆_ℎ≥𝑐, liminf_𝑡→∞𝐼_ℎ≥𝑐, and liminf_𝑡→∞𝐼_𝑣≥𝑐.

Let 𝑋 be a locally compact metric space with metric 𝑑 and let Ω be a closed nonempty subset of 𝑋 with boundary Ω and interior Ω^∘. Clearly, Ω^∘ is a closed subset of Ω. Let Φ_𝑡 be a dynamical system defined on Ω. A set 𝐵 in 𝑋 is said to be invariant if Φ(𝐵,𝑡) =𝐵. Define 𝑀_𝜕∶={𝑥∈Ω∶Φ(𝑡,𝑥)∈Ω, for all 𝑡≥0}.

Lemma 6 (see [14]). Assume that (a)Φ_𝑡 has a global attractor;(b) there exists 𝑀={𝑀₁,…,𝑀_𝑘} of pairwise disjoint, compact, and isolated invariant set on 𝜕Ω such that(1)⋃_{𝑥∈𝑀𝜕}𝜔(𝑥)⊆⋃^𝑘_𝑗=1𝑀_𝑗;(2)no subsets of 𝑀 form a cycle on 𝜕Ω;(3)each 𝑀_𝑗 is also isolated in Ω;(4)𝑊^𝑠(𝑀_𝑗)∩Ω^∘=𝜙 for each 1≤𝑗≤𝑘, where 𝑊^𝑠(𝑀_𝑗) is stable manifold of 𝑀_𝑗. Then Φ_𝑡 is uniformly persistent with respect to Ω.

Proof. We have Ω={(𝑆_ℎ,𝐼_ℎ,𝐼_𝑣)∶Π_ℎ/𝜇_ℎ≤𝑆_ℎ+𝐼_ℎ≤Π_ℎ/𝜔,0≤𝐼_𝑣≤Π_𝑣/𝜇_𝑣}, Ω^∘={(𝑆_ℎ,𝐼_ℎ,𝐼_𝑣)∶𝑆_ℎ,𝐼_ℎ,𝐼_𝑣>0}, 𝜕Ω=Ω/Ω^∘. Obviously, 𝑀_𝜕=𝜕Ω. Since Ω is bounded and positively invariant there exists a compact set 𝑀 in which all solutions of system (3) initiated in Ω ultimately enter and remain forever. On 𝑆_ℎ-axis we have 𝑆^′_ℎ=Π_ℎ−𝜇_ℎ𝑆_ℎ which means 𝑆_ℎ→Π_ℎ/𝜇_ℎ as 𝑡→∞. Thus 𝐸₀ is the only omega limit point on 𝜕Ω; that is, 𝜔(𝑥)=𝐸₀ for all 𝑥∈𝑀_𝜕. Furthermore 𝑀=𝐸₀ is a covering of Ω=⋃_{𝑥∈𝑀𝜕}𝜔(𝑥) because all solutions initiated on the 𝑆_ℎ-axis converge to 𝐸₀. Also 𝐸₀ is isolated and acyclic. This verifies that hypotheses (1) and (3) hold. When 𝑅₀>1, the “disease-free” equilibrium (DFE) 𝐸₀ is unstable from theorem (3) and also 𝑊^𝑠(𝑀)=𝜕Ω. Hypotheses (4) and (5) hold. There always admits a global attractor due to ultimate boundedness of solutions.

The boundedness of Ω and the above lemma imply that (3) has a compact absorbing set 𝐾⊂Ω [15]. Now we will prove that the quantity 𝑞₂<0. We choose a suitable vector norm |⋅| in 𝑅³ and a 3×3 matrix valued function 𝑃(𝑥)=0𝐼100_ℎ𝐼_𝑣0𝐼00^ℎ𝐼_𝑣. (14) Obviously 𝑃 is 𝐶¹ and nonsingular in the interior of Ω. Linearizing system (3) about an endemic equilibrium 𝐸^∗ gives the following Jacobian matrix: 𝐽=^{−𝜇ℎ−𝐼𝑣𝛿11+𝛼1𝐼𝑣−𝐼𝑣𝜃𝛿21+𝛼1𝐼𝑣0−𝑆ℎ𝛿11+𝛼1𝐼𝑣2−𝑆ℎ𝜃𝛿21+𝛼1𝐼𝑣2𝐼𝑣𝛿1}_{1+𝛼1𝐼𝑣}^{+𝐼𝑣𝜃𝛿2}_{1+𝛼1𝐼𝑣}^{−𝜔𝑆ℎ𝛿1}_{1+𝛼1𝐼𝑣2}^{+𝑆ℎ𝜃𝛿2}_{1+𝛼1𝐼𝑣20𝛽1𝜇𝑣1+𝛼2𝐼ℎ2Π𝑣−𝐼𝑣𝜇𝑣𝛽2𝜇𝑣Π𝑣−2𝐼𝑣𝜇𝑣−𝜇𝑣−𝐼ℎ𝛽11+𝛼2𝐼ℎ}. (15)

The second additive compound matrix of 𝐽(𝐸^∗) is given by𝐽^[2]=^{𝑏11𝑆ℎ𝛿11+𝛼1𝐼𝑣2+𝑆ℎ𝜃𝛿21+𝛼1𝐼𝑣2𝑆ℎ𝛿11+𝛼1𝐼𝑣2+𝑆ℎ𝜃𝛿21+𝛼1𝐼𝑣2𝛽1}_{𝜇𝑣1+𝛼2𝐼ℎ2}^Π_𝑣^−𝐼_𝑣^𝜇_𝑣^𝑏₂₂⁰_{0𝐼𝑣𝛿11+𝛼1𝐼𝑣+𝐼𝑣𝜃𝛿21+𝛼1𝐼𝑣𝑏33}, (16)where 𝑏₁₁=−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣𝑏−𝜔,₂₂=−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣+𝛽₂𝜇_𝑣Π_𝑣−2𝐼_𝑣𝜇_𝑣−𝜇_𝑣−𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ,𝑏₃₃𝛽=−𝜔+²𝜇_𝑣Π_𝑣−2𝐼_𝑣𝜇_𝑣−𝜇_𝑣−𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ. (17) The matrix 𝐵=𝑃_𝑓𝑃⁻¹+𝑃𝐽^[2]𝑃⁻¹ can be written in block form as 𝐵=^𝐵11𝐵12_𝐵21𝐵22, with 𝐵₁₁=−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣𝐵−𝜔,₁₂=𝐼_𝑣𝐼_ℎ𝑆_ℎ𝛿₁1+𝛼₁𝐼_𝑣²+𝑆_ℎ𝜃𝛿₂1+𝛼₁𝐼_𝑣²,𝐼_𝑣𝐼_ℎ𝑆_ℎ𝛿₁1+𝛼₁𝐼_𝑣²+𝑆_ℎ𝜃𝛿₂1+𝛼₁𝐼_𝑣²,𝐵₂₁=𝐼_ℎ𝐼_𝑣𝛽₁𝜇_𝑣1+𝛼₂𝐼_ℎ²Π_𝑣−𝐼_𝑣𝜇_𝑣0,𝐵₂₂=𝑀₂₂+𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣0𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣+𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣𝑀₃₃+𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣, (18) where 𝑀₂₂=−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣+𝛽₂𝜇_𝑣Π_𝑣−2𝐼_𝑣𝜇_𝑣−𝜇_𝑣−𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ,𝑀₃₃𝛽=−𝜔+²𝜇_𝑣Π_𝑣−2𝐼_𝑣𝜇_𝑣−𝜇_𝑣−𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ. (19) Consider the norm in 𝑅³ as |(𝑢,𝑣,𝑤)|=max(|𝑢|,|𝑣|+|𝑤|), where (𝑢,𝑣,𝑤) denotes the vector in 𝑅³. The Lozinskiĭ measure with respect to this norm is defined as 𝜇(𝐵)≤sup(𝑔₁,𝑔₂), where 𝑔₁=𝜇₁𝐵₁₁+𝐵₁₂,𝑔₂=𝜇₁𝐵₂₂+𝐵₂₁. (20) From system (3) we can write 𝐼^′_ℎ𝐼_ℎ=𝐼_𝑣𝐼_ℎ𝑆_ℎ𝛿₁1+𝛼₁𝐼_𝑣+𝑆_ℎ𝜃𝛿₂1+𝛼₁𝐼_𝑣𝐼−𝜔,_′𝑣𝐼_𝑣=𝐼_ℎ𝐼_𝑣𝛽₁𝜇_𝑣1+𝛼₂𝐼_ℎΠ_𝑣−𝜇_𝑣𝐼_𝑣+𝛽₂𝜇_𝑣Π_𝑣−𝜇_𝑣𝐼_𝑣−𝜇_𝑣. (21)

Since 𝐵₁₁ is a scalar, its Lozinskiĭ measure with respect to any vector norm in 𝑅¹ will be equal to 𝐵₁₁. Thus 𝐵₁₁=−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣𝐵−𝜔,₁₂=𝐼_𝑣𝐼_ℎ𝑆_ℎ𝛿₁1+𝛼₁𝐼_𝑣²+𝑆_ℎ𝜃𝛿₂1+𝛼₁𝐼_𝑣², (22) and 𝑔₁ will become 𝑔₁=−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣𝐼−𝜔+^𝑣𝐼_ℎ𝑆_ℎ𝛿₁1+𝛼₁𝐼_𝑣²+𝑆_ℎ𝜃𝛿₂1+𝛼₁𝐼_𝑣²𝑔₁=𝐼^′_ℎ𝐼_ℎ−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣. (23) Also |𝐵₂₁|=(𝐼_ℎ/𝐼_𝑣)(𝛽₁/𝜇_𝑣(1+𝛼₂𝐼_ℎ)²)(Π_𝑣−𝐼_𝑣𝜇_𝑣), |𝐵₁₂| and |𝐵₂₁| are the operator norms of 𝐵₁₂ and 𝐵₂₁ which are mapping from 𝑅² to 𝑅 and from 𝑅 to 𝑅², respectively, and 𝑅² is endowed with the 𝑙₁ norm. 𝜇₁(𝐵₂₂) is the Lozinskiĭ measure of 2×2 matrix 𝐵₂₂ with respect to 𝑙₁ norm in 𝑅². Consider 𝜇₁𝐵₂₂=sup𝑀₂₂+𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣+𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣+𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣,𝑀₃₃+𝐼^′_ℎ𝐼_ℎ−𝐼^′_𝑣𝐼_𝑣=sup−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣+𝛽₂𝜇_𝑣Π_𝑣−2𝐼_𝑣𝜇_𝑣−𝜇_𝑣−𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ+𝐼^′_ℎ𝐼_ℎ−𝐼_ℎ𝐼_𝑣𝛽₁𝜇_𝑣1+𝛼₂𝐼_ℎΠ_𝑣−𝜇_𝑣𝐼_𝑣−𝛽₂𝜇_𝑣Π_𝑣−𝜇_𝑣𝐼_𝑣+𝜇_𝑣+𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣+𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣𝛽−𝜔+²𝜇_𝑣Π_𝑣−2𝐼_𝑣𝜇_𝑣−𝜇_𝑣−𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ+𝐼^′_ℎ𝐼_ℎ−𝐼_ℎ𝐼_𝑣𝛽₁𝜇_𝑣1+𝛼₂𝐼_ℎΠ_𝑣−𝜇_𝑣𝐼_𝑣−𝛽₂𝜇_𝑣Π_𝑣−𝜇_𝑣𝐼_𝑣+𝜇_𝑣,𝜇₁𝐵₂₂=𝐼^′_ℎ𝐼_ℎ−𝛽₂𝜇_𝑣𝐼_𝑣𝜇_𝑣−𝐼_ℎ𝐼_𝑣𝛽₁𝜇_𝑣1+𝛼₂𝐼_ℎΠ_𝑣−𝜇_𝑣𝐼_𝑣−𝛽₁, (24) where 𝛽₁=min𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ+𝜇_ℎ,𝜔+𝐼_ℎ𝛽₁1+𝛼₂𝐼_ℎ. (25) Hence 𝑔₂≤(𝐼^′_ℎ/𝐼_ℎ)−(𝛽₂/𝜇_𝑣)(𝐼_𝑣𝜇_𝑣)−𝛽₁.

Thus 𝜇(𝐵)=sup𝑔₁,𝑔₂𝐼≤sup^′ℎ𝐼_ℎ−𝜇_ℎ−𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣−𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣,𝐼^′_ℎ𝐼_ℎ−𝛽₂𝜇_𝑣𝐼_𝑣𝜇_𝑣−𝛽₁,𝜇(𝐵)≤𝐼^′_ℎ𝐼_ℎ−𝛽₂, (26) where 𝛽₂=min𝜇_ℎ+𝐼_𝑣𝛿₁1+𝛼₁𝐼_𝑣+𝐼_𝑣𝜃𝛿₂1+𝛼₁𝐼_𝑣,𝛽₂𝜇_𝑣𝐼_𝑣𝜇_𝑣+𝛽₁. (27)Since (3) is uniformly persistent when 𝑅₀>1, so for 𝑇>0 such that 𝑡>𝑇 implies 𝐼_ℎ(𝑡)≥𝑐, 𝐼_𝑣(𝑡)≥𝑐 and (1/𝑡)log𝐼_ℎ(𝑡)<𝛽₂/2 for all (𝑆_ℎ(0),𝐼_ℎ(0),𝐼_𝑣(0))∈𝐾.

Thus 1𝑡^𝑡₀𝜇(𝐵)𝑑𝑡<log𝐼_ℎ(𝑡)𝑡−𝛽₂𝛽<−²2, (28) for all (𝑆_ℎ(0),𝐼_ℎ(0),𝐼_𝑣(0))∈𝐾, which further implies that 𝑞₂<0. Therefore all the conditions of Theorem 4 are satisfied. Hence unique endemic equilibrium 𝐸^∗ is globally stable in Ω.

6. Discussions and Simulations

In this paper, pine wilt disease transmission model with nonlinear incidence rates and horizontal transmission is proposed and analyzed. The basic reproduction number, 𝑅₀, of the model is obtained and with the help of this reproduction number the asymptotic behavior of the model is discussed. The variation of infected hosts and infected vectors is shown in Figure 1. It is not meaningful to consider the saturation level when transmission occurred during mating. Thus bilinear incidence has been considered. By simple calculation we see that 𝜕𝑅₀𝜕𝜇_𝑣2=−𝜔𝜇_ℎ𝜇³_𝑣Π_𝑣𝜔𝛽₂𝜇_ℎ+𝛽₁𝛿₁Π_ℎ+𝜃𝛽₁𝛿₂Π_ℎ<0,𝜕𝑅₀1𝜕𝜔=−𝜔²𝛽₁𝜇_ℎΠ_ℎ𝜇²_𝑣Π_𝑣𝛿₁+𝜃𝛿₂<0, (29) which means that 𝑅₀ is a decreasing function of 𝜇_𝑣 and 𝜔. The question arises, which parameter is more crucial in order to decrease the reproductive number? By using the definition given in [8] and parameter values Π_ℎ=100, Π_𝑣=400, 𝛽₂=0.00034, 𝜇_𝑣=0.00054, 𝛽₁=0.4, 𝜇_ℎ=0.000274, 𝜔=0.00137, 𝛿₁=0.01, 𝜃=0.00304, and 𝛿₂=0.01, we see that the sensitivity index of the reproductive number with respect to 𝜇_𝑣 is −2 and with respect to 𝜔 is −1. It means that the most sensitive parameter for 𝑅₀ is 𝜇_𝑣. Increasing the mortality rate of Monochamus alternatus by 10% decreases 𝑅₀ by 20%. Also, increasing the exploitation rate of infected pines by 10%  decreases  𝑅₀ by 10%. Thus control strategies, for example, setting out beetle traps, setting vertical wood traps, and using chemicals to kill sawyer beetles, by cutting down dead pine trees and disposing of them before the emergence of beetles can be useful for eradicating the disease.

Figure 1 

The infected population for Πℎ=0.22, Π𝑣=0.32, 𝜇ℎ=0.00002, 𝜔=0.003, 𝜇𝑣=0.004, 𝛽1=0.0004, 𝛽2=0.00034, 𝛿1=0.0016, 𝛿2=0.00016, 𝜃=0.00301, 𝛼1=0.001, and 𝛼2=0.001.

The above mentioned measures are very effective to control pine wilt disease but 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 due to which most owners of forests hesitate to use these measures.

However, we can decrease the endemic level of the disease by increasing the parameters 𝜇_𝑣 and 𝜔. Figures 2, 3, 4, and 5 show different endemic levels of 𝐼_𝑣 and 𝐼_ℎ with respect to the parameters 𝜇_𝑣 and 𝜔. We see that by increasing these parameters the infective levels of pine trees and vectors decrease.

Figure 2 

The effect of 𝜇𝑣 on infected pine trees for Πℎ=0.22, Π𝑣=0.32, 𝜇ℎ=0.02, 𝜔=0.03, 𝜇𝑣=0.004,0.0058,0.09,𝛽1=0.0004, 𝛽2=0.00034, 𝛿1=0.0016, 𝛿2=0.00016, 𝜃=0.00301, 𝛼1=0.001, and 𝛼2=0.001.

Figure 3 

The effect of 𝜔 on infected pine trees for Πℎ=0.22, Π𝑣=0.32, 𝜇ℎ=0.002, 𝜔=0.03,0.04,0.05,0.06, 𝜇𝑣=0.004, 𝛽1=0.0004, 𝛽2=0.00034, 𝛿1=0.0016, 𝛿2=0.00016, 𝜃=0.00301, 𝛼1=0.001, and 𝛼2=0.001.

Figure 4 

The effect of 𝜇𝑣 on infected vectors for Πℎ=0.22, Π𝑣=0.32, 𝜇ℎ=0.02, 𝜔=0.03, 𝜇𝑣=0.004,0.0058,0.09,𝛽1=0.0004, 𝛽2=0.00034, 𝛿1=0.0016, 𝛿2=0.00016, 𝜃=0.00301, 𝛼1=0.001, and 𝛼2=0.001.

Figure 5 

The effect of 𝜔 on infected vectors for Πℎ=0.22, Π𝑣=0.32, 𝜇ℎ=0.002, 𝜔=0.003,0.03,0.07,0.15, 𝜇𝑣=0.004, 𝛽1=0.0004, 𝛽2=0.00034, 𝛿1=0.0016, 𝛿2=0.00016, 𝜃=0.00301, 𝛼1=0.001, and 𝛼2=0.001.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work has been fully supported by the University Research Fund of the National University of Sciences and Technology. The author is grateful to a referee for helpful comments which led to improvement of the paper.