Backward bifurcation of a plant virus dynamics model with nonlinear continuous and impulsive control

1.   Introduction

Plant virus diseases are transmitted between plant individuals through vectors, including the tomato leaf roll virus and cassava mosaic virus. Plant viruses seriously affect crop yields, which can result in great economic losses and increased poverty, especially in developing countries ^[1]. According to ^[2], the annual wheat rust can cause yield loss of up to 80 percent, posing a great threat to global wheat production. Moreover, it is estimated that 12.5 million hectares were disturbed annually from 2003 to 2012 by plant disease, mostly Asia and Europe ^[3]. Therefore, the control of plant virus transmission remains a major economic and agricultural issue. The effective methods to control plant viruses include gene control, chemical control, and biological control. It is extremely effective to improve the antiviral ability of plants by changing genes, but it cannot be greatly used due to the high cost. Chemical control aims to control the vector by spraying insecticide. However, long-term use will not only pollute the environment, but also lead to drug resistance and increase the incidence of plant viruses. Biological control is widely used to control the spread of diseases by manually removing infected plants or introducing vector natural enemies, having significant effects and causing no pollution to the environment. Thus, biological control is increasingly used as an efficient method to control plant viruses.

Although biological control of plant viruses provides an effective method for controlling the spread of plant viruses, a successful control plan needs to understand the transmission mechanism of plant viruses in plants and vectors as well as the interaction between biological controllers and plants and vectors. Mathematical modeling has exerted a main and probably greater role in the epidemiology of plant viruses. Mathematical models enable us to comprehend the observed transmission routes and disease control of plant viruses, as well as grasp a deeper understanding of the mechanism of plant disease transmission dynamics. At present, there are numerous research achievements on biological control of plant viruses. As mentioned in ^[4], the authors consider the evolution of within-plant virus titer as the response to the application of various disease control methods. Moreover, they demonstrate that the process of new and improved disease control methods for viral diseases of vegetatively propagated staple food crops need to consider the evolutionary responses of the virus. According to ^[5], the authors concentrate on the contribution made by developing mathematical models with a scope of techniques that is extensive and their use for investigating plant virus disease epidemics. Their focus is on the extent to which models can help answer biological questions and raised questions in association with the epidemiology and ecology of plant viruses and the caused diseases. Based on ^[6,7], the authors applied optimization theory to investigate the plant vector virus model with continuous replanting and roguing, aiming to maximize the harvest of healthy plants and examine the best strategy to fight plant virus disease. In ^[8], the authors concentrate on a plant disease model with a latent period and nonautonomous phenomenon. Their focus is to explore the long-run behavior of the epidemic model.

In reality, considering that people control plant viruses at discrete time points, the establishment of impulse differential equations to reveal the propagation dynamics and control strategies of plant viruses has been favored by numerous scholars in recent years. In ^[9], the authors proposed and analyzed a plant virus disease model under periodic environment and pulse roguing. Based on their research, with the high infection rate, it can probably not be possible to eliminate the disease by easily roguing the infectious plant, and elevating the replanting rate is not conducive to the control of disease. In ^[10], to detect and design suitable plant disease control strategies, the authors proposed and investigated the dynamics of plant disease models by adopting continuous and impulsive cultural control strategies. Under a certain parameter space, it can be demonstrated that a nontrivial periodic solution occurs through a supercritical bifurcation. In ^[11], the authors formulated an epidemiological model for mosaic diseases considering plant and vector populations. They discovered that roguing is the most cost effective and beneficial management for mosaic disease eradication of plants if used at an appropriate rate and interval. In ^[12,13], with the purpose of eradicating plant diseases or keeping the number of infected plants below the economic threshold, the authors proposed and analyzed the plant disease models, including nonlinear impulsive functions and cultural control strategies. In ^[14], based on the purpose of minimizing losses and maximizing returns, the authors introduce a discontinuous plant disease model, which incorporates a threshold policy control. Their results suggest that we can adopt the proper replanting and roguing rates for designing the threshold policy, and thus the number of infected plants is within an acceptable level. Recently, pulse control methods can be applied in almost every field of applied science. The theoretical research and application analysis of pulse control can be found in many studies ^{[15,16,17,18]}.

However, most of the above studies consider only (1) Continuous removal of infected plants or replanting of healthy plants, without considering removal of infected plants at discrete time points. (2) Pulse removes the infected plants without destroying the vector on the plants. Therefore, considering the shortcomings of the above literature, based on the above two points, we aim to develop a plant virus disease transmission model with nonlinear continuous and pulse control, analyze its dynamic behavior, and discuss the impact on plant virus disease control when removing infected plants and vectors on plants.

2.   A model for continuous control

2.1.   System description

According to the model proposed in ^[19], we develop a vector-borne plant virus disease model with roguing control. The basic idea is that removing the infected plants will also eliminate the vectors on the plants. For plant populations, we classified plants into healthy plants $X(t),$ infected plants $Y(t),$ and the total number of plant populations at time $t$ with $N(t) = X(t)+Y(t).$ The recruitment rate of plant population is constant $r.$ The newly recruited plants enter healthy plant populations $X(t).$ The plant harvesting or death rate is denoted as $g.$ Healthy plants are infected through contact with infection vectors, and the infection rate is $\frac{k_{1}V(t)}{N(t)},$ where $k_{1}$ is the probability of transmission from infected vectors to healthy plant. The expression of infection rate $\frac{k_{1}V(t)}{N(t)}$ is acquired as follows. The probability that a vector selects a particular plant is assumed to be $\frac{1}{N(t)},$ which is the probability of each plant selected by total $N(t).$ Therefore, a plant receives in average $\frac{M(t)}{N(t)}$ contacts per unit of time. Later, the infection rate per healthy plant is offered by $\frac{M(t)}{N(t)}\frac{k_{1}V(t)}{M(t)}.$ Considering using removal control on infected plants, the removal rate of diseased plants is $\alpha.$

