Modelling the dynamics of Pine Wilt Disease with asymptomatic carriers and optimal control.

Pine wilt disease is a lethal tree disease caused by nematodes carried by pine sawyer beetles. Once affected, the trees are destroyed within a few months, resulting in significant environmental and economic losses. The role of asymptomatic carrier trees in the disease dynamics remains unclear. We developed a mathematical model to investigate the effect of asymptomatic carriers on the long-term outcome of the disease. We performed a stability and sensitivity analysis to identify key parameters and used optimal control to examine several intervention options. Our model shows that, with the application of suitable controls, the disease can be eliminated in the vector population and all tree populations except for asymptomatic carriers. Of the possible controls (tree injection, elimination of infected trees, insecticide spraying), we determined that elimination of infected trees is crucial. However, if the costs of insecticide spraying increase, it can be supplemented (although not replaced entirely) by tree injection, so long as some spraying is still undertaken.

Introduction

Among the vector-borne diseases of trees and plants, the most destructive are a red ring disease of palms and Pine Wilt Disease (PWD), whose causative agent is pine wood nematodes (PWNs)[1]. The vector for PWD is the pine sawyer beetle, which transfers nematodes to healthy host pine trees and usually kills host trees within a few months of infection. The lack of resin exudation of bark wounds become visible as a first symptom. The foliage become pale green in the second stage, yellow in the third stage and finally become reddish brown when the trees fail to resist against the disease. It is well-established that PWD has three different transmission paths: the first happens during maturation feeding[2]; the second during oviposition of the mature female on recently cut, dying, or dead pine trees through the oviposition wounds[3]; and the third is horizontal transmission, which happens during mating[4].

A number of epidemiological studies have been carried out to investigate the transmission dynamics of pine wilt disease[5-8]. These models investigating the spread and control of PWD are used to describe the host–vector interaction between nematode-carrying pine sawyers and pine trees. Lee[9] presented an epidemiological model of PWD and developed optimal-control strategies for the prevention of PWD. Khan et al.[10] introduced a dynamical model of PWD and investigates the stability of the disease with saturated incidence rate. They classified the total host tree size into three states: susceptible, exposed and infected host pine trees, while the vector size was also classified into three similar states. Ozair[11] included horizontal transmission and nonlinear incidence. The global stability of PWD in a model with nonlinear incidence rates was analyzed by Lee[5]. Optimal control has been used to study a variety of infectious disease[12-16], including plant diseases[17-19].

Asymptomatic carrier cases can play a critical role in the subsequent spread of PWD[20]. Asymptomatic infection increases the density of infected vectors, which further increases the level of infection in the host. Studies on asymptomatic infection in pine trees show that asymptomatic infected trees may remained infected for up to a year and may ultimately die[21]. Mathematical models that address pine tree dynamics with asymptomatic infections have previously been considered[22,23]. The effect of asymptomatic infection on neighboring trees has also been studied[20].

Here, we develop a dynamic model of PWD incorporating an asymptomatic carrier class and examine control policies that minimize implementation costs while protecting forests from the disease. To the best of our knowledge, none of the previous mathematical studies used optimal control to explore the transmission dynamics of the PWD in the presence of the asymptomatic carriers.

Model formulation

The total host (pine wood trees) and vector (beetles) are represented by N(t) and N(t), respectively. N(t) is further classified into four epidemiological classes: susceptible pine trees S(t), exposed pine trees E(t), asymptomatic carrier pine trees A(t) and infected pine trees I(t). N(t) is classified into three epidemiological classes: susceptible beetles S(t), exposed beetles E(t) and infected beetles I(t).

The recruitment rates of host trees and beetles are represented, respectively, by Λ and Λ, while the natural death rates of host pine trees and vector beetles are denoted by γ1 and γ2, and the disease mortality rate of host pine trees is represented by μ. Here, m and η are the respective rates of progression from the exposed class to the infected class in the host and vector populations. The term β1ψSI denotes the incidence rate, where β1 is the rate of transmission and ψ is the average number of daily contacts with vector adult beetles during maturation. β2 is the rate at which an infected beetle transmits a nematode through oviposition, with the average number of oviposition contacts per day denoted by θ. The termination of oleoresin exudation in susceptible trees without infection of nematode is denoted by α. We thus interpret β2θα as the transmission through oviposition, and hence β2θαSI represents the number of new infections. A fraction ω (0 ≤ ω ≤ 1) of the exposed tree class generates symptomatic infection, while the remaining fraction (1 − ω) generates asymptomatic infection. The vector incidence rate is given by the term KIS[15]. The schematic diagram for the PWD model is shown in Fig. 1.

Figure 1

Flow chart for the transmission of PWD. The short dashed arrows indicate the natural and the disease-specific death rates in each compartment. The long dashed arrows represent the interaction between the vector and pine trees. The long solid arrows represent the transition between compartments due to disease. The short solid arrows represent the recruitment.

