Dynamics of a delayed plant disease model with Beddington-DeAngelis disease transmission

1.   Introduction 

Plant has to rely solely on cellular innate immunity to deal with infections as it does not possess any form of mobile defence and therefore it exhibits many plant-specific characteristics ^[1]. Plants showing recovery phenotype to a specific viral infection initially become affected but later experience new growth that is progressively more resistant to the virus until they finally produce new virus-free leaves with complete immunity ^[2]. Generally, plant gains hypersensitive resistance by triggering the self-destruction of infected cells, with necrotic tissue forming at and around the infection site. The cells surrounding these necrotic lesions are usually found in the antiviral state. Although some of them may contain traces of the virus, the virus is unable to replicate. This can be explained by the post-transcriptional gene silencing (PTGS), or RNA interference ^[3]. Besides disease resistance, time delay in disease transmission, known as incubation period, is an important factor that arises from disease resistance. The length of this delay depends on the plant's resistance ^[4].

Vector-borne plant viral diseases have attracted the interest of researchers who are working with modeling approach ^[5,6]. Mathematical modeling of vector-borne plant disease is an effective approach to understand the transmission of disease in host plants. With this approach, some effective measures can be done to make appropriate disease control strategy ^[7].

For the modeling purposes, usually, the rate of infection in most of the models for viral disease is assumed to be bilinear in the carrier vector $v$ and susceptible host $x$. However, the actual incidence rate is probably not linear over the entire range of $v$ and $x$. Therefore, it is reasonable to assume the infection rate of plant-vector model as Beddington-DeAngelis functional response, $\frac{\lambda x v}{1+ax+cv}$, where $a, \; c, \; \lambda \geq 0$ are constants. The functional response was introduced by Beddington ^[8] and DeAngelis et al. ^[9]. It includes nonlinear incidence rate and bilinear functional response, Holling type II functional response and saturation response ^[10,11]. Mathematical models using Beddington-DeAngelis type infection rate have been proposed and analysed for viral disease dynamics ^[12,13].

In plant epidemiology, Holt et al. have studied the dynamics of mosaic disease of Cassava which is also carried through whitefly ^[14]. Venturino et al. have proposed a mathematical model for studying the dynamics of Jatropha curcas mosaic disease vectored by whitefly ^[15] and observed oscillations in the system due to high infection rate. In ^[16], authors have proposed a model for the dynamics of soil-borne plant disease with host demography and shown the limit cycle behaviour. In the above mentioned articles authors have used bilinear disease transmission in modeling process. On the other hand, in ^[17,18], authors have taken nonlinear incidence rate for vector-borne plant disease propagation.

Another important factor in a biological system is the time delay. Vector-borne plant disease models with delay can show stability switches, periodic oscillations, transcritical bifurcation etc. ^[14,15,16]. However, such an analysis is crucial to assess the persistence of infection and disease control.

Van der Plank ^[19] first incorporated delay differential equation with the single delay in plant epidemics assuming the time evolution of the density of infected host tissue. Despite many studies that followed Van der Plank's model and that adopted discrete time approximations of it, delay differential equations are rarely used in theoretical studies of plant disease. In ^[20], the authors introduced plant incubation period as a delay in a plant disease model proposed by Meng et al. ^[21] to notice the changes in the dynamics of the model. Also in ^[16], authors have proposed a model for the dynamics of soil-borne plant disease taking a time delay due to latent period of inoculum/vectors. In ^[22], the authors modified the model by Shi et al. ^[23] for vector-borne disease in plants by including multiple delays to account for the incubation periods of plants and latent period in vectors. They have also noticed the changes in the dynamics of the system, specifically changes in the solutions trajectories. They have derived the threshold value for delay induced destabilization. But they have missed important dynamics that delay can induce misleading the conclusions that is established in ^[24,25]. Large delay can stabilise a system or make the system disease-free ^[26]. Therefore, construction of vector-borne plant disease models with time delays requires a clear understanding of the relevant biological mechanism.

In this article, motivated by the above discussions, we have derived the mathematical model for plants viral disease including plants resistance and incubation delay (section 2). Then some basic properties of the model have been analysed (section 3). The effect of delay on the stability of equilibria has been focused (section 4, 5). Also, Hopf bifurcation due to the stability changes has been analysed and the outcomes are discussed and shown in figures (section 6). Finally, a discussion of the main results concludes the paper (section 7).