Vector populations are divided into non-infective $U(t)$ and infective vectors $V(t),$ and the total number of vector populations at time $t$ denote $M(t) = U(t)+V(t).$ The recruitment rate of vector population is constant $\mu.$ The newly recruited vectors enter non-infective vector populations $U(t).$ The vector death rate is $c.$ Uninfected vectors become infected through contact with infected plants, and the infection rate is $\frac{k_{2}Y(t)}{N(t)},$ where $k_{2}$ is the probability of transmission from infected plants to non-infective vectors.

The elimination of vectors occurs when the infected plants are removed with elimination rate of $\frac{\rho\alpha Y(t)}{N(t)}$ which is explained as follows. The average number of vector on each plant is $\frac{M(t)}{N(t)},$ and parameter $\rho$ represents the elimination ratio of vectors on infected plants while removing each infected plant. Since the number of infected plants removed is $\alpha Y(t),$ with the removal of the infected plants, the total number of eliminated vector is $\frac{\rho\alpha Y(t)M(t)}{N(t)}.$ Thus, the elimination rate of vectors is $\frac{\rho\alpha Y(t)M(t)}{N(t)}\frac{1}{M(t)}.$ Therefore, our basic plant disease model reads as

We suppose that all parameter values are strictly positive.

Next, we indicate the boundedness of system (2.1), and it can follow from system (2.1) that:

which suggests:

Similarly, we have

and thus, we obtain:

Therefore, the total plant population $N(t)$ and vector population $M(t)$ remain uniformly bounded. Moreover, the region

is positively-invariant, indicating that this study concentrates on the dynamics of system (2.1) on the set $\Omega$ in the following.

2.2.   Disease-free equilibrium and stability analysis

With no disease in both of the populations(i.e., $Y(t) = 0,$ $V(t) = 0$), the model (2.1) has one disease-free equilibrium $E_{0} = (N_{0}, 0, M_{0}, 0) = (\frac{r}{g}, 0, \frac{\mu}{c}, 0).$ By ^[20], we define the basic reproductive number

We note that

where

refers to the expected number of plants that one vector infects via its infectious life time. $R_{01}$ can be provided by the product of the infection rate of infectious vectors $k_{1}$ and the average duration in the infectious stage $\frac{1}{c}.$

Similarly, we have

which is the expected number of vectors that one plant infects via its infectious life time. $R_{02}$ is given by the product of the infection rate of infectious plants $k_{2}\frac{M_{0}}{N_{0}}$ and the average duration in the infectious stage $\frac{1}{g+\alpha}.$

The basic reproduction number equals to the geometric mean of $R_{01}$ and $R_{02}$ since infection from plant to plant goes via one generation of vectors. In addition, the local asymptotic stability result of equilibrium $E_{0}$ is provided:

Theorem 2.1. When $R_{0} < 1,$ the disease-free equilibrium $E_{0} = (\frac{r}{g}, 0, \frac{\mu}{c}, 0)$ remains locally asymptotically stable in $\Omega.$

Proof. The Jacobian matrix of system (2.1) at disease-free equilibrium $E_{0}$ can be provided by:

The characteristic polynomial of $J_{E_{0}}$ can be expressed by:

where $\varphi(\lambda) = \lambda^{2}+a_{1}\lambda+a_{2},$ with $a_{1} = c(g+\alpha)$ and $a_{2} = c(g+\alpha)(1-R_{0}^{2}).$

The roots of $P(\lambda)$ are $\lambda_{1} = -g,$ $\lambda_{2} = -c$ and the others roots are the roots of $\varphi(\lambda).$ Because $R_{0} < 1,$ all coefficients of $\varphi(\lambda)$ are always positive. Now, we just have to verify that the Routh-Hurwitz criterion holds for polynomial $\varphi(\lambda)$. Therefore, setting $H_{1} = a_{1},$ $H_{2} = \left| \begin{array}{ccc} a_{1} & 1 \\ 0 & a_{2} \\ \end{array} \right| = a_{1}a_{2}$. We have $H_{1} > 0$ and $H_{2} > 0$ if $R_{0} < 1.$ Thus, by the Routh-Hurwitz criterion the trivial equilibrium $E_{0}$ is locally asymptotically stable when $R_{0} < 1.$

In general, the biological implication of Theorem 2.1 is that with the basic reproduction number $R_{0}$ being less than $1$, a small number of infected vectors are introduced into the plant population, which will not cause disease outbreak, and the disease disappears in time. Nevertheless, in the following subsection, we present that the disease can persist even with $R_{0} < 1.$

2.3.   The existence of endemic equilibria and backward bifurcation

With the purpose of determining whether there is an endemic equilibrium, we need to find the solution of the algebraic system of equations acquired through equating the right sides of system (2.1) to zero.

For the sake of easier readability, we express the quantities as follows,

By addressing the equations in the system (2.4) in terms of $\lambda_{p}^{*}$ and $\lambda_{v}^{*}$, we can obtain:

and

By substituting (2.5) and (2.6) into the expression of $\lambda_{p}^{*}$ and $\lambda_{v}^{*}$, we have

where

We aimed to investigate the existence of endemic equilibria in the following cases:

(1) There is a unique endemic equilibrium if

($C < 0$) or ($B < 0$ and $C = 0$) or ($B < 0$ and $C > 0$ and $B^{2}-4AC = 0$);

(2) There are two endemic equilibria if

$B < 0$ and $C > 0$ and $B^{2}-4AC > 0;$

(3) There are no endemic equilibria otherwise.

Hence, we present the result in the theorem below.

Theorem 2.2. For system (2.1),

(1) If $R_{0} > 1,$ there is a unique endemic equilibrium $E^{*}.$

(2) If $R_{0} = 1$ and $B < 0$ there is a unique endemic equilibrium $E^{*}.$ Otherwise, there exists no endemic equilibrium.

(3) If $R_{0} < 1,$ and

(a) $B < 0$ and $B^{2}-4AC > 0,$ there are two equilibria $E_{1}$ and $E_{2}.$

(b) $B < 0$ and $B^{2}-4AC = 0,$ there exists a unique endemic equilibrium.

(c) There are no endemic equilibria otherwise.

Obviously, case $(3)$ (item(a)) of Theorem 2.2 indicates the possibility of backward bifurcation in the model (2.1). Next, by applying the method of literature ^[21] to model (2.1), we present a rigorous proof that model (2.1) experiences a backward bifurcation.

Theorem 2.3. If $\frac{c\alpha}{g} > r+k_{2}+2\rho\alpha,$ the direction of the bifurcation of system (2.1) at $R_{0} = 1$ is backward.

Proof. Let $k_{1}$ be the bifurcation parameter. To apply the method in reference^[21], we introduce the notation $x_{1} = X,$ $x_{2} = Y,$ $x_{3} = U,$ $x_{4} = V,$ the system (2.1) becomes

with $R_{0} = 1$ corresponding to $k_{1} = k_{1}^{*} = \frac{c^{2}r(g+\alpha)}{k_{2}\mu g}.$ The disease-free equilibrium is $E_{0} = (\frac{r}{g}, 0, \frac{\mu}{c}, 0).$ The linearization matrix of system (2.8) around the disease-free equilibrium when $k_{1} = k_{1}^{*}$ is

where $f = (f_{1}, f_{2}, f_{3}, f_{4})'.$ Clearly, $0$ indicates a simple eigenvalue of $D_{E_{0}}f.$ A right eigenvector related to $0$ eigenvalue is $\omega = [-\frac{cr(g+\alpha)}{k_{2}\mu g^{2}}, \frac{cr}{k_{2}\mu g}, -\frac{k_{2}+\rho\alpha}{ck_{2}}, \frac{1}{c}]',$ and the left eigenvector $\upsilon$ meeting $\omega\cdot\upsilon = 1$ is $\upsilon = [0, \frac{k_{2}\mu g}{r(c+g+\alpha)}, 0, \frac{c(g+\alpha)}{(c+g+\alpha)}].$ Algebraic calculations demonstrate that

The remaining of the second derivatives appear in the formula for $a$ and $b$ are all zero. The $a$ and $b$ presented in Theorem 4.1 of ^[21] are

Hence,

By ^[21], if $\frac{c\alpha}{g} > r+k_{2}+2\rho\alpha,$ then $a > 0,$ and the direction of the bifurcation of system (2.1) at $R_{0} = 1$ is backward.

The backward bifurcation phenomenon indicates that with the basic reproductive number $R_{0} < 1,$ the disease-free equilibrium is locally stable. Therefore, it is insufficient to only reduce the basic reproductive number $R_{0} < 1$ in order to eliminate the disease. Figure 1 depicts the associated bifurcation diagram. This clearly demonstrates the coexistence of two locally-asymptotically stable equilibria with $R_{0} < 1,$ proving that the model (2.1) reveals the phenomenon of backward bifurcation. Therefore, to eradicate plant diseases, it is of necessity to increase $\alpha,$ and thus $R_{0}$ is less than $R_{0}(\alpha_{0}),$ where $R_{0}(\alpha_{0})$ is the threshold for the appearance of two endemic equilibria.

Figure 2 presents the occurrence of the backward bifurcation. In this study, $R_{0} < 1,$ while according to the initial condition, the solution of the model (2.1) can approach either the endemic equilibrium point or the disease-free equilibrium point.

3.   A model for impulsive control

The current section extends model (2.1) by replacing the continuous infected plants with a periodic pulse roguing control strategy, which is shown to be more realistic. Therefore, our impulsive control model of plant virus diseases reads as

where $T$ indicates a fixed positive constant and suggests the periodic of the impulsive effect, where $n\in\mathcal{N}$ represents the positive integer set. The parameter $\alpha$ suggests the proportion of the infected plants, which is rogued at each pulse perturbation. Considering the inevitable impact of control measures on healthy plants, we use $p$ representing the proportion of the healthy plants, which is rogued at each pulse perturbation. Here, we assume that $p < \alpha.$ Similar to the continuous model, the ratio of removing uninfected and infected vectors while removing plants is $\frac{\rho\alpha Y(nT)}{N(nT)}.$ As there are only uninfected vectors on healthy plants, only uninfected vectors can be eliminated while removing healthy plants. The elimination ratio is $\frac{\rho p X(nT)}{N(nT)},$ where $\rho$ represents the vectors elimination rate of vectors, which probably means that some of the uninfected vectors and infected vectors will be occasionally eliminated due to the rogued plants.

3.1.   The existence and stability of the disease-free periodic solution

When $Y(t) = 0,$ $V(t) = 0,$ then model (3.1) is the subsystem below:

which presents the dynamics of the system in the absence of the infected plants and vectors. For system (3.2), we can address it in any impulsive interval $(nT, (n+1)T]$ and obtain:

Denote $X_{n} = X(nT^{+})$ and $U_{n} = U(nT^{+}),$ then:

There is a steady state $(X^{*}, U^{*}),$ which suggests that system (3.2) has a positive periodic solution $(X^{*}(t), U^{*}(t)),$ where $X^{*}(t) = (X^{*}-\frac{r}{g})\exp(-g(t-nT))+\frac{r}{g},$ $U^{*}(t) = (U^{*}-\frac{\mu}{c})\exp(-c(t-nT))+\frac{\mu}{c}$ with $X^{*} = \frac{r(1-p)(1-\exp(-gT))}{g(1-(1-p)\exp(-gT))}$ and $U^{*} = \frac{\mu(1-p\rho)(1-\exp(-cT))}{c(1-(1-p\rho)\exp(-cT))}.$ Hence, the following lemma can be obtained.

Lemma 3.1. System (3.2) has a positive periodic solution $(X^{*}(t), U^{*}(t)).$ For any solution of (3.2), this study obtains $X(t)\rightarrow X^{*}(t),$ $U(t)\rightarrow U^{*}(t)$ as $t\rightarrow \infty.$

Proof. From (3.3) and (3.4), we obtain the solution of (3.2) as

which indicates that $X(t)\rightarrow X^{*}(t)$ and $U(t)\rightarrow U^{*}(t)$ as $n\rightarrow \infty$ and $t\in(nT, (n+1)T].$

Therefore, system (3.1) has a disease-free periodic solution $(X^{*}(t), 0, U^{*}(t), 0)$ and its stability has been solved.

