Stability and bifurcation analysis of a multi-delay model for mosaic disease transmission

1.   Introduction

Mosaic disease is a plant disease that is caused by the viruses of the Begomovirus family ^[1]. It affects the productivity of economically important plants such as Cassava, Jatropha spp., Tomato etc ^[2]. This disease manifests itself in the substantial leaf damage, such as yellowing of the leaves and sap drainage, and it attacks the fruits, thus significantly reducing the yield of seeds ^[3,4]. Whitefly Bemisia tabaci (Gennadius) is the vector that carries the virus and transmits it to Jatropha plants ^[5,6], so the spread of the mosaic disease is largely determined by the distribution of whitefly vectors and the density of the host plants ^[7]. When an infected vector feeds on a healthy plant leaf or stem, the plant becomes infected after an incubation period. Similarly, a non-infected vector becomes infected upon feeding on an infected portion of the plant after the latest period ^[3,8].

The process of viral infection of plants is often characterized by a non-negligible duration, which represents various processes associated with production of virus particles by the vectors, feeding on the plants, entry of the virus into plant cells etc. ^[9]. When considering plant disease, it is difficult in practice to distinguish between different stages of the plant infectious state. Thus, identifying appropriate latency and incubation periods from observing the plant pathology status is a challenging problem ^[9,10,11]. From a mathematical perspective, delays in the disease transmission and development of symptoms can be effectively modeled using the formalism of delay differential equations (DDEs) ^[12,13].

Mosaic disease and leaf curl disease are two of the most common vector-borne viral diseases of major agricultural crops. They are caused by aphid vectors, such as whitefly (Bemisia sp.), which cause viral diseases of cassava, Jatropha, cotton, tomato, tobacco and other plants. Since whitefly mostly transmits viral disease in a persistent manner, there is a certain associated latency period ^[14]. Unfortunately, there are scant reliable data on the duration of the incubation and latency periods for various persistently-transmitted diseases. Both of these periods vary according to viral agents and host plant species. Moreover, even for the same virus and host plants, there is a variation in latency for different whitefly species. In most cases, the incubation period in plants is relatively long (from days to weeks), while the latency period in vectors is relatively short (from a few hour to days), though it can still be significant when compared to the duration of the vector life cycle. For example, in the case of ACMV, the latent period is $6$ hours, and the incubation period is $3$ to $5$ weeks ^[15]. Including time lags associated with disease development in vectors and plants in corresponding mathematical models is essential when considering the development of disease control policies. However, one should be mindful of the fact that precise estimates of those latency and incubation periods may be hard to obtain in experimental settings ^[16].

A number of mathematical models have considered the dynamics of plant diseases, starting with the seminal work of Van der Plank ^[17] who looked at the possibility of predicting whether an epidemic outbreak can occur, and what would be its size. Jeger et al. ^[18] provide a nice review of more recent work that has looked at various aspects of interactions between plants and disease-carrying vectors. Jackson and Chen-Charpentier ^[19,20] and Al Basir et al. ^[16] have recently studied the propagation of plant viruses while accounting for two time delays, one representing the incubation period of the plant, and the other being a shorter delay due to the incubation period of vector. While retaining this setting, here we focus on the specific interactions between whitefly and Jatropha or Cassava plants, using specific parameter values for these species in the simulations. Time-delayed models have also proved effective for the analysis of within-plant dynamics of immune responses to viral infections ^[21]. Several papers have investigated the effects of vector maturation delay in the context of vector-borne animal diseases ^[22,23,24]. These papers revealed that time delays can destabilize host-vector dynamics and generate periodic solutions through a Hopf bifurcation, more complex dynamics and chaos. Banerjee and Takeuchi ^[25], Al Basir et al. ^[16,26] have provided strong biological background for the suitable incorporation of time delay in the system, to avoid reaching misleading conclusions. As a further result, it has also been shown that delays can have a stabilizing role, whereas in many of the previous models it was observed that discrete-time delays often destabilize the co-existence steady-state.