The model is thus given bywith initial conditions

The total population sizes of host and vector are given by

For biological realism, we study the behaviour of the system (1) in the closed set

Nonnegative solutions of system (1) can be easily verified for appropriate initial values. The first four equations of (1) imply thatBy comparison theorem presented in[24], there exists t1 > 0, such that

Similarly, adding the last three equations of the system (1), we get

Using the comparison theorem again, there exists t2 > t1, such that

Hence, the solutions of the system (1) are bounded.

In the Supplementary Material, we determine and prove that the disease-free equilibrium (DFE) is globally asymptotically stable, which also rules out the possibility of a backward bifurcation. We also show that the endemic equilibrium is globally asymptotically stable, under certain conditions.

Sensitivity analysis of threshold dynamic

Due to uncertainties in experimental data, determining accurate outcomes from an epidemiological system is difficult[25]. To compensate for these uncertainties, we use partial rank correlation coefficients (PRCCs) to identify the impact of all parameters on . This technique measures the degree of the relationship between inputs and output of the system. Positive PRCCs indicate parameters that increase when they are increased, while negative PRCCs indicate parameters that decrease when they are increased. Parameters with PRCCs values greater than 0.4 in magnitude have a significant effect on the outcome.

Figure 2 illustrates the effect of parameter variations on for all fourteen parameters. Clearly, is most sensitive to γ1 and γ2, the natural death rates of pine trees and beetles, respectively; the latter can be controlled using insecticide (u3), while the former can be partially controlled by eliminating infected trees (u2). is also sensitive to the birth rates of pine trees and beetles, the latter of which can be controlled using insecticide (u3). The transmission rate K is also a sensitive parameter, which can be controlled by nematicide-injection and vaccination (u1).

Figure 2

Sensitivity of R0 to all input parameters.

Optimal control strategies

In this section, we introduce u1, u2 and u3 as three control measures that can affect PWD. The force of infection in the pine-tree population is reduced by (1 − u1), where precautionary measures efforts are denoted by u1; for example nematicide injection and vaccination. To keep the host tree population safe and to prevent infection, the nematicide-injection preventative control measure is used. We use the control variable u2 to describe elimination of infected host trees. Supplementary infections are extremely reduced by demolition and elimination of infected host trees. The removal of these infected trees guarantees that eggs, larvae and pupa that are occupying the host pines are devasted. Our third control variable represents spraying of insecticide and larvacide to kill adult insects and reduce the vector birth rate.

Model (1) is modified for optimal control as follows:with nonnegative initial conditions. The control functions u(t) = (u1, u2, u3) ∈ U associated to the variables S, E, A, I, S, E and I satisfy

The constants b0 and b1 are removal-rate constants whose inverses correspond to the average time spent in the relevant compartment. Since it is unlikely that infected trees will be removed within one day of infection, we set b1 = 1; hence the range 0 ≤ u2 ≤ 1 corresponds to a removal time between 1 day and infinite time. The objective functional for the optimal-control problem issubject to the control system (2). The constants L1, L2, L3, L4, B1, B2 and B3 are the weight or balancing constants, which measure the relative cost of interventions over the interval [0, T]. We seek optimal controls , such that

Clearly, the equations in the control system (2) are bounded above, and thus we can apply the results in[26] to model (2). Moreover, the set of control variables and the state variables is nonempty, and the set of control variables denoted by U is closed and convex. In the control problem (2), the right-hand side is continuous and bounded above by state variables and a sum of the bounded control, and can be expressed as a linear function of U having state- and time-dependent coefficients. Hence there exists constants m > 1 and l1, l2 > 0 such that the integrand L(y, u, t) of the objective functional J is convex and satisfies

We apply the results presented in[27] to justify the existence of (2) and to obey the above conditions. Clearly, the set of control and state variables are bounded and nonempty. The solutions are bounded and convex. Therefore the system is bilinear in the control variables. We verify the last condition:where L1, L2, L3, L4, B1, B2, B3, l1, l2 > 0 and m > 1. We have thus proved the following theorem.

Theorem 1. For the objective functional (4) and the control set (3) subject to the control system (2), there exists an optimal control such that .

In order to get the solution of the control problem, it is necessary to obtain the Lagrangian and the Hamiltonian of (2). The Lagrangian L is expressed as

By choosing X = (S, E, I, S, E, I), U = (u1, u2, u3) and λ = (λ1, λ2, λ3, λ4, λ5, λ6, λ7), the Hamiltonian can be written

We use Pontryagin’s Maximum Principle[28] to obtain the optimal solution of the control system (2). Since and are solutions to the control problem (2), there exist adjoint variables λ (i = 1, 2, 3, 4, 5, 6, 7) satisfying the following conditions:

Theorem 2. For the optimal-control measures and the state solutions of system (2), there exist adjoint variables λ (i = 1, 2, 3, 4, 5, 6, 7) such thatwith the transversally conditions

Furthermore, the controls are given by

Proof. To determine the required adjoint system (8) and the transversality conditions mentioned in (9), we utilize the Hamiltonian in (6). By applying the third condition of (7), we get (8). Applying the second condition of (7), we get (9).

Numerical Results

Unless mentioned otherwise, we use the fourth-order Runge–Kutta method over a timescale of 100 days. The input parameters for our simulations are L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B1 = 0.10, B2 = B3 = 10, c = 0.001241, b0 = 0.21; all other parameter values are shown in Table 1.

Table 1

Parameter interpretations and their sample values used in numerical simulations.

Parameter	Interpretation	Values	References
β1	transmission probability during maturation	0.0016600 day−1	[29]
β2	transmission probability of nematode through oviposition	0.0004000 day−1	[29]
ψ	contacts averagely made during maturation per day	0.2000000 day−1	[30]
γ1	natural death rate of host pine trees	0.0000301 day−1	[31]
γ2	natural death rate of vector beetles	0.0011764 day−1	[32]
θ	contacts averagely made during oviposition per day	0.0023000 day−1	assumed
m	progression rate of pine trees from EH to IH	0.0133000 day−1	assumed
ΛV	recruitment rate of susceptible vector	0.0132652 day−1	assumed
K	the rate at which the adult beetles carry PWN when they escape from dead trees	0.00305 day−1	[33]
η	progression rate of vectors from EV to IV	0.0100000 day−1	assumed
α	probability that host susceptible cease oleoresin exudation without infected by the nematode	0.0032000 day−1	assumed
ΛH	recruitment rate of host trees	0.0020210 day−1	assumed
ω	rate of symptomatic cases	0.1000000	assumed
μ	transfer rate from EV to IV	0.0022000 day−1	assumed

Elimination of infected trees (u2) and spraying of insecticides (u3)

We considered two controls: the elimination of infected trees (u2) and the spraying of insecticides (u3) in the absence of tree injection and vaccination. Figures 3 and 4 show the outcomes in both the absence and presence of control. Figure 3 shows the dynamics of the pine-tree population, while Fig. 4 shows the dynamics of the vector population. With these controls, we see a rapid increase in the population of susceptible trees (Fig. 3(a)) and eventual elimination of exposed and infected trees (Fig. 3(b,d)), with only the asymptomatic carriers remaining in the infected classes (Fig. 3(c)). The vector population is eventually depleted (Fig. 4(a–c)) in the presence of these two controls. The two control profiles u2 and u3 are bounded up to 0.4 and 0.8 (Fig. 4(d)). Biologically, u2 is the additional elimination rate of only infected trees, while u3 acts to simultaneously increase the removal rate of all vectors, while also decreasing the birth rate. Since all interventions range between 0 (no control) and 1 (complete control), this suggests that our objective can be achieved with only partial controls. Hence if infected trees are removed 2.5 days or later after infection or if insecticides/larvacides are up to 80% effective, the infection can be controlled.

Figure 3

The behaviour of the pine-tree population for the controls u2 and u3; (a) Susceptible pine trees, (b) Exposed pine trees, (c) Asymptomatic pine trees, (d) Infected pine trees.

Figure 4

The behaviour of the vector (beetles) population for the controls u2 and u3; (a) Susceptible beetles, (b) Exposed beetles, (c) Infected beetles, (d) Control profile.

Tree injection (u1) and spraying of insecticides (u3)

We next examine the combination of tree injection (u1) and insecticide spraying (u3). The results are shown in Figs. 5 and 6. With these two controls, there is a significant increase in the population of susceptible and exposed pine trees, while the population of asymptomatic carriers and infected pine trees are reduced but not eliminated (Fig. 5). This suggests that the elimination of infected pine trees has a significant impact on the disease. Note that the vector population is eliminated using these controls (Fig. 6).

Figure 5

The behaviour of the pine-tree population for the controls u1 and u3; (a) Susceptible pine trees, (b) Exposed pine trees, (c) Asymptomatic pine trees, (d) Infected pine trees.

Figure 6

The behaviour of the vector (beetles) population for the controls u1 and u3; (a) Susceptible beetles, (b) Exposed beetles, (c) Infected beetles, (d) Control profile.

Tree injection (u1) and elimination of infected trees (u2)

Considering u1 and u2 in combination, Figs. 7 and 8 illustrate that, without insecticide spraying, the control (minimization and/or elimination) of infection in the pine trees is not possible. While the population of susceptible pine trees has a slower decline with these control (Fig. 7(a)), the infection eventually takes over. Likewise, although the susceptible beetle population is recovered using these controls, the infection nevertheless eventually dominates (Fig. 8). It follows that, without insecticide spraying, the control of infection is not possible.

Figure 7

The behaviour of the pine-tree population for the controls u1 and u2; (a) Susceptible pine trees, (b) Exposed pine trees, (c) Asymptomatic pine trees, (d) Infected pine trees.

Figure 8

The behaviour of the vector (beetles) population for the controls u1 and u2; (a) Susceptible beetles, (b) Exposed beetles, (c) Infected beetles, (d) Control profile.

Complete control

We now apply all three controls in order to determine the ideal outcome (Figs. 9 and 10). Comparing Fig. 9 to Fig. 3, we see that susceptible pine trees recover faster and the disease is eliminated quicker, except for asymptomatic carriers. We thus see that the most effective strategy is to apply all three controls, although similar results can be achieved by applying only two controls: elimination of the infected pine trees (u2) and the spraying of insecticides (u3).

Figure 9

The behaviour of the pine-tree population for the controls u1, u2 and u3; (a) Susceptible pine trees, (b) Exposed pine trees, (c) Asymptomatic pine trees, (d) Infected pine trees.

Figure 10

The behaviour of the vector (beetles) population for the controls u1, u2 and u3; (a) Susceptible beetles, (b) Exposed beetles, (c) Infected beetles, (d) Control profile.

Temporal variation of control profiles

Next, we investigate the control profiles and their relationships to the weight constants. In Fig. 11, we fix the weight constants L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B2 = B3 = 10 and allow B1 to vary. In Fig. 12, we fix the weight constants L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B1 = 0.1, B3 = 10 and allow B2 to vary. In Fig. 13, we fix L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B1 = 0.1, B2 = 10 and allow B3 to vary. These variations represent fluctuating costs of implementing our controls.

Figure 11

Temporal variation of the control profile for L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B2 = B3 = 10; (a) B1 = 0.10, (b) B1 = 1, (c) B1 = 10, (d) B1 = 100.

Figure 12

Temporal variation of the control profile for L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B1 = 0.1, B3 = 10; (a) B2 = 0.010, (b) B2 = 0.10, (c) B2 = 1, (d) B2 = 10.

Figure 13

Temporal variation of the control profile for L1 = 0.01, L2 = 0.002, L3 = 0.0020, L4 = 0.003, B1 = 0.1, B2 = 10; (a) B3 = 0.010, (b) B3 = 0.1, (c) B3 = 1, (d) B3 = 10.

From Fig. 11, we see that, as the cost of u1 increases, the control profile is dominated by u3. That is, if tree injection becomes prohibitively expensive, the procedure can be replaced by increased insecticide spraying.

Figure 12 shows little variation in the control profiles as the cost B2 increases unless the cost is prohibitive. This suggests that the control u2 is worth implementing, even at high cost. The combination of u1 and u3 alone does not eliminate infection, so it follows that elimination of infected trees is essential to disease control. This may hinder disease eradication if the costs of elimination become prohibitively expensive.

Figure 13 shows that if the cost of insecticide spraying increases, the control profile is dominated by tree injection. Interestingly, while the combination of u2 and u3 produced superior results to the combination of u1 and u2, the latter combination can still produce effective results if supplemented by a small amount of insecticide spraying.

Discussion

We developed a mathematical model to examine the effect of asymptomatic carriers of Pine Wilt Disease (PWD) on the long-term course of disease. We showed that the disease-free equilibrium was globally asymptotically stable and that the endemic equilibrium was globally asymptotically stable under some conditions. A sensitivity analysis identified key parameters: natural death rates in trees and beetles; birth rates in both trees and beetles; and transmission rates from trees to beetles.

We applied several controls to our system: tree injection, insecticide spraying and elimination of infected trees. These were chosen in conjunction with the most sensitive parameters except for the natural birth and death rates of trees, since our ultimate goal is the preservation of trees. We showed that the disease can be eliminated using suitable controls, except for the asymptomatic carriers. By including this class, our model showed that the disease may remain endemic, requiring permanent control, even in the best-case scenario.

Examining the controls in detail, we found that elimination of infected trees is critical, especially when used in conjunction with insecticide spraying. If the cost of insecticide spraying becomes prohibitive, it can be partially replaced by tree injection. However, if the costs of elimination of infected trees becomes prohibitive, there is no simple replacement, which may result in runaway costs.

It follows that we can control the disease using suitable interventions, but it will not be eliminated due to the presence asymptomatic carriers. The presence of infection in these carriers suggests that infection can restart in nearby healthy trees. It follows that our control measures must be undertaken continually unless such asymptomatic carriers can be identified and removed. This has long-term implications for disease management and economic investment.

Supplementary Information.