2.   The mathematical model 

Plant and vector populations are considered without explicitly including the mosaic virus. Here, $x(t)$, $y(t)$ denote the abundances of susceptible and infected plants ($m^{-2}$) respectively. Number of infected vector (virus carrier) is denoted by $v(t)$ (plant$^{-1}$) at any time $t$.

Logistic growth for susceptible plant due to the finite area of the plantation is specified by net growth rate $r$ and carrying capacity $k$. A successful contact between infected vectors with a healthy plant makes it infected ^[27], $m$ is the natural death rate, $d$ is disease related death rate of infected plants.

Let $\lambda$ be the maximum contact rate from an infected vector to a susceptible plant. Usually the rate of infection in most of the models for viral disease is assumed to be bilinear in the carrier vector $v$ and susceptible host $x$. However, the actual incidence rate is probably not linear over the entire range of $v$ and $x$. Thus, from the above discussion, it is reasonable to assume a Beddington-DeAngelis type infection rate i.e., $\frac{\lambda x v}{1+ax+cv}$, where $a$ is the plant resistance rate, $c$ is rate for the crowding effect of vectors ^[12,13].

Insect vector is recruited in the system at a rate $b$ and is proportional to infected plant population. The infective insects do not get sick from the virus, they are just carriers ^[28]. Let $\mu$ be the mortality rate of infected vector.

The effect of time delay in terms of incubation period of plant is modeled via the term $\lambda e^{-m\tau}\frac{x(t-\tau)v(t-\tau)}{1+ax(t-\tau)+cv(t-\tau)}$. The term $e^{-m\tau}$ represents the survival probability of plant through the incubation time $\tau$.

Based on the above assumptions, we have the following mathematical model:

Let $C$ denote the Banach space of continuous functions $\phi : [-\tau, 0] \rightarrow \mathbb{R}^3$ equipped with the sup-norm,

where, $\phi = (\phi_ 1, \phi_2, \phi_3)\in C([-\tau, 0], \mathbb{R}^3)$. For biological reasons, populations always have nonnegative values, therefore, the initial function for model (2.1) is taken as below:

3.   Basic properties of the model 

In this section, some basic properties such as non-negativity, boundedness of solutions and possible equilibria are analysed.

3.1.   Non-negativity and boundedness 

Lemma 1. All solutions of system (2.1) are non-negative on $[0, +\infty)$ under the initial conditions (2.1) and are ultimately bounded.

Proof. We rewrite the first equation of (2.1) as

Let $f(x, y, v) = r\left(1-\frac{x+y}{k}\right) -\frac{\lambda v}{1+ax+cv}$, then from above equation, we have

Since, $x(0) = \phi_1(0) > 0$, then $x(t) > 0$ for $t \geq 0.$

Using the method of steps, we show the non-negativity of $y$ and $v$.

The second equation of system (2.1) gives,

for some constant $T > 0$ and for all $t \in (0, T]$. Using the standard comparison principle, we have $y(t)\geq 0$ and $v(t) \geq 0$ for all $t \in (0, T]$.

The above argument can be repeated to deduce the non-negativity of $y$ and $v$ on the interval $t\in (T, 2T]$ and so on the successive intervals $t \in (nT, (n + 1)T ], n = 2, 3, \dots$ to include all positive times.

To ensure the model biologically plausible, both plant and vector populations have to remain bounded during their time evolution.

To prove the boundedness of $x(t)$ and $y(t)$, let

Then,

Now, $(r+q)x-\frac{rx^2}{k}$ is a quadratic in $x$ and its maximum value is $\frac{k(r+q)^2}{4r}$ and $e^{-m\tau}\leq1$. Using this, we can write

This implies

and hence, the $x(t)$ and $y(t)$ are bounded. As $y$ is bounded, using the well-known comparison principle, from the third equation of (2.1), $v(t)$ is also bounded.

3.2.   Basic reproduction number 

We follow the method established in the paper by Hefferman et al. ^[29] for calculating $\mathcal{R}_0$.