In the particular context of mosaic disease affecting the population of J. curcas plant, several models have been proposed that have looked at the dynamics of control of disease ^{[27,28,29,30,31,32]}. Venturino et al. ^[31] have proposed a mathematical model for the dynamics of Jatropha curcas plant's mosaic disease, though they focused on a constant disease transmission rate. These earlier models of mosaic disease have assumed a constant disease transmission without including a period of latency, i.e., once an infected vector has fed on healthy biomass, it immediately becomes infected. To make that approach more realistic, in this paper we propose an epidemic model for the mosaic disease that explicitly accounts for disease latency and incubation as delayed processes.

In ^[16], authors have studied a model for vector borne disease and studied the effects of delays. In this paper, we have analyzed a particular vector (whitefly) and particular disease (mosaic disease) in a particular plant (Cassava/Jatropha). We have additionally discussed the transcritical bifurcation analytically and numerically. Also, we have studied how the periods of oscillation depends on the model parameters.

The outline of the paper is as follows. In the next section, we derive a time-delayed model of mosaic disease transmission and discuss its main properties. Section 3 is devoted to analysis of feasibility and stability of different steady states of the model. In Section 4, numerical techniques are used to perform and to illustrate the behavior of the model in different dynamical regimes. The paper concludes in Section 5 with a discussion of results and future outlook.

2.   Model derivation

The following assumptions are made in order to establish the mathematical model for mosaic disease transmission.

To analyze the spread of mosaic disease on a given farm, we consider the overall plant biomass, which includes stem, leaves etc. To simplify the model, we do not explicitly include a separate compartment for mosaic virus, but instead concentrate on plant and vector populations. To this end, we introduce $S(t)$ and $I(t)$ as the amounts of healthy and infected plant biomass, respectively, while uninfected and infected whitefly populations are denoted by $U(t)$ and $V(t)$.

Due to the fact that any plantation has a finite area on which plants can grow, we will assume a logistic growth for healthy plant biomass, with the growth rate $r$ and the carrying capacity $K$. Infected vectors transmit the disease to susceptible plant at rate $\lambda$ from infected plant.

Thus the equations for plant population become:

For the dynamics of vector population, we consider $b$ to be the net growth rate of healthy vectors, and $a$ to be the maximum number of vectors that can survive on a plant. Then, $a(S+I) > 0$ is the overall carrying capacity of the vectors on all plants. We will denote, by $m$, the rate of removal of infected plant biomass, and by $\mu$ the sum of the natural and the virus-related mortality rates for infected vectors ^[31,33]. Uninfective vectors themselves pick up the disease from the infected plants at rate $\beta$ while feeding on the infected plants ^[31,33]. Thus the equation for vectors read as:

To better understand the biological reason for introducing a delayed system, we now discuss some features of the whitefly-cassava ecological interactions. Whiteflies need to feed for at least 4 hours on the young leaves of a cassava plant with mosaic before they acquire the virus. Whiteflies that acquired virus in 4–6 hours, then require another 4 hours to become viruliferous. Once viruliferous, they can infect healthy plants during a feeding period of 15 minutes, but longer periods give more infections. Adult whiteflies can continue to infect healthy plants for up to 48 hours after initial acquisition of the virus. Cassava plants fed upon even by a single viruliferous fly sometimes became infected ^[34]. Symptoms appear after a 3–5 week latent period ^[35,36,37]. Thus, transmission of mosaic disease through whitefly vector is a delayed process with two distinct delays: one describes plant incubation, and another is the delay in whitefly becoming infectious i.e., the latent period. Incubation in the infection of a plant by a vector refers to a delayed process, where a healthy biomass takes some time after an infected whitefly has fed on it to actually become infectious itself.

To model this mathematically, we represent transmission of infection from vectors to plants at time $t$ by a term $\lambda e^{-m\tau_1}S(t-\tau_1)V(t-\tau_1)$, where $m$ and $\lambda$ are positive constants. This term represents healthy plants that became infected at time $t-\tau_1$ and have survived from natural death for the duration of incubation period $\tau_1$. Similarly, introducing a latency period $\tau_2\in R^+$, the term describing the transmission of disease from infected plants to healthy vectors is taken to be $\beta e^{-\mu\tau_2}U(t-\tau_2)I(t-\tau_2)$, which represents healthy vectors infected at time $t-\tau_2$ that have survived for the duration of latency period $\tau_2$.

