1 Introduction

Mosaic disease of plants, caused by begomoviruses, affects a variety of important agricultural crops around the world, including tomatos, peppers, cucumbers, aubergines, with estimated losses reaching 40–70% of annual harvest. In the tropical regions, mosaic disease is, perhaps, the biggest challenge to growing cassava and Jatropha. With cassava being a major staple food crop in many parts of the tropics, and Jatropha curcas emerging recently as one of the most promising sources of biofuel, the mosaic disease of these crops is having a profound negative effect on communities growing them. Mosaic disease is obligately transmitted to plants by insect vectors, primarily the whitefly B. tabaci, and results in mottled and distorted leaves, veinal netting, as well as stunting of plants, often causing them to produce virtually no yield.

In terms of control of mosaic disease, a number of different strategies have been proposed, some focusing on plants, such as roguing (removing infected plant biomass) and replanting (replacing the infected plants with healthy ones), with others targeting disease vectors, e.g. using insecticides, biological control, and developing vector-resistant plant varieties [1, 2]. A number of mathematical models have analysed the effectiveness of these approaches, and also studied specific optimal strategies for disease control. One of the earliest models, proposed by Chan and Jeger [3], looked into whether roguing can be used to eradicate mosaic disease. Since that model did not properly account for various characteristics of the process of disease transmission, it was subsequently improved in Jeger et al. [4] to include the details of virus development in vectors. More advanced models of roguing have been analysed by Jeger et al. [1], Gao et al. [5], Luo et al. [6], and Basir and Roy [7]. Spatial aspects of various strategies for disease management associated with roguing and replanting have been analysed by Sisterson and Stenger [8]. These were also included in the bio-economic analysis of disease control measures by Atallah et al. [9], while Venturino et al. [10] considered optimal control problem of using insecticide for control of mosaic virus disease in J. curcas plantations. Jeger et al. [2] provide a nice overview of modelling work on control of plant virus disease in different settings.

In recent years, several promising strategies have been put forward for control of plant disease and protection of crops against parasites, while avoiding potential dangers to humans and environment, associated with the use of chemical pesticides. Among these strategies, an important role is played by developing tools based on RNA interference (RNAi), which is a fundamental biological process, through which eukaryotic cells are able to post-transcriptionally control expression of specific genes [11,12,13,14]. In the specific context of plant protection against parasites, RNAi can be used in two different ways. It can protect plants from infection by stopping the expression of plant genes that are essentially involved in the process of infection, or by halting the expression of viral genes within plant cells [15,16,17,18]. Alternatively, it can be used to target specific parasites that consume RNAi complexes when feeding on plants, and this then results in the reduction of their fecundity and increased mortality through silencing of essential parasite housekeeping genes [13, 19, 20]. Both of these approaches rely on supplying plants with appropriate dsRNA (double-stranded RNA) that would trigger RNAi process, and this dsRNA can be delivered into plants through various means, including development of transgenic plants [21, 22], root soaking [13, 23], as well as topical application (spraying) [24, 25]. Recent work has demonstrated the feasibility of RNAi for control of B. tabaci whitefly [26,27,28,29,30], while mathematical models have elucidated a number of important aspects of RNAi in plants [31,32,33].

An alternative but related approach to improving performance of crops and protection from disease is based on the use of biopesticides, or natural microbial biostimulants (MBs) obtained from metabolites of soil micro-organisms [34,35,36,37,38]. Among various sources of natural biostimulants are soil actinobacteria of Streptomyces and Saccharopolyspora genera [39,40,41,42]. One prominent example is Streptomyces avermitilis, which produces avermectins, known to have nematicidal and insecticidal properties [43,44,45,46], and a similar effect has been recently demonstrated for other streptomycetes [39, 40, 47]. Recent work has shown that MBs produced as combinations of metabolites and cell culture supernatant of such soil streptomycetes as S. netropsis [48], and S. violaceus [49] provide a wide range of plants with strong antagonistic activity against various phytopathogenic bacteria and micromycetes. These MBs contain a wide range of antibacterial and plant stimulating compounds, as well as amino acids, lipids and plant hormones, which together significantly contribute towards disease resistance and improvement of plant growth and development [50, 51]. This positive effect of biostimulants on plant growth has been demonstrated through improved callogenesis and organogenesis in wheat [52, 53]. In terms of protection against parasites, extensive experimental studies have comprehensively demonstrated that these MBs trigger within-plant production of RNAi products complementary with mRNA (messenger RNA) of plant parasitic nematodes, and upon consumption by nematodes, these products trigger their death by silencing essential nematode housekeeping genes. More specifically, nematicidal effects of these MBs on root-knot nematode M. incognita have been shown in vitro [54], and specific RNAi-mediated bioprotective effects of MBs have been demonstrated in wheat against the cereal cyst nematode H. avenae [52, 53], in rapeseed [55, 56], Brassica rapa subs. pekinensis (Chinese cabbage) [57], and sugar beet [58] against the cyst nematode H. schachtii, and in cucumber and sugar beet against both of these nematodes [59]. A natural biostimulant Regoplant has been shown to provide effective RNAi-based protection of common horse chestnut against the horse-chestnut leaf miner Cameraria ohridella [60].