We consider the next generation matrix $G$ comprised of two parts namely $F$ and $V$, $F_i$ are the new infections, while the $V_i$ transfers of infections from one compartment to another.

and

where $E_1(k, 0, 0)$ is the disease-free equilibrium and indices $i$, $j$ correspond, respectively, to $y$ and $v$. The basic reproduction number $\mathcal{R}_0$ is the dominant eigenvalue of the matrix $G = FV^{-1}$, i.e., Eq (3.3).

Remark 1. It is important to note that $\mathcal{R}_0$ does not depend on $c$ (crowding effect) and hence it is independent that how efficiently vectors interact with healthy plants. At the same time, taking the limit of $k\to\infty$ gives

In light of the fact that $\mathcal{R}_0$ is monotonically decreasing with increasing $a$, this suggests that eradication of disease, as represented by a stable disease-free steady state $E_1$ is possible if $\mathcal{R}_0^{\infty} < 1$. The available means to achieve this is by increasing the resistance rate or incubation period. This can be done using nutrients, biocontrolling agents, genetic engineering, chemicals etc. ^{[31,32,33,34]}.

3.3.   Equilibria 

The model (2.1) has three equilibria namely (ⅰ) the trivial equilibrium $E_0(0, 0, 0)$, (ⅱ) the disease-free equilibrium, $E_1(k, 0, 0)$, and (ⅲ) the unique coexistence equilibrium, here $E^*(x^*, y^*, v^*)$ is determined as follows.

From the third equation of (2.1), we can write $v^*$ as

Considering the summation of the first and second equation times $e^{m\tau}$ of (2.1), we get $y^*$ in term of $x^*$ as

Putting the value of $v^*$ and $y^*$ in the first equation of (2.1), we get $x^*$ as the positive root of the following equation:

where,

and $q = (m+d)$.

We can check that $L(0) > 0$ always and $L(k) < 0$ when $\mathcal{R}_0 > 1$. Hence there exists $x^*$ satisfying $0 < x^* < k$ if $\mathcal{R}_0 > 1$. Note that $\mathcal{R}_0 > 1$ implies that $b\lambda e^{-m\tau} > a\mu q$, which shows $L_0 < 0$ when $\mathcal{R}_0 > 1$. Therefore the solution $x^*$ of $L(x^*) = 0$ is unique and satisfies $0 < x^* < k$ which ensures $y^* > 0$, from (3.5).

Thus we have the following result for the feasibility of endemic steady state $E^*$.

Proposition 1. Let us assume $\mathcal{R}_0 > 1$. Then there exists a unique endemic equilibrium $E^*$ of system (2.1).

3.4.   Characteristic equation 

In this section, the local stability of the delayed system (2.1) is studied around the equilibria. For this, we linearize the system (2.1) about any equilibrium $E(x, y, v)$ and we get

Here $P, Q$ are $3\times3$ matrices, given below:

The characteristic equation of the delay system (2.1), at any equilibrium point, is given by

This gives

Here,

The trivial equilibrium $E_0(0, 0, 0)$ is always unstable in view of the fact that one eigenvalue ($\xi = r$) is positive. Now we study the stability of disease-free and endemic equilibria.

4.   Stability of disease-free equilibrium 

For the disease-free equilibrium $E_1$,

where $p_{22} = -(d+m)$. Now, $\mid {\xi}{I}-P-e^{-{\xi}{\tau}}Q\mid = 0$ gives the following characteristic equation at $E_1$,

The following two cases are raised: (ⅰ) $\tau = 0$ and (ⅱ) $\tau > 0$.

(ⅰ) For $\tau = 0$, the characteristic equation is

Thus, at the disease-free equilibrium $E_1(k, 0, 0)$, one eigenvalue is $-r < 0$, and the rest of the eigenvalues will be negative or with negative real part if $(m+d)\mu-\frac{b \lambda k}{1+ak} > 0$, which is equivalent to $\mathcal{R}_0 < 1$. Thus for the stability of $E_1$, we have the following theorem.

Theorem 1. For $\tau = 0$, the disease-free equilibrium $E_1$ of the model (2.1) is stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

(ⅱ) For $\tau > 0$, one eigenvalue is $-r < 0$, and the remaining roots satisfy the quadratic equation