Note that in contrast to ^[16], here the abundance of insects depends on the availability of plants on which to feed. Mathematically this is expressed by a logistic term in the vector equation, of Leslie-Gower type, i.e. the carrying capacity depends on the plant population.

With the above assumptions, the model for the dynamics of the mosaic disease takes the following form

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

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

Biologically, this initial condition together with $\phi_2(0)+\phi_4(0) > 0$ means that at time $t = 0$, at least some number of plants and vectors are already infected. While this condition may be quite natural for the vectors, as without them carrying the disease, the model would make no sense, for plants this would also very quickly be satisfied, once the vectors start transmitting the infection. It can be shown that the solution $(S(t), I(t), U(t), V(t))^T$ of model (2.3) with the initial condition (2.4) exists and is unique on $[0, +\infty)$ ^[13,38].

Before proceeding with the analysis of model (2.3), we have to established some of its well-posedness properties.

Theorem 1. All solutions of model (2.3) with the initial condition (2.4) remain non-negative for $t\geq 0$.

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

Let $\widetilde f(S, I, V) = r\left(1-\frac{S+I}{K}\right) -{\lambda K}$, then from above equation, we have

Since, $S(0) = \phi_1(0)\geq 0$ then $S(t) \geq 0$ for $t \geq 0.$ From the second, third and fourth equation of (2.3), we can write

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

We can repeat the above argument to deduce the non-negativity of $I$, $U$ and $V$ on the interval $t\in (t_c, 2t_c]$ and so on the successive intervals $t \in (nt_c, (n + 1)t_c ], n = 2, 3, \dots$ to include all positive times. □

To ensure the model remains biologically plausible, both plant and vector populations have to remain bounded during their time evolution. If we denote by $M(t)$ the total plant biomass, i.e., $M(t) = S(t)+I(t)$, then it satisfies the equation

which implies

where $B = \max[M(0), K]$, and hence, the total plant biomass is bounded. If we denote by $N(t)$ the total plant biomass, i.e., $N(t) = U(t)+ V(t)$. Non-negativity of solutions implies that the third and fourth equation of model (2.3) can be rewritten as follows,

which implies

Thus, the region

where $B = \max[S(0)+I(0), K]$, is positively invariant and attracting, in view of the fact that the solutions exceeding the $\limsup$ will ultimately decrease, with all solutions of the model (2.3) with initial conditions (2.4) being contained inside this region for sufficiently large $t$.

3.   Steady states and their stability

The system (2.3) has up to three possible equilibria, which include

(a) the healthy plant-only equilibrium, $E_1 = (K, 0, 0, 0)$,

(b) a disease-free equilibrium, $E_2 = (K, 0, aK, 0)$,

(c) a co-existence equilibrium, $E^* = (S^*, I^*, U^*, V^*)$,

where

and $S^*$ satisfies the quartic equation

where

To ensure biological feasibility of $E^*$, i.e., positivity of all state variables, one has to require $S^* < K$.

Descartes's rule of signs ^[39,40] states that the number of real positive roots of polynomials with real coefficients is equal to the number of sign changes in those coefficients, or is less than the number of sign changes by an even number. In the case where there is a single sign change in polynomial coefficients, this theorem immediately implies the existence of a single real positive root. In the particular case of the quartic equation (3.1), since $\ell_0 > 0$ and $\ell_1 > 0$, we have the following result.

Theorem 2. If $\ell_4 < 0$ and either $\ell_3\cdot \ell_2 > 0$, or $\ell_3 < 0$ and $\ell_2 > 0$, then Eq (3.1) has a single real positive root.

Remark 1. Since the third equation of the system (2.3) contains a denominator that tends to zero as $S+I\to 0$, in order to explore whether $(0, 0, 0, 0)$ can be considered a steady state of the system, we need to investigate if the singularity at the origin can be removed. This can be done by checking whether the limit of the right-hand side of our system as $(S, I, U, V)\rightarrow (0, 0, 0, 0)$ exists and is equal to $(0, 0, 0, 0)$, in which case the system can be continuously extended to $(0, 0, 0, 0)$. If we start on a hyperplane $U = V = 0$ with $S > 0$ and $I > 0$, then the right-hand sides of the third and fourth equations are zero, and as $S\rightarrow 0$ and $I\rightarrow 0$, the right-hand sides of the first two equations also tend to zero. Hence, $\lim_{(S, I, U, V)\rightarrow (0, 0, 0, 0)}(S, I, U, V) = (0, 0, 0, 0)$ when computed on this hyperplane. On the other hand, if we start on a hyperplane $S = I = 0$ with $U > 0$ and $V > 0$, the right-hand sides of the first two equations are zero, but for the third equation we would actually be starting at $\infty$, and this will continue as $U\rightarrow0$ and $V\rightarrow0$, which indicates that the limit as $\lim_{(S, I, U, V)\rightarrow (0, 0, 0, 0)}(S, I, U, V)$ does not exist, and therefore, $(0, 0, 0, 0)$ cannot be considered a feasible steady state.

