# Stability and Sensitivity Analysis of a Plant Disease Model with Continuous Cultural Control Strategy

Authors: Zhang Zhonghua; Suo Yaohong
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

In this paper, a plant disease model with continuous cultural control strategy and time delay is formulated. Then, how the time delay affects the overall disease progression and, mathematically, how the delay affects the dynamics of the model are investigated. By analyzing the transendental characteristic equation, stability conditions related to the time delay are derived for the disease-free equilibrium. Specially, whenR0=1, the Jacobi matrix of the model at the disease-free equilibrium always has a simple zero eigenvalue for all τ≥0. The center manifold reduction and the normal form theory are used to discuss the stability and the steady-state bifurcations of the model near the nonhyperbolic disease-free equilibrium. Then, the sensitivity analysis of the threshold parameter R0 and the positive equilibrium E* is carried out in order to determine the relative importance of different factors responsible for disease transmission. Finally, numerical simulations are employed to support the qualitative results.


---

## Body

## 1. Introduction

Viral disease is a key constraint on the production of staple food crop in the lesser developed countries. Diseases caused by plant viruses in cassava (Manihotesculenta), sweet potato (Ipomoeabatatas) and plantain (musa spp.) are among the main constraints on sustainable production of these vegetatively propagated staple food crops, see Rybicki and Pietersen [1], Dahal et al. [2], Gibson and Aritua [3], and Thresh and Cooter [4]. Furthermore, new viral strains frequently emerge, some of which bring devastating consequences, such as the current pandemic of the virus causing cassava mosaic disease (CMD) in Africa, see Gibson et al. [5]. Governments and farmers have been evolving practices for combatting with the plagues suffered by crops. The growing understanding of the interactions between pathogen and host has enabled us to develop a wide array of measures for the control of specific plant diseases. Large investments are underway to alleviate poverty and malnutrition by developing new or more effective control strategies, which include plant breeding for resistance to the virus, control of vectors as well as crop sanitation through removal of diseased plants, and improved selection of planting material for these vegetatively propagated crops. Such experiences have led to the development of integrated management concepts for virus diseases that combine available host resistance, cultural, chemical, and biological control measures. Examples of how epidemiological information can be used to develop effective integrated disease management (IDM) strategies for diverse situations have been described in [6–9]. IDM involves the selection and application of a wide range of control strategies that minimize losses and maximize returns. A cultural control strategy including replanting or roguing diseased plants is a widely accepted treatment for plant epidemics which involves periodic inspections followed by removal of the infected plants.Mathematical models of plant-virus epidemics were developed to provide detailed explanation on how to describe, analyze, and predict epidemics of plant disease for the ultimate purposes of developing and testing control strategies and tactics for crop protection, see van den Bosch et al. [10], Chan and Jeger [11], Fishman et al. [12], and the references cited therein. A simple model for plant disease with a continuous cultural control strategy, such as replanting and roguing or removing, is as follows:
(1)dS(t)dt=σϕ-βS(t)I(t)-ηS(t),dS(t)dt=σ(1-ϕ)+βS(t)I(t)-(η+ω)I(t),
where S and I denote the respective number of the susceptible and infected plants and β is the transmission rate. Crop is planted fromin vitro propagated virus free material, from cuttings taken from a previous crop, or from a combination of these methods. New plants enter the system by continual replanting at rate σ with a proportion ϕ for the susceptible plants and 1-ϕ for the infected ones. Removal occurs for sanitation at rate ω or for death at rate η, and 1/η denotes either harvest time or the end of reproductive lifetime of plants. Neglecting the biological meaning of the parameters, (1) can be seen as an SI human disease model with immigration of the infective individuals and its extended forms have been extensively studied, for example, van den Bosch et al. [10], Brauer and van den Driessche [13], Tang et al. [14], and the references cited therein.In [12], Fishman et al. divided the time axis (0,∞) into equal periods of length T and developed a model with periodic control strategy as follows:
(2)dxi(t)dt=Axi(t)(1-xi(t)),0≤t≤T,
where A is the infection rate and xi(t) is the ratio between the number ni(t) of the infected trees at instant t of period i and the total number Ni of trees at the beginning of the period. Equation (2) characterizes the temporal spread of an epidemic in a closed plant population with periodic removal of infected plants and has an application to the spread ofcitrus tristeza virus disease. Economic evaluations and comparisons between two policies, eradication or no eradication, were given, and further simulations of the model and sensitivity analysis for a wide range of parameters were presented. They concluded that the discovery-eradication program is economically justified and superior to allowing the disease to progress unchecked. The results were helpful in evaluating policies of controlling the disease, and (2) could be modified to simulate other plant epidemics with periodic treatments.Based on [12], Tang et al. in [14] proposed periodic pulse replanting and roguing strategies and changed (1) into the following:
(3)dS(t)dt=-βS(t)I(t)-ηS(t),t≠nT,dI(t)dt=βS(t)I(t)-ηI(t),t≠nT,S(nT+)=S(nT)+σϕ,t=nT,I(nT+)=(1-ω)I(nT)+σ(1-ϕ)ϕ,t=nT,
where n=1,2,…, and T is a fixed positive constant and denotes the period of the impulsive effect. The sufficient conditions under which the infected plant free periodic solution with fixed moments is globally stable were obtained.It can be seen that the reversion of the infected plants is not considered in the above models. In [10], to study the effect of the reversion of infected plants on the transmission of disease, van den Bosch et al. constructed the following model:
(4)dS(t)dt=σϕ+σ(1-ϕ)r(1-p)I(t)+S(t)(1-p)I(t)+S(t)-βS(t)I(t)-ηS(t),dI(t)dt=σ(1-ϕ)(1-r)(1-p)I(t)(1-p)I(t)+S(t)+βS(t)I(t)-(η+ω)I(t).
Some cuttings from infected plants may be healthy due to reversion with probability r. Cuttings are selected visually or using diagnostic methods and discarded with probability p. They found a threshold parameter R0 and proved that (4) owns a unique positive equilibrium when R0>1, and its stability was numerically checked. They have shown that the development of new and improved disease control methods for viral diseases of vegetatively propagated staple food crops ought to take evolutionary responses of the virus into consideration. Not doing so leads to a risk of failure, which can result in considerable economic losses and increased poverty. Specificallyin vitro propagation, diagnostics, and breeding methods carry a risk of failure due to the selection for virus strains that build up a high within-plant virus titre. For vegetatively propagated crops, sanitation by roguing has a low risk of failure owing to its combination of selecting for low virus titre strains as well as increasing healthy crop density.Latent infection of plants by pathogens has been recognized for many years and it is often considered one of the highest levels of parasitism, since the host and the parasite coexist with minimal damage to the host, see Sinclair [15]. Latent infection is important in the epidemiology, the control of the plant diseases, and also in breeding for resistance or tolerance to a pathogen; see Chan and Jeger [11]. An understanding of latent infection contributes to development of effective control measures, as does an understanding of penetration, colonization, disease expression, and yield losses. Many plant diseases, such as chlorotic leaf distortion of sweet potato, citrus black spot, Colletotrichumgloeosporioides, Alternariaalternata, possess latent period, that is, the time elapsed between exposure to a pathogenic organism and when symptoms and signs are first apparent; see Ames et al. [16], Miles et al. [17], Chakraborty [18], and Karaoglanidis et al. [19]. The earliest plant disease model is due to van der Plank [20], and in which he used a delay differential equation to represent I(t), the density of host tissue first infected at or before time t,
(5)dI(t)dt=R(I(t-p)-I(t-p-i))(1-I(t)).
The latent period p and the infectious period i are constant, and the parameter R is the corrected basic rate of infection. The key threshold is the progeny-parent ratio, iR:avalue>1 indicates that the density of infected tissue will increase. Unfortunately delay differential equations are difficult to analyze; see Madden [21] and Murray [22] and, despite widespread adoption of discrete time approximations to van der Plank’s model in early simulations of plant disease, see Teng [23], Zadoks [24] and a number of often subtle mathematical analyses that followed by Jeger [25, 26], Kushalappa and Ludwig [27], Waggoner [28], the model is now rarely used in theoretical studies. However, it is so influential and still of significant historical interest; see Cunniffe et al. [29].To our knowledge, how the latent period affects the dynamics of (4) remains unknown. Motivated by this, using a time delay τ to denote the latent period of plant disease and replacing the incidence rate of the disease by βS(t-τ)I(t-τ) we formulate the following model:
(6)dS(t)dt=σϕ+σ(1-ϕ)r(1-p)I(t)+S(t)(1-p)I(t)+S(t)-βS(t-τ)I(t-τ)-ηS(t),dI(t)dt=σ(1-ϕ)(1-r)(1-p)I(t)(1-p)I(t)+S(t)+βS(t-τ)I(t-τ)-(η+ω)I(t).The remainder of the paper is arranged as follows. Section2 discusses the stability of the disease-free equilibrium under the conditions of R0≠1, or R0=1, τ=0; Section 3 computes the normal forms on the center manifold and investigates the dynamical behaviors of system (6) near the disease-free equilibrium when R0=1, τ>0; Section 4 is the sensitivity analysis of R0 and the positive equilibrium to all of the parameters in system (6); Section 5 performs numerical simulations to illustrate the qualitative results; Section 6 makes some conclusions.


## 2. Stability of the Disease-Free Equilibrium

Let(S,I) be an equilibrium of system (6). Then, (S,I) must solve the following algebraic equations:
(7)σϕ+σ(1-ϕ)r(1-p)I+S(1-p)I+S-βSI-ηS=0,σ(1-ϕ)(1-r)(1-p)I(1-p)I+S+βSI-(η+ω)I=0.For convenience, we introduce the following threshold parameter:(8)R0=(1-ϕ)(1-r)(1-p)η2+βση(η+ω).Theorem 1.
IfR0≤1, system (6) only has a disease-free equilibrium E0 and if R0>1, system (6) has a unique positive equilibrium E* except for the disease-free equilibrium E0, where E0=(σ/η,0), E*=(S*,I*), S*=(-a1-a12-4a0a2)/2a2, I*=(σ(1-p)(1-r)(1-ϕ)-S*(η+ω-βS*))/(η+ω-βS*)(1-p), a0=σ(p-1)(ϕ(r-1)-r)(η+ω), a1=η2p-βσ+ηωp+βσp+ω2+ηω, and a2=-β(ηp+ω).By the transformationx1(t)=S(t)-(σ/η), x2(t)=I(t), system (6) can be rewritten as the following equivalent form:
(9)dx(t)dt=Mx(t)+Nx(t-τ)+f(x(t)),
where(10)M=(-η-(1-ϕ)(p-1)(r-1)η0(1-ϕ)(1-r)(1-p)η-(η+ω)),N=(0-βση0βση),x(t)=(x1(t)x2(t)),f(x1(t),x2(t))=(-βx1(t-τ)x2(t-τ)+η2(r-1)(p-1)(1-ϕ)x2(t)((p-1)x2(t)-x1(t))η(p-1)x2(t)-ηx1(t)-σβx1(t-τ)x2(t-τ)-η2(r-1)(p-1)(1-ϕ)x2(t)((p-1)x2(t)-x1(t))η(p-1)x2(t)-ηx1(t)-σ).It can be calculated that the characteristic equation of system (9) at the origin owns the following form:
(11)(λ+η)(λ-η(1-ϕ)(1-r)(1-p)+η+ω-βσηe-λτ)=0.
Then, neglecting the time delay τ, we have the following.Theorem 2.
IfR0>1, the disease-free equilibrium E0 is unstable, while it is asymptotically stable if R0<1.Clearly, the origin is a nonhyperbolic equilibrium of system (9) if R0=1 and τ=0. To discuss its stability, we use the transformation as x1(t)=u(t)-((η+ω)/η)v(t), x2(t)=v(t). Then, system (9) is changed into
(12)du(t)dt=-ηu(t)+f1(u(t),v(t)),dv(t)dt=f2(u(t),v(t)),
where
(13)f1(u(t),v(t))=-ωv(t)ησ(σ-(ηp+ω)v(t)+ηu(t))×((t)2ησ(2η(r+(1-r)(ϕ+p(1-ϕ)))+ω-η)u(t)hhh+σ(η2(1+(1-r)(1-ϕ)(p2+1))+ω2+2ωηgghhhhhh×((1-r)(p(1-ϕ)+ϕ)+r)η2)v(t)ghh-η2(η((1-r)(p(1-ϕ)+ϕ)+r)+ω)u(t)2ghh+η(η(1+p)+2ω)ghh×(η((r+p(1-r))(1-ϕ)+ϕ)+ω)u(t)v(t)ghh-(η+ω)(ω+ηp)ghh×(η((r+p(1-r))(1-ϕ)+ϕ)+ω)v(t)2),f2(u(t),v(t))=-v(t)σ((ηp+ω)v(t)-ηu(t)-σ)×((t)2ση(2η(ϕ(1-p)(1-r)+p(1-ϕ)+r)ghhhhh-η+ω)u(t)g+σ(ω2+η2((p2(1-r)+r)(1-ϕ)+ϕ)ghhghhh+2ηω((1-ϕ)·(p(1-r)+r)+ϕ))v(t)g-η2(η((p(1-ϕ)+ϕ)(1-r)+r)+ω)u(t)2g+η(η(1+p)+2ω)g×(η((1-ϕ)(p(1-r)+r)+ϕ)+ω)u(t)v(t)g-(η+ω)(ω+ηp)g×(η((1-ϕ)(p(1-r)+r)+ϕ)+ω)v(t)2).By the existence theorem in the center manifold theory, see Wiggins [30] for details, there exists a center manifold for system (12), which can be locally expressed as follows:
(14)Wc(0)={(u,v)∈ℝ2∣u=h(v),∥v∥<δ,h(0)=0,Dh(0)=0∣δ>0ℝ2},
where δ is sufficiently small, Dh is the derivative of h with respect to v.Now, the first task is to compute the center manifoldWc(0). For the purpose, we assume h(v) has the form
(15)u=h(v)=h2v2+h3v3+h4v4+h5v5+⋯,
where hi, i=2,…,5, are constants to be determined in the following. By the invariance of Wc(0) under the dynamics of (12), h(·) satisfies
(16)Dh·f2(h,v)+ηh-f1(h,v)=0.
Substituting (11) into (12), and then equating coefficients on each power of v to zero, yields
(17)h2=-ω2ση2(2ωη((1-r)(ϕ(1-p)+p)+r)ghh+η2((1-ϕ)(p2(1-r)+r)+ϕ)+ω2),h3=-ω4σ2η3(2ϕ2(r-1)2(p-1)2(p+1)2+2ϕ(1-r)ghhhhh×(1-p)(2r(1-p)(1+p)2+p2·(2p+3))ghhhhh-2(p2r(1-p)(1+p)2ghhhhhhhhh+p(p2+p-1)(rp-p-r))η4ghhhhh+(10ϕ2ω(1-r)2·(1+p)(1-p)2+ϕ(1-r)ghghhhhhh×(1-p)(20p2+13p-1-20r(1-p2))ghghhhhhh+ω(10p3(r-1)2+(1-r)(10r+3)p2ghffhhfffhhhhhh+2(1-r)(5r-2)p+r(10r-1)))η3ghfffffhh+(12ϕ2ω2(1-r)2·(1-p)2+ϕω2ghffffffghghh×(1-r)(1-p)(29p-24rp+24r-5)ghffffffghghh+ω2((r-1)(12r-17)p2)ghffffffghghh+24r(1-r)p+(3r+2)(4r-1)ϕ2)η2ghhhhh+(12ω2(1-r)(1-p)ϕghfffffhhhfh-ω3(12(1-r)·(1-p)-11)ω2)η+3ω4ϕ2(r-1)2),h4=-12σ3η((2h2p2σ((1-r)(1-ϕ)(1-p)-pr)ghhhhhh-(1-r)(1-ϕ)(1-p)p3ω)η3ghhhh+((1-p)6h2σpω(1-r)(1-p)(1-ϕ)gghhhhhg+2h22σ2(1-2(1-r)(1-p)(1-ϕ))ghhhhhhg-3h3σ2·((1-ϕ)(p2(1-r)+r)+r)ghhhhhhg-3p2ω2(1-r)(1-ϕ)(1-p))η2ghhhh-(h3σ2ω(8(1-r)·(1-p)(1-ϕ)-7)ghhhhhhg-4h2σω2(1-r)(1-p)(1-ϕ)ghhhhhhg+2h22σ2ω2-3pω3(1-r)·(1-p)(1-ϕ))ηghhhh-4h3σ2ω2-(1-r)(1-r)(1-ϕ)(2h2p2σ((1-r)(1-ϕ)(1-p)-pr)).On substituting (18) into the second equation of system (12), we obtain the following equation on the center manifold Wc(0):
(18)dv(t)dt=d2v(t)2+d3v(t)3+⋯,
where
(19)d2=-2σ×(((p2-1)(1-ϕ)(1-r)-1)η2ghh+2ω(1-(1-p)(1-r)(1-ϕ))η+ω2(p2-1)),d3=-3ησ2(2p2(r-1)(ϕ-1)(p-1)η4ghh+ω(2ϕ2(r-1)2(p-1)2(p+1)+ϕ(1-p)ghhhhhh·(4r2(p-1)(p+1)-r(8p2+5p-5)ghhhhhhhhh+4p2+5p-1)ghhhhhh+p(2p2+3p-1)-r·(4p2+p-6)ghhhhhh+2r2(1-p)2(1+p))η3ghh+ω2(ϕ2((pr+4)4(1-r)2+p(p-2)ghhhhhhhhhh+p2r(r-2)+pr(pr+4))ghhhhhhh+ϕ(1-p)ghhhhhhh×(9p+1-r(17p-7)-8r2(1-p))ghhhhhhh+5p2-2+(1-p)(1+9p)ghhhhhhh+4r2(1-p)2)η2ghh+ω3(3-4(1-p)(1-r)(1-ϕ))η+ω4).Therefore, by using the center manifold theorem [30], we have the following results.Theorem 3.
IfR0=1, the disease-free equilibrium E0 is unstable since d2<0.It can be seen that (11) always has a negative root -η. To investigate the stability of the disease-free equilibrium when τ>0, we reduce (11) to the following
(20)λ+p1+q1e-λτ=0,
where p1=η+ω-η(1-ϕ)(1-r)(1-p), q1=-(βσ/η).Suppose thatλ=iy, y>0 is a root of (20), then y satisfies
(21)p1+q1cos⁡(yτ)=0,y-q1sin(yτ)=0.
Eliminating the trigonometric functions yields
(22)y2q12=q12-p12q12.It can be computed that(23)q12-p12=-η+ωq12((1-R0)2+2βση(η+ω)(1-R0)).
Then (22) leads to a contradiction if R0<1, and we have the following results.Theorem 4.
IfR0<1, the disease-free equilibrium E0 is asymptotically stable for any τ>0. If R0>1+(βσ/η(η+ω)), E0 is unstable for any τ>0.It can be seen thatq12-p12>0 if 1<R0<1+(βσ/η(η+ω)), which means that (22) has a unique positive root y0=q12-p12. That is, there is a single pair of purely imaginary roots ±iy0 of (20). Let τk=(1/y0)arccos⁡(-p1/q1)+2kπ/y0, k=1,2,3,….Theorem 5.
If1<R0<1+(βσ/η(η+ω)), the disease-free equilibrium E0 is unstable for any 0<τ<τ0.Theorem 6.
Ifτ=τk and 1<R0<1+(βσ/η(η+ω)), system (6) undergoes Hopf bifurcations at the disease-free equilibrium E0.Proof.
Differentiating (20) with respect to τ yields
(24)(1-q1τe-λτ)dλdτ=λq1e-λτ,
which gives
(25)(dλdτ)-1=-τ-λ-1λ(λ+p1).
When λ=±iy0 and τ=τk, we obtain
(26)(dℜ𝔢{λ}dτ)-1|λ=±iy0,τ=τk=ℜ𝔢(dλdτ)-1|λ=±iy0,τ=τk=y02y04+p12y02>0.
The transversality condition for Hopf bifurcations holds. Then system (6) undergoes Hopf bifurcations at the disease-free equilibrium E0 when τ=τk.


## 3. Normal Forms on the Center Manifold for a Simple Zero Eigenvalue

It can be obtained that zero is a simple eigenvalue ofR0=1 and τ≠1/q1. In this section, we refer the reader to Hale and Lunel [31] for notation and general results on the theory of retarded functional differential equations (RFDES). To determine the dynamic properties of the disease-free equilibrium E0 with R0=1, τ>0, we have to compute the normal forms on the center manifold. The method we use is based on the center manifold reduction and normal form theory; see Faria and Magalhaes [32, 33]. In the following, we shall compute the normal form of model (6) associated with the zero eigenvalue.For convenience, we rescale the time byt→t/τ to normalize the delay of system (9) and obtain
(27)du1(t)dt=τ[×(η((p-1)u2(t)-u1(t))-σ)-1-ηu1(t)-η(1-ϕ)(p-1)(r-1)u2(t)-βσηu2(t-1)-βu1(t-1)u2(t-1)+(η2(r-1)(p-1)(1-ϕ)ghh×(u2(t)(p-1)-u1(t))u2(t)η2)×(η((p-1)u2(t)-u1(t))-σ)-1],du2(t)dt=τ[(η((p-1)u2(t)-u1(t))-σ)-1(η(1-ϕ)(1-r)(1-p)-(η+ω))u2(t)+βσηu2(t-1)+βu1(t-1)u2(t-1)-(η2(1-ϕ)(1-r)(1-p)ghh×(u2(t)(p-1)-u1(t))u2(t)η2)×(η((p-1)u2(t)-u1(t))-σ)-1].Let𝒞=𝒞([-1,0],ℝ2) be the Banach space of continuous functions from [-1,0] into ℝ2 with supremum norm. We define zt∈𝒞, as zt(θ)=z(t+θ), θ∈[-1,0]. Using the perturbation: of β=β0+ϵ, where β0=(η/σ)(η+w-(1-ϕ)(1-r)(1-p)η), system (27) can be written as the functional differential equation
(28)z˙(t)=L(ϵ)(zt)+F(zt,ϵ),
where ϵ∈V is considered as a parameter, V is a neighborhood of zero in space of real numbers, L:𝒞×V→ℝ2 is a parameterized family of bounded linear operators, and F:𝒞×V→ℝ2 is a function with F(0,ϵ)=0, ∂F(0,ϵ)/∂z=0 for all ϵ∈ℝ, and they have the following respective forms:
(29)L(ϵ)(φ)=L0φ+L1(ϵ)φ,F(φ,ϵ)=F2(φ,ϵ)+F3(φ,ϵ)+h.o.t,
where h.o.t stands for the higher-order terms, φ=(φ1,φ2)T, and(30)F2(φ,ϵ)=τ(η2σ(r-1)(p-1)(1-ϕ)((p-1)φ22(0)-φ1(0)φ2(0))-β0φ1(-1)φ2(-1)η2σ(1-r)(1-p)(ϕ-1)((p-1)φ22(0)-φ1(0)φ2(0))+β0φ1(-1)φ2(-1)),F3(φ,ϵ)=τ(-ϵ0φ1(-1)φ2(-1)+η3σ2(r-1)(p-1)(1-ϕ)φ2(0)·(φ1(0)-(p-1)φ2(0))((p-1)φ2(0)-φ1(0))ϵ0φ1(-1)φ2(-1)+η3σ2(r-1)(p-1)(ϕ-1)φ2(0)·(φ1(0)-(p-1)φ2(0))((p-1)φ2(0)-φ1(0))),L0φ=τ(-η(1-ϕ)(1-p)(r-1)φ2(0)-β0σηφ2(-1)-ηφ1(0)(η(-ϕ)(1-r)(1-p)-η-ω)φ2(0)+β0σηφ2(-1)),L1(ϵ)φ=τ(-σηϵφ2(-1)σηϵφ2(-1)).From the Riesz representation theorem the linear mapL can be expressed in integral form as follows:
(31)L(ϵ)(φ)=∫-10‍dηϵ(θ)φ(θ),
where ηϵ(·) is a bounded variation function on [-1,0] and can be defined as
(32)ηϵ(θ)=τMδ(θ)+τN(ϵ)δ(θ+1),δ(·) is the Dirac delta function, and N(ϵ)=(0-((β0+ϵ)σ/η)0(β0+ϵ)σ/η).Letℝ2* be the 2-dimensional vector space of row vectors and denote 𝒞*=𝒞([-1,0],ℝ2*). We define the adjoint bilinear form on 𝒞*×𝒞 as follows:
(33)〈ψ(s),φ(θ)〉=ψ(0)φ(0)-∫-10‍∫0θ‍ψ(ξ-θ)dηϵ(θ)φ(ξ)dξ,
where ψ=(ψ1,ψ2)∈𝒞* and φ=(φ1,φ2)T∈𝒞.LetA(ϵ) be the infinitesimal generator of the flow for the linear system z˙(t)=L(ϵ)zt, with spectrum σ[A(ϵ)]. The adjoint operator A*(ϵ) is defined as the infinitesimal generator for the solution operator of the adjoint equation in 𝒞*(34)w˙(t)=-∫-10‍w(t-θ)dηϵ(θ).It is well known that the eigenvalues ofA(ϵ) with zero real parts play an important role in the bifurcation theory of RFDES. Let A0=A(0) and denote all of its singular eigenvalues by the set Λ0; that is,
(35)Λ0={λ∈σ(A0)∣ℜ𝔢λ=0}.
According to the discussion in Section 2, we have Λ0={0}.Using the formal adjoint theory for FDEs in Hale and Lunel [31], the phase space 𝒞 can be decomposed by Λ0 as 𝒞=P⊕Q, where P is the generalized eigenspace associated with the eigenvalues in Λ0, Q={φ∈𝒞|〈ψ,φ〉=0forallψ∈P*}, and the dual space P* is the generalized eigenspace for A*(0) associated with the eigenvalues in Λ0. In particular, we consider bases for P and P* denoted by Φ and Ψ, respectively, and satisfying 〈Ψ(s),Φ(θ)〉=1. Then, we can choose Φ and Ψ as follows:
(36)Φ(θ)=(η+ω,-η)T,-1≤θ≤0,Ψ(s)=(0,k),0≤s≤1,wherek=-1η+β0στ.LetB=0. Then, the following equations satisfy simultaneously
(37)Φ˙=ΦB,Ψ˙=-BΨ.As shown in Faria and Magalhaes [32, 33], an appropriate phase space for considering normal forms of (28) is the Banach space BC of functions from [-1,0] into ℝ2 which are uniformly continuous on [-1,0) with a jump discontinuity at 0. Then, the elements of BC have the form φ+X0ρ, where φ∈𝒞, ρ∈ℝ2, and
(38)X0(θ)={I,θ=0,0,-1≤θ<0,
so that BC is identified with 𝒞×ℝ2 with the norm |φ+X0ρ|=|φ|C+|ρ|ℝ2.Letπ:BC→P denote the projection
(39)π(φ+X0ρ)=Φ(〈Ψ,φ〉+Ψ(0)ρ),φ∈𝒞,ρ∈ℝ2,
and then the decomposition 𝒞=P⊕Q yields a decomposition of BC by Λ0 as the topological direct sum BC=P⊕Ker⁡π with the property Q⊂Ker⁡π, where Q is an infinite-dimensional complementary subspace of P and 𝒞 as shown above. Now, we decompose zt∈𝒞1 in (28) as zt=Φx(t)+y, where x(t)∈ℝ and y∈Q1=Q∩𝒞1, and 𝒞1 is the subset of 𝒞 consisting of continuously differentiable functions. We rewrite system (28) as
(40)z˙(t)=L0zt+(L(ϵ)-L0)zt+F(zt,ϵ),
and then, under the composition zt=Φx(t)+y, system (28) can be decomposed as a system of ODEs in ℝ×Ker⁡π as follows:
(41)x˙=Bx+Ψ(0)((L(ϵ)-L0)(Φx+y)+F(Φx+y,ϵ)),y˙=AQ1y+(I-π)X0×((L(ϵ)-L0)(Φx+y)+F(Φx+y,ϵ)),
where AQ1=y˙+X0(L(y)-y˙(0)) is an operator from Q1 into Ker⁡π.Considering the Taylor expansion of the functions to the right of (41), we have
(42)x˙=Bx+∑j≥2‍1j!fj1(x,y,ϵ),y˙=AQ1y+∑j≥2‍1j!fj2(x,y,ϵ),
where
(43)1j!fj1(x,y,ϵ)=Ψ(0)[1(j-1)!Lj-1(ϵ)(Φx+y)+1j!Fj(Φx+y,ϵ)],1j!fj2(x,y,ϵ)=(I-π)X0[1(j-1)!Lj-1(ϵ)(Φx+y)+1j!Fj(Φx+y,ϵ)].
It can be obtained that(44)12!f21(x,y,ϵ)=kτ[ϵση(y2(-1)-ηx)+η2σ(1-ϕ)(1-r)(1-p)ghg×((p-1)(ηx-y2(-1))2+(ηx-y2(0))((η+ω)x+y1(0)))-β0((η+ω)x+y1(-1))(ηx+y2(-1))εση],12!f22(x,y,ϵ)=(I-π)X0·(-ϵστη(y2(-1)-ηx)-β0τ((η+τ)x+y1(-1))(y2(-1)-ηx)-η2τσ(1-ϕ)(1-r)(1-p)[(p-1)(y2(0)-ηx)2-((η+ω)x+y1(0))(y2(0)-ηx)]ϵστη(y2(-1)-ηx)+β0τ((η+τ)x+y1(-1))(y2(-1)-ηx)-η2τσ(1-ϕ)(1-r)(1-p)[(p-1)(y2(0)-ηx)2-((η+ω)x+y1(0))(y2(0)-ηx)]).As for autonomous ODEs inℝ2, the normal forms are obtained by a recursive process of changes of variables. At a step j, the terms of order j≥2 are computed from the terms of the same order and from the terms of lower orders already computed in previous steps. Assume that steps of orders 2,3,…,j-1 have already been performed; this leads to
(45)x˙=Bx+∑l≥2j-1‍1l!gl1(x,y,ϵ)+1j!f~j1(x,y,ϵ)+h.o.t,y˙=AQ1y+∑l≥2j-1‍1l!gl2(x,y,ϵ)+1j!f~j2(x,y,ϵ)+h.o.t.
Following the algorithm of Faria and Magalhaes [32, 33] at step j, using a change of variables of the form
(46)(x,y)=(x^,y^)+Uj(x^,ϵ)≡(x^,y^)+[Uj1(x^,ϵ),Uj2(x^,ϵ)],
where x,x^∈ℝ, y,y^∈Q1, and Uj1:ℝ2→ℝ, Uj2:ℝ2→Q1 are homogeneous polynomials of degree j in x^ and ϵ, after dropping the hats for simplification of notations, system (42) can be put into the normal form
(47)x˙=Bx+∑j≥2‍1j!gj1(x,y,ϵ),y˙=AQ1y+∑j≥2‍1j!gj2(x,y,ϵ),
where
(48)gj1(x,y,ϵ)=f~j1(x,y,ϵ)-[DxUj1(x,ϵ)Bx-BUj1(x)],gj2(x,y,ϵ)=f~j2(x,y,ϵ)-[DxUj2(x,ϵ)Bx-AQ1Uj2(x)].
It can be verified that system (28) satisfies nonresonance conditions since Λ0={0}; see Faria and Magalhaes [32]. Then, the locally invariant manifold of system (28) tangent to P at zero must be y=0 and the flow on this manifold is given by 1-dimensional ODE as follows:
(49)x˙=Bx+12!g21(x,0,ϵ)+13!g31(x,0,ϵ)+h.o.t.
The nonlinear terms in (47) are in normal form in the classical sense with respect to matrix B. In applications, gj1(x,0,ϵ) usually can be determined by the following procedure.Definition 7.
Forj≥2, let Mj denote the operator defined in Vj(ℝ2×Kerπ), with values in the same place, by
(50)Mj(h1,h2)=(Mj1h1,Mj2h2),(Mj1h1)(x,ϵ)=Dxh1(x,ϵ)Bx-Bh1(x,ϵ),(Mj2)(x,ϵ)=Dxh2(x,ϵ)Bx-AQ1[h2(x,ϵ)],
with domain D(Mj)=Vj2(ℝ2)×Vj2(Q1). Here, we use the notation Vj2(Y) to denote the space of homogeneous polynomials of degree j in 2 variables (x,ϵ)∈ℝ2, with coefficients in a Banach space Y.According to Faria and Magalhaes [32, 33], we have
(51)Uj(x)=Mj-1PI,jf~j(x,0,ϵ)∈Ker⁡(Mj)c,gj1(x,0,ϵ)=(I-PI,j)f~j1(x,0,ϵ)∈Im⁡(Mj1)c,
where PI,j=(PI,j1,PI,j2) is the projection of Vj2(ℝ2)×Vj2(Ker⁡π) on Im⁡(Mj1)×Im⁡(Mj2).SinceB=0, it can be checked that [Im⁡(M21)]2=span⁡{x2,xϵ,ϵ2}. Then, we have
(52)12g21(x,0,ϵ)=kτ(-ϵσx+η2σ(1-ϕ)(1-r)(1-p)ghhh×((p-1)η2x2+η(η+ω)x2-β0η(η+ω)x2)η2σ)=kτ(-ϵσx+k1x2),
where k1=(η3/σ)(1-ϕ)(1-r)(1-p)((pη+ω)-β0(η+ω)). Further the normal form of (41) on the invariant local center manifold y=0 is given by
(53)x˙=kτ(-ϵσx+k1x2)+h.o.t.Furthermore, ifk1=0, we have to compute g31(x,0,ϵ). It can be obtained that Ker⁡(M21)c=0, and
(54)U2(x,ϵ)=M2-1PI,2(f21(x,0,ϵ)f22(x,0,ϵ))=(0h2(θ)(x,ϵ)),
where h2(θ)(x,ϵ)=(h2(1)(θ)(x,ϵ),h2(2)(θ)(x,ϵ))T is the unique solution in U22(Q1) of the equation
(55)M22h2(θ)(x,ϵ)=(I-π)X0(2ϵστx-2ϵστx).
Let
(56)h2(i)(θ)(x,ϵ)=∑|q|=2‍h2,q(i)(θ)(x,ϵ)q=h220(i)(θ)x2+h211(i)(θ)xϵ+h202(i)(θ)ϵ2,ghhhhhhhhhhhhhhhhhhhhhhhi=1,2.
By (55), we get
(57)Mh2(0)(xϵ)+N(0)h2(-1)(x,ϵ)-h˙2(0)(x,ϵ)=(2στϵx-2ϵτϵx),(h˙220(1)(θ)x2+h˙211(1)(θ)xϵ+h˙201(1)(θ)ϵ2h˙220(2)(θ)x2+h˙211(2)(θ)xϵ+h˙201(2)(θ)ϵ2)+(2k(η+ω)τϵx-2kηστϵx)=0,
which leads to
(58)-ητh220(1)(0)-τ(η+ω-β0ση)h220(2)(0)-β0στηh220(2)(-1)-h˙220(1)(0)=0,-ητh211(1)(0)-τ(η+ω-β0ση)h211(2)(0)-β0στηh211(2)(-1)-h˙211(1)(0)=2στ,-ητh202(1)(0)-τ(η+ω-β0ση)h202(2)(0)-β0στηh202(2)(-1)-h˙202(1)(0)=0,-β0στηh220(2)(0)+β0στηh220(2)(-1)-h˙220(1)(0)=0,-β0στηh211(2)(0)+β0στηh211(2)(-1)-h˙211(1)(0)=-2στ,-β0στηh202(2)(0)+β0στηh202(2)(-1)-h˙202(1)(0)=0.After computation, we arrive at(59)h220(θ)=(00),h202(θ)=(00),h211(θ)=(c1+aθc2+bθ),
where c1=-(β0τσ(η+ω)b/2η(η+β0στ))+(2kσω/η), c2=β0στb/2(η+β0στ), a=-2k(η+ω)στ, and b=2kηστ.Then, we have(60)U2(x,ϵ)=(0aθ+c1bθ+c2)xϵ.
Substituting the change of variables (x,y)=(x^,y^)+(1/2)U2(x^,ϵ) into (48) and dropping the hats, we have(61)13!f~31(x,0,ϵ)=στϵηΨ(0)((b2-c2)xϵ(c2-b)xϵ)+τΨ(0)(η5σ2(r-1)(p-1)(1-ϕ)(ω+pη)2x3+η(η+ω)ϵx2-η5σ2(r-1)(p-1)(1-ϕ)(ω+pη)2x3-η(η+ω)ϵx2)=kτ((c2-b)σηϵ2x-η(η+ω)ϵx2-η5σ2(r-1)(p-1)(1-ϕ)(ω+pη)2x3).Since [Im⁡(M31)]c=span⁡(x3,x2ϵ,xϵ2,ϵ3), we have
(62)13!g31(x,0,ϵ)=Proj[Im⁡(M31)]c13!f~31(x,0,ϵ)=kτσ(c2-b)ηϵ2x-kτη(η+ω)ϵx2-kτη5σ2(r-1)(p-1)(1-ϕ)(ω+pη)2x3,
which, together with (53), implies that when k1=0, the normal form on the invariant local center manifold y=0 is given by
(63)x˙=(-kστ+kτσ(c2-b)ηϵ)ϵx-kτη(η+ω)ϵx2-kτη5σ2(r-1)(p-1)(1-ϕ)(ω+pη)2x3+h.o.t.By the center manifold theorem and the bifurcation theorem, see Carr [34], Guckenheimer and Holmes [35], Hale and Lunel [31], Wiggins [30], and the references cited therein, the dynamics of the delayed differential (9) is topologically equivalent to that of (63) at the sufficiently small neighborhood of ϵ=0. Therefore, by the normal forms on the center manifold equations (9) and (63), the following results can be obtained immediately.Theorem 8.
Assume thatR0=1; that is, ϵ=0. (1)
Ifk1≠0, the zero solution of (9) is unstable.(2)
Ifk1=0, the zero solution of (9) is stable due to (kτη5/σ2)(r-1)(p-1)(1-ϕ)(ω+pη)2>0.Theorem 9.
Assume thatβ=β0+ϵ and |ϵ| is a sufficiently small positive number.(1)
Ifk1≠0, (9) undergoes a transcritical bifurcation at the zero solution.
(2)
Ifk1=0, (9) undergoes a generic pitchfork bifurcation at the zero solution.


## 4. Sensitive Analysis

Sensitivity indices allow us to measure the relative change in a variable when a parameter changes. The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives. Here, one adopts the following definition as described by Chitnis et al. [36].Definition 10.
The normalized forward sensitivity index of a variable,u, that depends differentiably on a parameter, p, is defined as γpu=(∂u/∂p)×(p/u).Generally speaking, initial disease transmission is directly related to the threshold parameterR0, and the disease prevalence is directly related to the positive equilibrium E*, specifically to the magnitude of I*. In the following, we calculate the sensitivity indexes of R0 and the positive equilibrium E* to the parameters in system (6). The values of parameters used in this section are chosen from the literatures Gibson and Aritua [3], Holt and Chancellor [37], and van den Bosch et al. [10], as shown in Table 1.Table 1
Description of variables and parameters in model (6).
Parameter
Description
Estimated values
Reference
ϕ
Fraction planted from virus free material
—
—
r
Ratio of reversion
0.08
van den Bosch et al. [10]
p
Probability of detecting an infected cutting
(
0,0.2316
)
van den Bosch et al. [10]
η
Harvest rate
2
×
1
0
-
3
Holt and Chancellor [37]
β
Transmission rate
6.4
×
1
0
-
3
van den Bosch et al. [10]
σ
Planting rate
1.5
×
1
0
-
3
Gibson and Aritua [3]
ω
Roguing rate
3.3
×
1
0
-
4
van den Bosch et al. [10]Table2 shows the sensitivity indices of R0 with respect to each of the seven parameters in system (6). The indices with positive signs increase the value of R0 when they are increased and those having negative signs decrease the value of R0 when they are increased. It can be seen from Table 2 that γϕR0∈[-0.338963,0], γrR0∈[-0.022013,0], γpR0∈[-0.033138,0], γηR0∈[-1.860169,-1.353863], and γβR0=γσR0∈[0.746847,1] as ϕ is varying in [0,1]. The most sensitive parameters to R0 are the transmission rate β and the planting rate σ, which means that increasing (or decreasing) β and σ by 10% increases (or decreases) R0 by more than 7.46%. Obviously, to maintain the reproduction of a plant, it is not practical to reduce R0 by decreasing the planting rate σ. We also observe that R0 is most sensitive to the harvest rate η in an inversely proportional way. That is to say R0 is a decreasing function of η, where increasing η will decrease R0 and decreasing η will increase R0. Furthermore R0 is a decreasing function of ϕ, r, p, and ω. These results suggest that intervention strategies should be targeted to transmission rate β, harvest time 1/η, fraction planted fromin vitro propagated, virus free, material ϕ, the ratio of reversion r, the probability of detecting an infected cutting p, and the roguing rate ω.Table 2
Sensitivity indices ofR0 to all parameters.
Parameterj
ϕ
r
p
η
β
σ
ω
γ
j
R
0
3
×
1
0
-
6
ϕ
l
1
(
ϕ
)
3
×
1
0
-
7
(
ϕ
-
1
)
l
1
(
ϕ
)
4
×
1
0
-
7
(
ϕ
-
1
)
l
1
(
ϕ
)
-
l
2
(
ϕ
)
l
1
(
ϕ
)
1
×
1
0
-
5
l
1
(
ϕ
)
1
×
1
0
-
5
l
1
(
ϕ
)
−0.139831
wherel1(ϕ)=1.3×10-5-3×10-6ϕ and l2(ϕ)=1.7×10-5+5×10-7ϕ.As the analytical expressions for the sensitivity indices of the positive equilibriumE* do not possess good structures, we therefore evaluate the sensitivity indices at the parameter values in Table 1, and then show the sensitivity indices of the state variables for low and high planting rates of virus free material (e.g., ϕ=0.2,0.8) in Tables 3 and 4, respectively.Table 3
Sensitivity indices ofE* with ϕ=0.2.
Parameter
S
*
I
*
N
*
ϕ
+
0.663446
-
0.110946
+
0.552500
r
+
0.230764
-
0.038590
+
0.192174
p
+
0.107352
-
0.010476
+
0.096876
Wη
+
1.006333
-
1.195680
-
0.189352
β
-
1.114529
+
0.163789
-
0.928150
σ
-
0.114529
+
1.186579
+
1.072050
ω
+
0.222730
-
0.177077
+
0.045653Table 4
Sensitivity indices ofE* with ϕ=0.8.
Parameter
S
*
I
*
N
*
ϕ
+
0.797581
-
0.543699
+
0.253882
r
+
0.017339
-
0.011820
+
0.005519
p
+
0.041648
-
0.008413
+
0.033235
η
+
0.888967
-
2.147856
-
1.258887
β
-
1.022566
+
0.697064
-
0.325502
σ
-
0.022566
+
1.697066
+
1.674500
ω
+
0.156165
-
0.246287
-
0.090121It can be seen from Tables3 and 4 that the most sensitive parameter to the population size N* is the planting rate σ for both the low and high proportions of virus free material, followed in an inverse way by the transmission rate β and the ratio of reversion r. The most sensitive parameter to the infected plant I* is the harvest rate η in an inverse way, followed by the planting rate σ. The most sensitive parameter to the susceptible plant S* is the transmission rate β in an inverse way, followed by the harvest rate η.


## 5. Number Simulations

Numerical simulation is an effective method to investigate the properties of system (6). In this section, we present numerical simulations to support the theoretical analysis given in Sections 2 and 3.Firstly, letτ=0, β=0.01, ϕ=0.1, p=0.01, r=0.01, ω=0.01, σ=0.01, and η=0.01. Then, it is easy to obtain that R0≈0.94<1. From the second case of Theorem 2, the disease-free equilibrium is asymptotically stable if R0<1 and τ=0. Figure 1 exhibits that all of the solutions approach the disease-free equilibrium E0.τ
=
0, R0<1. The solutions of system (6) approach to the disease-free equilibrium E0=(1,0).
(a)
(b)
(c)Secondly, letτ=0.5, ϕ=0.1, p=0.1, r=0.1, η=0.1, ω=0.001, σ=0.01, and β=0.01. Then, it is easy to obtain that R0≈0.73<1. From Theorem 4, the disease-free equilibrium is asymptotically stable if R0<1 and τ>0. Figure 2 exhibits that all of the solutions approach the disease-free equilibrium E0.τ
>
0, R0<1. The solutions of system (6) tends to the disease-free equilibrium E0=(0.1,0).
(a)
(b)
(c)Finally, letτ=0.5, ϕ=0.1, p=0.1, r=0.1, η=0.1, ω=0.01, σ≈0.0204, and β≈0.1818. Then, we have R0=1 and k1=0. From the second case of Theorem 9, the disease-free equilibrium is stable if R0=1 and k1=0. Figure 3 exhibits that all of the solutions approach the disease-free equilibrium E0.τ
>
0, R0=1. The solutions of system (6) tends to the disease-free equilibrium E0=(0.20405,0).
(a)
(b)
(c)


## 6. Conclusions

It is still an interesting topic to determine how the time delay affects overall disease progression and, mathematically, how the delay affects the dynamics of systems. In this paper, a delayed disease model with continuous cultural control strategy is formulated. By analyzing the transendental characteristic equation, we derive some conditions on the stability of the disease-free equilibrium. The disease-free equilibrium is stable, ifR0<1, and unstable if R0>1, 0≤τ<τ0. Then, by virtue of the methods in Guckenheimer and Holmes [35] and Wiggins [30], the stability of the nonhyperbolic disease-free equilibrium is discussed if R0=1 and τ=0. Furthermore, if R0=1, τ>0 and ≠1/q1, zero is a simple root of the characteristic equation (20), and the delayed model (6) exhibits rich dynamics. Employing the center manifold reduction and normal form theory due to Faria and Magalhaes [32, 33], we compute the normal form of system (6) associated with zero eigenvalue. By analyzing the normal form, we get that the disease-free equilibrium is unstable if k1≠0 and stable if k1=0. From the view of bifurcation, system (6) undergoes transcritical bifurcation if k1≠0 and generic pitchfork bifurcation if k1=0 at the origin with small ϵ.Then, the sensitivity analysis ofR0 and the positive equilibrium E* is carried out in order to determine the relative importance of different factors responsible for disease transmission. The results show that the most sensitive parameters to R0 are the transmission rate β and the planting rate σ, the most sensitive parameter to the state variable I* is the harvest rate η followed by the planting rate σ, the most sensitive parameter to the state variable S* is the transmission rate β followed by the harvest rate η, and the most sensitive parameter to the population size N* is the planting rate σ for both the low and high proportions of virus free material ϕ, followed by the transmission rate β and the ratio of reversion r.At last, numeric simulations are used to verify the qualitative results obtained in the present paper.



---

*Source: 207959-2014-03-16.xml*