and we have the following theorem.

Theorem 2. Disease-free equilibrium $E_1$ of the model (2.1) for $\tau > 0$ is asymptotically stable for $\mathcal{R}_0 < 1$, and unstable for $\mathcal{R}_0 > 1$.

Proof. Assume that $\mathcal{R}_0 > 1$, then $H(0, \tau) = (m+d)\mu - \frac{b k\lambda e^{-m\tau}}{1+ak} = (m+d)\mu(1-\mathcal{R}_0) < 0$. Since $\lim_{\xi\rightarrow\infty}H(\xi, \tau) = \infty$, there exists at least one real positive root of the characteristic equation (4.3), implying that $E_1$ is unstable.

Next, let us assume $\mathcal{R}_0 < 1$. We note that in this case for $\tau = 0$ all characteristic roots of Eq (4.3) have negative real part (Theorem 1) and therefore, our goal is to show that as $\tau > 0$, none of the characteristic roots can reach the imaginary axis.

Let us assume by contradiction that for some $\tau > 0$, $\xi = i\kappa$ is a root of (4.3). Putting $\xi = i\kappa$ in (4.3) and separating real and imaginary parts, we have

Here $M_1 = (m+d+\mu)$, $M_2 = (m+d)\mu$, and $M_3 = \frac{b k\lambda e^{-m\tau}}{1+ak}$. Squaring and then adding the above equations, we finally obtain

Note that $\mathcal{R}_0 < 1$ implies $M_3 < M_2$. Since this equation has no real roots for $\kappa$, it means that the characteristic equation (4.3) cannot have purely imaginary roots. Thus, for $\mathcal{R}_0 < 1$, the steady state $E_1$ is asymptotically stable for all $\tau\geq 0$.

5.   Stability of Endemic equilibrium 

For the stability of endemic equilibrium $E^*$, again we consider two cases, $\tau = 0$ and $\tau > 0$.

When $\tau = 0,$ the characteristic at $E^*$ equation becomes

Denoting by $\sigma_1 = m_1, \; \sigma_2 = (m_2+m_4), \; \sigma_3 = (m_3+m_5)$ and using the Routh-Hurwitz criterion ^[30], for $\tau = 0$, we can get the following theorem.

Theorem 3. For $\tau = 0$, the coexistence equilibrium point $E^*(x^*, y^*, v^*)$ is stable if the conditions

are satisfied.

Now, we study the local Hopf bifurcation of $E^*$. Any of the parameters of the model may be a bifurcation parameter. We assume, $\alpha$ as the generic bifurcating parameter of the system. The following theorem characterises the possible occurrence of Hopf bifurcation for $\tau = 0$.

Theorem 4. The system (2.1) with $\tau = 0$, undergoes Hopf bifurcation around the endemic equilibrium point $E^*$ at $\alpha = \alpha^*$, where the critical parameter value is contained in the following domain

Primes denote the derivative with respect to $\alpha$.

Proof. Using the condition, $\sigma_1 \sigma_2 -\sigma_3 = 0$, the characteristic equation (5.1) becomes

which has three roots namely $\xi_1 = +i\sqrt{\sigma_2}, \; \xi_2 = -i\sqrt{\sigma_2}$ and $\xi_3 = -\sigma_1$. Therefore, a pair of purely imaginary eigenvalues exists for $\sigma_1 \sigma_2 -\sigma_3 = 0$. Now we verify the transversality condition.

Differentiating the characteristic equation (5.1) with respect to $\alpha$, we have

Therefore,

and thus Hopf bifurcation occurs at $\alpha = \alpha^*$.

For $\tau > 0$, the characteristic equation at $E^*$ becomes

which is a transcendental equation in $\xi$ with infinitely many roots. It is known that $E^*$ is locally stable (or unstable) if all the roots of the corresponding characteristic equation have negative real parts (or at least one root has a positive real part). Also, Hopf bifurcating periodic solution exists if purely imaginary roots exist. We shall check the possible occurrence of Hopf bifurcation.

For this, suppose $\xi = i w(\tau)$ is a root of the Eq (5.6), then we have

Separating real and imaginary parts, we obtain the following equations,

Squaring and adding the real and imaginary parts, we get