The characteristic equation for eigenvalues $\rho$ of linearisation near any steady point $E(\bar{S}, \bar{I}, \bar{U}, \bar{V})$ has the form

where the matrices $A$ and $D$ have the following nonzero entries

and

The characteristic Eq (3.2) has the explicit form

where the coefficients are given below,

The steady state $E_1$ is unstable for any parameter values, as one of its characteristic eigenvalues is equal to $b$. At the steady state $E_2$, two eigenvalues of the characteristic equation are $-b$ and $-r$, and the remaining roots satisfy the transcendental equation

If we define the basic reproduction number (a short description on how it is determined is given in Appendix A) as

then we have the following result concerning stability of the steady state $E_2$.

Theorem 3. Disease-free equilibrium $E_2$ of the model (2.3) is asymptotically stable for $R_0 < 1$, linearly neutrally stable for $R_0 = 1$, and unstable for $R_0 > 1$.

Proof. Assume $R_0 > 1$, then $H(0) = m\mu- e^{-m\tau_1}e^{-\mu\tau_2}aK^2 \beta\lambda = m\mu(1-R_0) < 0$. Since $\lim_{\rho\rightarrow \infty}H(\rho) = \infty$, there exists at least one real positive root of the characteristic equation (3.8), implying that $E_2$ is unstable.

If $R_0 = 1$, then $\rho = 0$ is a simple characteristic root of (3.8). Let $\rho = \eta+i\kappa$ by any other characteristic root, then the Eq (3.8) turns into

with $M_1 = m+\mu$, $M_2 = m\mu$, and $M_3 = aK^2\beta e^{-\mu\tau_2}\lambda e^{-m\tau_1}$. Separating real and imaginary parts gives

where $\tau_1+\tau_2 = \tau$. Squaring and adding these two equations, we obtain

which can only be satisfied for $\eta < 0$. Thus, for $R_0 = 1$, the steady state $E_2$ is linearly neutrally stable.

Finally, let us assume $R_0 < 1$, which implies $aK^2\beta e^{-\mu\tau_2} \lambda e^{-m\tau_1} < m\mu$, i.e. $M_3 < M_2$. We note that in this case for $\tau = 0$ all characteristic roots of Eq (3.8) have negative real part, and, therefore, our goal is to show that as $\tau_1 > 0$, none of the characteristic roots can reach the imaginary axis. Let us assume by contradiction that for some $\tau > 0$, $\rho = i\kappa$ is a root of (3.8). Substituting this into (3.11) yields

Since this equation has no positive real roots for $\kappa$, it means that the characteristic equation (3.8) cannot have purely imaginary roots, and hence, $R_0 < 1$ the steady state $E_2$ is linearly asymptotically stable for all $\tau_1, \tau_2\geq 0$. □

Now that we have established that the steady state $E_2$ undergoes a bifurcation as $R_0$ crosses unity, it is instructive to investigate whether this change in stability of $E_2$ is associated with the emergence of another equilibrium, in which case it is called a supercritical bifurcation ^[41,42]. While Theorem 3 applies for $\tau_1, \tau_2\geq 0$, before we investigate the stability of the co-existence equilibrium we focus our attention at the special case of $\tau_1 = 0 = \tau_2$, which will provide some initial insights.

Theorem 4. When $\tau_i = 0, \; i = 1, 2$, the disease-free steady state $E_2$ undergoes a supercritical (or forward) transcritical bifurcation at $R_0 = 1$.

Proof. To prove this result, we use the bifurcation theory methodology based on the analysis of dynamics of the system dynamics on the centre manifold ^[11,43].