In this article, we propose and analyse a mathematical model for control of mosaic disease using RNAi-mediated bioprotective effects of MBs. In the next section we derive the model and discuss its basic properties. Section 3 contains analysis of feasibility and stability of different equilibria. In Sect. 4, we perform numerical stability analysis and simulations of the model. The paper concludes in Sect. 5 with the discussion of results and future work.

2 Model derivation

To analyse how the mosaic disease can be controlled using MBs, we consider a population of plants (this can be J. curcas or cassava) that can be exposed to mosaic disease, which is spread by the B. tabaci whitefly. We assume that in the absence of mosaic disease, healthy plants, whose population is denoted by X(t), reproduce logistically with the linear growth rate r and the carrying capacity K. Once whiteflies infect a healthy plant, it becomes infected and moves into a population of infected plants Y(t). Following Holt et al. [61], we do not distinguish between individual parts of the plant, such as stem, leaves etc. that can become infected, but rather consider the overall plant biomass. To simplify the model, we do not explicitly include a separate compartment for the virus population, but focus instead on their whitefly vectors. Since the vectors are relying on plants for their food source, in the model we assume that the uninfected vectors U(t) grow logistically with a growth rate b and a carrying capacity given by the total plant population \(X+Y\) multiplied by the maximum vector abundance m.

The virus that causes the mosaic disease is assumed to be transmitted from infected plants to uninfected vectors at rate \(\beta \), and from infected vectors to uninfected plants at rate \(\lambda \). Begomoviruses that cause mosaic disease are circulative-persistent viruses [62], which means that once the whitefly vectors become infected, they will remain infectious for the rest of their lifetime [1, 61], as is common with vector-borne infections. The reason for this is that when whiteflies feed on infected plants, they ingest the virus contained in the plant sap with their stylets, and subsequently the virus crosses the filter chamber and the midgut to be then translocated into the primary salivary glands [62, 63]. When these vectors then feed on healthy plants, virus particles circulating in the whitefly saliva will enter these plants and start infection in them. We will denote by a 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.

Finally, MBs are applied to plants at rate B, and they will have an effect on the plant growth rate r and the rate of disease transmission from vectors to plants \(\lambda \), and they will also cause additional mortality of vectors at rate c.

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

$$\begin{aligned} \begin{array}{l} \displaystyle {\frac{dX}{dt} = r(B)X\left[ 1-\frac{X+Y}{K}\right] -\lambda (B) XV,}\\ \\ \displaystyle {\frac{dY}{dt} = \lambda (B)XV-aY,}\\ \\ \displaystyle {\frac{dU}{dt} = b (U+V)\left[ 1-\frac{U+V}{m(X+Y)}\right] -\beta U Y-c(B)(X+Y)U,}\\ \\ \displaystyle {\frac{dV}{dt} = \beta UY - \mu V-c(B)(X+Y)V.} \end{array} \end{aligned}$$
(1)

Experimental evidence suggests that streptomycete-derived MBs stimulate growth and viability of crops, as has been recently demonstrated during in vitro, in vivo and greenhouse experiments on wheat [52, 53, 64], cucumbers [56], rapeseed [59], and Chinese cabbage [57]. Thus, we replace a constant growth rate of uninfected plants by a growing saturated function of the amount of applied MBs

$$\begin{aligned} r(B)=\frac{r_0(1+B)}{1+sB}, \end{aligned}$$
(2)

with \(r(0) = r_0\) and \(r(\infty ) = r_0/s\). Molecular genetic analyses show that MBs trigger a significant increase in the within-plant synthesis of si/miRNA homologous with plant genes associated with proteins that are important for a successful infection, as well as with viral genes, whose expression is needed for viral replication in plants [19, 53, 55, 56, 58, 59]. To model this effect, we take the rate of disease transmission from infected vectors to plants in the form of a decreasing function of the amount of applied MBs as follows,

$$\begin{aligned} \lambda (B)=\frac{\lambda _0}{1+pB}. \end{aligned}$$
(3)

Finally, once whiteflies feed on plants that have been treated with MBs (irrespective of whether the plants themselves are healthy or infected), the RNAi products they ingest will result in vector mortality through silencing essential vector genes [13, 19, 20]. Formally, this is represented in the model as an additional death term for vector population, with the mortality rate being proportional to the total plant biomass, and the proportionality rate c(B) being given by a growing function of the amount of applied MBs,

$$\begin{aligned} c(B)=\frac{c_0B}{1+c_1B}. \end{aligned}$$
(4)

In the case where no MBs are used, we have \(c(0)=0\), and for an increasing amount of biostimulants we have a saturation at \(c(\infty )=c_0/c_1\).

It is straightforward to show that, subject to non-negative initial conditions, all state variables will remain non-negative and bounded for all \(t\ge 0\), which implies that the model is well-posed.

3 Steady states and their stability

For any parameter values, the model (1) has an axial equilibrium \(E_1{=}(K,0,0,0)\) characterised by the absence of vector population and the presence of only healthy plants. Provided \(b{>}cK\), it also has a disease-free steady state \(\displaystyle {E_2{=}(K,0,mK(b{-}cK){/}b,0)}\). Finally, the model (1) can also have one or more endemic equilibria \(E^*=(X^*,Y^*,U^*,V^*)\) with

$$\begin{aligned} Y^* = \frac{rX^*(K-X^*)}{aK+rX^*},\quad V^* = \frac{aY^*}{\lambda X^*}, \quad U^*=\frac{a(\mu +cX^*+cY^*)}{\lambda \beta X^*}, \end{aligned}$$

where \(X^*\) is the positive root of the following equation

$$\begin{aligned} b(U^*{+}V^*)\Big [m(X^*{+}Y^*)-(U^*{+}V^*)\Big ]-m U^* (X^*{+}Y^*)\Big [\beta Y^*{+}c(X^*{+}Y^*)\Big ]{=}0.\nonumber \\ \end{aligned}$$
(5)

Explicitly, this is a quintic equation

$$\begin{aligned} \omega _0(X^*)^5+\omega _1(X^*)^4+\omega _2(X^*)^3+\omega _3(X^*)^2+\omega _4X^*+\omega _5=0, \end{aligned}$$
(6)

with the coefficients

$$\begin{aligned} \omega _0= & {} brL_1-m\lambda \beta K(a+r) F_1,\quad \omega _1=abKL_1+brL_2-m\lambda \beta K(a+r) F_2, \\ \omega _2= & {} abKL_2+br L_3-m\lambda \beta K(a+r) F_3, \quad \omega _3=abKL_3+br L_4, \\ \omega _4= & {} abKL_4+brL_5,~~~~\omega _5=abKL_5, \end{aligned}$$

where

$$\begin{aligned} L_1= & {} A_1B_1,\quad L_2=A_1B_2+A_2B_1,\quad L_3=A_1B_3+A_3B_1+A_2B_2,~~\\ L_4= & {} A_2B_3+A_2B_3,~~L_5=A_3B_3,\quad F_1 = -\beta r(\mu r+caK+crK),\\ F_2= & {} K(\mu r+caK+crK)(ca+cr+\beta r),\quad F_3=\mu a K^2(ca+cr+\beta r),\\ A_1= & {} -\beta r,~~~ A_2=\mu r+caK+crK+\beta cK,\quad A_3=\mu a K,\\ B_1= & {} a \beta r+\lambda m \beta K(a+r),\quad B_2{=}-a[\mu r+caK+crK+\beta r K],\quad B_3{=}-a^2\mu K. \end{aligned}$$

The endemic equilibrium exists only if \(K-X^*>0\).

The Jacobian of the system (1) at any steady state \({\bar{E}}({\bar{X}},{\bar{Y}},{\bar{U}},{\bar{V}})\) is given by

$$\begin{aligned} J=[J_{ij}]=\left[ \begin{array}{cccc} J_{11} &{} ~~-\frac{r{\bar{X}}}{K} &{} 0 &{} -\lambda {\bar{X}}\\ \ \\ \lambda {\bar{V}} &{} ~~-a &{} 0 &{} \lambda {\bar{X}}\\ \ \\ \frac{b({\bar{U}}+{\bar{V}})^2}{m({\bar{X}}+{\bar{Y}})^2} -c{\bar{U}}&{} ~~\frac{b({\bar{U}}+{\bar{V}})^2}{m({\bar{X}}+{\bar{Y}})^2}-\beta {\bar{U}}-c{\bar{U}}&{} ~~~J_{33}&{} ~~~ b-\frac{2b({\bar{U}}+{\bar{V}})}{m({\bar{X}}+{\bar{Y}})}\\ \ \\ -c{\bar{V}}&{} \beta {\bar{U}}-c{\bar{V}}&{} \beta {\bar{Y}}&{} -\mu -c({\bar{X}}+{\bar{Y}})\\ \end{array} \right] , \end{aligned}$$

where,

$$\begin{aligned} J_{11}=r\left[ 1-\frac{(2{\bar{X}}+{\bar{Y}})}{K}\right] -\lambda {\bar{V}}, \text{ and } J_{33}=b-\beta {\bar{Y}}-\frac{2b({\bar{U}}+{\bar{V}})}{m({\bar{X}}+{\bar{Y}})} -c({\bar{X}}+{\bar{Y}}). \end{aligned}$$

At the axial equilibrium \(E_1(K,0,0,0)\), the eigenvalues of the Jacobian matrix are \(-\,r\), \(-\,a\), \(b-c(B)K\), and \(-\,\mu -cK\), suggesting that this steady state is stable for \(b<c(B)K\) and unstable for \(b>c(B)K\). At \(b=c(B)K\), the steady state \(E_1\) undergoes a transcritical bifurcation, at which point it loses stability, and the disease-free steady state \(E_2\) emerges.

If we define

$$\begin{aligned} R_0=\frac{mK^2 \beta \lambda (B)(b-c(B)K)}{ab(\mu +c(B)K)}, \end{aligned}$$
(7)

as the basic reproduction number, then we have the following result.

Proposition 1

Suppose the condition \(b>c(B)K\) holds. The disease-free equilibrium \(E_2\) is stable for \(R_0<1\), unstable for \(R_0>1\), and undergoes a steady-state bifurcation at \(R_0=1\).

Proof

At the steady state \(E_2\), two characteristic eigenvalues can be easily found as \(-b\) and \(-r\), while the remaining eigenvalues satisfy a quadratic equation

$$\begin{aligned} \xi ^2+(a+\mu +c(B)K)\xi +\frac{1}{b}\left[ ab(\mu +c(B)K)-mK^2 \beta \lambda (B)(b-c(B)K)\right] =0.\nonumber \\ \end{aligned}$$
(8)

Since \(a+\mu +c(B)K>0\), the condition for stability of \(E_2\) reduces to

$$\begin{aligned} ab\left( \mu +c(B)K\right) -mK^2\beta \lambda (B)\left( b-c(B)K\right) >0\quad \Longleftrightarrow \quad R_0<1. \end{aligned}$$

Thus, for \(R_0<1\) the steady state \(E_2\) is stable, and for \(R_0>1\) it is unstable. As \(R_0\) increases from \(R_0<1\), at the moment it passes through 1, one of the two roots of the quadratic equation goes through zero from left to right along the real axis, making the steady state \(E_2\) unstable, which implies that at \(R_0=1\), \(E_2\) undergoes a steady-state bifurcation. \(\square \)

The characteristic equation at the endemic equilibrium \(E^*\) has the form

$$\begin{aligned} \xi ^4 + \alpha _1 \xi ^3 + \alpha _2\xi ^2 + \alpha _3\xi + \alpha _4 = 0, \end{aligned}$$
(9)

where,

$$\begin{aligned} \alpha _1= & {} - J_{11} - J_{22} - J_{33} - J_{44},\\\\ \alpha _2= & {} J_{11}(J_{22}+J_{33}+J_{44}-J_{21})-J_{14}(J_{41}-J_{42})+J_{33}(J_{22}+J_{44})-J_{34}J_{43},\\\\ \alpha _3= & {} J_{11}[(J_{33}+J_{44})(J_{21}-J_{22})+J_{24}(J_{42}-J_{41})]+(J_{11}+J_{22})(J_{34}J_{43}-J_{33}J_{44})\\&-\,J_{14}[J_{41}(J_{22}+J_{33})-J_{42}(J_{33}+J_{21})+J_{43}(J_{32}-J_{31})],\\\\ \alpha _4= & {} J_{14}[(J_{11}+J_{22})(J_{33}J_{41}+J_{31}J_{43})+(J_{11}+J_{21})(J_{33}J_{42}-J_{32}J_{43})]\\&+\,J_{11}(J_{21}-J_{22})(J_{34}J_{43}-J_{33}J_{44}), \end{aligned}$$

and we have used the fact that \(J_{12}=J_{11}\) and \(J_{24}=-J_{14}\). Since the values of all parameters and state variables at the endemic equilibrium \(E^*\) are positive, we have

$$\begin{aligned} J_{22}=-a<0,\quad J_{44}=-\mu -c(X^*+Y^*)<0, \quad J_{11}=-\frac{rX^*}{K}<0, \end{aligned}$$

and

$$\begin{aligned} \displaystyle {J_{33}=-bV^*\left[ 1-\frac{U^*+V^*}{m(X^*+Y^*)}\right] -\frac{bU^*(U^*+V^*)}{m(X^*+Y^*)}}. \end{aligned}$$

From equation (5) it follows that \(U^*+V^*<m(X^*+Y^*)\), which implies that \(J_{33}<0\). Taken together, we then find that

$$\begin{aligned} \alpha _1=- J_{11} - J_{22} - J_{33} - J_{44}>0, \end{aligned}$$

for any parameter values, for which the endemic steady state \(E^*\) is biologically feasible. Considering \(\gamma \) as a generic bifurcation parameter (i.e. any of the parameters of the model), it is possible to prove the following result [7, 10].

Theorem 1

  1. (a)

    The endemic equilibrium \(E^*\) is stable, if the following conditions are satisfied

    $$\begin{aligned} \nonumber&\alpha _2> 0,\quad \alpha _3> 0,\quad \alpha _4> 0,\\&\alpha _1 \alpha _2 - \alpha _3> 0,\quad \alpha _1\alpha _2\alpha _3-\alpha _3^2-\alpha _4\alpha _1^2>0. \end{aligned}$$
    (10)
  2. (b)

    The endemic equilibrium \(E^*\) undergoes a Hopf bifurcation at \(\gamma =\gamma ^*\in (0,\infty )\) if and only if

    $$\begin{aligned} \begin{array}{l} \alpha _2(\gamma ^*)> 0,\quad \alpha _3(\gamma ^*)> 0,\quad \alpha _4(\gamma ^*)>0, \quad \alpha _1(\gamma ^*)\alpha _2(\gamma ^*)- \alpha _3(\gamma ^*)>0,\\ \\ \alpha _1(\gamma ^*)\alpha _2(\gamma ^*)\alpha _3(\gamma ^*)-\alpha _3^2(\gamma ^*) -\alpha _4(\gamma ^*)\alpha _1^2(\gamma ^*)=0,\\ \\ \alpha _1^3(\gamma ^*)\alpha _2'(\gamma ^*)\alpha _3(\gamma ^*)\left[ \alpha _1(\gamma ^*)-3\alpha _3(\gamma ^*)\right] \ne \left[ \alpha _2(\gamma ^*)\alpha _1^2(\gamma ^*)-2\alpha _3^2(\gamma ^*)\right] \\ \times \left[ \alpha _3'(\gamma ^*)\alpha _1^2(\gamma ^*)-\alpha _1'(\gamma ^*)\alpha _3^2(\gamma ^*)\right] , \end{array} \end{aligned}$$
    (11)

    where primes denote differentiation with respect to \(\gamma \).

4 Numerical results

To better understand how different parameters affect the dynamics of the model, in this section we explore stability of different steady states numerically, as well as solve the model directly to illustrate its behaviour in different dynamical regimes. Baseline values of parameters are taken from Holt et al. [61], Venturino et al. [10], and Al Basir and Roy [7]. Regarding carrying capacity of plants, we consider a plantation with a maximum of 50 plants, which, with a recommended spacing of around 3 m \(\times \) 3 m for cultivating Jatropha [65] would correspond to a plot of 450 \(\hbox {m}^2\), which is consistent with average-size fields and farm plots used for growing Jatropha [66, 67] and cassava [68]. It should be noted that although we have chosen to fix the maximum number of plants, an alternative is to rather treat this carrying capacity as plant density per some pre-defined plot area [1, 61], so that one could then consider a large plantation or a field consisting of a multitude of such plots. Yet another alternative is to instead interpret the carrying capacity as a characteristic of cumulative weight of a single plant, and the biomass of entire plantation could be obtained by multiplying this weight by the overall number of plants. Similarly to carrying capacity, there is significant variability in reported/observed numbers of whiteflies feeding on a single plant, ranging from 0 to 250 [69, 70] vectors per plant, and we have chosen the maximum vector abundance per plant to have the baseline value of \(m=80\). In Fig. 1 we fix the values of all parameters and allow the rate \(\lambda _0\) of disease transmission from infected vectors to plants to vary, which allows us to explore how the steady-state values of the infected biomass and vector population change with the basic reproduction number. For \(R_0<1\), the endemic steady state is not feasible, and, in accordance with Proposition 1, the disease-free steady state \(E_2\) is feasible and stable. As the basic reproduction number increases and crosses unity, the endemic steady state \(E^*\) becomes biologically feasible, and for small values of \(R_0\) it is stable, with the values of \(Y^*\) and \(V^*\) increasing with \(R_0\). For some value of the basic reproduction number, the equilibrium values of \(Y^*\) and \(V^*\) reach their maxima, and subsequent increase of \(R_0\) results in the reduction of these values. Furthermore, for sufficiently high value of \(R_0\), the steady state \(E^*\) loses its stability via a Hopf bifurcation in agreement with Theorem 1.

Table 1 Values of the parameters used in the numerical simulations [7, 10, 61]
Full size table
Fig. 1
figure 1

Steady state values of the infected plant biomass \(Y^*\) and the infected vector \(V^*\) population depending on \(R_0\). The parameter values are \(c_0=0.05,c_1=5,\) \(p =0.05,~s =0.01,~ B=0.2\), and the remaining parameter values are taken from Table 1. Solid (dashed) line denotes stable (unstable) steady state

Full size image

Figure 2 illustrates the behaviour of the model for different values of \(\lambda _0\), with other parameters chosen in such a way that it ensures biological feasibility of the endemic equilibrium \(E^*\). For sufficiently small \(\lambda _0<\lambda _0^*\), the equilibrium \(E^*\) is stable, and for \(\lambda >\lambda _0^*\), the Hopf bifurcation has taken place, and the system exhibits sustained periodic oscillations around \(E^*\). Further details of this transition to instability are provided in Fig. 3, which shows a bifurcation diagram for the steady state \(E^*\) depending on the transmission rate \(\lambda _0\). One observes that higher values of the disease transmission rate not only cause instability of the equilibrium \(E^*\), but also result in a higher amplitude of oscillations around this steady state. Figure 4 illustrates similar behaviour when all parameters are fixed, and the baseline rate \(c_0\), at which MBs cause death of vectors through RNAi, is varied. Counter intuitively, we observe that increasing this rate, i.e. making MBs more efficient, results in the suppression of oscillations and stabilisation of the endemic equilibrium \(E^*\).

Fig. 2
figure 2

Numerical solution of the model (1) for parameter values as in Fig. 1 and \(\lambda _0=0.0005\) (solid), \(\lambda _0=0.0008\) (dashed), and \(\lambda _0=0.00101\) (dotted)

Full size image
Fig. 3
figure 3

Bifurcation diagram of the endemic steady state \(E^*\) depending on the transmission rate \(\lambda _0\), with other parameter values as in Fig. 1. Solid line shows a stable \(E^*\), while dots indicate minimum and maximum values of the periodic solution around \(E^*\) where it is unstable

Full size image
Fig. 4
figure 4

Numerical solution of the model (1) for parameter values as in Fig. 1 with \(\lambda _0=0.0025\), and \(c_0=0.1\) (solid), \(c_0=0.08\) (dashed), and \(c_0=0.05\) (dotted)

Full size image

Figure 5 demonstrates different types of dynamics that can be exhibited by the system depending on various parameters. Figure 5a suggests that for a sufficiently small value of the product of two disease transmission rates \(\lambda \) and \(\beta \), the disease-free steady state \(E_2\) is stable, and as the value of this product increases, this steady state loses its stability via a steady-state bifurcation as discussed in Proposition 1, giving rise to a stable endemic steady state \(E^*\). Surprisingly, a further increase of the product of two disease transmission rates eventually destabilises \(E^*\) via Hopf bifurcation, as discussed in Theorem 1, leading to the onset of stable periodic oscillations around this state. Figure 5b shows that for a higher rate of disease transmission, a larger amount of biostimulants B is required to eradicate the disease. Parameter values in this plot satisfy the condition \(b>c(B)K\) for any positive value of B, which means that the axial equilibrium \(E_1\) is always unstable, while the disease-free equilibrium \(E_2\) is feasible everywhere. Figure 5c indicates that increasing the rate of RNAi-induced mortality of vectors \(c_0\) first results in the stabilisation of the endemic steady state, and then it eliminates the disease by making the endemic steady state infeasible and stabilising the disease-free steady state, and finally stabilises the axial equilibrium \(E_1\), characterised by a complete eradication of vectors. Figure 5d shows that the total amount of biostimulants B plays a complementary role in disease control in that increasing either of \(c_0\) or B leads to a sequence of transitions between stable endemic, disease-free, and axial steady states. For smaller values of B and \(c_0\), the endemic steady state \(E^*\) exists and is unstable, and increasing either of those parameters leads to stabilisation of this steady state. A further increase in those parameters, which is equivalent to decreasing the value of the transmission rate \(\lambda (B)\) defined in (3) and increasing the amount of applied biostimulants c(B) defined in (4), results in reducing the value of the basic reproduction number \(R_0\) as given in (7). As soon as the value of \(R_0\) crosses the threshold value of 1, in accordance with Proposition 1, the disease-free steady state \(E_2\) becomes stable. A further increase in B and/or \(c_0\) eventually results in violating the condition \(b>c(B)K\), resulting in making the disease-free steady state \(E_2\) biologically infeasible, and the axial equilibrium \(E_1\) stable.

Fig. 5
figure 5

Stability of steady states of the system (1). Disease-free equilibrium \(E_2\) is stable in red region; endemic equilibrium \(E^*\) is stable in violet region and unstable in yellow region; axial equilibrium \(E_1\) is stable in light blue region. Parameter values are taken from Table 1, except for \( p=1.2, s=0.01, c_1=2.5\). a \(B=5.2\), \(c_0 =0.01\), b \(\beta =0.0025\), \(c_0 =0.01\), c \(c_1=5, \beta =0.0012, B=5.2\), d \(c_1=5, \beta =0.0012, \lambda _0=0.0005\)

Full size image

To obtain a better insight into how the biostimulant affects the endemic steady state \(E^*\), in Fig. 6 we have plotted the bifurcation diagram of this steady state with B being a free parameter. This figure illustrates how increasing the amount of biostimulant results in the stabilisation of \(E^*\) as B crosses some critical value \(B^*\), in agreement with Theorem 1. One also observes that when \(E^*\) is stable, for sufficiently small values of the applied biostimulant B, all components of the endemic steady state are growing with B, suggesting a counter-intuitive result that through enhancement of plant growth, the biostimulant also effectively sustains the presence of infection. For much higher values of B, the populations of infected plants and vectors decrease, until eventually the endemic steady state becomes infeasible, and the system tends to a stable disease-free steady state, as shown earlier in Fig. 5.