Denote $z = w^2$ and the Eq (5.10) becomes

where,

By Descartes rule of signs, we can say that if the conditions $S_1 \geq0, \; S_2\geq0, \; S_3\geq0$ are satisfied, then it is easy to show that the Eq (5.11) has no positive real roots.

Lemma 2. If the conditions (i) $\sigma_1 > 0, \; \sigma_3 > 0$ and $\sigma_1\sigma_2-\sigma_3 > 0,$ and (ii) $S_1\geq0, \; S_3\geq 0$ and $S_2\geq0$ are satisfied then the infected steady state $E^*$ is asymptotically stable for all $\tau\geq0$.

Remark 2. We can say that if the model parameters satisfy the conditions given in Lemma 2, the steady state $E^*$ of the delay model (2.1) is asymptotically stable for all $\tau > 0$ i.e., the stability of the system is independent of delay.

If $S_3 < 0$, then Eq (5.11) satisfies $F(0) < 0$ and $\lim_{z\rightarrow\infty} F(z) = \infty$. Thus Eq (5.11) has at least one positive root, $r_*$ for $S_3 < 0$, and consequently the characteristic equation (5.6) has a pair of purely imaginary roots $\pm i \sqrt{r_*} = iw,$ say.

We have the following lemmas.

Lemma 3. If $S_3 < 0$, then there exists at least one positive root of Eq (5.11).

Let us assume (5.11) has three positive roots denoted as $r_k, \; k = 1, 2, 3$. Then (5.10) has three positive roots $w_r = \sqrt{r_k}$. Let us define

For $w(\tau_0) = w_0$, the Eqs (5.8) and (5.9) give,

if $(m_1m_4-m_3)w^3_0-(m_2m_3-m_4m_5)w_0^2 > 0,$ and

if $(m_1m_4-m_3)w^3_0-(m_2m_3-m_4m_5)w_0^2\leq0$. Thus we have the following theorem.

Theorem 5. Suppose that the interior equilibrium point $E^*$ exists and is locally asymptotically stable for $\tau = 0$ i.e., $\sigma_1 > 0, \; \sigma_3 > 0$ and $\sigma_1\sigma_2-\sigma_3 > 0$ hold. If $S_3 < 0$, then $E^*$ is asymptotically stable for $\tau < \tau_0$ and unstable for $\tau > \tau_0$, where

When $\tau = \tau_0$, Hopf bifurcation occurs, i.e., a family of periodic solution bifurcates at $E^*$ as $\tau$ passes through the critical value $\tau_0$ provided the transversality condition,

is fulfilled.

Proof. Differentiating (5.8) and (5.9) with respect to $\tau$, we have

where,

Therefore, $\frac{d}{d\tau}[Re \{\lambda(\tau)\}]_{\tau = \tau_0, w = w_0}\neq0$ if $P(w_0)R(w_0)-Q(w_0)S(w_0) \neq 0.$ Thus, the transversality condition is satisfied and hence Hopf bifurcation occurs at $\tau = \tau_0$.

Remark 3. Theorem 5 only holds if (5.11) has at least one positive real root. Lemma 3 guarantees the occurrence of a positive root. However, when $S_3\geq0$, it is impractical to consider the analytical distribution of roots.

6.   Numerical simulations 

In this section, the analytical results of the previous sections are corroborated by numerical simulations. We have varied the important parameters, within their feasible range, to observe their impact on both the solution trajectories and equilibria of the model (2.1). We have plotted Figures 2–5 using $ode45$ solver in MATLAB. Figure 6 is plotted using the sign of eigenvalues at an equilibrium point. The values of parameters used in the numerical simulations are given in Table 1.

Figure 1 shows how the steady state values of the infected plant change with the infection transmission rate $\lambda$. This figure shows that the endemic steady state $E^*$ becomes feasible for some minimum value of $\lambda$ that corresponds to $R_0 = 1$ as discussed in Proposition 1. For values of $\lambda$ just above its critical value $\lambda^*$, the steady state $E^*$ is unstable.

In Figure 2, numerical solution of the system (2.1) is plotted using different values of $\lambda$ (the disease transmission rate from vectors to plant) for $\tau = 0$. The endemic equilibrium, $E^*$, is stable for some lower value of $\lambda$ and system has a periodic orbit for $\lambda = 0.4$. Oscillation increases the value of $\lambda$. In particular, there exists a critical value of $\lambda$ below which $E^*$ is stable but if $\lambda$ crosses this value then the system undergoes bifurcation and a stable periodic solution appears which is clearly presented in the figure. From Figure 3 and using Theorem 4, we observe that Hopf bifurcation occurs at $\lambda = 0.3667$. The higher transmission rate causes larger amplitude of oscillations around this steady state $E^*$.

In Figure 4, numerical solution of the system is plotted for three different value of delay $\tau$ for a fixed value of $\lambda = 0.35$ i.e., the $E^*$ is stable for $\tau = 0$ at that value of $\lambda$. The system (2.1) has periodic orbits when $\tau$ lies in the range $\tau\in(0, 6)$. That is endemic equilibrium $E^*$ is stable for higher values of delay. It is clear from Figure 5 that the stability switches of the endemic equilibrium $E^*$ in terms of time delay $\tau$ at two points, nearly at $\tau \approx1.012$ days and $\tau\approx5.56$ days from stable to unstable and from unstable to stable, respectively.

Figure 6 illustrates the stability of equilibrium points in different parameter planes. It is observed that the stability of equilibria depends on time delay $\tau$, the plant infection rate $\lambda$ and vector growth rate $b$ and immune response rate $a$. This figure shows that in the parameter region with $\mathcal{R}_0 < 1$, the system is free of disease, as signified by the stable disease-free equilibrium $E_1$, and in the region with $\mathcal{R}_0 > 1$, the disease is present i.e., the endemic steady state $E^*$ is feasible there.

Figure 6(a) shows that disease-free state is stable when the product of infection rate $\lambda$ and infective vector growth rate $b$ is below some threshold value. It also illustrates that when the product of transmission rate $\lambda$ and $b$ crosses some threshold value, the endemic steady state becomes feasible and stable but it loses its stability for higher values of both $\lambda$ and $b$. Figure 6(b), (c) suggest that disease transmission significantly depends on the resistance of plant and the incubation delay. Critical value of $\lambda$ for Hopf bifurcation depends on these two parameter. Region of stability of equilibria widens as both of them are increased. Figure 6(d) shows that $E^*$ is unstable for small $a$ and then becomes stable for larger time delay $\tau$. It is also clear that disease resistance rate $a$ and delay parameter $\tau$ have the ability to transform the unstable endemic into a stable disease-free state.

7.   Discussions and conclusions 

In this article, we have formulated a mathematical model using Beddington-DeAngelis functional response type infection rate and studied the dynamics of vector-borne plant disease. We have derived the basic reproduction number, $\mathcal{R}_0$, that helps us in understanding the dynamics of the model depending on the time delay and disease resistance. In the case where $\mathcal{R}_0 < 1$, the disease-free state, $E_1$, is stable. This is also corroborated by the numerical calculations. When $\mathcal{R}_0 > 1$, however, the disease-free state becomes unstable and the endemic steady state becomes feasible. Oscillations in the solution of the susceptible plant, infected plant and infective vectors have been observed around the endemic equilibrium point. Resistance rate of the vector has been chosen as the continuation parameter, since it is one value that can be enhanced using nutrients which prevent the disease from occurring ^[28,34].

Numerical simulations have shown that the model can support resistant and recovery-type behaviours. We also found that the delay preserves the stability property of disease-free equilibrium by failing to induce stability switches when the equilibrium is stable. It is interesting to see that when endemic equilibrium is unstable and a periodic solution emerges at the endemic equilibrium (without delay), large delay or the resistance rate leads to the disappearance of the periodic solutions and later to disease-free situation ^[35]. This can be done using different control strategies (Remark 1).

In conclusion, Beddington-DeAngelis functional response type disease transmission rate has been successfully fitted in the model. The predictions from the model concerning particular dynamical regimes depending on parameters can provide the managers with effective methods for disease monitoring, as well as for developing proper techniques for successful disease management.

Acknowledgement 

YT's research is supported by Japan Society for the Promotion of Science "Grant-in-Aid 20K03755".

Conflict of interest 

All authors declare no conflicts of interest in this paper.