For a dynamical system

with a non-hyperbolic steady state $x_0$ at a particular value of parameter $\phi = \phi_0$, this approach provides conditions for a supercritical (also named forward) or subcritical (also named backward) bifurcation depending on the signs of two coefficients, namely,

and

where ${\bf v}$ and ${\bf w}$ denote, respectively, the left and right eigenvectors corresponding to the zero eigenvalue of the Jacobian matrix of system evaluated at $x = x_0$.

For $\tau = 0$, assuming $R_0 = 1$ implies $aK^2\beta \lambda = m\mu$. Considering without loss of generality $\lambda$ as the bifurcation parameter, the criticality $R_0 = 1$ corresponds to $\lambda^* = m\mu/aK^2\beta$. At the steady state $E_2$, two eigenvalues of the characteristic equation are $-b$ and $-r$, and the remaining roots satisfy the quadratic equation

For $R_0 = 1$, one eigenvalue is zero and another is $-(m+\mu)$. Thus, for $R_0 = 1$ the disease-free equilibrium $E_2$ is a non-hyperbolic equilibrium.

The right eigenvector of $A(E_2)$ is

Similarly, the left eigenvector is

Substituting these expressions into (3.13) and (3.14) gives the following values for the coefficients $B_1$ and $B_2$:

which, in light of Theorem A in ^[43], implies that at $R_0 = 1$, the steady state $E_2$ does indeed undergo a supercritical transcritical bifurcation, i.e. as $R_0$ crosses the value of 1 from $R_0 < 1$ to $R_0 > 1$, the steady state $E_2$ loses its stability, and another positive and stable steady state arises, which, in the case of model (2.3), is the co-existence steady state $E^*$. □

Having just demonstrated that for $\tau_i = 0, \; i = 1, 2$, the steady state $E_2$ undergoes a supercritical bifurcation at $R_0 = 1$, giving rise to a stable co-existence equilibrium $E^*$, the next question is how the stability of that equilibrium changes when parameters are varied. Staying with the case $\tau_1 = 0 = \tau_2$, one can prove the following result ^[27,31].

Theorem 5. Denoting differentiation with respect to $\lambda$ by prime and letting

for $\tau_i = 0, \; i = 1, 2$, the co-existence equilibrium $E^*$ has the following properties.

(a) It is linearly asymptotically stable if and only if the following conditions are satisfied

(b) It undergoes a Hopf bifurcation at $\lambda = \lambda^*\in(0, \infty)$ if and only if

To investigate whether increasing time delay $\tau_i, \; i = 1, 2$ can affect stability of the co-existence steady state $E^*$, we look at the characteristic equation (3.3) at the steady state $E^*$. Stability changes of $E^*$ can only occur if the characteristic equation (3.3) has purely imaginary solutions.

Four cases may occur which are discussed below.

Case I: $\tau_1 > 0, \; \tau_2 = 0$

In this case the characteristic equation becomes

where,

Substituting $\rho = i\omega$ into this equation and separating real and imaginary parts gives

Squaring and adding these two equations yields the following equation for the Hopf frequency $\omega$:

where $z = \omega^2$, and

Provided the Routh-Hurwitz criterion for Eq (3.21) holds, all of its roots will have a negative real part, and hence, there will be no purely imaginary roots of the characteristic equation (3.19). Thus, we have the following result

Proposition 1. For $\tau_2 = 0$, the co-existence steady state $E^*$ is linearly asymptotically stable for all $\tau_1\geq 0$, if the conditions of (3.17) hold, and the following conditions are satisfied for all $\tau_1 > 0$,

Instead, if $\delta_4 < 0$, Eq (3.21) possesses at least one real root $z_0 > 0$, and then, $\rho = \pm i\omega_0$ with $\omega_0 = \sqrt{z_0}$ will be two roots of the characteristic equation (3.19).

From the system (3.20) one can determine the value of time delays at which this pair of purely imaginary roots occurs,

Of course, $\delta_4 < 0$ is just one possibility where the characteristic Eq (3.19) has a pair of purely complex eigenvalues, and this can also happen for $\delta_4 > 0$, provided some of $\delta_1$, $\delta_2$, or $\delta_3$ are sufficiently negative to ensure that equation (3.21) has at least one real positive root. Whenever this happens, we have the following result.

Theorem 6. For $\tau_2 = 0$, the steady state $E^*$ undergoes a Hopf bifurcation at $\tau_1 = \tau_1^{*}$, provided the following condition holds

where

Proof. In light of the above analysis, at $\tau_1 = \tau_1^{*}$, the characteristic equation (3.19) has a pair of purely imaginary eigenvalues. Hence, to complete the proof of the theorem, it remains to prove the transversality condition, i.e., that the characteristic eigenvalues cross the imaginary axis. To this end, we differentiate characteristic equation (3.19) with respect to $\tau_1$ to obtain

Evaluating this at $\tau_1 = \tau_1^{*}$ and using relations (3.20), we find

Since the denominator of this expression is always positive, if the condition of the theorem holds, this means that the transversality condition

is satisfied, and thus, the steady state $E^*$ undergoes a Hopf bifurcation at $\tau_1 = \tau_1^*$. □

Case II: When $\tau_1 = 0, \tau_2 > 0$

For $\tau_1 = 0$, the characteristic equation in this case takes the following form,

Here, $A_1 = l_1, \; \; A_2 = (l_2+b_1), \; \; A_3 = (l_3+b_2), \; \; A_4 = (l_4+b_3)$, $C_1 = a_1+q_1, \; \; C_2 = a_2+q_2, \; \; C_3 = a_3+q_3$.

The analysis is similar to Case I, therefore we will not show it.

Case III: $\tau_1$ is fixed in the interval $(0, \tau_1^*)$ and $\tau_2 > 0$

When $\tau_1 > 0, \; \tau_2 > 0$, the characteristic equation is:

The effects of both delays can be analyzed using the following theorem.

Theorem 7. Suppose that the endemic equilibrium is asymptotically stable for $\tau_1\in(0, \tau_1^*)$. Now, if $\delta_4 < 0$ holds, then there exists $\tau_2^0$, for which a stability switch occurs at $E^*$, when $\tau_2$ crosses the critical value of $\tau_2^0$. Furthermore, $E^*$ will undergo a Hopf-bifurcation when $\tau_2 = \tau_2^0$, provided that

A formal proof of the Theorem 7 can be found in the article by Misra et al. ^[44] (see Theorems 5 and 6 of ^[44]). A similar result can be provided for $\tau_1$ when $\tau_2\in(0, \tau_2^*)$.

Remark 2. At the endemic equilibrium $E^*$, the characteristic equation (3.26) has delay-dependent coefficients (i.e. for $\tau_i, i = 1, 2$) and it is quite involved. Therefore, it is difficult to analytically obtain information on the nature of the eigenvalues and on the conditions for occurrence of stability switches. But the nature of the eigenvalues can be investigated at the endemic state through numerical simulations.

4.   Numerical stability analysis and simulations

In this section, we investigate different types of behaviors of the system and analyze how various system parameters affect stability of the steady states and the dynamics of the model. We follow the work of Holt et al. (1997) ^[33] to set baseline values of parameters. The death rate $c$ of vectors is taken the same value of $c = 0.12$, the growth rate of plants is chosen as $r = 0.1$ which is in the range of $0.025-0.2$, and the carrying capacity of plants $K = 1$ is chosen to be at the higher value of the range $0.01-1$ explored in ^[33]. Vector abundance $a$ was investigated in ^[33] in the range $0-2500$, with the baseline value of $a = 500$, and we chose this to be $a = 300$, while plant death/roguing rate $m$ is chosen to be $m = 0.1$ within the range of $0–0.033$ explored in ^[33]. Finally, the transmission rates $\lambda$ and $\beta$ from vectors to plants, and from plants to vectors, were varied in the range $0–0.03$, respectively, $0-0.03$, against the range $0.002-0.032$ studied in ^[33].

Figure 1 shows how the values of the basic reproduction number $R_0$ varies with the time delay $\tau$, the plant infection rate $\lambda$, and vector infection rate $\beta$ in accordance with expression (3.9). This figure demonstrates a parameter region with $R_0 < 1$ where the disease is eradicated, as signified by the stable disease-free equilibrium $E_2$, and a region with $R_0 > 1$, where the disease is present at some constant level, i.e., the endemic steady state $E^*$ is feasible, and, as already mentioned in the previous section, for parameter values given in Table 1, it is also unique. At $R_0 = 1$, a transcritical bifurcation of the disease-free steady state occurs, as described in Theorem 4, at which point this steady state loses its stability, and the disease becomes endemic, i.e. the co-existence steady state becomes biologically feasible. One observes that for higher values of the time delay, a higher infection rate is required for the disease to persist, as indicated by the feasibility of the endemic steady state. Also, for a fixed duration of time delay, it is not so much an individual rate of disease transmission, but rather the product of transmission rates $\lambda$ from vector to plants, and $\beta$ from plants to vectors that determines whether the disease-free steady state is stable or the disease is endemic, as transpires from the expression for the basic reproduction number (3.9).

Figure 2 illustrates how the steady state values of the infected plant density and infective vectors change $R_0$ for $\tau_i = 0, \; i = 1, 2$. This figure shows that the co-existence steady state $E^*$ only becomes biologically feasible starting with some minimum value of $\lambda$ (transmission rate from vectors to plants) that corresponds to $R_0 = 1$ as discussed in Theorem 4, and for values of $R_0$ just above this critical value, the steady state $E^*$ is stable. Interestingly, starting from some threshold, further increasing the disease transmission rate $\lambda$ actually results in the reduction of the steady state value of the numbers of infected plants and vectors. Moreover, for sufficiently high values of $\lambda > \lambda^* = 0.01489$, the endemic steady state loses its stability via Hopf bifurcation, illustrating the results of Theorem 5.

To demonstrate the system's dynamics in different parameter regimes, we show in Figure 3 the behavior of the model with $\tau_i = 0, \; i = 1, 2$ and other parameters chosen in such a way that the endemic equilibrium $E^*$ is biologically feasible. When $\lambda < \lambda^* = 0.01489$, the endemic steady state $E^*$ is stable, and for $\lambda > \lambda^*$, the Hopf bifurcation has taken place, and, as a result, the system exhibits sustained periodic oscillations shown in Figure 4, in accordance with Theorem 3.4. Once can observe that the higher transmission rate not only causes the instability of this equilibrium, but also results in higher amplitude of oscillations around this steady state.

We compute the periods for different values of $\lambda$ using MatCont (See Figure 5). At the point where a limit cycle emerges, the period is 132.45 days, which happens at $\lambda^* = 0.01489$. Period increases with the value of $\lambda$, and at $\lambda = 0.025$, the period is 175 days, i.e., almost six months.

Bifurcation diagram of endemic equilibrium $E^*$ is shown in Figure 6, taking $\tau_1$ as the main parameter. We observed that $E^*$ becomes stable for large value of the delay parameter $\tau_1$. Figure 7 illustrates stability of the endemic steady state $E^*$ depending on different parameters, indicating that $E^*$ is unstable for small $\tau_1$ (which corresponds to high $R_0$) and stable for larger $\tau_1$. One can also note that for very small rate $\lambda$ of disease transmission from vectors to plants, the co-existence steady state $E^*$ is stable for all values of $\tau_1\geq 0$, for which it is feasible, as discussed in Proposition 1. We also observe that the product of transmission rate $\lambda$ from vectors to plants and $\beta$ from plants to vectors has to exceed some threshold value for the endemic steady state to be feasible, and for smaller values of this product the endemic state is stable but then also loses its stability via a Hop bifurcation, giving rise to stable periodic oscillations. Figure 7(d) shows the region of stability of $E^*$ in $\tau_1-\tau_2$ parameter plane. Large delays seems to imply that the coexistence equilibrium ${E}^*$ is stable, in spite of the high infection rate. When the delay crosses the threshold value, the endemic equilibrium becomes infeasible and the disease-free equilibrium becomes stable.

While we have focused our attention on the effects of various parameters and the time delay on the stability of the endemic steady state, it is also worth mentioning that the delay does not have such a pronounced effect on the stability of the disease-free equilibrium $E_{2}$ beyond a single stability change via the transcritical bifurcation at $R_0 = 1$. According to expression (3.9), increasing time delays $\tau_1$ and $\tau_2$ reduces the basic reproduction number, thus making the disease-free equilibrium stable, and making the solutions approach it more rapidly.