Fig. 6
figure 6

Bifurcation diagram of the endemic steady state \(E^*\), depending on the amount of biostimulants, B, with parameter values taken from Table 1, except for \(c_0=0.05, c_1=5, p =0.05, s =0.01, \lambda _0=0.0025\). Solid line shows a stable \(E^*\), while dots indicate minimum and maximum values of the periodic solution around \(E^*\) where it is unstable

Full size image

5 Discussion

Due to the wide spread of mosaic disease in different geographical regions, and its major negative economic and societal effects associated with a substantial reduction in yields of major agricultural crops, it is of paramount importance to establish effective means for control of this disease. In this paper we have analysed a mathematical model for the use of natural microbial biostimulants to improve plant performance and protect them against mosaic disease. The protective effect of MBs based on RNA interference was taken to be in the form of reducing the level of transmission of mosaic virus from infected whitefly to healthy plants, as well as in causing mortality of whitefly by silencing their essential housekeeping genes. Analysis of stability of the model’s steady states has revealed how different parameters characterising the performance of MBs, such as their ability to reduce disease transmission from whitefly to plants, or their effect on the whitefly mortality can result in eradication of the disease, maintaining it at some steady level, or causing periodic oscillations in the level of disease. One interesting observation is that increasing the rate of disease transmission from either vectors to plants, or from plants to vectors, actually causes an instability of the endemic steady state via a Hopf bifurcation, which is associated with the emergence of stable oscillations in the populations of plants and vectors. For the fixed parameters characterising disease transmission, increasing the amount of MBs being used, or increasing their effectiveness in causing vector mortality, first results in suppressing oscillations and stabilising disease at some steady level, then eradicating it and moving the system to a stable disease-free steady state, and eventually in stabilising a boundary equilibrium, which represents a complete elimination of the vector population. Again, one should note a rather counterintuitive result that for intermediate values of these two parameters increasing their values leads to a stabilisation of endemic steady state. Our analysis establishes critical values of parameters required for disease eradication, while also showing various alternative mechanisms through which basic reproduction number can be brought to the value below unity, which provides several alternative practical routes for control of mosaic disease in specific farming circumstances.

There are several directions, in which the work presented in this paper can be extended to make it more realistic, as well as more practically oriented. In the current version of the model, we only included the effect of MBs on reducing disease transmission from infected vectors to healthy plants, and this restriction can be removed by also considering an analogous effect of reducing disease transmission from plants to vector, based on the idea that successful within-plant RNAi reduces the amount of viable free virus particles that are produced within plant cells and which can be consumed by vector when feeding on plants. One can also include in the model interactions between the effects of MBs and other available strategies for control of mosaic disease, such as roguing [7]. Since the use of MBs largely relies on farmers’ awareness about the mosaic disease, one can also consider how different types of disease awareness and corresponding use of MBs affect progression and control of mosaic disease, as has been recently explored in the context of using nutrients and insecticides [71]. It should be noted that there are different methods of applying MBs to plants, including soaking of seeds prior to sowing, applying them to the medium, in which the plants are growing, or spraying them directly on growing plants [52, 54, 64]. Each of these approaches has its own benefits and limitations in terms of costs, ease of implementation, and the efficiency of inducing within-plant production of relevant RNAi compounds. Laboratory and field experiments have shown the effectiveness of MBs in protecting cereal crops and vegetables against insects and plant parasitic nematodes [50, 52], clearly indicating the feasibility of using MBs to protect plants against parasites and pathogens, in a manner similar to existing approaches using RNAi for protecting crops against plant viruses [13, 25, 72]. This paper has identified main effects of MBs on the dynamics of mosaic disease from the perspective of disease control and eradication. At the same time, to ensure that MBs deliver their desired effect in an optimal way in terms of cost-effectiveness and minimum negative effect on the environment, one can extend the analysis of our model by considering it from the perspective of optimal control theory in a manner similar to Venturino et al. [10]. Combined with more detailed experimental data on transmission of begomovirus between plants and vectors, this would then guide the development of practical recommendations on the best use of microbial biostimulants subject to various logistical and financial constraints, which would greatly facilitate the implementation of control of mosaic disease by farmers in different settings.