Impact of ontogenic changes on the dynamics of a fungal crop disease model motivated by coffee leaf rust

Abstract

Ontogenic resistance has been described for many plant-pathogen systems. Conversely, coffee leaf rust, a major fungal disease that drastically reduces coffee production, exhibits a form of ontogenic susceptibility, with a higher infection risk for mature leaves. To take into account stage-dependent crop response to phytopathogenic fungi, we developed an SEIR-U epidemiological model, where U stands for spores, which differentiates between young and mature leaves. Based on this model, we also explored the impact of ontogenic resistance on the sporulation rate. We computed the basic reproduction number $\mathcal {R}_0$, which classically determines the stability of the disease-free equilibrium. We identified forward and backward bifurcation cases. The backward bifurcation is generated by the high sporulation of young leaves compared to mature ones. In this case, when the basic reproduction number is less than one, the disease can persist. These results provide useful insights on the disease dynamics and its control. In particular, ontogenic resistance may require higher control efforts to eradicate the disease.

Appendices

Appendix A Bifurcation analysis

Theorem 6

(Castillo-Chavez and Song 2004) Consider the following ordinary differential equations, with a parameter $\psi$:

$$\begin{aligned} \frac{dx}{dt}= f(x, \psi ), \quad f:\mathbb {R}^n \times \mathbb {R} \rightarrow \mathbb {R}^n \quad \text {and} \quad f \in C^2(\mathbb {R}^n \times \mathbb {R}). \end{aligned}$$

(A1)

Without loss of generality, it is assumed that 0 is an equilibrium for system (A1) for all values of the parameter $\psi$, that is $f(0, \psi ) \equiv 0$ for all $\psi$. Assume

$A_1$::

$A = D_x f(0, 0) = (\frac{\partial f_i}{\partial x_j}(0,0))$ is the linearization matrix of system (A1) around the equilibrium 0 with $\psi$ evaluated at 0. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts;

$A_2$::

Matrix A has a nonnegative right eigenvector u and a left eigenvector v corresponding to the zero eigenvalue. Let $f_k$ be the $k^{th}$ component of f and

$$\begin{aligned} a=\sum _{k,i,j=1}^{n}v_k u_i u_j \frac{\partial ^2 f_k}{\partial x_i\partial x_j}(0,0) \quad \text {and} \quad b=\sum _{k,i=1}^{n}v_k u_i \frac{\partial ^2 f_k}{\partial x_i\partial \psi }(0,0). \end{aligned}$$

The local dynamics of (A1) around 0 are totally determined by a and b.

• 1. $a > 0$, $b > 0$. When $\psi < 0$ with $\Vert \psi \Vert \ll 1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when $0 < \psi \ll 1$, 0 is unstable and there exists locally asymptotically stable equilibrium;

• 2. $a < 0$, $b < 0$. When $\psi < 0$ with $\Vert \psi \Vert \ll 1$, 0 is unstable; when $0 < \psi \ll 1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium;

• 3. $a > 0$, $b < 0$. When $\psi < 0$ with $\Vert \psi \Vert \ll 1$, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0 < \psi \ll 1$, 0 is stable, and a positive unstable equilibrium appears;

• 4. $a < 0$, $b > 0$. When $\psi$ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

The linearisation of system (2) in the neighbourhood of DFE $x^{1}$ and around $\mathcal {R}_0=1$, using expression in Eq. (21) setting $\gamma _1$ to $\zeta$, yields

$$\begin{aligned} \frac{dx(t)}{dt}=J(x^1)x(t), \end{aligned}$$

(A2)

where

$$\begin{aligned} J(x^1)=\begin{pmatrix} -\frac{r_0 \beta }{\mu _M} &{}\quad 0 &{}\quad 0 &{}\quad \frac{\mu _M (\beta +\mu _J)}{\beta } &{}\quad \frac{\mu _M (\beta +\mu _J)}{\beta } &{}\quad \mathcal {Q} &{}\quad \mathcal {Q} &{}\quad - \frac{\omega _1\nu \mu _M}{\beta +\mu _M}\\ 0 &{}\quad -k_1 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \frac{\omega _1\nu \mu _M}{\beta +\mu _M} \\ 0 &{}\quad \theta &{}\quad -k_2 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0\\ \beta &{}\quad 0 &{}\quad 0 &{}\quad -\mu _M &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad - \frac{\omega _2\nu \beta }{\beta +\mu _M} \\ 0 &{}\quad \beta &{}\quad 0 &{}\quad 0 &{}\quad -k_3 &{}\quad 0 &{}\quad 0 &{}\quad \frac{\omega _2\nu \beta }{\beta +\mu _M} \\ 0 &{}\quad 0 &{}\quad \beta &{}\quad 0 &{}\quad \theta &{}\quad -k_4 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \alpha &{}\quad -k_5 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad \zeta &{}\quad 0 &{}\quad 0 &{}\quad \eta \zeta &{}\quad 0 &{}\quad -k_6 \end{pmatrix}, \end{aligned}$$

where $\mathcal {Q}=\frac{\mu _M (\beta +\mu _J)+r_0(\varepsilon -1)\beta }{\beta }$.

The matrix $J(x^1)$ admits $\rho =0$ as eigenvalue, the other eigenvalues still having a negative real part. Thereby, assumption $A_1$ of Theorem 6 is verified.

Let us now verify assumption $A_2$. We need to compute the left and right eigenvectors of matrix $J(x^1)$ associated to the eigenvalue $\rho =0$. The left eigenvector, denoted by $v=(v_1, v_2, v_3, v_4, v_5, v_6,v_7, v_8)$, satisfies $v J(x^1) =\textbf{0}$

Solving the above equation, yields

$$\begin{aligned} \left\{ \begin{aligned} v_1&=0,\\ v_2&= \frac{ \zeta \theta [\beta \eta (k_2+k_3)+k_3k_4]}{k_1k_2k_3k_4} v_8,\\ v_3&= \frac{\zeta (\beta \eta +k_4)}{k_2k_4} v_8,\\ v_4&=0,\\ v_5&=\frac{\eta \zeta \theta }{k_3k_4}v_8,\\ v_6&=\frac{\eta \zeta }{k_4}v_8,\\ v_7&=0,\; v_8>0. \end{aligned} \right. \end{aligned}$$

Similarly, the right eigenvector of matrix $J(x^1)$, denoted by $u=(u_1, u_2, u_3, u_4, u_5,u_6,u_7,u_8)^T$, satisfies $J(x^1)u =\textbf{0}$,

Solving the above equation, we obtain

$$\begin{aligned} \left\{ \begin{aligned} u_1&= \frac{u_8\nu \mu _M \mathcal {X}}{(\beta +\mu _M)((\beta +\mu _M)\mu _M-r_0\beta )k_1k_2k_3k_4k_5},\\ u_2&=\frac{ \omega _1\nu \mu _M}{(\beta +\mu _M)k_1} u_8,\\ u_3&= \frac{\theta \omega _1\nu \mu _M}{(\beta +\mu _M)k_1k_2} u_8,\\ u_4&= \frac{\beta (-\nu \omega _2u_8+(\beta +\mu _M)u_1)}{(\beta +\mu _M)\mu _M} \\ u_5&=\frac{\beta \nu (\mu _M\omega _1+\omega _2k_1)}{(\beta +\mu _M)k_1k_3} u_8,\\ u_6&= \frac{ \beta \nu \theta (\mu _M\omega _1(k_2+k_3)+\omega _2k_1k_2)}{(\beta +\mu _M)k_1k_2k_3k_4} u_8,\\ u_7&= \frac{ \nu \theta \alpha \beta (\mu _M\omega _1(k_2+k_3)+\omega _2k_1k_2)}{(\beta +\mu _M)k_1k_2k_3k_4k_5} u_8,\\ u_8&>0. \end{aligned} \right. \end{aligned}$$

where

$$\begin{aligned} \left\{ \begin{aligned} \mathcal {X}&=\mathcal {X}_1+\mathcal {X}_2+\mathcal {X}_3+\mathcal {X}_4,\\ \mathcal {X}_1&= (\beta +\mu _J)\omega _1((k_2+k_3)(\alpha +k_5)\theta +k_2k_4k_5)\mu _M^2,\\ \mathcal {X}_2&= ((r_0(\varepsilon -1)\omega _1+k_1\omega _2)k_2+\omega _1 k_3 r_0(\varepsilon -1))\beta +\mu _J\omega _2k_1k_2,\\ \mathcal {X}_3&= ((\alpha +k_5)\theta +k_1k_2k_4k_5(\beta \omega _2-k_3\omega _1+\mu _J\omega _2))\mu _M,\\ \mathcal {X}_4&= k_1\omega _2k_2(\beta r_0(\varepsilon -1)(\alpha +k_5)\theta -k_3k_4k_5(\beta +\mu _J)). \end{aligned} \right. \end{aligned}$$

The expression b defined in Assumption $A_2$ of Theorem 6 for system (2) is rewritten by

$$\begin{aligned} b =\sum _{k,i=1}^{8}v_k u_i \frac{\partial ^2 F_k}{\partial x_i\partial \gamma _1}(x^1, \zeta ), \end{aligned}$$

where $F_k=\dot{x}_k$ with $k=1,2,\ldots ,8$ corresponding to system (2).

The only term of partial derivative that is non-null corresponds to

$$\begin{aligned} v_8\frac{\partial ^2 F_8}{\partial x_3 \partial \gamma _1}(x^1, \zeta )=v_8 \; \text { and }\; v_8\frac{\partial ^2 F_8}{\partial x_6 \partial \gamma _1}(x^1, \zeta )=v_8 \eta . \end{aligned}$$

Substituting the above expression into the expression of b, yields

$$\begin{aligned} b&=v_8 u_3\frac{\partial ^2 F_8}{\partial x_3 \partial \gamma _1}(x^1, \zeta )+ v_8 u_6 \frac{\partial ^2 F_8}{\partial x_6 \partial \gamma _1}(x^1, \zeta )=v_8(u_3 + u_6 \eta ), \end{aligned}$$

Then

$$\begin{aligned} b=v_8(u_3 + u_6 \eta )>0 \end{aligned}$$

(A3)

Let us compute the expression of a, defined in Assumption $A_2$ of Theorem 6 for system (2), which is

$$\begin{aligned} a =\sum _{i,j,k=1}^{8}v_k u_i u_j \frac{\partial ^2 F_k}{\partial x_i\partial x_j}(x^1, \zeta ). \end{aligned}$$

The only terms that are non-null correspond to

$$\begin{aligned} v_2\frac{\partial ^2 F_2}{\partial x_1\partial x_8}(x^1, \zeta )&=v_2\frac{\partial ^2 F_2}{\partial x_8 \partial x_1}(x^1, \zeta ) = \frac{\omega _1\nu M_S^1}{(J_S^1+M_S^1)^2}v_2\\ v_2 \frac{\partial ^2 F_2}{\partial x_i\partial x_8}(x^1, \zeta )&=v_2\frac{\partial ^2 F_2}{\partial x_8 \partial x_j}(x^1, \zeta )= -\frac{\omega _1\nu J_S^1}{(J_S^1+M_S^1)^2}v_2\quad \text {for } i,j=2,\ldots ,7.\\ v_5 \frac{\partial ^2 F_5}{\partial x_4\partial x_8}(x^1, \zeta )&=v_5\frac{\partial ^2 F_5}{\partial x_8 \partial x_4}(x^1, \zeta )= \frac{\omega _2\nu J_S^1}{(J_S^1+M_S^1)^2}v_5\\ v_2 \frac{\partial ^2 F_5}{\partial x_i\partial x_8}(x^1, \zeta )&=v_5\frac{\partial ^2 F_5}{\partial x_8 \partial x_j}(x^1, \zeta )= -\frac{\omega _2\nu M_S^1}{(J_S^1+M_S^1)^2}v_5\quad \text {for } i,j=2,\ldots ,7. \end{aligned}$$

Substituting these above expressions into the expression of a, yields

$$\begin{aligned} a&=2 u_8\left[ v_2 \sum _{j=1}^{7} u_j\frac{\partial ^2 F_2}{\partial x_j\partial x_8}(x^1, \zeta )+v_5\sum _{j=1}^{7} u_j\frac{\partial ^2 F_5}{\partial x_j\partial x_8}(x^1, \zeta )\right] ,\\&= 2u_8v_2\left( \frac{\omega _1 \nu M_S^1}{(J_S^1+M_S^1)^2} u_1 -\frac{\omega _1 \nu J_S^1}{(J_S^1+M_S^1)^2}(u_2+u_3+u_4+u_5+u_6+u_7)\right) \\&\quad +2u_8v_5\left( \frac{\omega _2\nu J_S^1}{(J_S^1+M_S^1)^2} u_1 -\frac{\omega _2 \nu M_S^1}{(J_S^1+M_S^1)^2}(u_2+u_3+u_4+u_5+u_6+u_7)\right) . \end{aligned}$$

Replace the values of $J_S^1$, $M_S^1$, $v_2$ and $v_5$, we obtain

$$\begin{aligned} \begin{aligned} a&= \frac{2\zeta \mu _M r_1\omega _1 \nu \beta \theta u_8v_8(\mu _M\sum _{i=2}^{7}u_i-\beta u_1)(k_4k_3+\eta \beta (k_2+k_3))}{k_1k_2k_3k_4(\beta +\mu _M)^2((\beta +\mu _J)\mu _M-\beta r_0))}\\&\quad +\frac{2 r_1 \theta \mu _M \zeta \beta \nu \omega _2 u_8v_8\eta (\beta (u_1+u_2+u_3+u_5+u_6+u_7)-\mu _M u_4) }{(\beta +\mu _M)^2((\beta +\mu _J)\mu _M-\beta r_0)k_3k_4}.\\ \end{aligned} \end{aligned}$$

(A4)

The expression b given in Eq. (A3) is always positive, then the direction of bifurcation depends on the sign of a given by (A4). This means that we are in the case 1 and 4 of Theorem 6.

Appendix B Other types of bifurcation

Beside the forward and backward bifurcation diagrams presented in section 4. Figure 6 shows two different bifurcation diagrams, where the endemic equilibrium disappears for large $\mathcal {R}_0$ values while the trivial equilibrium becomes stable. Subplot (a) presents a forward bifurcation, as in Fig. 2; subplot (b) presents a backward bifurcation, as in Fig. 4. In both cases, these large $\mathcal {R}_0$ values lead to the plantation destruction. This phenomenon is commonly associated with frequency-dependence (Boots and Sasaki 2003; De Castro and Bolker 2005; Ryder et al. 2007; Best et al. 2012). Boots and Sasaki (2003) explored diseases transmitted via frequency-dependent methods, like sexually transmitted and vector-borne infections, which can evolve to push host populations towards extinction by maximizing their basic reproduction number. Authors proved that, independently from the form of the relationship between transmission and virulence, there is always the possibility of evolution to extinction when transmission is frequency-dependent. De Castro and Bolker (2005) conducted a review identifying the mechanisms of disease-induced extinction. Based on various studies, they concluded that frequency-dependence is a theoretical mechanism driving host extinction. Furthermore, their findings suggested that there is little empirical evidence indicating that diseases transmitted in a frequency-dependent manner pose a significant threat. Ryder et al. (2007) demonstrated that in host-parasite systems normally exhibiting density-dependent transmission, parasite-induced extinction could occur when adding a little frequency-dependent transmission. The findings of Best et al. (2012) underscore the crucial role of local processes in extinctions driven by parasites. They show that disease extinction may occur both for local density-dependent and frequency-dependent transmission, but that is more likely for the latter.

Nevertheless, we justified our frequency-dependent transmission by a dense plantation. In the process of plantation destruction, this high density assumption is not satisfied anymore, which limits the relevance of this result in the CLR case.

Fig. 6

Subplots present the ($\mathcal {R}_0, J_S$)-bifurcation diagram with stable (blue curve) and unstable (red curve) equilibria of system (1).The bifurcation parameter is $\gamma _1$, $\mathcal {R}_0 = \zeta \gamma _1$ and $\gamma _2 = \eta \gamma _1$, with $\zeta$ defined in Eq. (21). Subplot (a) shows a forward bifurcation using parameter values in Table 4 except $\beta =0.0003$/day and $\alpha =0.95$ leading to $\zeta =26.37$ spores/leaf.day. Subplot (b) shows the backward bifurcation using parameter values in Table 4 except $d=0.8961$ spores/leaf.day leading to $\zeta =7.77$ spores/leaf.day (colour figure online)

Appendix C Analysis of the model without stage structure

When we remove the stage structure, we obtain a system with five equations, with state variables (S, L, I, U, B):

$$\begin{aligned} \left\{ \begin{aligned}&\dot{S}= r_0 N_1-r_1 N S-\frac{\nu \omega U S}{N}-\mu S,\\&\dot{L}= \frac{\nu \omega U S}{N}-(\mu +\theta )L,\\&\dot{I}= \theta L-(\mu +\alpha )I,\\&\dot{R}=\alpha I-(\mu +d)R,\\&\dot{U}= \gamma I-(\nu +\mu _U)U. \end{aligned} \right. \end{aligned}$$

(C5)

System (C5) is defined in domain

$$\begin{aligned} \Gamma _{n}= \mathbb {R}_+^5 \setminus \{(0,0,0,0,U) \;\text {with} \;U\ge 0\}. \end{aligned}$$

The DFE of system (C5) is $E_{n} = (S_{n}, 0, 0, 0, 0)$, where $S_{n} =\frac{r_0-\mu }{r_1}$ exists when $r_0-\mu >0$. If the condition does not hold, the coffee plantation cannot survive.

The basic reproduction number of system (C5) is given by

$$\begin{aligned} \mathcal {R}_{0,n} =\frac{\nu \theta \gamma }{(\mu +\theta )(\mu +\alpha )(\nu +\mu _U)}. \end{aligned}$$

Lemma 7

According to Van den Driessche and Watmough (2002), the DFE $E_{n}$ of system (2) is locally asymptotically stable (LAS) in $\Gamma _{n}$ when $\mathcal {R}_{0,n}< 1$, and unstable if $\mathcal {R}_{0,n}>1$.

Theorem 8

If $\mathcal {R}_{0,n}<1$, the DFE $E_{n}$ is globally asymptotically stable in region $\Gamma _{n}$ and unstable otherwise.

Proof

Consider the cooperative linear subsystem

$$\begin{aligned} \left\{ \begin{aligned}&\dot{L}= \nu \omega U -(\mu +\theta )L,\\&\dot{I}= \theta L-(\mu +\alpha )I,\\&\dot{R}=\alpha I-(\mu +d)R,\\&\dot{U}= \gamma I-(\nu +\mu _U)U, \end{aligned} \right. \end{aligned}$$

(C6)

which is an upper-bound of system (C5) without the first equation, as $\frac{S}{N}\le 1$. One can easily show that (0, 0, 0, 0) is asymptotically stable if $\mathcal {R}_{0,n}<1$. Therefore, (0, 0, 0, 0) is GAS in the (L, I, R, U)-subsystem of system (C5).

Substituting $L = 0, I = 0, R=0$ and $U = 0$ into the first equation of system (C5) yields:

$$\begin{aligned} \dot{S}= (r_0-\mu )S\left( 1- \frac{r_1S}{(r_0-\mu )}\right) . \end{aligned}$$

The above logistic equation has two equilibria $S_0=0$ and $S_{n} =\frac{r_0-\mu }{r_1}$, which are unstable and GAS, respectively. Hence, since the solutions of system (C5) are bounded, $S(t) \longrightarrow S_{n}$. Thus, one can conclude that the DFE $E_{n}$ is GAS in region $\Gamma _{n}$ for system (C5). $\square$

This analysis shows that there cannot be a backward bifurcation in the system without stage structure (C5).