5.   Conclusions

In this research, a time-delayed model for the dynamics of the transmission of mosaic disease. Analytical and numerical analysis have provided conditions for feasibility and stability of various steady states of the model depending on parameters and the time delay. Increasing the rate of disease transmission from vectors to plants results in the destabilization of the disease-free equilibrium through a transcritical bifurcation, and the emergence of a stable endemic equilibrium. Interestingly, further increase of this transmission rate can lead to a destabilization of the endemic steady state, and the emergence of stable periodic oscillations, whose amplitude is itself growing with the transmission rate.

The time delay can play a dual role: As is often found in time-delayed models, it can by itself lead to a destabilization of the endemic state, but interestingly, it can also provide a mechanism for suppression of oscillations and recovery of stability for the endemic steady state that was otherwise unstable in the absence of the time delays. Numerical solution of the transcendental characteristic equation has provided further insights into how stability of the endemic steady state depends on parameters, and it has shown that the time delay plays a destabilizing role for smaller disease transmission rates from vectors to plants, and a largely stabilizing role for higher values of this transmission rate. We have also seen that increasing either of the transmission rates, i.e., from vectors to plants, or from plants to vectors, can result in the loss of stability by the endemic steady state.

Another important aspect of our model is elucidating the effects of the latency period on the type of dynamics that occurs between plants and their vectors, whereas long latency favors maintenance of a steady level of disease, shorter periods of latency result in the sustained oscillations in the populations of infected plants and vectors. Practical importance of this observation is two-fold. Since the latency period itself cannot be changed, and the fact that the critical value of the time delay in the model depends on several parameters, it is possible to switch the dynamics in a specific agricultural setting between a steady state and oscillations by controlling those other parameters, such as plant vector abundance per plant, or disease transmission rate.

In ^[16], constant growth was considered for the vector population. Impact of different parameters (such as delay, rate of infection etc.) on non-infective vectors was not observed in ^[16]. In this research, we have taken Leslie–Gower type growth for the vector population. We have considered a specific vector (whitefly) and the disease (mosaic disease) transmitted by the vector to the plants like Jatropha curcas, Cassava, etc. Thus we have used the parameter values that are used for Cassava plants by Holt et al. 1997, ^[33] and for Jatropha plant by Venturino et al. 2016, ^[31]. Moreover, we have studied the Forward bifurcation analytically and numerically (Figure 2). Also, we have shown the dependence of the bifurcation parameter on the period of the limit cycle (Figure 5).

The results of this research provide several practically important insights for control and mitigation of the effects of mosaic disease. The basic reproduction number that can be computed using values for fundamental parameters describing the specific farming situation, provides an effective cumulative characteristic of the disease severity that can be used to compare different outbreaks. While it is not possible to change the underlying biological properties of plants and vectors, it is possible to modify their interactions by administering various sprays and insecticides. These would act to reduce the effective infection rate of disease transmission from vectors to plants, and could ultimately lead to eradication of the disease, as represented by stabilization of the disease-free steady state. Furthermore, when the system exhibits sustained oscillations in the level of infection, having a detailed information about how the period and amplitude of those oscillations depend on other parameters, it is possible to optimize administering various interventions, such as fertilizers or insecticides, by exploiting the troughs of the oscillations.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

EV is a member of the INDAM research group GNCS. He has partially been supported by the project "Metodi numerici nelle scienze applicate" of the Dipartimento di Matematica "Giuseppe Peano" of the Università di Torino.

Conflict of interest

Ezio Venturino is an editorial boardmember for AIMS Mathematics and was not involved in the editorial review or the decision to publish this article.

The authors declare that there is no conflict of interest.

Appendix A

Basic reproduction number

We follow the method established in the paper by Hefferman et al. ^[45] for calculating $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_2(K, 0, aK, 0)$ is the disease-free equilibrium and indices $i$, $j$ correspond, respectively, to $I$ and $V$.

The basic reproduction number $R_0$ is the dominant eigenvalue of the matrix $G = FV^{-1}$, i.e., Eq (5.1).

and it is the dominant eigenvalue of the matrix ${\bf G} = {\bf F}{\bf V}^{-1}$.