To determine the stability of disease-free periodic solution $(X^{*}(t), 0, U^{*}(t), 0)$ of system (3.1), we first calculate the basic reproduction number for the impulsive model (3.1) with the application of the next infection operator for the piecewise continuous periodic system being proposed in ^[22]. $A(t)$ can be denoted as a $n\times n$ matrix, $\Phi_{A(.)}(t)$ is suggested as the fundamental solution matrix of the linear ordinary differential system $\chi^{'} = A(t)\chi,$ and $r(\Phi_{A(.)}(t))$ is indicated as the spectral radius of $\Phi_{A(.)}(t).$ Define $X(t) = x(t)+X^{*}(t),$ $Y(t) = y(t),$ $U(t) = u(t)+U^{*}(t),$ $V(t) = v(t),$ $\chi(t) = (x(t), u(t), y(t), v(t)).$ The corresponding linear system (3.1) reads as

where

with

and

Let $\Phi_{G}(t) = \Phi_{ij},$ $1\leq i, j\leq 2,$ the fundamental solution matrix of system (3.6). Then, we have $\Phi^{'}_{G}(t) = G(t)\Phi_{G}(t)$ with the initial value $\Phi_{G}(0) = I_{4}.$ Based on the further computation, it is suggested that:

then we have

It is obvious that $r(H_{1}\exp(G_{1}T)) < 1,$ by ^[22], the basic reproduction number for system (3.6) is provided as follow

Using Floquet theory, if $R_{1} < 1,$ we obtain the following result.

Theorem 3.2. The disease-free periodic solution $(X^{*}(t), 0, U^{*}(t), 0)$ of model (3.1) is locally asymptotically stable when $R_{1} < 1.$

3.2.   Persistence of the disease

According to the persistence of the system, it can be found that plants (susceptible and infected) and vectors (susceptible and infected) can coexist, indicating that some conditions are satisfied, and the disease will not disappear. If we aim to eliminate the disease, the persistence indicates that control strategies are inefficient. Moreover, the permanent condition acquired from the exploration of the system can offer scientific support for determining the key factors causing failure and the effectiveness of control strategies, which later benefit us in establishing a good control program.

Theorem 3.3. When $R_{2} = r(H^{'}_{3}\Phi_{F-V_{r}}(T)) > 1,$ the disease is uniform persistence, i.e., there exists $\eta > 0$ such that $\lim\limits_{t\rightarrow \infty}\inf{Y(t)}\geq \eta > 0, \lim\limits_{t\rightarrow \infty}\inf{V(t)}\geq \eta > 0.$ Here,

Proof. At first, we claim that there exists $\eta > 0,$ and thus:

Next, there is a $t_{1} > 0$ such that $Y(t) < \eta, V(t) < \eta$ for all $t\geq t_{1}.$ According to the first and third equations in system (3.1), the following can be obtained:

and focus on the auxiliary system

Since system (3.16) refers to a quasimonotone increasing system, based on the comparison theorem ^[23], the following can be obtained:

Using the solution of system (3.2), we acquire that model (3.16) admits a globally asymptotically stable positive periodic solution $z^{*}(t) = (z^{*}_{1}(t), z^{*}_{2}(t)),$ and $\lim\limits_{\eta\rightarrow 0}z^{*}(t) = (X^{*}(t), U^{*}(t)).$ Then, there exists $\eta_{1}$ small enough and $\varepsilon_{1} > 0$ such that $z^{*}_{1}(t)\geq X^{*}(t)-\varepsilon_{1}$ and $z^{*}_{2}(t)\geq U^{*}(t)-\varepsilon_{1}$for $\eta < \eta_{1}.$ Therefore, by the comparison theorem, there is $t_{2}\geq t_{1}$ and $\varepsilon_{2} > 0,$ and thus

By substituting the above inequalities into the second and fourth equations of system (3.1), we acquire:

Supplementing the equations of plant population and vector population respectively provides:

Denoting $(N_{1}(t), M_{1}(t))$ is the solution of system (3.20) and $(X_{1}(t), U_{1}(t))$ as that of system (3.2). Comparing the equations in system (3.20) and (3.2) yields $N_{1}(t)\leq X_{1}(t), M_{1}(t)\leq U_{1}(t).$ Since $X_{1}(t)\leq X^{*}(t),$ $U_{1}(t)\leq U^{*}(t),$ then we have $N_{1}(t)\leq X^{*}(t),$ $M_{1}(t)\leq U^{*}(t),$ and system (3.19) can be modified as

Further, we concentrate on an auxiliary system

It can be re-expressed as

The solution of system (3.23) can be shown as $(z_{3}(t), z_{4}(t))^{T} = \Phi_{F-V_{r}}(t-nT)(z_{3}(nT^{+}), z_{4}(nT^{+}))^{T}.$ Then, $(z_{3}((n+1)T^{+}), z_{4}((n+1)T^{+}))^{T} = H^{'}_{3}\Phi_{F-V_{r}}(T)(z_{3}(nT^{+}), z_{4}(nT^{+}))^{T}.$ While $R_{2} > 1,$ $z_{3}(t)\rightarrow \infty$ and $z_{4}(t)\rightarrow \infty$ as $t\rightarrow \infty.$ Consequently, $\lim\limits_{t\rightarrow \infty}Y(t) = \infty$ and $\lim\limits_{t\rightarrow \infty}V(t) = \infty.$ Moreover, this is a contradiction with the $Y(t)$ and $V(t).$ Thus, the claim is demonstrated, i.e.,

We obtain the two possibilities.

(ⅰ) $Y(t)\geq\eta$ and $V(t)\geq\eta$ for all large $t$;

(ⅱ) $Y(t)$ and $V(t)$ oscillations concerning $\eta$ for all large $t$.

When case (ⅰ) is true, our proof is completed. Next, we will consider case (ⅱ). Because $\lim\limits_{t\rightarrow \infty}\sup{Y(t)}\geq \eta,$ $\lim\limits_{t\rightarrow \infty}\sup{V(t)}\geq \eta,$ we can select a $t_{1}\in (n_{1}T, (n_{1}+1)T]$ such that $Y(t_{1})\geq \eta$ and $V(t_{1})\geq \eta.$ If the oscillation exists, there must be another $t_{2}\in (n_{2}T, (n_{2}+1)T]$ such that $Y(t_{2})\geq \eta$ and $V(t_{2})\geq \eta,$ in which $n_{2}-n_{1}\geq 0$ is finite. Next, this study will focus on the solution of system (3.1) in the interval $[t_{1}, t_{2}]:$

Which suggests that

Moreover, it follows from

Then we obtain:

Let $\delta_{1} = \min\{\eta(1-\alpha)^{n_{2}-n_{1}}\exp(-g(n_{2}-n_{1})T), \eta(1-\rho\alpha)^{n_{2}-n_{1}}\exp(-c(n_{2}-n_{1})T)\},$ then $\delta_{1} > 0$ cannot be infinitely small due to the $n_{2}-n_{1}\geq 0$ is finite, leading to $Y(t)\geq \delta_{1} > 0$ and $V(t)\geq \delta_{1} > 0.$

When $t > t_{2},$ we take the same steps and get another non-infinitesimal positive $\delta_{2}.$ Therefore, the sequence $\{\delta_{i}\}(i = 1, 2, 3, ...),$ in which $\delta_{i} = \min\{\eta(1-\alpha)^{n_{i+1}-n_{i}}\exp(-g(n_{i+1}-n_{i})T), \eta(1-\rho\alpha)^{n_{i+1}-n_{i}}\exp(-c(n_{i+1}-n_{i})T) \}$ is non-infinitesimal since $n_{i+1}-n_{i}\geq 0$ is finite. The solution of system (3.1) $Y(t)\geq \delta_{k} > 0$ and $V(t)\geq \delta_{k} > 0$ is true in the time interval $[t_{k}, t_{k+1}], t_{k}\in(n_{k}T, (n_{k}+1)T], t_{k+1}\in(n_{k+1}T, (n_{k+1}+1)T].$ Let $\eta = \min\{\delta_{1}, \delta_{2}, ... \},$ then $\eta\in \{\delta_{1}, \delta_{2}, ...\}$ and it is shown that $Y(t)\geq \eta$ and $V(t)\geq\eta$ for all $t\geq t_{1}.$ The proof is complete.

Remark: It becomes easy to observe that $R_{2}$ is less than $R_{1}.$ Therefore, when $R_{1} > R_{2} > 1,$ the disease persists. Nevertheless, due to the complexity of the model, we obtain only the stronger condition for the persistence in the (3.3), that is $R_{2} > 1,$ than $R_{1} > 1.$

3.3.   Bifurcation analysis of nontrivial endemic periodic solution

In the current subsection, to study the influence of the removal rate of infected plants on the model, we choose $\alpha$ as the bifurcation parameter to investigate the possible behaviors of nontrivial periodic solution. Now, we proceed to explore bifurcation based on the bifurcation theorem of ^[24] and ^[25]. Let $x_{1}(t) = X(t), x_{2}(t) = U(t), x_{3}(t) = Y(t), x_{4}(t) = V(t).$ Then, we employ the following notations in model (3.1)

First, the following notations are introduced: The solution vector $\xi(t): = (x_{1}(t), x_{2}(t), x_{3}(t), x_{4}(t)),$ the mapping $F(\xi(t)) = (F_{1}(\xi(t)), F_{2}(\xi(t)), F_{3}(\xi(t)), F_{4}(\xi(t))):$ $R^{4}\rightarrow R^{4}$ by the right hand side of the first four equations of system (3.29), and the mapping with impulse to be

Let $\Phi(t, \xi_{0})$ the solution of the system which is consisted of the first four equations of system (3.29), where $\xi_{0} = \xi(0).$ Next, $\xi(T) = \Phi(T, \xi_{0}) = :\Phi(\xi_{0})$ and $\xi(T^{+}) = \Theta(\alpha, \Phi(\xi_{0})).$ The operator $\Psi$ is defined by:

and denote by $D_{\xi}\Psi$ the derivative of $\Psi$ with respect to $\xi.$ Later, $\xi$ indicates a periodic solution of period $T$ for system (3.29) if and only if $\Psi(\alpha, \xi_{0}) = \xi_{0}.$ Therefore, to obtain the nontrivial periodic solution of system (3.29), we are required to demonstrate the presence of nontrivial fixed point of $\Psi.$

All parameters are fixed except the infective plants removal rate $\alpha.$ Denote that $\alpha_{0}$ is the critical remove rate, corresponding to $R_{1} = r(H_{3}\Phi_{F-V_{R}}(T)) = 1.$ Let $\overline{\xi} = (x^{*}_{1}, x^{*}_{2}, 0, 0)$ be the disease-free periodic solution of system (3.29).

Denote $\alpha = \alpha_{0}+\overline{\alpha},$ $\xi = \xi_{0}+\overline{\xi},$ with $\xi_{0}$ is the starting point for the disease-free periodic solution with the removal rate $\alpha_{0},$ and let

Next, the fixed point problem can be expressed:

The derivative of $N(\overline{\alpha}, \overline{\xi})$ with respect to $\xi$ provides:

where $E_{4}$ refers to the identity matrix. Based on system (3.29), the equality can be obtained:

with $D_{\xi}\Phi(0, \xi_{0}) = E_{4}$ and $\Phi(t, \xi_{0}) = (\Phi_{1}(t, \xi_{0}), \Phi_{2}(t, \xi_{0}), 0, 0),$ then Eq (3.33) takes the form

The following can be obtained:

where $\textbf{O} = (0, 0, 0, 0).$ The essential condition for the bifurcation of nontrivial zeros of function $N$ refers to that the determinant of the Jacobian matrix $D_{\xi}N(0, \textbf{O})$ is shown to be equal to zeros, for example,

It is not difficult to observe from (3.35) that $\det(E_{2}-H_{1}e^{G_{1}T})\neq 0.$ Then $\det(D_{\xi}N(0, \textbf{O})) = 0$ can reduce to $\det(E_{2}-H_{3}\Phi_{F-V_{r}}(T)) = 0.$ If $R_{1} = r(H_{3}\Phi_{F-V_{r}}(T)) = 1,$ then $\det(E_{2}-H_{3}\Phi_{F-V_{r}}(T)) = 0.$ Now we investigate the sufficient conditions for the existence of bifurcation nontrivial period solutions. With the consideration of convenience, we denote the elements in matrix $D_{\xi}N(0, \textbf{O})$ as

from the calculation in Appendix A, we can obtain the expression of each element in the above matrix as follows:

Then $\det(E_{2}-H_{3}\Phi_{F-V_{r}}(T)) = 0$ suggests that there is a constant $k\neq 0$ such that $c_{0} = ka_{0}$ and $d_{0} = kb_{0}.$ Moreover, since $\det(D_{\xi}N(0, \textbf{O})) = 0,$ we cannot employ the Implicit Function Theorem for giving variable $\xi$ as a function of $\alpha.$ By ^[24] and ^[25], we perform a Lyapunov-Schmide reduction to obtain a system of equations, where the Implicit Function Theorem can be applied. It is easy to observe that dimKer$(D_{\xi}N(0, \textbf{O})) = 1$ and a basis in Ker$(D_{\xi}N(0, \textbf{O}))$ is

Moreover, $Y_{2} = (1, 0, 0, 0),$ $Y_{3} = (0, 1, 0, 0)$ and $Y_{4} = (0, 0, 1, 0)$ compose a basis in Im$(D_{\xi}N(0, \textbf{O})).$ Later, for $\overline{\xi}\in R^{4},$ based on the decomposition $R^{4} = Ker(D_{\xi}N(0, \textbf{O}))\oplus Im(D_{\xi}N(0, \textbf{O})),$ the following can be obtained:

where $\alpha_{i}\in R(i = 1, 2, 3, 4)$ are special. Thus, (3.31) is equivalent to

From the first three equations of (3.39), the following can be obtained:

Thus, from the Implicit Function Theorem, one can deal with equations (3.39) as $i = 2, 3, 4$ regarding $(0, \textbf{O})$ with respect to $\alpha_{i},$ $i = 2, 3, 4$ as function of $\overline{\alpha}$ and $\alpha_{1},$ and find $\tilde{\alpha_{i}} = \tilde{\alpha_{i}}(\overline{\alpha}, \alpha_{1})$ and thus $\tilde{\alpha_{i}}(0, 0) = 0,$ $i = 2, 3, 4,$ and

$i = 1, 2, 3.$ Next, only $N(\bar{\alpha}, \bar{\xi}) = 0$ if and only if

It can be known that $f(\bar{\alpha}, \alpha_{1})$ vanishes at $(0, 0).$ Therefore, it becomes essential to compute higher order derivatives of $f(\bar{\alpha}, \alpha_{1})$ up to the order $i$ for which $D^{i}f(0, 0)\neq 0.$ Similar to the calculation in ^[24] and ^[25], the first partial derivatives of $f$ about $\bar{\alpha}$ and $\alpha_{1}$ satisfy

(See Appendix B). Therefore, we obtain $Df(0, 0) = (0, 0).$ Then, it becomes essential to calculate the second partial derivative of $f.$

Let $A = \frac{\partial^{2}f(\bar{\alpha}, \alpha_{1})}{\partial\bar{\alpha}^{2}},$ $B = \frac{\partial^{2}f(\bar{\alpha}, \alpha_{1})}{\partial\bar{\alpha}\alpha_{1}}$ and $C = \frac{\partial^{2}f(\bar{\alpha}, \alpha_{1})}{\partial\alpha^{2}_1{}}.$ From the calculation in Appendix C, it can be found that $A = 0,$

and

Therefore, we have

where

Furthermore,

Thus, for $B\neq 0,$ we can employ the Implicit Function Theorem, implying that $\bar{\alpha} = \varphi_{1}(\alpha_{1}),$ such that for all $\alpha_{1}$ near $0,$ $\tilde{f}(\varphi_{1}(\alpha_{1}), \alpha_{1}) = 0.$ Furthermore, for $C\neq 0,$ we can also be $\alpha_{1} = \varphi_{1}(\bar{\alpha}),$ such that for all $\bar{\alpha}$ near $0,$ $\tilde{f}(\varphi_{2}(\bar{\alpha}), \bar{\alpha}) = 0.$ Then, if $BC\neq 0,$ we have $\frac{\alpha_{1}}{\bar{\alpha}}\simeq \frac{-2B}{C}.$ There is a supercritical bifurcation to a nontrivial periodic solution near the fixed point $\xi_{0},$ if $BC < 0,$ and it is subcritical, if $BC > 0.$ We know that the threshold $R_{1}$ lowers with the increasing $\alpha,$ then a supercritical bifurcation in the $\bar{\alpha}-\alpha_{1}$ plane indicates a backward bifurcation in the model. Additionally, the subcritical bifurcation equated to a forward bifurcation. Finally, the following theorem below can be concluded.

Theorem 3.4. If $BC < 0,$ system (3.29) admits a backward bifurcation as $\alpha$ passes via the critical value $\alpha_{0}.$ Then, a forward bifurcation will occur $BC > 0.$

4.   Numerical simulations

First, the analysis result of model (2.1) suggests that the system has backward bifurcation with the basic reproduction number $R_{0} = 1.$ As shown in Figure 1, we fixed other parameters to let $\alpha$ change, and gave the backward bifurcation diagram of the system (2.1). It can be observed from the Figure 1 that when $R_{0}(\alpha_{0})\leq R_{0} < 1$, the system has a stable positive equilibrium point and an unstable positive equilibrium point. The disease can be eradicated only when $R_{0} < R_{0}(\alpha_{0}).$ Thus, when the infected plants are removed, the disease can be eradicated only if the removal rate $\alpha > \alpha_{0}.$ As presented in Figure 2, we choose the parameter such that $R_{0} = 0.9975 < 1$ changes the initial value of the system. With the initial value being $(7, 0.5,750, 30),$ the system has a disease free equilibrium point that remains globally asymptotically stable. When the initial value is $(70, 0.5,750, 30),$ the system possesses a stable positive solution. Therefore, it is necessary to consider the initial infection when preventing and controlling plant viruses.

Second, the disease control threshold condition $R_{1}$ is obtained by analyzing the model (3.1). However, owing to the complexity of the model, we cannot obtain the explicit expression of the threshold $R_{1}.$ To describe the effect of parameters on the threshold conditions, the method in ^[26] is used to give a sensitivity analysis of the threshold $R_{1}.$ Owing to the lack of information, we choose uniform distribution for all parameters and possible ranges of parameters in Table 1. It can be observed from Figure 3 that among the four control parameters, the increasing $\alpha$ and $\rho$ can decrease $R_{1},$ while increasing $p$ and $T$ leads to an increase in $R_{1}.$ As presented in Figure 4, we fixed other parameters except the removal rate $p$ of healthy plants and applied numerical simulation to explore the impact of the removal rate $p$ of healthy plants on the threshold $R_1$. The results demonstrated that the threshold $R_1$ increased with the increasing removal rate $p$ of healthy plants. Therefore, when controlling the spread of plant diseases, we should attempt to avoid removing healthy plants. Moreover, because $R_{1}$ depends on $\alpha$ and $\rho$, we use a three-dimensional surface (Figure 5) to illustrate the dependence of $R_{1}$ on these two parameters. The numerical simulation indicates that high values of $\rho$ and $\alpha$ will guarantee $R_{1} < 1.$ As a result, when the control strategy is implemented, not only should the infected plants be removed, but also the vectors should be eliminated together, aiming to ensure that the plant virus disease can be controlled.

Finally, it can be found from Theorem 3.4 that system (3.1) may have backward bifurcation. In Figure 6, we select the parameter to make $R_{1} = 0.9775 < 1.$ For different initial values, it is found that (1) the system with initial value $(2, 0.1, 1000, 15)$ has a stable disease eradication periodic solution; and (2) with the initial value being $(2, 0.1, 1000,150),$ the system has a stable positive periodic solution.

5.   Discussion

Roguing control is a very vital and effective control strategy for plant virus transmission. In ^[27], the author considered the removal of infected plants and did not consider the elimination of vectors on plants when removing infected plants. We consider two kinds of plant virus transmissions control strategies to remove infected plants and vectors in a continuous and pulse manner.

For control strategies with continuous removal of infected plants and vectors, through the analysis of system (2.1), we obtain that when the threshold $R_{0} < 1$, the disease-free equilibrium is locally asymptotically stable, and the system may have two positive equilibrium points. Using the conclusion in ^[21], we demonstrate the sufficient condition for the existence of backward bifurcation in system (2.1). Based on the results, the disease can be eliminated only when $R_{0} < R_{0}(\alpha_{0}),$ not $R_{0} < 1.$

For control strategies with an impulse to remove infected plants and vectors, through the analysis of system (3.1), we find that when the threshold $R_{1} < 1$, the disease-free periodic solution of system (3.1) is locally asymptotically stable. Due to the complexity of the model, we use a strong threshold condition threshold $R_{2} > 1$ ($R_{2} < R_{1}$) to prove the persistence of the disease. Using the nonlinear fixed point theory, we give the conditions for forward and backward bifurcation of impulsive control systems. The results show that the disease can be eradicated only when $R_{1} < R_{1}(\alpha_{0}),$ not $R_{1} < 1.$

Backward bifurcation of the system indicates that based on different initial conditions, the system may have stable disease free periodic solutions or stable positive periodic solutions, which undoubtedly brings great difficulties to plant disease control. Our results suggest that there is a backward bifurcation phenomenon in continuous systems, and there is also a backward bifurcation phenomenon in impulsive control systems. Therefore, when controlling plant virus diseases, if the initial infection scale is large, we should kill a large number of plants and vectors to control the disease.

Finally, we note that the introduction of natural enemies will control the number of vectors. This is undoubtedly beneficial to the spread of plant viruses. Therefore, how to establish a plant virus transmission model with the mixed control of natural enemies and the removal of infected plants will be our focus in the later stage.

Use of AI tools declaration

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

Acknowledgments

The current work was supported by the National Natural Science Foundation of Qinghai Province (No.2022-ZJ-T02).

Conflict of interest

We declare that we have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

A.   The Expression of each element of matrix $D_{\xi}(0, \textbf{O})$

From

at $\xi_{0} = (\xi_{10}, \xi_{20}, 0, 0)$ ($\xi_{10} = x^{*}_{1}(0), \xi_{10} = x^{*}_{1}(0)$), we have

From the variational equation related to the first four equation of system (3.29),

with the initial condition $D_{\xi}\Phi(0, \xi_{0}) = E_{4}$ and $\Phi(t, \xi_{0}) = (\Phi_{1}(t, \xi_{0}), \Phi_{2}(t, \xi_{0}), 0, 0),$ we obtain

Then, from (3.30) and (A.3), we get

where

B.   The first partial derivatives of $f$

We obtain:

Since

then

Consider Eq (3.41) as $i = 1,$ we have

One can similarly acquire from Eq (3.41) as $i = 2, 3$ that

It can be deduced from (B.4) and (B.5) that

Thus $\frac{{\partial f(0, \textbf{O})}}{{\partial \bar \alpha }} = 0.$ Additionally,

because $Y_{1}$ is a basis in Ker($D_{\xi}(0, \textbf{O})$, that is:

Therefore, we can deduce from (12.3) and (B.7) that

Similarly, as $i = 2, 3,$ we can obtain that

From (B.9) and (B.10) we obtain that

Since

submitting (B.11) into (B.12), we get $\frac{\partial f(0, \textbf{O})}{\partial\alpha_{1}}.$

C.   The second-order derivatives of $f$

Let $A = \frac{\partial^{2} f(0, \textbf{O})}{\partial\bar{\alpha}^{2}},$ $B = \frac{\partial^{2} f(0, \textbf{O})}{\partial\bar{\alpha}\partial\alpha_{1}},$ $C = \frac{\partial^{2} f(0, \textbf{O})}{\partial\alpha^{2}_{1}}.$

Calculation of $A.$

submitting (B.2) and (B.6) into (C.1), it can be deduced that

Consider Eq (3.41) as $i = 3,$ we obtains that

Similarly, by submitting (B.2) and (B.6) into (C.3), we can obtain that

then $\frac{\partial^{2}\tilde{\alpha}_{4}(0, \textbf{O})}{\partial\bar{\alpha}^{2}} = 0.$ By substitute the result of the above formula into (C.2), we have that $A = \frac{{{\partial ^2}f(0, O)}}{{\partial {{\bar \alpha }^2}}} = 0.$

Calculation of $B.$

First, we calculate the value of $\frac{\partial^{2}\tilde{\alpha_{4}}(0, \textbf{O})}{\partial\bar{\alpha}\partial{\alpha_{1}}},$

submitting (B.2) and (B.6) into (C.5), we can therefore deduce that

It can be calculated that

submitting (B.2), (B.6) and (C.6) into (C.7), we have

Calculation of $C.$

Consider Eq (3.41) as $i = 1,$ we obtains that

then we have that

Similarly, we can obtain from Eq (3.41) as $i = 2, 3$ that

By solving (C.11) and (C.12), we can get the values of $\frac{\partial^{2}\tilde{\alpha}_{i}(0, \textbf{O})}{\partial\alpha^{2}_{1}},$ $i = 2, 3, 4,$ and submit it as $i = 4$ into (C.9), one can obtain: