We therefore introduce the linear space $\mathbf{H} = (\mathbf{L}^2(\Omega))^5$ endowed with the usual inner product:

$$ \|(y, z)\|_{\mathbf{H}} = \sum_{i=1}^{5} (y_i, z_i)_{\mathbf{L}^2(\Omega)}, \quad y = (y_1, y_2, y_3, y_4, y_5), \quad z = (z_1, z_2, z_3, z_4, z_5); $$

$\mathbf{H}$ is a Hilbert space. Now define on $\mathbf{H}$ the linear operator $A$ such that:

$$ Ay = (\mu_p y_1, \mu_p y_2, \mu_p y_3, \mu_p y_4, -\varepsilon \Delta y_5), \quad (4.11) $$

with domain $D(A) = (\mathbf{L}^2(\Omega))^4 \times \mathbf{H}^2(\Omega)$, where

$$ \mathbf{H}_{\eta}^{2}(\Omega) := \{u \in \mathbf{H}^{1}(\Omega), \Delta u \in \mathbf{L}^{2}(\Omega), \frac{\partial U}{\partial \eta} = 0\} $$

**Remark 4.2.** For $u \in \mathbf{H}_{\eta}^2(\Omega)$ the boundary condition $\frac{\partial U}{\partial \eta} = 0$ implies that

$$ \forall v \in \mathbf{H}^1(\Omega), \quad \int_{\Omega} -\Delta u v \, dx = \int_{\Omega} \nabla u \cdot \nabla v \, dx. \quad (4.12) $$

We will need the following result.

**Lemma 4.3 ([121]).** Let $\varepsilon$ be a positive real number. Let $f \in L^2(\Omega)$ and $g$ the trace on $\partial\Omega$ of an element of $\mathbf{H}^1(\Omega)$. Assume that $\Omega$ is a bounded subset of $\mathbb{R}^N$, $N = 1, 2$ with smooth boundary and consider the following stationary boundary value problem:

$$ \left\{ \begin{array}{ll} u - \varepsilon \Delta u = f & \text{a.e. on } \Omega, \\ \frac{\partial u}{\partial \eta} = g & \text{a.e. on } \partial\Omega. \end{array} \right. \quad (4.13) $$

Then, the problem (4.13) has a unique solution $u \in \mathbf{H}^1(\Omega)$. Moreover $u$ belongs to $\mathbf{H}^2(\Omega)$ and there exists a universal constant $M > 0$ that only depends on $\varepsilon$, $N$ and $\Omega$ such that

$$ \|u\|_{H^1(\Omega)} \leq M(\|f\|_{L^2(\Omega)} + \|g\|_{L^2(\partial\Omega)}). $$

We have the following result about the existence of solution of subsystems (4.5).

**Lemma 4.4.** Assume that the initial conditions (4.3) holds, then there exists a unique local solution of problem (4.5) defined on $[0, \mathcal{T}_{\max}) \times \Omega$. More precisely:

$$ (S, L, I, J, U) \in C((0, \mathcal{T}_{\max}); D(A)) \cap C([0, \mathcal{T}_{\max}]; \mathbf{H}) \cap C^1((0, \mathcal{T}_{\max}]; \mathbf{H}); $$

Moreover,

$$ (S, L, I, J, U) \in [C([0, \mathcal{T}_{\max}) \times \bar{\Omega}) \cap C^{1,2}((0, \mathcal{T}_{\max}) \times \Omega)]^5 \text{ if } (\phi_S, \phi_L, \phi_I, \phi_J, \phi_U) \in (C(\Omega))^4 \times L^2(\Omega). $$

**Proof of Lemma 4.4.** System (4.5)-(4.3) can be written abstractly in the Hilbert space $\mathbf{H}$ in the following form:

$$ \left\{ \begin{array}{ll} y'(t) = Ay(t) + F(y(t)), & t \ge 0, \\ y(0) = y_0 \in D(A), & \end{array} \right. \quad (4.14) $$
where $y = (S, L, I, J, U)$, $y_0 = (\phi_S, \phi_L, \phi_I, \phi_J, \phi_U)$, $Ay$ is defined as in Equation (4.11) and:

$$F(y) = \left( \Lambda_p - \frac{\omega_p \nu U}{N} S, \frac{\omega_p \nu U}{N} S - \theta_p L, \theta_p L - (\alpha_p + d)I, \alpha_p I, \gamma_p I - (\nu + \mu_U)U \right).$$

Note that, since $\Omega$ is bounded, for $y = (S, L, I, J, U) \in \mathbf{H}$, one has $F(y) \in \mathbf{H}$. Now let us show that the linear operator $A$ is a maximal monotone operator on $\mathbf{H}$.

* $\forall v \in D(A), (Av, v) = \sum_{i=1}^{4} \|v_i\|_{L^2}^2 - \varepsilon \int_{\Omega} \Delta u_5 u_5 dx$. Since $u_5 \in H_{\eta}^2(\Omega)$, the identity (4.12) is satisfied with $v = u_5$ and one obtains that:

$$(Av, v) = \sum_{i=1}^{4} \|v_i\|_{L^2}^2 + \varepsilon \int_{\Omega} |\nabla u_5|^2 dx \geq 0,$$

which implies that $A$ is monotone.

* Now, let us show that $\forall v \in D(A), \exists u \in \mathbf{H}, u + Au = v$. Let $v = (v_1, v_2, v_3, v_4, v_5) \in D(A)$. We are looking for $u = (u_1, u_2, u_3, u_4, u_5) \in D(A)$, such that $u + Au = v$. With this is mind, one has:

$$u + Au = v \iff (1 + \mu_p)u_i = v_i, \quad i = 1, \dots, 4 \quad \text{and} \quad \begin{cases} u_5 - \varepsilon \Delta u_5 = v_5 & \text{in } \Omega, \\ u_5 \in H_{\eta}^2(\Omega). \end{cases}$$

Therefore, for $i = 1, \dots, 4$, $u_i = \frac{1}{1+\mu_p}v_i \in L^2(\Omega)$ and one obtains $u_5$ by solving the following elliptic problem:

$$\left\{ \begin{array}{ll} u_5 - \varepsilon \Delta u_5 = v_5 & \text{in } \Omega, \\ \frac{\partial u_5}{\partial \nu} = 0 & \text{on } \partial \Omega. \end{array} \right. \qquad (4.15)$$

According to Lemma 4.3, there exists a unique $u_5 \in H_{\eta}^2(\Omega)$ solution of (4.15) and $A$ is maximal. By the Hille-Yosida theorem [121], we conclude that $-A$ is the infinitesimal generator of a $C_0$ semigroup of contractions on $\mathbf{H}$.

* We now show that $F$ is Lipschitz continuous in both variables. Let $y = (S, L, I, J, U)$, $\tilde{y} = (\tilde{S}, \tilde{L}, \tilde{I}, \tilde{J}, \tilde{U}) \in \mathbf{H}$. Therefore, one has:

$$F(y) - F(\tilde{y}) = \left( \omega_p \nu \left( -\frac{U}{N} S + \frac{\tilde{U}}{\tilde{N}} \tilde{S} \right), \omega_p \nu \left( \frac{U}{N} S - \frac{\tilde{U}}{\tilde{N}} \tilde{S} \right) - \theta_p (L - \tilde{L}), \right. \\ \left. \theta_p (L - \tilde{L}) - (\alpha_p + d)(I - \tilde{I}), \alpha_p (I - \tilde{I}), \gamma_p (I - \tilde{I}) - (\nu + \mu_U)(U - \tilde{U}) \right).$$

Since $\dot{N} = \Lambda_p - \mu_p N - dI \geq \Lambda_p - \mu_p N - dN$, we have $N(t, x) \geq \min(\frac{\Lambda_p}{\mu_p+d}, \tilde{\phi}(x))$. Moreover, according to Lemma 4.2, $\tilde{U}$ is bounded. This implies that $\frac{\tilde{U}}{N}$ is bounded. As $\frac{S}{N} \leq 1$, from the identity $\frac{U}{N}S - \frac{\tilde{U}}{\tilde{N}}\tilde{S} = \frac{U-\tilde{U}}{N}S - \frac{\tilde{U}}{\tilde{N}}(\tilde{S}-S) + \frac{\tilde{U}S}{N\tilde{N}}(\tilde{N}-N)$, we obtain:

$$\|F(y) - F(\tilde{y})\|_{\mathbf{H}} \leq K \|y - \tilde{y}\|_{\mathbf{H}},$$

where $K$ is a constant that depends on the constant $\omega_p, \nu, \theta_p, \alpha_p, d, \gamma_p$ and $\mu_U$. Thus $F$ is uniformly Lipschitz continuous on $\mathbf{H}$. Using the fact that the solution are positive, we can
now apply Theorem 1.6 page 189 of [121] and conclude that system (4.14) has a unique local strong solution on $[0, T_{\max})$ in the sense that $(S, L, I, J, U) \in C((0, T_{\max}); D(A)) \cap C([0, T_{\max}); H) \cap C^1((0, T_{\max}); H)$, we can also use the method present in [179] to conclude prove. This completes the proof.

□

**Theorem 4.1.** *Subsystems (4.5), with initial conditions in $H$ that verify $0 < \phi_S + \phi_L + \phi_I + \phi_J \le N_m$, $0 \le \phi_U \le U_m$ for all $x \in \bar{\Omega}$, admit a unique global solution in $C((0, \infty); D(A)) \cap C([0, \infty); H) \cap C^1((0, \infty); H)$ such that, for all $(t, x) \in \mathbb{R}_+ \times \Omega$*

$$S(t, x) > 0, \quad L(t, x) > 0, \quad I(t, x) > 0, \quad J(t, x) > 0, \quad 0 < N(t, x) \le N_m, \quad 0 < U(t, x) \le U_m$$

$$\text{where } N_m = \max\left\{\frac{\Lambda_p}{\mu_p}, \|\phi\|_{\infty}\right\} \text{ and } U_m = \max\left\{\max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\}\frac{\gamma_p}{(\nu+\mu_U)}, \|\phi_U\|_{\infty}\right\}$$

*Proof of Theorem 4.1.* Using Lemma 4.1, Lemma 4.2 and Lemma 4.4 we conclude the existence and uniqueness of bounded positive solution of subsystems (4.5) defined on $[0, \infty) \times \Omega$. we can conclude the proof of this theorem

□

### 4.3.2 Equilibria and their stability

Here, we compute the equilibria of subsystems (4.5) with $p = D$ or $R$ and study their stability. Let

$$\mathcal{R}_0^{(p)} = \frac{\gamma_p}{(\nu + \mu_U)} \frac{\nu\omega_p}{(\theta_p + \mu_p)} \frac{\theta_p}{(\alpha_p + \mu_p + d)}$$

The expression $\mathcal{R}_0^{(p)}$ is derived in (4.38) in the proof of the local stability of the disease free equilibrium. It corresponds to the basic reproduction number in the dry ($p = D$) or rainy ($p = R$) season, since (i) $\frac{\gamma_p}{(\nu+\mu_U)}$ represents the mean number of new Uredospores generated by a single infectious branch and (ii) $\frac{\nu\omega_p}{(\theta_p+\mu_p)} \frac{\theta_p}{(\alpha_p+\mu_p+d)}$ measures the average number of new infectious branches generated by a single uredospore introduced in a completely susceptible field. We prove the following results for the stability of equilibria of subsystems (4.5).

**Lemma 4.5.**

1. If $\mathcal{R}_0^{(p)} < 1$ the disease-free equilibrium (DFE) $Q_p^0 = (\frac{\Lambda_p}{\mu_p}, 0, 0, 0, 0)$ is locally asymptotically stable (LAS).

2. If $\mathcal{R}_0^{(p)} > 1$ and close to one, the DFE is unstable and there exists a unique positive endemic equilibrium $Q_p^* = (S_p^*, L_p^*, I_p^*, J_p^*, U_p^*)$, where

$$I_p^* = \frac{\Lambda_p(\nu + \mu_U)(\mathcal{R}_0^{(p)} - 1)}{\gamma_p\omega_p\nu - d(\nu + \mu_U)}, \quad S_p^* = \frac{\Lambda_p - dI_p^*}{\mu_p\mathcal{R}_0^{(p)}}, \qquad (4.16)$$

$$L_p^* = \frac{(\alpha_p + \mu_p + d)}{\theta_p} I_p^*, \quad J_p^* = \frac{\alpha_p}{\mu_p} I_p^* \quad \text{and} \quad U_p^* = \frac{\gamma_p}{(\nu + \mu_U)} I_p^*$$

which is locally asymptotically stable.

*Proof of Lemma 4.5.* The proof is given in Appendix (Section 4.6).

□
Using Lyapunov theory and LaSalle's principle [88], we prove the global stability of the DFE, which implies that the CLR will dwindle until extinction, whatever the initial number of Uredospores and infectious branches.

**Theorem 4.2.** If $\mathcal{R}_0^{(p)} < 1$, the DFE $Q_p^0$ of subsystems (4.5) is globally asymptotically stable (GAS).

*Proof of Theorem 4.2.* Consider the function

$$V_1 = a_1 L + b_1 I + c_1 J + d_1 U. \qquad (4.17)$$

where $a_1, b_1, c_1$ and $d_1$ are positive constants to be chosen later. Consider then the following Lyapunov function candidate:

$$W_1 = \int_{\Omega} V_1 dx = \int_{\Omega} (a_1 L + b_1 I + c_1 J + d_1 U) dx, \qquad (4.18)$$

where $\Omega$ is the spatial domain defined above.

Let us consider the set $\mathcal{V} = C((0, \infty); (\mathbf{L}^2(\Omega))^3 \times \mathbf{H}_\eta^2(\Omega)) \cap C([0, \infty); (\mathbf{L}^2(\Omega))^4) \cap C^1((0, \infty); (\mathbf{L}^2(\Omega))^4)$. For any solution $(L, I, J, U) \in \mathcal{V}$ of subsystems (4.5) with positive initial condition $(\phi_L(x), \dots, \phi_U(x))$, $W_1(L, \dots, U)$ is positive. Also, $W_1 = 0$ if and only if $(L, I, J, U) = (0, 0, 0, 0)$ on $\Omega$.

The time derivative of the Lyapunov function $W_1$ along the trajectories of subsystems (4.5) satisfies:

$$\frac{dW_1}{dt} = \int_{\Omega} \dot{V}_1(L, I, J, U) dx. \qquad (4.19)$$

From equation (4.17), one has that

$$\begin{aligned} \dot{V}_1 &= a_1 \frac{\partial L}{\partial t} + b_1 \frac{\partial I}{\partial t} + c_1 \frac{\partial J}{\partial t} + d_1 \frac{\partial U}{\partial t}, \\ &= a_1 \left( \frac{\omega_p \nu U}{N} S - (\theta_p + \mu_p)L \right) + b_1 (\theta_p L - (\alpha_p + \mu_p + d)I) + c_1 (\alpha_p I - \mu_M J) \qquad (4.20) \\ &\quad + d_1 (\gamma_p I - (\nu + \mu_U)U) + d_1 \varepsilon \Delta U. \end{aligned}$$

Now, using the fact that $\frac{S}{N} \le 1$, equation (4.20) becomes:

$$\dot{V}_1 \le (-a_1(\theta_p+\mu_p)+b_1\theta_p)L+(-b_1(\alpha_p+\mu_p+d)+c_1\alpha_p+d_1\gamma_p)I+(-d(\nu+\mu_U)+a_1\omega_p\nu)U-c_1\mu_pJ+d_1\varepsilon\Delta U. \qquad (4.21)$$

Then, the positive constants $a_1, b_1, c_1$ and $d_1$ are chosen such that:

$$\left\{ \begin{array}{l} -a_1(\theta_p + \mu_p) + b_1\theta_p = 0, \\ -b_1(\alpha_p + \mu_p + d) + c_1\alpha_p + d_1\gamma_p = 0, \\ -d_1(\nu + \mu_U) + a_1\omega_p\nu = 0. \end{array} \right. \qquad (4.22)$$

Solving the above system (4.22) yields:

$$a_1 = \theta_p, \quad b_1 = \theta_p + \mu_p, \quad d_1 = \frac{\theta_p \omega_p \nu}{\nu + \mu_U} \quad \text{and} \quad c_1 = \frac{(\theta_p + \mu_p)(\alpha_p + \mu_p + d)}{\alpha_p} (1 - \mathcal{R}_0^{(p)}).$$

Note that $\mathcal{R}_0^{(p)} < 1$ ensures that $c_1 > 0$. Replacing the above expressions of $a_1, b_1, c_1$ and $d_1$
into equation (4.21), yields:

$$ \dot{V}_1 \le -c_1\mu_p J + d_1\varepsilon\Delta U. \qquad (4.23) $$

Then equation (4.19) becomes:

$$ \frac{dW_1}{dt} \le -\int_{\Omega} c_1\mu_p J dx + d_1\varepsilon \int_{\Omega} \Delta U dx. $$

From the Neumann bounded condition $\frac{\partial U}{\partial \eta} = 0$, one has that:

$$ \int_{\Omega} \Delta U dx = \int_{\partial \Omega} \frac{\partial U}{\partial \eta} = 0. \qquad (4.24) $$

Thus, if $\mathcal{R}_0^{(p)} < 1$ one has:

$$ \frac{dW_1}{dt} \le -\int_{\Omega} c_1\mu_p J dx. \qquad (4.25) $$

Thus, $\mathcal{R}_0^{(p)} < 1$ ensures that $\frac{dW_1}{dt} \le 0$ for all $L, I, J, U \ge 0$.

The domain $\Gamma_\rho = \{(L, I, J, U) \in \mathcal{V} : W_1(L, I, J, U) \le \rho\}, \rho > 0$, is compact and includes the origin. Moreover, it is positively $\mathcal{V}$-invariant with respect the last four equations of subsystems (4.5). In fact, if the initial condition $(\phi_L, \phi_I, \phi_J, \phi_U) \in \Gamma_\rho$, we have

$$ \frac{dW_1(L, I, J, U)}{dt} \le 0 \Rightarrow W_1(L, I, J, U) \le W_1(\phi_L, \phi_I, \phi_J, \phi_U) \le \rho $$

Also, $\frac{dW_1}{dt} = 0$ if and only if $J = 0$, which from subsystems (4.5) is satisfied over a time interval of non-zero length only if $L = 0, I = 0$ and $U = 0$. Then it is easy to see that the largest invariant set defined over $\mathcal{E} = \{(L, I, J, U) \in \mathcal{V} : \frac{dW_1(L, I, J, U)}{dt} = 0\}$ is the singleton $\{(0, 0, 0, 0)\}$. Therefore, using LaSalle-Krasowski invariant principle [[54], Theorem 2] and the LAS in Lemma 4.5, one can conclude that $(L, I, J, U) = (0, 0, 0, 0)$ is GAS when $\mathcal{R}_0^{(p)} < 1$. Similar proofs adapting Lasalle's principle for PDE can be found in [152, 174].

Replacing of the value $U = 0$ into the first equation of subsystems (4.5) yields:

$$ \dot{S} = \Lambda_p - \mu_p S. \qquad (4.26) $$

The above equation (4.26) has a unique equilibrium $S_p^0 = \frac{\Lambda_p}{\mu_p}$, which is GAS. Hence, since the solutions of subsystems (4.5) are bounded, $S(x, t) \to S_p^0$. Thus, one can conclude that the DFE $Q_p^0 = (S_p^0, 0, 0, 0, 0)$ is globally asymptotically stable for subsystems (4.5). This concludes the proof. $\square$

### 4.3.3 Comparison of the dynamics in the dry and rainy seasons

The objective of this section is to compare the dynamics of subsystems (4.5) during the rainy and dry seasons. To this end, we first compare the basic reproduction numbers $\mathcal{R}_0^{(D)}$ and $\mathcal{R}_0^{(R)}$, and the number of infectious branches at the endemic equilibrium $I_R^*$ and $I_D^*$. The former will indicate if and when CLR can persist within the coffee plantation, and the latter the intensity of the infection.

Simple biological hypotheses give the relation between the parameter values during the dry and rainy seasons. They mainly rely on the fact that some mechanisms require humid conditions:
coffee tree growth is faster during the rainy season ($\Lambda_D \le \Lambda_R$), germination occurs mostly during the rainy season ($\omega_D \le \omega_R$), the lesion progresses faster during the rainy season ($\theta_D \le \theta_R$) and the mortality rate of the branches is higher during the dry season (as a consequence of harvest, which always damages some branches) i.e ($\mu_R \le \mu_D$). Hence, we have:

$$ \Lambda_D \le \Lambda_R, \quad \omega_D \le \omega_R, \quad \theta_D \le \theta_R, \quad \mu_R \le \mu_D, \quad \gamma_D = \gamma_R \quad \text{and} \quad \alpha_D = \alpha_R. \qquad (4.27) $$

### 4.3.3.1 Comparison of the basic reproduction numbers

Herein we compare the basic reproduction numbers of the subsystems:

$$ \mathcal{R}_0^{(R)} = \frac{\gamma_R}{(\nu + \mu_U)} \frac{\nu \theta_R \omega_R}{(\theta_R + \mu_R)(\alpha_R + \mu_R + d)} \quad \text{and} \quad \mathcal{R}_0^{(D)} = \frac{\gamma_D}{(\nu + \mu_U)} \frac{\nu \theta_D \omega_D}{(\theta_D + \mu_D)(\alpha_D + \mu_D + d)}. $$

Using (4.27), one has:

$$ \mathcal{R}_0^{(D)} \le \mathcal{R}_0^{(R)}. \qquad (4.28) $$

This means that the GAS of the DFE $Q_R^0$ during the rainy season implies the GAS of the DFE $Q_D^0$ during the dry season and the existence and local asymptotic stability of the endemic equilibrium $Q_D^*$ during the dry season implies the existence and local asymptotic stability of the endemic equilibrium $Q_R^*$ during the rainy season. Biologically speaking, if the disease is present during the dry season, it will be also present during the rainy season while, if the disease is absent during the rainy season, it will be also absent during the dry season.

### 4.3.3.2 Comparison of the infectious branches at the endemic level

To compare the disease at the endemic level, both basic reproduction numbers need to be greater than one, $1 < \mathcal{R}_0^{(p)}$ for $p = D$ and $R$. The expression of $I_p^*$ depends on the parameters identified in (4.27). The partial derivatives of $I_p^*$ with respect to parameters $\Lambda_p, \omega_p, \theta_p$ and $\mu_p$ are:

$$ \begin{aligned} \frac{\partial I_p^*}{\partial \Lambda_p} &= \frac{(\nu + \mu_U)(\mathcal{R}_0^{(p)} - 1)}{\gamma_p \omega_p \nu - d(\nu + \mu_U)} > 0, & \frac{\partial I_p^*}{\partial \omega_p} &= \frac{\gamma_p \nu (\theta_p (\alpha_p + \mu_p) + \mu_p (\alpha_p + \mu_p + d))}{[\gamma_p \omega_p \nu - d(\nu + \mu_U)] (\theta_p + \mu_p) (\alpha_p + \mu_p + d)} > 0, \\ \frac{\partial I_p^*}{\partial \theta_p} &= \frac{\gamma_p \nu \omega_p \mu_p}{(\omega_p + \mu_p)} > 0, & \frac{\partial I_p^*}{\partial \mu_p} &= -\frac{\gamma_p \nu \omega_p \theta_p (\theta_p + 2\mu_p + \alpha_p d)}{[\gamma_p \omega_p \nu - d(\nu + \mu_U)] (\theta_p + \mu_p) (\alpha_p + \mu_p + d)} < 0. \end{aligned} $$

This shows that $I_p^*$ is an increasing function if $\Lambda_p, \omega_p, \theta_p$ and a decreasing function of $\mu_p$. Now, using (4.27), we can deduce that:

$$ I_R^* > I_D^*. \qquad (4.29) $$

Therefore, if it reached equilibrium, the CLR would be more severe during the rainy season than during the dry season.

## 4.4 Numerical simulations

Herein, we present the results of numerical simulations of system (4.1) using a non standard finite difference method [108] presented in Appendix (Section 4.7). We take $T = 365$ days and $\tau = 120$ days which correspond to the durations of the year and the dry season in Cameroon, respectively. We set $\Omega = [0, 100]$ meters. We consider that initially all branches are healthy, that
<figure><img src="image_18.png" /><figcaption>Figure 4.2: Spatio-temporal simulation of system (4.1,4.4,4.30) when $\gamma_D = \gamma_R = 1.5$, which leads to $\mathcal{R}_0^{(D)} = 0.52$ and $\mathcal{R}_0^{(R)} = 0.95$. All other parameter values are given in Table 4.1. Subplots represent healthy branches $S$, latent branches $L$, infectious branches $I$, leafless branches $J$, Uredospores $U$ and berries $B$.</figcaption></figure>

Uredospores are concentrated in the middle of the plantation and that there are no berries, so the initial conditions are, $\forall x \in \Omega$:

$$ \left\{ \begin{array}{l} S_0(x) = 100, \quad L_0(x) = I_0(x) = J_0(x) = 0, \\ U_0(x) = 2000 \sin^{14}\left(\frac{\pi x}{100}\right) \\ B_0(x) = 0. \end{array} \right. \qquad (4.30) $$

Table 4.1 summarizes the parameter values for system (4.1) used for numerical simulations. We suppose that $\gamma_p$ and $\alpha_p$ have the same value during the dry and rainy seasons.

### 4.4.1 CLR dynamics

We first choose $\gamma_p = 1.5$ for $p = D$ and $R$, so that $\mathcal{R}_0^{(D)} = 0.52$ and $\mathcal{R}_0^{(R)} = 0.95$. Figure 4.2 presents the spatio-temporal evolution of system (4.1) with these parameters in domain $\Omega$ during time interval $[0, 1500]$ days. Uredospores are initially present in the middle of the $\Omega$ interval, which induces local CLR infection during the first year. Subsequently, the infection dies out, as $L, I, J$ and $U$ go to 0. The $S$ and $B$ variables regain their healthy levels, corresponding to the stationary periodic disease-free solution. This illustrates Theorem 4.2: since $\mathcal{R}_0^{(D)}$ and $\mathcal{R}_0^{(R)}$ are smaller than one, the disease should vanish during both the dry and rainy seasons; we see here that connecting both seasons preserves this property.

Then we choose $\gamma_p = 8$ for $p = D$ and $R$, so that $\mathcal{R}_0^{(R)} = 5.09$ and $\mathcal{R}_0^{(D)} = 2.77$. In this
case, we have shown in Lemma 4.5 that the endemic equilibria exist and is locally asymptotically
stable during both seasons. Numerical results are shown in Figure 4.3. One can observe that
the peaks of the infectious branches and Uredospores flatten very quickly compared to the other
variables. The convergence towards the endemic stationary solution is clear at the beginning
of the third year (after 1000 days approximately). Oscillations are observed, as the endemic
equilibrium is higher during the rainy season (subsystem (4.5) with $p = R$) than the dry season
(subsystem (4.5) with $p = D$). Note that Uredospores hardly reach the edges of the plantation
during the first year, so that healthy branches are produced in large quantity, contrary to the
middle of the plantation where Uredospores directly infect branches, which limits growth. Hence,
when the Uredospores reach the edges during the second year, there is a large infection peak at
the edges, which contrasts with the middle of the plantation.

Figure 4.4 is a temporal representation of Figure 4.3 for three different values of *x* taken in the centre of the domain (*x* = 50m, yellow curve), on the border (*x* = 1m, blue curve) and in-between (*x* = 25m, red curve). Convergence towards the stationary endemic solution can also be observed here, as all curves overlap for the three *x* values at the end of the third year and on. At the centre of the domain, the system very rapidly converges towards this stationary solution, but it takes a year for *x* = 25m. On the border of the domain, branches remain susceptible during almost two years. Then the number of latent *L* and infectious *I* branches peak, as Uredospores *U*, initially in the centre of the domain, reach its borders. These peaks occur during the second rainy season. However, no such peak is observed for leafless branches *J*. Indeed, the number of leafless branches starts to increase at the end of the second rainy season (around time *t* = 600 days), but the arrival of the dry season reduces the number of branches before it can peak very high. Solutions then oscillate between the endemic equilibria of the dry and rainy seasons (purple and green dashed lines, respectively). Infection has a negative impact on the production, which is reduced to less than a third of its value without disease.

These simulations allow to conclude that the dynamics of system (4.1), can be inferred from
the dynamics of subsystems (4.5) during the dry and rainy seasons.

4.4.2 Control of CLR

Integrated pest management currently used for the control of CLR [3] relies on the following
methods:

* cultural practices, including stumping of diseased plants, shading, spore traps [5];

• chemical control by copper, triazole or dithiocarbamate fungicides [25];

• genetic control based on wild rust-resistant *Coffea* species [5];

* biological control based on *Hemileia vastatrix* natural enemies, e.g. mycoparasites such as *Lecanicillium lecanii* [26, 168], or antagonistic bacteria such as *Pseudomonas putida* P286 and *Bacillus thuringiensis* B157 [61, 62]; it can also rely on natural allies of the crop such as benefic endophytes [147, 146].

Biological control is a potentially powerful tool for managing coffee leaf rust that allows organic certification. Therefore, we introduced in this study a mycoparasite (a parasitic fungus whose host is another fungus), *Lecanicillium lecanii*, known to hamper the reproduction of *H. vastatrix*. In our model, this biocontrol agent reduces the production of Uredospores by infectious
<figure><img src="image_19.png" /><figcaption>Figure 4.3: Spatio-temporal simulation of system (4.1,4.4,4.30) when $\gamma_D = \gamma_R = 8$, which leads to $\mathcal{R}_0^{(R)} = 5.09$ and $\mathcal{R}_0^{(D)} = 2.77$. All other parameter values are given in Table 4.1. Subplots represent healthy branches $S$, latent branches $L$, infectious branches $I$, leafless branches $J$, Uredospores $U$ and berries $B$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 4.4: Temporal simulation, for three values of $x$, of system (4.1,4.4,4.30) when $\gamma_D = \gamma_R = 8$, which leads to $\mathcal{R}_0^{(R)} = 5.09$ and $\mathcal{R}_0^{(D)} = 2.77$ (as in Figure 4.3). All other parameter values are given in Table 4.1. $x$ is chosen in the centre of the domain (yellow curve), on the border (blue curve) and in-between (red curve). Subplots represent healthy branches $S$, latent branches $L$, infectious branches $I$, leafless branches $J$, Uredospores $U$ and berries $B$.</figcaption></figure>
branches with an efficiency $q$. The production rate of Uredospores $\gamma_p$ hence becomes $(1 - q)\gamma_p$ and system (4.5) becomes:

$$ 
\left\{ 
\begin{aligned}
\partial_t S &= \Lambda_p - \frac{\omega_p \nu U}{N} S - \mu_p S, \\
\partial_t L &= \frac{\omega_p \nu U}{N} S - (\theta_p + \mu_p)L, \\
\partial_t I &= \theta_p L - (\alpha_p + \mu_p + d)I, \\
\partial_t J &= \alpha_p I - \mu_p J, \\
\partial_t U &= \varepsilon \Delta U + (1 - q)\gamma_p I - (\nu + \mu_U)U.
\end{aligned}
\right. 
\qquad (4.31) 
$$

The basic reproduction numbers of system (4.31), for $p = D$ and $R$, are then:

$$ \mathcal{R}_0^{(p)} = \frac{(1-q)\gamma_p \nu}{(\nu + \mu_U)} \frac{\theta_p \omega_p}{(\theta_p + \mu_P)(\alpha_p + \mu_p + d)}. \qquad (4.32) $$

The higher the biocontrol efficiency $q$, the lower the reproduction numbers. So high enough values of $q$ should allow to control the disease.

We simulated system (4.31) for different values of biocontrol efficiency $q$ over 8 years. Figure 4.5 presents the temporal evolution of the state variables ($S, L, I, J, U, B$) integrated over domain $\Omega$. Table 4.2 provides the reproduction numbers defined in equation (4.32), as well as the berry production and production loss during the 8th year, the latter being computed by reference to the disease-free case $q = 1$.

* $q = 0$ corresponds to no control (red curves) and $\mathcal{R}_0^{(R)} > \mathcal{R}_0^{(D)} > 1$, as in previous Section 4.4.1. The yield loss is as high as 80.1%.

* For a 50% efficiency (cyan curves), which still corresponds to $\mathcal{R}_0^{(R)} > \mathcal{R}_0^{(D)} > 1$, the number of spores $U$ drops as expected, as the mycoparasite hampers the production of spores. So the number of healthy branches $S$ increases. However, the number of latent $L$, infectious $I$ and leafless $J$ branches, barely change. Indeed, the production of latent branches depends on the product $SU$. The production of berries $B$ then increases, as it is mostly driven by the healthy branches, but the yield loss is still high at 68.2%.

* For a 70% efficiency (magenta curves), corresponding to $\mathcal{R}_0^{(R)} > 1 > \mathcal{R}_0^{(D)}$, these observations are amplified: spores seriously drop; healthy branches and berries largely increase; latent, infectious and leafless branches decrease. The yield loss is 33.5%.

* For a 75% efficiency (black curves), also corresponding to $\mathcal{R}_0^{(R)} > 1 > \mathcal{R}_0^{(D)}$, infection remains very low, leading to an acceptable 6.0% yield loss.

* An 76% efficiency (blue curves), still corresponding to $\mathcal{R}_0^{(R)} > 1 > \mathcal{R}_0^{(D)}$, almost yields the same results as a perfect efficiency ($q = 1$). Infection is almost negligible, so the yield loss is as low as 1.8%.

* $q = 1$ (green curves) corresponds to a perfect (and unrealistic) efficiency, with $\mathcal{R}_0^{(R)} = \mathcal{R}_0^{(D)} = 0$, leading to the disease extinction and no yield loss.

Hence a 75% biocontrol efficiency is enough to sustain the berry production in the plantation, with a negligible yield loss. Higher efficiencies, moreover, achieve disease eradication.
<figure><img src="image_21.png" /><figcaption>Figure 4.5: Temporal simulations, for various efficiencies $q$ of the biocontrol mycoparasite, of controlled system (4.31,4.4,4.30) when $\gamma_D = \gamma_R = 8$. All other parameter values are given in Table 4.1. Subplots represent healthy branches $S$, latent branches $L$, infectious branches $I$, leafless branches $J$, Uredospores $U$ and berries $B$, these variables being integrated over domain $\Omega$.</figcaption></figure>

<table><caption>Table 4.2: Impact of biocontrol efficiency $q$ on reproduction numbers, defined in (4.32), and on berry production during the 8th and last year.</caption><thead><tr><th>Biocontrol efficiency</th><th>q = 0</th><th>q = 0.5</th><th>q = 0.7</th><th>q = 0.75</th><th>q = 0.76</th><th>q = 1</th></tr></thead><tbody><tr><td>R<sub>0</sub><sup>(D)</sup></td><td>2.77</td><td>1.38</td><td>0.83</td><td>0.69</td><td>0.66</td><td>0</td></tr><tr><td>R<sub>0</sub><sup>(R)</sup></td><td>5.09</td><td>2.54</td><td>1.52</td><td>1.27</td><td>1.22</td><td>0</td></tr><tr><td>Number of berries ×10<sup>5</sup></td><td>33.2</td><td>53.1</td><td>111</td><td>157</td><td>164</td><td>167</td></tr><tr><td>Yield loss (%)</td><td>80.1</td><td>68.2</td><td>33.5</td><td>6.0</td><td>1.8</td><td>0</td></tr></tbody></table>
## 4.5 Conclusion

In this paper, we have proposed and analysed a PDE model that describes the dispersal of CLR in a coffee plantation during the rainy and dry seasons and its behaviour over time. Furthermore, we computed the disease-free and endemic equilibria of the two subsystems defined during the rainy and dry seasons. We showed that the basic reproduction numbers during the two seasons can determine the dynamics of global model: when the basic reproduction number is less than one during the rainy season, then CLR globally decreases till extinction; when it is greater than one for the dry season, then CLR persists.

We implemented a biocontrol in our model, corresponding to a mycoparasite such as *Lecani-cillium lecanii*, which hampers CLR reproduction at all times. This solution was tested in Mexico [26] but is still under development. A rather high biocontrol efficiency (75% at least) is necessary in our model to control the disease, but lower efficiencies still improve coffee production notably. Moreover, the mycoparasite is applied all year round, so it is not easily implemented in practice and it involves important costs. It would be interesting to study when to deploy the mycopara- site in a cost-efficient way. An ideal mycoparasite should sustain the dry season and efficiently control CLR, so that the coffee plantation during the rainy season would suffer reasonable yield losses. In further work, we will also include cultural management and other biocontrol agents, in particular natural endophytes which are affordable for growers starting a new plantation.

Several extensions to this work are considered: (i) adding a stage structure on the coffee branches; (ii) simplifying the model, using an impulsive formalism for the dry season as in [96, 182], in order to obtain analytical results on the global model behaviour; (iii) solving an optimal control problem, consisting in maximising coffee production while minimising the control costs.

## 4.6 Appendix: Mathematical analysis

### Proof of Lemma 4.5

1. Local stability of the disease-free equilibrium

Let $(S, L, I, J, U)$ be a solution of subsystems (4.5). Then, according to Kiehöfer [81], this solution can be written in the following form:

$$ (S, L, I, R, U) = Q_p^0 + (W_1, W_2, W_3, W_4, W_5), \quad x \in \Omega, \quad t > 0. \quad (4.33) $$

With this in mind, subsystems (4.5) can be written in the following compact form:

$$ \frac{\partial W}{\partial t} = \mathcal{D}\Delta W + F(W), \quad (4.34) $$

where $\mathcal{D} = \text{diag}(0, 0, 0, 0, 0, \varepsilon)$. The linearisation of subsystems (4.34) in the neighbourhood of $Q_p^0$ is:

$$ \frac{\partial W}{\partial t} = \mathcal{D}\Delta W + L(Q_p^0)W, \quad (4.35) $$

where $L(Q_p^0)$ is the Jacobian matrix at the DFE $Q_p^0$ of subsystems (4.5) in the absence of
diffusion, that is:

$$L(Q_p^0) = \begin{pmatrix} -\mu_p & 0 & 0 & 0 & -\omega_p\nu \\ 0 & -(\theta_p + \mu_p) & 0 & 0 & \omega_p\nu \\ 0 & \theta_p & -(\alpha_p + \mu_p + d) & 0 & 0 \\ 0 & 0 & \alpha_p & -\mu_p & 0 \\ 0 & 0 & \gamma_p & 0 & -(\nu + \mu_U) \end{pmatrix}.$$

Let $g_j, j \in \mathbb{N}$, be the $j^{th}$ eigenfunction of operator $-\Delta$ with Neumann boundary conditions, so that:

$$\left\{ \begin{array}{ll} -\Delta g_j = \lambda_j g_j & \text{in } \Omega, \\ \frac{\partial g_j}{\partial \eta} \bigg|_{\partial \Omega} = 0 & \text{on } \partial \Omega, \end{array} \right. \qquad (4.36)$$

where $\lambda_j$ are the associated eigenvalues verifying $0 = \lambda_0 < \lambda_1 < \lambda_2 < \dots$

According to [81], the expanded expression of $W$ in equation (4.35) can be written as:

$$W = \sum_{j=0}^{\infty} Y_j(t) g_j(x),$$

where each $Y_j(t) \in \mathbb{R}^5$. Substituting this expression into equation (4.35) yields:

$$\frac{dY_j}{dt} = (L(Q_p^0) - \lambda_j D)Y_j. \qquad (4.37)$$

Thus, the DFE $Q_p^0$ is stable if and only if each $Y_j(t) \to 0$ when $t \to \infty$, that is, if and only if all the eigenvalues of matrix $H_j = L(Q_p^0) - \lambda_j D$ have negative real parts. Matrix $H_j$ can be written in the following form:

$$H_j = \begin{pmatrix} -\mu_p & 0 & 0 & 0 & -\omega_p\nu \\ 0 & -(\theta_p + \mu_p) & 0 & 0 & \omega_p\nu \\ 0 & \theta_p & -(\alpha_p + \mu_p + d) & 0 & 0 \\ 0 & 0 & \alpha_p & -\mu_p & 0 \\ 0 & 0 & \gamma_p & 0 & -(\nu + \mu_U + \varepsilon\lambda_j) \end{pmatrix}.$$

The characteristic polynomial of matrix $H_j$ is given by

$$P_1(X) = (X + \mu_p)^2(X^3 + a_2 X^2 + a_1 X + a_0)$$

where

$$\left\{ \begin{array}{l} a_2 = \theta_p + 2\mu_p + \alpha_p + d + \nu + \mu_U + \varepsilon\lambda_j \\ a_1 = (\theta_p + \mu_p)(\alpha_p + \mu_p + d + \nu + \mu_U + \varepsilon\lambda_j) + (\alpha_p + \mu_p + d)(\nu + \mu_U + \varepsilon\lambda_j) \\ a_0 = (\theta_p + \mu_p)(\alpha_p + \mu_p + d)(\nu + \mu_U + \varepsilon\lambda_j) - \gamma_p\theta_p\omega_p\nu \end{array} \right.$$

Firstly one can observe that $a_2 > 0, a_1 > 0$ and

$$a_0 = (\theta_p + \mu_p)(\alpha_p + \mu_p + d)(\nu + \mu_U + \varepsilon\lambda_j) \left( 1 - \frac{\gamma_p \theta_p \omega_p \nu}{(\theta_p + \mu_p)(\alpha_p + \mu_p + d)(\nu + \mu_U + \varepsilon\lambda_j)} \right)$$
Hence $a_0 > 0$ if $\frac{\gamma_p \theta_p \omega_p \nu}{(\theta_p + \mu_p)(\alpha_p + \mu_p + d)(\nu + \mu_U + \varepsilon \lambda_j)} < 1$, which is satisfied for all (non negative) $\lambda_j$ if it is satisfied for $\lambda_0 = 0$. The condition $a_0 > 0$ then leads to:

$$ \mathcal{R}_0^{(p)} = \frac{\gamma_p \theta_p \omega_p \nu}{(\theta_p + \mu_p)(\alpha_p + \mu_p + d)(\nu + \mu_U)} < 1. \quad (4.38) $$

Secondly, the expression

$$ a_2 a_1 - a_0 = (\alpha_p + \mu_p + d + \nu + \mu_U + \varepsilon \lambda_j)(\theta_p + \mu_p)^2 + (\theta_p + \mu_p)(\alpha_p + \mu_p + d + \nu + \mu_U + \varepsilon \lambda_j) \\ + (\alpha_p + \mu_p + d)(\nu + \mu_U + \varepsilon \lambda_j)) + \gamma_p \theta_p \omega_p \nu > 0 $$

Using the fact that $a_2 > 0$, $a_1 > 0$, $a_0 > 0$ and $a_2 a_1 > a_0$, the Routh-Hurwitz stability criterion indicates that all the eigenvalues of matrix $H_j$ have negative real parts. Hence the DFE $Q_p^0$ of subsystems (4.5) is LAS if and only $\mathcal{R}_0^{(p)} < 1$.

## 2. Existence and local stability of the endemic equilibrium $Q^*$

Suppose that $\mathcal{R}_0^{(p)} > 1$. The expression (4.16) of the endemic equilibrium can easily be established as the solution of a set of linear equations, derived from equating the right-hand side of equation (4.5) to zero. We first prove that $I_p^*$ and $S_p^*$ are positive. Recall that:

$$ I_p^* = \frac{\Lambda_p(\nu + \mu_U)(\mathcal{R}_0^{(p)} - 1)}{\gamma_p \omega_p \nu - d(\nu + \mu_U)} \quad \text{and} \quad S_p^* = \frac{\Lambda_p(\gamma_p \omega_p \nu - d(\nu + \mu_U)\mathcal{R}_0^{(p)})}{\mu_p \mathcal{R}_0^{(p)}(\gamma_p \omega_p \nu - d(\nu + \mu_U))}. $$

Note that:

$$ \begin{aligned} \mathcal{R}_0^{(p)} > 1 &\Rightarrow \frac{\gamma_p \omega_p \nu \theta_p}{(\nu + \mu_U)(\theta_p + \mu_p)(\alpha + \mu_p + d)} > 1, \\ &\Rightarrow d < \frac{\gamma_p \omega_p \nu \theta_p}{(\nu + \mu_U)(\theta_p + \mu_p)} - \alpha_p - \mu_p, \\ &\Rightarrow d < \frac{\gamma_p \omega_p \nu \theta_p}{(\nu + \mu_U)(\theta_p + \mu_p)}, \\ &\Rightarrow d(\nu + \mu_U) < \frac{\gamma_p \omega_p \nu \theta_p}{(\theta_p + \mu_p)} < \gamma_p \omega_p \nu \\ &\Rightarrow \gamma_p \omega_p \nu - d(\nu + \mu_U) > 0, \\ &\Rightarrow I_p^* > 0. \end{aligned} $$

Now using the fact that $\frac{\theta_p d}{(\theta_p + \mu_p)(\alpha_p + \mu_p + d)} < 1$, one has:

$$ \begin{aligned} \frac{\gamma_p \omega_p \nu \theta_p d(\nu + \mu_U)}{(\nu + \mu_U)(\theta_p + \mu_p)(\alpha_p + \mu_p + d)} &< \gamma_p \omega_p \nu \Rightarrow (\nu + \mu_U)\mathcal{R}_0^{(p)} < \gamma_p \omega_p \nu, \\ &\Rightarrow \gamma_p \omega_p \nu - d(\nu + \mu_U)\mathcal{R}_0^{(p)} > 0 \\ &\Rightarrow S_p^* > 0 \end{aligned} $$

This concludes the existence of the endemic equilibrium.

Second, we investigate the local stability of the endemic equilibrium, using the following theorem.
**Theorem 4.3** (Castillo-Chavez and Song [27]). Consider the following ordinary differential equations, with a parameter $\psi$:

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

Without loss of generality, it is assumed that 0 is an equilibrium for system (4.39) 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 (4.39) 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

$$ a = \sum_{k,i,j=1}^{n} v_k u_i u_j \frac{\partial^2 f_k}{\partial z_i \partial z_j}(0,0), \quad \text{and} \quad b = \sum_{k,i=1}^{n} v_k u_i \frac{\partial^2 f_k}{\partial z_i \partial \psi}(0,0). $$

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

(a) $a > 0, b > 0$. When $\psi < 0$ with $||\psi|| \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 a negative and locally asymptotically stable equilibrium;

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

(c) $a > 0, b < 0$. When $\psi < 0$ with $||\psi|| \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;

(d) $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.

In order to apply this theorem, subsystems (4.5) can be written as follows, with $z = (S, L, I, J, U)$:

$$ \left\{ \begin{aligned} \frac{\partial t z_1}{\partial t} &= F_1 = \Lambda_p - \mu_p z_1 - \frac{\omega_p \nu z_5}{z_1 + z_2 + z_3 + z_4}, \\ \frac{\partial t z_2}{\partial t} &= F_2 = \frac{\omega_p \nu z_5}{z_1 + z_2 + z_3 + z_4} - (\theta_p + \mu_p) z_2, \\ \frac{\partial t z_3}{\partial t} &= F_3 = \theta_p z_2 - (\alpha + \mu_p + d) z_3, \\ \frac{\partial t z_4}{\partial t} &= F_4 = \alpha_p z_3 - \mu_p z_4, \\ \frac{\partial t z_5}{\partial t} - \varepsilon \Delta z_5 &= F_5 = \gamma_p z_3 - (\nu + \mu_U) z_5, \end{aligned} \right. \qquad (4.40) $$

Solving $\mathcal{R}_0^{(p)} = 1$, we obtain the following bifurcation value for parameter $\omega_p$:

$$ \omega_p^* = \frac{(\mu_p + \theta)(\alpha_p + \mu_p + d)(\nu + \mu_U)}{\nu \gamma_p \theta_p}. $$
We linearise this system at the DFE $Q_p^0$, as previously in equation (4.37), setting $\omega_p$ to $\omega_p^*$. We need to determine the eigenvalues of matrix $L_{\omega_p^*} - \lambda_j D$, where $\lambda_j$ is an eigenvalue of the Laplacian operator $-\Delta$ (simplified notation $L_{\omega_p^*}$ is used instead of $L_{\omega_p^*}(Q_p^0)$). With $\lambda_j = \lambda_0 = 0$, the matrix admits $\beta_0 = 0$ as eigenvalue, the other eigenvalues still having a negative real part. For the other $\lambda_j$, all eigenvalues have negative real parts. So assumption $A_1$ of theorem 4.3 is verified.

Let us now verify assumption $A_2$. We need to compute the left and right eigenvectors of matrix $L_{\omega_p^*} - \lambda_j D$ associated with the eigenvalue $\beta$. The left eigenvector, denoted by $v = (v_1, v_2, v_3, v_4, v_5)$, satisfies the following equation:

$$v(L_{\omega_p^*} - \lambda_j D - \beta I) = 0,$$

where $I$ and $0$ are the identity matrix and null vector of dimension 5, respectively. For $\beta = \beta_0 = 0$ one has:

$$ \left( v_1 \ v_2 \ v_3 \ v_4 \ v_5 \right) \begin{pmatrix} -\mu_p & 0 & 0 & 0 & -\omega_p^* \nu \\ 0 & -(\theta_p + \mu_p) & 0 & 0 & \omega_p^* \nu \\ 0 & \theta_p & -(\alpha_p + \mu_p + d) & 0 & 0 \\ 0 & 0 & \alpha_p & -\mu_p & 0 \\ 0 & 0 & \gamma_p & 0 & -(\nu + \mu_U + \varepsilon \lambda_j) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

which gives:

$$ \left\{ \begin{array}{l} v_1 = 0, \\ v_2 = \frac{\gamma_p \theta_p}{(\mu_p + \theta_p)(\alpha_p + \mu_p + d)} v_5, \\ v_3 = \frac{\gamma_p}{(\alpha_p + \mu_p + d)} v_5, \\ v_4 = 0, \\ v_5 > 0. \end{array} \right. $$

Similarly the right eigenvector of matrix $L_{\omega_p^*} - \lambda_j D$, denoted by $u = (u_1, u_2, u_3, u_4, u_5)^T$, satisfies the following equation:

$$(L_{\omega_p^*} - \lambda_j D - \beta I)u = 0.$$

For $\beta = \beta_0 = 0$, one has:

$$ \begin{pmatrix} -\mu_p & 0 & 0 & 0 & -\omega_p^* \nu \\ 0 & -(\theta_p + \mu_p) & 0 & 0 & \omega_p^* \nu \\ 0 & \theta_p & -(\alpha_p + \mu_p + d) & 0 & 0 \\ 0 & 0 & \alpha_p & -\mu_p & 0 \\ 0 & 0 & \gamma_p & 0 & -(\nu + \mu_U + \varepsilon \lambda_j) \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, $$
which gives:

$$
\left\{
\begin{aligned}
u_1 &= -\frac{\omega_p^* \nu}{\mu} u_5, \\
u_2 &= \frac{\omega_p^* \nu}{\mu + \theta} u_5, \\
u_3 &= \frac{\omega_p^* \nu \theta}{(\mu + \theta)(\alpha_p + \mu_p + d)} u_5, \\
u_4 &= \frac{\alpha_p \omega_p^* \nu \theta_p}{(\mu_p + \theta_p)(\alpha_p + \mu_p + d)} u_5, \\
u_5 &> 0.
\end{aligned}
\right.
$$

Let us now compute $a$, defined in assumption $A_2$ of theorem 4.3, for system (4.40):

$$a = \sum_{i,j,k=1}^{5} v_k u_i u_j \frac{\partial^2 F_k}{\partial z_i \partial z_j} (Q_p^0, \omega_p^*).$$

The only terms that are non null correspond to:

$$\frac{\partial^2 F_2}{\partial z_i \partial z_5} (Q_p^0, \omega_p^*) = v_2 \frac{\partial^2 F_2}{\partial z_5 \partial z_j} (Q_p^0, \omega_p^*) = -\frac{\omega_p \nu^*}{S^0} \quad \text{for } i,j = 2,3,4.$$

Substituting these terms into the expression of $a$, one obtains:

$$
\begin{aligned}
a &= v_2 \left[ 2u_2 u_5 \frac{\partial^2 F_2}{\partial z_2 \partial z_5} (Q_p^0, \omega_p^*) + 2u_3 u_5 \frac{\partial^2 F_2}{\partial z_3 \partial z_5} (Q_p^0, \omega_p^*) + 2u_4 u_5 \frac{\partial^2 F_2}{\partial z_4 \partial z_5} (Q_p^0, \omega_p^*) \right], \\
&= -v_2 u_5 (u_2 + u_3 + u_4) \frac{\omega_p^* \nu}{S^0} < 0.
\end{aligned}
$$

Let us now compute $b$, defined in assumption $A_2$ of theorem 4.3, for system (4.40):

$$b = \sum_{k,i=1}^{5} v_k u_i \frac{\partial^2 F_k}{\partial z_i \partial \omega_p} (Q_p^0, \omega_p^*).$$

The only term that is non null corresponds to:

$$\frac{\partial^2 F_2}{\partial z_5 \partial \omega_p} (Q_p^0, \omega_p^*) = \frac{\nu}{S_0}.$$

Substituting this term into the expression of $b$, one obtains:

$$
\begin{aligned}
b &= v_2 u_5 \frac{\partial^2 F_2}{\partial z_5 \partial \omega_p} (Q_p^0, \omega_p^*) \\
&= v_2 u_5 \frac{\nu}{S_0} > 0.
\end{aligned}
$$

Thus, $a < 0$ and $b > 0$. Using Theorem 4.3, we conclude that when bifurcation parameter changes from $\omega_p < \omega_p^*$ to $\omega_p^* < \omega_p$, i.e. the basic reproduction number changes from $\mathcal{R}_0^{(p)} < 1$ to $\mathcal{R}_0^{(p)} > 1$, the DFE changes from GAS to unstable. Moreover, when the basic reproduction number is close to one, the endemic equilibrium $Q_p^*$ appears and is locally asymptotically stable. This completes the proof.
## 4.7 Appendix: Numerical scheme

### 4.7.1 Discretization of model system

Herein, we develop a non-standard finite difference scheme that is consistent with the dynamics of the continuous PDE model system (3.1). To do so, model system (3.1) can be written in the following compact form

$$ \frac{\partial P}{\partial t} = D \frac{\partial^2 P}{\partial^2 x} - A_0(P)P + A_1, \qquad (4.41) $$

where

$$ P = (S, L, I, R, B_e, U)^T, \quad D = \mathrm{diag}(0, 0, 0, 0, 0, \varepsilon), \quad A_1 = (\Lambda, 0, 0, 0, 0, 0)^T, $$

and

$$ A_0(P) = \begin{pmatrix} \mu_M + \mu(t) + \frac{\omega\nu U}{N} & 0 & 0 & 0 & 0 & 0 \\ -\frac{\omega\nu U}{N} & (1-q)\theta + \mu_M + \mu(t) & 0 & 0 & 0 & 0 \\ 0 & -(1-q)\theta & \alpha + \mu_M + d + \mu(t) & 0 & 0 & 0 \\ 0 & 0 & -\alpha & \mu_M + \mu(t) & 0 & 0 \\ -\delta_1(t) & -\delta_2(t) & -\delta_3(t) & -\delta_4(t) & \mu_e & 0 \\ 0 & 0 & -\gamma & 0 & 0 & (\nu + \mu_U) \end{pmatrix}. $$

For the temporal discretization, we replace the continuous time variable $t \in [0, \infty)$ by the discrete variables $t_n = ndt$, $n \in \mathbb{Z}$ where $dt$ is the temporal step size. Also, for the spatial discretization, the continuous space variable $x \in \Omega$ is replaced by the discrete variables, $x_j = jdx$, $j \in \mathbb{Z}$ where $dx$ is the spatial step size. With this in mind, one has

$$ \frac{\partial P}{\partial t} = \frac{P_j^{n+1} - P_j^n}{\psi(dt)} + o(\psi(dt)), \qquad (4.42) $$

and

$$ \frac{\partial^2 P}{\partial^2 t} = \frac{P_{j+1}^n - 2P_j^n + P_{j-1}^n}{(dx)^2} + o(dx)^2. \qquad (4.43) $$

Now, plugging Eqs. (4.42) and (4.43) in Eq. (4.41), gives

$$ \frac{P_j^{n+1} - P_j^n}{\psi(dt)} = D \frac{P_{j+1}^n - 2P_j^n + P_{j-1}^n}{(dx)^2} - A_0(P_j^n)P_j^{n+1} + A_1 \qquad (4.44) $$

Which can be rewritten as follows:

$$ P_j^{n+1} = A_1^{-1} A_2 (P_{j+1}^n + P_{j-1}^n) + A_1^{-1} A_3 P_j^n + A_1^{-1} A_1, \qquad (4.45) $$
where

$$ \mathcal{A}_1 = \begin{pmatrix} C_1 & 0 & 0 & 0 & 0 & 0 \\ O_1 & C_2 & 0 & 0 & 0 & 0 \\ 0 & O_2 & C_3 & 0 & 0 & 0 \\ 0 & 0 & O_3 & C_4 & 0 & 0 \\ O_4 & O_5 & O_6 & 0 & C_5 & 0 \\ 0 & 0 & O_8 & 0 & 0 & C_6 \end{pmatrix}, \quad \mathcal{A}_2 = \mathrm{diag}(0, 0, 0, 0, 0, d_U) $$

and

$$ \mathcal{A}_3 = \mathrm{diag}(1, 1, 1, 1, 1, 1 - 2d_U), $$

with

$$ \left\{ \begin{array}{l} C_1 = 1 + \psi(dt)[\mu_M + \mu(t) + \frac{\omega\nu U}{N}], \quad C_2 = 1 + \psi(dt)((1-q)\theta + \mu_M + \mu(t)), \\ \\ C_3 = 1 + \psi(dt)(\alpha + \mu_M + d + \mu(t)), \quad C_4 = 1 + \psi(dt)(\mu_M + \mu(t)), \quad d_U = \frac{\varepsilon\psi(dt)}{(dx)^2}, \\ \\ C_5 = 1 + \psi(dt)\mu_e, \quad C_6 = 1 + \psi(dt)(\nu + \mu_U), \\ \\ O_1 = -\psi(dt)\frac{\omega\nu U}{N}, \quad O_2 = -\psi(dt)(1-q)\theta, \quad O_3 = -\psi(dt)\alpha, \quad O_4 = \psi(dt)\delta_1(t), \\ \\ O_5 = \psi(dt)\delta_2(t), \quad O_6 = \psi(dt)\delta_3(t), \quad O_7 = \psi(dt)\delta_4(t) \quad \text{and} \quad O_8 = \psi(dt)\gamma \end{array} \right. $$

Now, using the notation of P defined in Eq.(4.41), and after the calculation of $\mathcal{A}^{-1}$, we obtain the following for model system (4.45) :

$$ \left\{ \begin{aligned} S_j^{n+1} &= \frac{1}{C_1} S_j^n + \frac{\Lambda}{C_1}, \\ L_j^{n+1} &= -\frac{O_1}{C_1 C_2} S_j^n + \frac{1}{C_2} L_j^n - \frac{\Lambda O_1}{C_1 C_2}, \\ I_j^{n+1} &= \frac{O_1 O_2}{C_1 C_2 C_3} S_j^n - \frac{O_2}{C_2 C_3} L_j^n + \frac{1}{C_3} I_j^n + \frac{\Lambda O_1 O_2}{C_1 C_2 C_3}, \\ R_j^{n+1} &= -\frac{O_1 O_2 O_3}{C_1 C_2 C_3 C_4} S_j^n + \frac{O_2 O_3}{C_2 C_3 C_4} L_j^n - \frac{O_3}{C_3 C_4} I_j^n + \frac{1}{C_4} R_j^n - \frac{\Lambda O_1 O_2 O_3}{C_1 C_2 C_3 C_4}, \\ B_{ej}^{n+1} &= \xi_1 S_j^n + \xi_2 L_j^n - \frac{(C_4 O_6 - O_3 O_7)}{C_3 C_4 C_5} I_j^n - \frac{O_7}{C_4 C_5} R_j^n + \frac{1}{C_5} B_{ej}^n + \xi_1 \Lambda, \\ U_j^{n+1} &= -\frac{O_1 O_2 O_8}{C_1 C_2 C_3 C_6} S_j^n + \frac{O_2 O_8}{C_2 C_3 C_6} L_j^n - \frac{O_8}{C_3 C_6} I_j^n + \frac{1 - 2d_U}{C_6} U_j^n + \frac{d}{C_6} (U_{j-1}^n + U_{j+1}^n) - \frac{\Lambda O_1 O_2 O_8}{C_1 C_2 C_3 C_6}, \end{aligned} \right. \quad (4.46) $$

where

$$ \xi_1 = \frac{-C_4 O_4 C_2 C_3 + C_4 O_1 O_5 C_3 - C_4 O_1 O_2 O_6 - O_1 O_2 O_3 O_7}{C_1 C_2 C_3 C_4 C_5} $$

and

$$ \xi_2 = \frac{-C_3 C_4 O_5 + C_4 O_2 O_6 - O_2 O_3 O_7}{C_2 C_3 C_4 C_5} $$
### 4.7.2 Discretization of boundary conditions

The discretization of the Neumann boundary conditions gives the following systems

$$ \left\{ \begin{array}{l} \displaystyle \frac{P_1^{n+1} - P_1^n}{\psi(dt)} = \frac{D}{(dx)^2}(P_2^n - 2P_1^n + P_0^n), \\[1.5em] P_2^n = P_0^n, \end{array} \right. \qquad (4.47) $$

and

$$ \left\{ \begin{array}{l} \displaystyle \frac{P_{Ne+1}^n - P_1^n}{\psi(dt)} = \frac{D}{(dx)^2}(P_{Ne+2}^n - 2P_{Ne+1}^n + P_{Ne}^n), \\[1.5em] P_{Ne+2}^n = P_{Ne}^n, \end{array} \right. \qquad (4.48) $$

Combining Eq. (4.47) and Eq. (4.48) yields

$$ \left\{ \begin{array}{l} P_1^{n+1} = (Id - 2A_2)P_1^n + 2A_2P_2^n, \\[1.5em] P_{Ne+1}^{n+1} = (Id - 2A_2)P_{Ne+1}^n + 2A_2P_{Ne}^n, \end{array} \right. \qquad (4.49) $$

Using the notation of $P$ defined in Eq. (4.41), Eq.4.49 can be writteen as follows:

$$ \left\{ \begin{array}{l} S_1^{n+1} = S_1^n, \quad L_1^{n+1} = L_1^n, \quad I_1^{n+1} = I_1^n, \quad R_1^{n+1} = R_1^n, \quad B_{e1}^{n+1} = B_{e1}^n \\ U_1^{n+1} = (1 - 2d_U)U_1^n + 2d_U U_2^n, \quad S_{Ne+1}^{n+1} = S_{Ne+1}^n, \\ L_{Ne+1}^{n+1} = L_{Ne+1}^n, \quad I_{Ne+1}^{n+1} = I_{Ne+1}^n, \quad R_{Ne+1}^{n+1} = R_{Ne+1}^n, \\ B_{eNe+1}^{n+1} = B_{eNe+1}^n, \quad U_{Ne+1}^{n+1} = (1 - 2d_U)U_{Ne+1}^n + 2d_U U_{Ne}^n, \end{array} \right. \qquad (4.50) $$

### 4.7.3 Consistence of scheme

The numerical scheme (4.46) is consistent if $d_U \le \frac{1}{2}$ where $d_U = \frac{\epsilon\psi(dt)}{(dx)^2}$

Taking $\psi(dt) = \frac{e^{Mt} - 1}{M}$, $M = \max\{\mu_M, \mu_M + (1 - q)\theta, \mu_M + \alpha + d, \mu_U + \nu, \mu_{Be}\}$, we have the result.
IMPULSIVE MODELLING OF RUST
DYNAMICS AND HYPERPARASITE
RELEASES FOR BIOCONTROL

In this chapter, the model is seasonal, but does not distinguish between young and mature leaves.
It includes the dynamics of a biocontrol agent with impulsive releases. The model is obtained
using the following hypotheses.

$H_1$: The model is based on a classical epidemiological model, in which the host population is subdivided by health status (healthy leaves $S$ and infected leaves $I$), considering the crop leaf as an individual. The model includes crop production (flowers $F$ and berries $B$) by leaves during the rainy season and harvest at the end of the season. The infection being mediated by fungus spores (uredospore $U$) released in the plantation, their dynamics are also included.

$H_2$: The recruitment of healthy leaves is constant.

$H_3$: The force of infection is density-dependent: it is proportional to the number of uredospores and inversely proportional to the total number of leaves.

$H_4$: The model alternates dry and rainy seasons. The dynamics are continuous during the rainy seasons and represented by discrete events during the dry seasons.

$H_5$: Space is not explicitly taken into account in this model, uredospores are homogeneously distributed on leaves.

$H_6$: The model uses an impulsive differential equations (IDE) formalism.

$H_7$: The model includes the dynamics of hyperparasites ($P$) which feed on uredospores. Hy-
perparasites are released one or several times per rainy season.

We then obtain an IDE model with five compartments. The aim of this last study is to
implement a biocontrol strategy of CLR based on hyperparasites.

The work presented in this chapter is based on preliminary results published in *IFAC-PapersOnLine*¹ [37] and presented in Appendix A.1. In this preliminary work the impact of control is represented by adding an extra negative term in the dynamics of uredospores. We then, use optimal control to maximise the yield and reduce the cost of control.
Contents

<table><tr><td>5.1</td><td>Introduction</td><td>91</td></tr><tr><td>5.2</td><td>Coffee leaf rust dynamics</td><td>92</td></tr><tr><td>5.2.1</td><td>Modelling of CLR</td><td>92</td></tr><tr><td>5.2.2</td><td>Mathematical analysis</td><td>94</td></tr><tr><td>5.2.3</td><td>Numerical simulations of CLR dynamics</td><td>99</td></tr><tr><td>5.3</td><td>Biocontrol of CLR using hyperparasites</td><td>102</td></tr><tr><td>5.3.1</td><td>Controlled periodic disease free solution and its stability</td><td>102</td></tr><tr><td>5.3.2</td><td>Semi-numerical analysis of the controlled model</td><td>104</td></tr><tr><td>5.3.3</td><td>Asymptotic behaviour for yearly releases</td><td>106</td></tr><tr><td>5.4</td><td>Multiple releases of hyperparasites per year</td><td>107</td></tr><tr><td>5.4.1</td><td>Multiple release controlled periodic disease free solution and its stability</td><td>107</td></tr><tr><td>5.4.2</td><td>Impact of release frequency on CLR control</td><td>112</td></tr><tr><td>5.5</td><td>Discussion</td><td>116</td></tr></table>

**Abstract:** Fungal diseases cause serious damages in crop worldwide. In particular, coffee leaf rust (CLR), caused by fungus *Hemileia vastatrix* attacks coffee leaves and reduces coffee yield. This paper presents a multi-seasonal model of the CLR development in the coffee plantation with continuous dynamics during the rainy season and a discrete event to represent the simpler dynamics during the dry season. Biological control using hyperparasites through one or more discrete introduction events over the year is then added. Analytical and semi-numerical studies are performed to identify how much and how frequently hyperparasites need to be introduced through the definition of a threshold value, as a function of various parameters. We show that the best strategy to efficiently control the disease depends on the hyperparasite mortality: low mortality parasites need be released only once a year, while high mortality parasites should be released more frequently to ensure their persistence in the plantation. This work hence provides qualitative and quantitative bases for the deployment of hyperparasite-based biocontrol, a promising alternative to fungicides for rust control.

**Highlights:**

* Original impulsive model of multi-seasonal fungal disease dynamics

* Explicit dynamics of hyperparasites, acting as predators of spores

* Hyperparasite-based biocontrol can drastically reduce disease impact on crop production

* Hyperparasite release frequency modulates the control efficiency

**Keywords:** crop protection, coffee leaf rust, hybrid model, seasonality, Floquet theory, stability
**5.1 Introduction**

Many of the most serious crop diseases are caused by fungi such as rusts [3]. Mathematical models of fungal diseases have received a lot of attention from researchers. For instance, Pivonia at al. studied the seasonal appearance of rusts in the United States, in particular soybean rust, thanks to a general disease model [124]. Rimbaud et al. investigated how the spatial deployment of resistant cultivars affects the resistance efficiency and durability, using a demogenetic model, for a seasonal crop infected by a fungal-like pathogen [134]. Fleming used a mathematical model to prove that a complex of polyphagous non-synchronised predators and parasites is likely to control only low-density cereal rust populations, but the author also proved that control at low rust density can delay epidemic development and thus substantially reduce yield losses [49]. Mammeri et al. studied the impact of spatial heterogeneities on the spread and control of grapevine powdery mildew [97], but only during a cropping season. There are fewer models of fungal diseases that target perennial hosts. For instance, Ravigné et al. looked at the impact of sexual and asexual reproduction on the epidemiological dynamics and showed that they could induce cyclic persistence of the disease, as it can observed for banana Sigatoka diseases [130]. Desprez-Loustau et al. developed a seasonal eco-evolutionary model of oak powdery mildew in Europe, based on a within-season and between-season transmission trade-off, which captures the main features of the disease, that is seasonality and pathogen species coexistence [36].

Plant growth and disease spread may be affected by seasonal patterns, which then need to be included in epidemiological models. It is the case of crops cultivated in temperate climates or tropical regions that alternate dry and rainy seasons. As the dynamics can substantially differ between the seasons with possibly rapid transitions, impulsive or semi-discrete models have been developed for several plant fungal diseases, as for oak powdery mildew in the study cited above. Among these models, one can also cite Tankam-Chedjou et al. who built a model to describe and control the dynamics of a banana soilborne pest in a multi-seasonal framework, optimising the fallow period durations between cropping seasons [155]. Periodic patterns are not necessarily linked to seasonality and can also be due to impulsive control strategies introduced in the epidemiological models. For instance, in the biological control framework, Nundloll et al. studied the periodic release of predators, natural enemies of the plant pest of interest, and determined the minimal predator rate required to eradicate the pest [115]. Xinzhu et al. formulated and analysed a model with continuous cultural control and with impulsive cultural control strategies such as replanting and/or removing diseased plants; they concluded that impulsive removing of diseased plants is more efficient and more economical than continuous removing [105]. Nembot et al. developed a model of cocoa black pod rot disease, caused by *Phytophthora megakarya*, and showed the impact of periodic impulsive sanitary harvests on the disease dynamics [113]. Semi-discrete models are hence particularly appropriate for seasonal plant and pathogen dynamics, as well as impulsive control strategies.

Control methods used to fight fungal diseases include chemical fungicide application, cultural practices, and the use of resistant cultivars [145]. These methods induce significant labour and/or material costs. Moreover, chemical fungicides are harmful to the environment [82], potentially also to farmers and consumers [70]; in addition, they may affect non-target organisms and induce pest resistance [76]. Therefore, researchers are currently investigating alternative control methods such as biocontrol [75, 6]. Several biological agents have been tested to control fungal plant diseases, among which various bacteria. Antagonist bacteria such as *Bacillus* species make
plants more resistant to fungal infections [144, 68, 148, 34]. Other investigations on biocontrol focused on hyperparasites, which are parasites whose host is often an insect pest but can also be a fungus. For example, several insects such as *Mycodiplosis* (Diptera) [67] and fungi such as *Lecanicillium lecanii* (previously called *Verticillium lecanii*) [104, 55] feed on rust spores.

As a foundation for our model, we considered coffee leaf rust (CLR), which is caused by a fungus, *Hemileia vastatrix*. It attacks the lower leaves of the coffee tree and causes premature defoliation, which reduces the photosynthetic capacity and weakens the tree. Leaf fall causes abortion of a large part of the flowers and fruits, as well as desiccation of shoots. It has direct and indirect economic impacts on coffee production. Direct impacts include decreased quantity and quality of yield. In some cases, more than 70% of the coffee production is lost [58, 8]. Indirect impacts include increased costs to control the disease.

The fungus *Hemileia vastatrix* is a basidiomycete, which, like most fungi, produces spores (used for reproduction). Its control is achieved similarly to that of other fungal diseases. For instance, the antagonist bacteria *Bacillus subtilis*, isolated from the rhizosphere of coffee crops, largely reduced the growth of CLR (up to 68%) under in vitro conditions [34]. The fungus *Cladosporium hemileiae*, inhibits the evolution of *H. vastatrix* mycelium [149]. The fungus *Lecanicillium lecanii* [168] and the insect *Mycodiplosis* [139, 63] can feed on spores. In all these biocontrol studies, the dose and time of application of the bacteria or hyperparasites affected their effectiveness in controlling coffee leaf rust development. In this work, we chose to focus on a *Mycodiplosis*-like hyperparasite, which was reported as a promising biocontrol agent and paid special attention to the quantity and frequency of the releases. To this aim, mathematical modelling can help identifying when and how to release a biocontrol agent.

Our goal in this paper is to control a fungal disease in the field, by means of hyperparasite releases, using a mathematical modelling approach. To achieve this goal, we extend an original impulsive CLR model that we previously developed [37]. In particular, we explicitly represent the hyperparasite dynamics, the latter acting as predators of spores. Moreover, we consider in this work discrete releases of the hyperparasite instead of a continuous control.

This paper is organised as follows. Section 5.2 is devoted to the formulation of the impulsive model of CLR, with continuous dynamics during the rainy seasons and discrete events for the dry seasons. It also presents the mathematical analysis of the model using Floquet theory and simulations of this impulsive model to illustrate the theoretical results. In Section 5.3, we introduce the biocontrol in the model, based on hyperparasite yearly releases. We semi-analytically study the stability of the control model to obtain the stability regions as functions of parameter values. Section 5.4 presents the impact of multiple yearly releases of hyperparasites, with contrasting results according to the value of the hyperparasite mortality rate. The last Section 5.5 concludes the paper with a discussion of the main results and possible perspectives for future work.

## 5.2 Coffee leaf rust dynamics

### 5.2.1 Modelling of CLR

In this section, we formulate a mathematical model for the CLR in the coffee plantation. To do so, we consider the dynamics of fungus during the production season, which corresponds to the rainy season for some countries, through ordinary differential equations; while the non-
<figure><img src="image_22.png" /><figcaption>Figure 5.1: Diagram of the CLR multi-seasonal model. Model equations are given in system (5.1). State variables are: healthy leaves (S), infectious leaves (I), uredospores (U), flowers (F) and berries (B). T corresponds to the length of rainy season and n the year number.</figcaption></figure>

production period, which is the dry season, is represented by impulses since the dynamics are simpler during that period. More precisely, we consider that in a coffee plantation we can find: susceptible leaves *S*, which are healthy leaves that have not (yet) been attacked, infected leaves *I*, uredospores *U*, which the fungus *H. vastatrix* uses for its asexual reproduction, flowers *F* and berries *B*.

During the production period with length $T > 0$ which represents the rainy season, the recruitment of healthy leaves occurs at rate $\Lambda$. Uredospores are deposited on leaves of all leaves at rate $\nu$. When the conditions are favourable, they germinate with efficiency rate $\omega$ and the susceptible leaves become infected leaves. All leaves undergo natural mortality with baseline rate $\mu$ and the infected leaves have an additional mortality rate $d$ due to the disease. Uredospores are produced by infected leaves at rate $\gamma$ and lose their ability to infect coffee leaves at constant rate $\mu_U$. Susceptible and infected leaves produce the flowers with constant rate $\delta_S$ and $\delta_I$ respectively and $\delta_S > \delta_I$. The flowers become berries at rate $\theta$. At the end of the production period, harvest occurs instantaneously and we consider that the dry season or non-production period can be summarised in a discrete time event due to the slow growth of leaves during the dry season. The impact of harvest and dry season reduces the number of susceptible leaves, infected leaves and uredospores with rates $\phi_S$, $\phi_I$ and $\phi_U$ respectively. Using the fact that uredospores lose their ability to infect quickly, we assume that $\phi_U$ is very close to 0. The representation of this evolution in the coffee plantation over a multi-year period, gives the multi-seasonal model, given in Figure 5.1, with the continuous dynamics during the rainy season and switching event for the dry season.

Using flowchart in Figure 5.1, we can write the following impulsive differential system:
$$
\text{Rainy season for } t \neq nT \text{:}
\begin{cases}
\dot{S} = \Lambda - \omega\nu U \frac{S}{N} - \mu S, \\
\dot{I} = \omega\nu U \frac{S}{N} - (\mu + d)I, \\
\dot{U} = \gamma I - (\nu + \mu_U)U, \\
\dot{F} = \delta_S S + \delta_I I - (\theta + \mu_F)F, \\
\dot{B} = \theta F - \mu_B B;
\end{cases}
\tag{5.1}
$$

Dry season:

$$
\left\{
\begin{aligned}
S(nT^+) &= \varphi_s S(nT), \\
I(nT^+) &= \varphi_I I(nT), \\
U(nT^+) &= \varphi_U U(nT), \\
F(nT^+) &= 0; \\
B(nT^+) &= 0.
\end{aligned}
\right.
$$

where $S(t)$, $I(t)$, $U(t)$, $F(t)$, $B(t)$ and $N(t) = S(t) + I(t)$ represent the number of healthy leaves, infected leaves, uredospores, flowers, berries and total number of leaves, respectively, at time $t$; with $0 < \varphi_S, \varphi_I, \varphi_U < 1$ and $\varphi_U$ close to 0.

5.2.2 Mathematical analysis

Herein, we first presented basic properties of the model (5.1). Secondly, the periodic disease free
solution is computed and its stability is proven. Finally, numerical simulations are performed to
illustrate our mathematical results.

5.2.2.1 Basic properties of the model

Let $\mathbb{R}_+^5 = \{X \in \mathbb{R}^5, X \ge 0\}$. Denote by $f = (f_1, f_2, f_3, f_4, f_5)$ the map given by the right-hand side of system (5.1) when $t \ne nT$, and consider the initial conditions $(S(0), I(0), U(0), F(0), B(0)) \in \mathbb{R}_+^* \times \mathbb{R}_+^4$.

The solutions of system (5.1) are non-negative. Indeed, suppose that one of the variables,
that we will denote $x$, is equal to 0 at some instant, with all the others being non-negative.
A quick analysis shows that $\dot{x} \ge 0$, so that $x$ cannot become negative. This implies that the
solutions are non-negative.

By adding the first and second equations of system (5.1), we obtain the dynamics of total
leaves which satisfies:

$$
\begin{cases}
\dot{N} = \Lambda - \mu N - dI, \quad t \neq nT; \\
N(nT^+) \leq \varphi_S N(nT).
\end{cases}
\tag{5.2}
$$

From equation (5.2), and using the non-negativity of variable $I$, we obtain $\Lambda - (\mu + d)N \le$
$\dot{N} \le \Lambda - \mu N$ and $N(nT^+) \le \varphi_S N(nT)$. It follows that $\min(\Lambda/(\mu + d), N(0)) \le N(t) \le$
$\max(\Lambda/\mu, N(0)) = \Gamma_N$, which implies that $S, I \le \Gamma_N$. For $\Lambda/(\mu + d) < N(0) < \Lambda/\mu$, replacing $S$
and $I$ by this upper bound value in the uredospore and flower dynamics, we obtain the following
systems:

$$ egin{cases} 
\dot{U} \leq \gamma \frac{\Lambda}{\mu} - (\nu + \mu_U)U, & t \neq nT; \\ 
U(nT^+) = \varphi_U U(nT); 
\end{cases} 
\quad (5.3) $$

$$ egin{cases} 
\dot{F} \leq \delta_S \frac{\Lambda}{\mu} + \delta_I \frac{\Lambda}{\mu} - (\theta + \mu_F)F, & t \neq nT; \\ 
F(nT^+) = 0. 
\end{cases} 
\quad (5.4) $$

Identifying the values where the upper bounds of (5.3) and (5.4) are equal to 0, one has $U(t) \leq \max(\frac{\gamma\Lambda}{\mu(\nu+\mu_U)}, U(0))$ and $F(t) \leq \max(\frac{\Lambda(\delta_S+\delta_I)}{\mu(\theta+\mu_F)}, F(0))$. Replace the upper bound of $F$ in the berry dynamics, we obtain

$$ egin{cases} 
\dot{B} \leq \theta \frac{\Lambda(\delta_S + \delta_I)}{\mu(\theta + \mu_F)} - \mu_B U, & t \neq nT; \\ 
B(nT^+) = 0. 
\end{cases} 
\quad (5.5) $$

Identifying the values where the upper bounds of (5.5) is equal to 0, one has $B(t) \leq \max(\frac{\Lambda\theta(\delta_S+\delta_I)}{\mu(\theta+\mu_F)\mu_B}, B(0))$. Hence, we have shown the boundedness of $S$, $I$, $U$, $F$ and $B$. Also the region given by

$$ G = \left\{ (S, I, U, F, B) \in \mathbb{R}_+^5 \mid \begin{array}{l} S(t) + I(t) \leq \frac{\Lambda}{\mu}, \quad U(t) \leq \frac{\gamma\Lambda}{\mu(\nu + \mu_U)}, \\ F(t) \leq \frac{\Lambda(\delta_S + \delta_I)}{\mu(\theta + \mu_F)}, \quad B(t) \leq \frac{\Lambda\theta(\delta_S + \delta_I)}{\mu(\theta + \mu_F)\mu_B} \end{array} \right\} $$

is positively invariant for impulsive system (5.1).

The solutions of system (5.1) are piecewise continuous functions $X: \mathbb{R}_+ \to \mathbb{R}_+^4$. The non-null constant parameter $\Lambda$, allows to the variable $N$ remain not close to zero and the smooth properties of the right side of system (5.1) guarantee the existence and uniqueness of positive solutions of system (5.1).

Since the berry state variables $F$ and $B$ is not present in the other equations of system (5.1), we do no need to consider $F$ and $B$ in the rest of the mathematical analysis. We then analyse this new system:

$$ egin{cases} 
\dot{S} = \Lambda - \frac{\omega\nu U}{N}S - \mu S, & t \neq nT; \\ 
\dot{I} = \frac{\omega\nu U}{N}S - (\mu + d)I, & t \neq nT; \\ 
\dot{U} = \gamma I - (\nu + \mu_U)U, & t \neq nT; \\ 
S(nT^+) = \varphi_S S(nT); \\ 
I(nT^+) = \varphi_I I(nT); \\ 
U(nT^+) = \varphi_U U(nT). 
\end{cases} 
\quad (5.6) $$
### 5.2.2.2 Periodic disease free solution and its stability

The periodic disease free solution (PDFS) occurs when $I = 0$ and $U = 0$. Replacing these values in system (5.6), we obtain

$$ egin{cases} \dot{S} = \Lambda - \mu S, & t \neq nT; \\ S(nT^+) = \varphi_S S(nT). \end{cases} $$

Solving the above equation for $t \in (nT, (n+1)T]$, yields

$$ S(t) = \frac{\Lambda}{\mu} + \left( S(nT^+) - \frac{\Lambda}{\mu} \right) e^{-\mu(t-nT)}. \quad (5.7) $$

This implies that for $t = (n+1)T$, one has

$$ S((n+1)T) = \frac{\Lambda}{\mu} + \left( S(nT^+) - \frac{\Lambda}{\mu} \right) e^{-\mu T}. $$

Using the impulsive condition $S((n+1)T^+) = \varphi_S S((n+1)T)$ yields

$$ S((n+1)T^+) = \varphi_S \left( \frac{\Lambda}{\mu} + (S(nT^+) - \frac{\Lambda}{\mu})e^{-\mu T} \right). \quad (5.8) $$

The fixed point of equation (5.8) is given by

$$ S(nT^+) = \frac{\Lambda \varphi_S (1 - e^{-\mu T})}{\mu (1 - \varphi_S e^{-\mu T})} > 0. \quad (5.9) $$

Substituting the value of $S(nT^+)$ into equation (5.7), for $t \in (nT, (n+1)T]$, the solution is

$$ S^T(t) = \frac{\Lambda}{\mu} \left[ 1 - \frac{(1 - \varphi_S) e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right]. \quad (5.10) $$

Finally, the PDFS is $X^T(t)=(S^T(t), 0, 0)$, where $S^T(t)$ is defined above.

For the stability, we first study the local stability of the PDFS $X^T(t)$ using small amplitude perturbation methods. For that, let us denote $\tilde{X}(t) = X(t) - X^T(t)$, where $X(t) = (S(t), I(t), U(t))^T$ and $\tilde{X}(t)$ is understood to be small amplitude perturbations. Substituting the expression of $\tilde{X}(t)$ in system (5.6) gives

$$ egin{cases} \dot{\tilde{S}} = - \frac{\omega \nu \tilde{U}}{\tilde{N} + S^T(t)} (\tilde{S} + S^T(t)) - \mu \tilde{S}, & t \neq nT; \\ \dot{\tilde{I}} = \frac{\omega \nu \tilde{U}}{\tilde{N} + S^T(t)} (\tilde{S} + S^T(t)) - (\mu + d) \tilde{I}, & t \neq nT, \\ \dot{\tilde{U}} = \gamma \tilde{I} - (\nu + \mu_U) \tilde{U}, & t \neq nT; \\ \tilde{S}(nT^+) = \varphi_S \tilde{S}(nT); \\ \tilde{I}(nT^+) = \varphi_I \tilde{I}(nT); \\ \tilde{U}(nT^+) = \varphi_U \tilde{U}(nT). \end{cases} \quad (5.11) $$

where $\tilde{N}(t) = \tilde{S}(t) + \tilde{I}(t)$ and $N^T(t) = S^T(t)$.
The linearization of system (5.11) in the neighbourhood of $0_{\mathbb{R}^3}$ is

$$ \begin{cases} \dot{\tilde{X}}(t) = A\tilde{X}(t), \ t \neq nT; \\ \tilde{X}(nT^+) = \text{diag}(\varphi_S, \varphi_I, \varphi_U)\tilde{X}(nT); \end{cases} \qquad (5.12) $$

where

$$ A = \begin{pmatrix} -\mu & 0 & -\omega\nu \\ 0 & -(\mu + d) & \omega\nu \\ 0 & \gamma & -(\nu + \mu_U) \end{pmatrix}. $$

Solving the first equation of system (5.12) for $t \in (0, T]$ yields

$$ \tilde{X}(t) = \Phi_A(t)\tilde{X}(0^+), \qquad (5.13) $$

where

$$ \Phi_A(t) = e^{At} = \begin{pmatrix} e^{-\mu t} & * & * \\ 0 & \psi_{22}(t) & \psi_{23}(t) \\ 0 & \psi_{32}(t) & \psi_{33}(t) \end{pmatrix} $$

and

$$ \left\{ \begin{aligned} \psi_{22}(t) &= \frac{1}{2\beta} \left[ (\beta + k_1 - k_2)e^{\lambda_1 t} + (\beta - k_1 + k_2)e^{\lambda_2 t} \right], \\ \psi_{23}(t) &= \frac{\nu\omega}{\beta} \left[ -e^{\lambda_1 t} + e^{\lambda_2 t} \right], \psi_{32}(t) = \frac{\gamma}{\beta} \left[ -e^{\lambda_1 t} + e^{\lambda_2 t} \right], \\ \psi_{33}(t) &= \frac{1}{2\beta} \left[ (\beta - k_1 + k_2)e^{\lambda_1 t} + (\beta + k_1 - k_2)e^{\lambda_2 t} \right], \end{aligned} \right. $$

with

$$ \left\{ \begin{array}{ll} k_1 = \mu + d, & k_2 = \nu + \mu_U, \\ \alpha = k_1 + k_2, & \beta = \sqrt{(k_1 - k_2)^2 + 4\gamma\omega\nu}, \\ \lambda_1 = -\frac{\alpha}{2} - \frac{\beta}{2} & \text{and } \lambda_2 = -\frac{\alpha}{2} + \frac{\beta}{2}. \end{array} \right. $$

Using the impulsive conditions $\tilde{X}(nT^+) = \text{diag}(\varphi_S, \varphi_I, \varphi_U)\tilde{X}(nT)$, the solution given by equation (5.13) becomes

$$ \tilde{X}(nT^+) = \text{diag}(\varphi_S, \varphi_I, \varphi_U)\Phi_A(T)\tilde{X}((n-1)T^+). \qquad (5.14) $$

From the above equation, the monodromy matrix $M = \text{diag}(\varphi_S, \varphi_I, \varphi_U)\Phi_A(T)$, i.e.

$$ M = \begin{pmatrix} \varphi_S e^{-\mu T} & * & * \\ 0 & \varphi_I \psi_{22}(T) & \varphi_I \psi_{23}(T) \\ 0 & \varphi_U \psi_{32}(T) & \varphi_U \psi_{33}(T) \end{pmatrix}. $$

Due to the block-triangular form of monodromy matrix $M$, there is no need to calculate the exact form of (*) for the following analysis.

The Floquet multipliers of the monodromy matrix $M$ are given by $\chi_1 = \varphi_S e^{-\mu T} < 1$ and the eigenvalues of the sub-matrix

$$ M_1 = \begin{pmatrix} \varphi_I \psi_{22}(T) & \varphi_I \psi_{23}(T) \\ \varphi_U \psi_{32}(T) & \varphi_U \psi_{33}(T) \end{pmatrix}. \qquad (5.15) $$
The PDFS $X^T(t)$ is locally asymptotically stable if the Floquet multipliers of the monodromy matrix $M$ stay inside the unit circle, which is equivalent to the Floquet multipliers of sub-matrix $M_1$ staying inside the unit circle. To obtain this result, the submonodromy matrix $M_1$ needs to satisfy the following Jury conditions [173]:

$$ \left\{ \begin{array}{l} -\mathrm{Tr}(M_1) - \det(M_1) < 1, \\ \det(M_1) < 1, \\ \mathrm{Tr}(M_1) - \det(M_1) < 1, \end{array} \right. \qquad (5.16) $$

where

$$ \left\{ \begin{array}{l} \mathrm{Tr}(M_1) = \varphi_I \psi_{22}(T) + \varphi_U \psi_{33}(T), \\ \det(M_1) = \varphi_I \varphi_U (\psi_{22}(T)\psi_{33}(T) - \psi_{23}(T)\psi_{32}(T)). \end{array} \right. $$

Using the hypothesis $\varphi_U \to 0$, conditions (5.16) hold if and only if $\mathcal{R} = \varphi_I \psi_{22}(T) < 1$. Finally, we obtain the following Lemma for local stability.

**Lemma 5.1.** The PDFS $X^T(t)$ is locally asymptotically stable provided that $\mathcal{R} < 1$ and unstable else, where

$$ \mathcal{R} = \frac{\varphi_I}{2\beta} \left[ (\beta + k_1 - k_2)e^{\lambda_1 T} + (\beta - k_1 + k_2)e^{\lambda_2 T} \right]. \qquad (5.17) $$

Moreover, the following result about the global asymptotic stability of the PDFS $X^T(t)$ of system (5.6) holds.

**Theorem 5.1.** The PDFS $X^T(t)$ of system (5.6) is globally asymptotically stable if $\mathcal{R} < 1$ and unstable else.

*Proof of Theorem 5.1.* Let us prove the global attractivity of the PDFS $X^T(t)$. Consider the following subsystem:

$$ \left\{ \begin{array}{l} \dot{I} = \frac{\omega\nu U}{N} S - (\mu + d)I, \quad t \neq nT; \\ \dot{U} = \gamma I - (\nu + \mu_U)U, \quad t \neq nT; \\ I(nT^+) = \varphi_I I(nT); \\ U(nT^+) = \varphi_U U(nT). \end{array} \right. \qquad (5.18) $$

The trivial solution of system (5.18) is $(0, 0)$ and the system is cooperative because $\partial I/\partial U = \frac{\omega\nu S}{N} > 0$, $\partial U/\partial I = \gamma > 0$. As $S(t)/N(t) < 1$, one has that system (5.18) is upper bounded by the following cooperative system:

$$ \left\{ \begin{array}{l} \dot{I}_1 = \omega\nu U_1 - (\mu + d)I_1, \quad t \neq nT; \\ \dot{U}_1 = \gamma I_1 - (\nu + \mu_U)U_1, \quad t \neq nT; \\ I_1(nT^+) = \varphi_I I_1(nT); \\ U_1(nT^+) = \varphi_U U_1(nT). \end{array} \right. \qquad (5.19) $$

Applying Kamke's theorem [33], one has that $I(t) \le I_1(t)$, $U(t) \le U_1(t)$ when $I_1(0) = I(0)$ and $U_1(0) = U(0)$. When $\mathcal{R} < 1$, using the previous result of local stability, $(I_1(t), U_1(t)) \to (0^+, 0^+)$. Then, one has

$$ (I_1(t), U_1(t)) \to (0^+, 0^+) \Rightarrow (I(t), U(t)) \to (0^+, 0^+). $$
Since $U \to 0^+$, for any arbitrary positive $\epsilon_U$ there exists $t_0 > 0$ such that for $t \ge t_0$, $U(t) \le \epsilon_U$. Using non-negativity of the solutions and the fact that $S(t)/N(t) \le 1$ into the first equation of system (5.6), one has:

$$ \Lambda - \omega\nu\epsilon_U - \mu S(t) \le \dot{S}(t) \le \Lambda - \mu S(t) $$

Applying the comparison principle on the above differential inequalities and using the impulsive condition, one has $S_1(t) \le S(t) \le S_2(t)$, with $S_1(t)$ and $S_2(t)$ the solutions of the following impulsive differential equations respectively:

$$ \begin{cases} \dot{S}_1 = \Lambda - \omega\nu\epsilon_U - \mu S_1, & t \ne nT; \\ S_1(nT^+) = \varphi_S S_1(nT), & t = nT; \end{cases} $$

and

$$ \begin{cases} \dot{S}_2 = \Lambda - \mu S_2, & t \ne nT; \\ S_2(nT^+) = \varphi_S S_2(nT). & t = nT; \end{cases} $$

One can observe that $S_1(t) \to S_{\epsilon_U}^T(t)$ and $S_2(t) \to S^T(t)$ asymptotically. This implies that

$$ S_{\epsilon_U}^T(t) \le S(t) \le S^T(t). \qquad (5.20) $$

Since we can write $\epsilon_U$ as small as we want using the fact that $U(t) \to 0$, we obtain $S_{\epsilon_U}^T(t) \to S^T(t)$ when $\epsilon_U \to 0$. Finally, using (5.20), one has that $S(t) \to S^T(t)$. Which implies that the PDFS is globally attractive. In other words, independently from the initial conditions ($S_0, I_0, U_0$), one has that

$$ (S(t), I(t), U(t)) \to (S^T(t), 0, 0) $$

This concludes the proof of global asymptotic stability. $\square$

### 5.2.3 Numerical simulations of CLR dynamics

Herein, we present the results of numerical simulations of system (5.1) using ode45 in Matlab to integrate the differential equations during the rainy seasons, with the impulses giving each year's initial conditions. We take $T = 250$ days, which is the length of rainy season in Cameroon [151] and also corresponds to the length of coffee production period. For the simulations, we suppose that all leaves are initially healthy, which means that $I(0) = 0$ leaf, with $S(0) = 500$ leaves. Infection is initiated by uredospores, with $U(0) = 3000$ spores. There are no flowers and berries initially, $F(0) = 0$ flower and $B(0) = 0$ berry, because simulations start at the beginning of the production period. The initial conditions are then

$$ (S(0), I(0), U(0), F(0), B(0)) = (500, 0, 3000, 0, 0). \qquad (5.21) $$

Table 5.1 summarises the parameter definitions and values used in this paper and that this calculation of the values of some parameters values.

- The maturation of flowers takes 180 days, which represents stage I + stage II + stage III presented by Torres Castillo et al [159], then the maturation rate is $\theta = 1/180 \approx 0.0055/\text{day}$.

- In the coffee tree, a node produces 30 flowers during the owering period ([31], page 59).
<table><caption>Table 5.1: Description and values of parameters for system (5.22).</caption><thead><tr><th>Symbol</th><th>Biological meaning</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>T</td><td>Rainy season duration</td><td>250 days</td><td>[151]</td></tr><tr><td>Λ</td><td>Recruitment rate of S</td><td>8 leaves/day</td><td>Assumed</td></tr><tr><td>ω</td><td>Germination rate</td><td>0.065 leaves/spore</td><td>[180, 131]</td></tr><tr><td>μ</td><td>Mortality rate of leaves</td><td>0.0034/day</td><td>Assumed</td></tr><tr><td>γ</td><td>Sporulation rate by I</td><td>2 spores/leaf.day</td><td>[19]</td></tr><tr><td>d</td><td>Mortality rate due to CLR</td><td>0.056/day</td><td>Assumed</td></tr><tr><td>ν</td><td>Deposition rate</td><td>0.09/day</td><td>[19]</td></tr><tr><td>μ<sub>U</sub></td><td>Mortality rate of U</td><td>0.015/day</td><td>Assumed</td></tr><tr><td>φ<sub>S</sub></td><td>Reduction rate of S</td><td>0.7</td><td>Assumed</td></tr><tr><td>φ<sub>I</sub></td><td>Reduction rate of I</td><td>0.4</td><td>Assumed</td></tr><tr><td>φ<sub>U</sub></td><td>Reduction rate of U</td><td>0.1</td><td>Assumed</td></tr><tr><td>δ<sub>S</sub></td><td>Flower production rate by S</td><td>0.08 flowers/leaf.day</td><td>[28]</td></tr><tr><td>δ<sub>I</sub></td><td>Flower production rate by I</td><td>0.04 flowers/leaf.day</td><td>[28]</td></tr><tr><td>μ<sub>F</sub></td><td>Mortality rate of F</td><td>2,4 × 10<sup>-3</sup>/day</td><td>[31]</td></tr><tr><td>θ</td><td>Maturation of F</td><td>0.0055/day</td><td>[159]</td></tr><tr><td>μ<sub>B</sub></td><td>Mortality rate of B</td><td>2,4 × 10<sup>-3</sup>/day</td><td>[31]</td></tr><tr><td>Λ<sub>P</sub></td><td>Yearly released quantity of P</td><td>variable predators</td><td>Assumed</td></tr><tr><td>a</td><td>Consumption rate</td><td>variable spores/predator.day</td><td>Assumed</td></tr><tr><td>e</td><td>Biomass transformation rate</td><td>0.7 predators/spore</td><td>Assumed</td></tr><tr><td>K</td><td>Saturation constant of P</td><td>100000 spores</td><td>Assumed</td></tr><tr><td>μ<sub>P</sub></td><td>Mortality rate of P</td><td>0.003 &amp; 0.1/day</td><td>Assumed</td></tr><tr><td>φ<sub>P</sub></td><td>Reduction rate of P</td><td>0.3</td><td>Assumed</td></tr></tbody></table>

Each node has 2 leaves, so a leaf corresponds to 15 flowers. The production of flowers in the model is hence $\delta_S = 0.08$ flowers/leaf.day, we suppose $\delta_I = 0.04$ flowers/leaf.day due to the disease which reduce the production of flowers by the infected leaves.

- In ([31], page 28), on 100 flowers 65 set fruit after pollination. This implies $e^{-180\mu_F} = 0.65$, computation gives $\mu_F = 2.4 \times 10^{-3}/\text{day}$. We suppose that the berries have the same mortality rate, then $\mu_B = 2.4 \times 10^{-3}/\text{day}$.

In Theorem 5.1, we show that when the spectral radius $\mathcal{R}$ is greater than one, the periodic disease free solution (PDFS) is unstable. In this case, CLR can establish itself, which is confirmed by the solid curve of Figure 5.2, drawn for $\mathcal{R} = 2.53 > 1$.

The dashed curve of Figure 5.2, drawn for $\mathcal{R} = 0.27 < 1$, shows that, for the initial conditions considered, the number of infected leaves and uredospores converges towards zero. The system then converges to the PDFS. We can observe in subplot (a) that, once the stationary regime is established, there are considerably less healthy leaves with endemic CLR (solid curve) than when CLR disappears (dashed curve). The same phenomenon can be observed in subplot (d), which represents the berry trajectories. At the end of the 6th year, the number of berries is 9844 without disease (dashed curve) and 4290 with disease (solid curve), indicating a yield loss larger than 56% due to CLR.
<figure><img src="image_23.png" /><figcaption>Figure 5.2: Impact of CLR on the production of coffee berries. The plots present the trajectories of impulsive differential system (5.6) when the PDFS is stable (dashed curve, $\gamma = 1.6$ spores/leaf.day $\Rightarrow \mathcal{R} = 0.27 < 1$) and when the PDFS is unstable (solid curve, default value $\gamma = 2$ spores/leaf.day $\Rightarrow \mathcal{R} = 2.53 > 1$). Subplots represent: (a) healthy leaves S, (b) infectious leaves I, (c) uredospores U and (d) berries B. Remaining parameter values are given in Table 5.1 and initial conditions by (5.21).</figcaption></figure>
## 5.3 Biocontrol of CLR using hyperparasites

To limit the impact of CLR on coffee berry production, we now introduce a biocontrol agent in the model, more specifically an hyperparasite such as *Mycodiplosis*, which consumes uredospores of *H. vastatrix*. The new model is then an extension of system (5.6), to which we add an equation for the hyperparasites (or predators) *P*, that prey on the uredospores with consumption rate *a*, biomass transformation rate *e* and mortality rate *μ*<sub>P</sub>. A proportion *φ*<sub>P</sub> survives the dry season. Our aim is to identify how many and when to release hyperparasites in order to control CLR. In this section, we assume that a quantity Λ<sub>P</sub> of hyperparasites is released once a year, at the start of the production season. This translates into a Λ<sub>P</sub> jump in the hyperparasite population at each switching moment. The model with biocontrol is then given by

Rainy season for $t \neq nT$:

$$ 
\begin{cases}
\dot{S} = \Lambda - \omega\nu U \frac{S}{N} - \mu S; \\
\dot{I} = \omega\nu U \frac{S}{N} - (\mu + d)I; \\
\dot{U} = \gamma I - (\nu + \mu_U)U - aP \frac{U}{K+U}; \\
\dot{P} = eaP \frac{U}{K+U} - \mu_P P; \\
\dot{F} = \delta_S S + \delta_I I - (\theta + \mu_F)F, \\
\dot{B} = \theta F - \mu_B B.
\end{cases}
$$

Dry season:

$$ 
\begin{cases}
S(nT^+) = \varphi_s S(nT); \\
I(nT^+) = \varphi_I I(nT); \\
U(nT^+) = \varphi_U U(nT); \\
P(nT^+) = \varphi_P P(nT) + \Lambda_P; \\
F(nT^+) = 0. \\
B(nT^+) = 0.
\end{cases}
\qquad (5.22) $$

System (5.22) has the same properties as system (5.6) for the existence, uniqueness, non-negativity and boundedness of solutions. We pursue the mathematical analysis and simulations of system (5.22) without the $\dot{F}$ and $\dot{B}$ equations, since berry variable *B* is not present in the other equations.

### 5.3.1 Controlled periodic disease free solution and its stability

The controlled periodic disease free solution (cPDFS) occurs when $I = 0$ and $U = 0$. The healthy leaves dynamics are unchanged compared to the PDFS and still satisfy (5.10) while the hyperparasite dynamics are

$$ 
\begin{cases}
\dot{P} = -\mu_P P, \quad t \neq nT; \\
P(nT^+) = \varphi_P P(nT) + \Lambda_P, \quad t = nT.
\end{cases}
\qquad (5.23) $$

Solving equation (5.23) yields

$$ P((n+1)T) = P(nT^+)e^{-\mu_P T}. $$

Combined, the above equation and the impulsive condition give

$$ P((n+1)T^+) = \varphi_P P(nT^+)e^{-\mu_P T} + \Lambda_P. $$

Then, solving for the fixed point of the above discrete equation and using the dynamics of hyperparasites at the cPDFS yields

$$ P^T(t) = \frac{\Lambda_P}{1 - \varphi_P e^{-\mu_P T}} e^{-\mu_P (t-nT)}. \qquad (5.24) $$
Thus, the cPDFS for system (5.22) is $Y^T(t) = (S^T(t), 0, 0, P^T(t))$, where $S^T(t)$ and $P^T(t)$ are given in equation (5.10) and (5.24), respectively. We have the following result for the global stability of cPDFS.

**Lemma 5.2.** The cPDFS $Y^T(t) = (S^T(t), 0, 0, P^T(t))$ of system (5.22) is globally asymptotically stable when $\mathcal{R} < 1$.

*Proof of Lemma 5.2.* Let $\tilde{Y}(t) = Y(t) - Y^T(t)$ where $Y(t) = (S(t), I(t), U(t), P(t))$ and $\tilde{Y}(t)$ is understood as to be small amplitude perturbations. Substituting the expression of $\tilde{Y}(t)$ in the system (5.6) gives

$$ 
\left\{ 
\begin{aligned}
\dot{\tilde{S}} &= -\frac{\omega\nu\tilde{U}}{\tilde{N} + S^T(t)}(\tilde{S} + S^T(t)) - \mu\tilde{S}, & t \neq nT; \\
\dot{\tilde{I}} &= \frac{\omega\nu\tilde{U}}{\tilde{N} + S^T(t)}(\tilde{S} + S^T(t)) - (\mu + d)\tilde{I}, & t \neq nT; \\
\dot{\tilde{U}} &= \gamma\tilde{I} - (\nu + \mu_U)\tilde{U} - a(\tilde{P} + P^T(t))\frac{\tilde{U}}{K + \tilde{U}}, & t \neq nT; \\
\dot{\tilde{P}} &= ea(\tilde{P} + P^T(t))\frac{\tilde{U}}{K + \tilde{U}} - \mu_P\tilde{P}, & t \neq nT; \\
\tilde{S}(nT^+) &= \varphi_S\tilde{S}(nT); \\
\tilde{I}(nT^+) &= \varphi_I\tilde{I}(nT); \\
\tilde{U}(nT^+) &= \varphi_U\tilde{U}(nT); \\
\tilde{P}(nT^+) &= \varphi_P\tilde{P}(nT); 
\end{aligned}
\right. 
\qquad (5.25) 
$$

where $\tilde{N}(t) = \tilde{S}(t) + \tilde{I}(t)$ and $N^T(t) = S^T(t)$.

The linearisation of system (5.25) around $0_{\mathbb{R}^4}$ is

$$ 
\left\{ 
\begin{aligned}
\dot{\tilde{Y}}(t) &= G(t)\tilde{Y}(t), & t \neq nT; \\
\tilde{Y}(nT^+) &= \mathrm{diag}(\varphi_S, \varphi_I, \varphi_U, \varphi_P)\tilde{Y}(nT); 
\end{aligned}
\right. 
\qquad (5.26) 
$$

where

$$ G(t) = \begin{pmatrix} -\mu & 0 & -\omega\nu & 0 \\ 0 & -(\mu + d) & \omega\nu & 0 \\ 0 & \gamma & -(\nu + \mu_U) - \frac{aP^T(t)}{K} & 0 \\ 0 & 0 & \frac{eaP^T(t)}{K} & -\mu_P \end{pmatrix}. $$

Solving system (5.26) for $t \in (0, T]$ yields $\tilde{Y}(t) = \Phi_G(t)\tilde{Y}(0^+)$, with $\Phi_G(t)$ the fundamental matrix that satisfies

$$ \frac{d\Phi_G(t)}{dt} = G(t)\Phi_G(t), $$

where $\Phi_G(0) = I$. Re-ordering the variables as follows, $(\tilde{S}, \tilde{P}, \tilde{I}, \tilde{U})$, we note that the system is block-triangular. The stability is then determined by considering separately $(\tilde{S}, \tilde{P})$ and the pair $(\tilde{I}, \tilde{U})$. The Floquet multipliers for the first two are $\varphi_S e^{-\mu T}$ and $\varphi_P e^{-\mu_P T}$, which are smaller than one. With this in mind, the stability of system (5.26) reduces to the stability of the following sub-system

$$ 
\left\{ 
\begin{aligned}
\dot{\tilde{Y}}_1(t) &= G_1(t)\tilde{Y}_1(t), & t \neq nT; \\
\tilde{Y}(nT^+) &= \mathrm{diag}(\varphi_I, \varphi_U)\tilde{Y}(nT); 
\end{aligned}
\right. 
\qquad (5.27) 
$$
where

$$Y_1(t) = (I(t), U(t)) \quad \text{and} \quad G_1(t) = \begin{pmatrix} -(\mu + d) & \omega\nu \\ \gamma & -(\nu + \mu_U) - \frac{aP^T(t)}{K} \end{pmatrix}.$$

One can observe that

$$G_1(t) < B = \begin{pmatrix} -(\mu + d) & \omega\nu \\ \gamma & -(\nu + \mu_U) \end{pmatrix}.$$

The monodromy matrix associated to matrix B is the matrix $M_1$ given in equation (5.15), which is stable if $\mathcal{R} < 1$, where $\mathcal{R}$ is given in equation (5.17). Using the fact that system (5.27) is cooperative, we can conclude that, when $\mathcal{R} < 1$, the cPDFS $Y^T(t)$ is globally asymptotically stable, which implies that CLR dwindles until extinction. □

### 5.3.2 Semi-numerical analysis of the controlled model

In order to analyse system (5.22) when $\mathcal{R} > 1$, we need to compute the monodromy matrix associated with matrix $G_1(t)$. This is done by numerically solving the linear system (5.27) using the initial conditions $\tilde{Y}_1(0^+) = (1, 0)$ and $\tilde{Y}_1(0^+) = (0, 1)$. The two solutions evaluated at time $T$ are put together to obtain the fundamental matrix. Then we use the impulsive condition to obtain the monodromy matrix.

The monodromy matrix associated with $\tilde{Y}_1$ is

$$M_{G_1} = \begin{pmatrix} \varphi_I & 0 \\ 0 & \varphi_U \end{pmatrix} \Phi_{G_1}(T).$$

We analyse the spectral radius $\mathcal{R}_c$ of this matrix and check that it is smaller than one for the local stability of the cPDFS $Y^T(t)$ of system (5.22). We term this method “semi-numerical” as it requires numerical computations, but not extensive simulations of the system [56].

In Figure 5.3, we plot the threshold level of the spectral radius $\mathcal{R}_c = 1$ of controlled system (5.22) for parameter pairs $(\Lambda_P, a)$ (subplot (a)) and $(\nu, \gamma)$ (subplot (b)), as well as the threshold stability level $\mathcal{R} = 1$ of uncontrolled system (5.1) for parameter pair $(\nu, \gamma)$. In sub-plot (a), the regions below and above the blue curve are, respectively, the unstable and stable regions of the cPDFS of controlled system (5.22). As expected, stability of the cPDFS is guaranteed for large values of the yearly quantity of predators released $\Lambda_P$, and of the predator consumption rate $a$. Also, when the predator capacity to consume uredospores is higher, less predators are needed for the extinction of the disease, and vice versa. The blue dot corresponds to parameter values ensuring a stable cPDFS, values that are used for the blue curve in subplot (b).

In subplot (b), as in sub-plot (a), the blue curve separates the unstable and stable regions of the cPDFS of the controlled system; moreover, the black curve separates these two regions for the uncontrolled system. The parameter region between these two curves hence represents the gain obtained by adding the biocontrol: in this region, CLR goes extinct with biocontrol, but persists without. Subplot (b) shows that, in order to eliminate CLR through the proposed biocontrol effort, the sporulation and deposition rate should not be too large.
<figure><img src="image_24.png" /><figcaption>Figure 5.3: Stability regions of periodic disease free solutions as functions of model parameters. Blue curves represent the threshold $R_C = 1$ separating the stable and unstable regions of the model with predator (5.22), for different values of parameter pair $(\Lambda_P, a)$ in subplot (a) and $(\nu, \gamma)$ in subplot (b); moreover, the black curve $R=1$ in subplot (b) separates the stable and unstable regions of the model without predator (5.1). The blue dot in subplot (a) corresponds to the $(\Lambda_P, a)$ values used for the blue curve in subplot (b). Similarly, the grey dot in subplot (b) corresponds to the default $(\nu, \gamma)$ values used for the blue curve in subplot (a). Remaining parameter values are given in Table 5.1.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 5.4: Impact of biocontrol on the system dynamics. 6-year simulations of model (5.1) without predator (black curve) and model (5.22) with predator (blue curve), using initial conditions (5.28): (a) healthy leaves S; (b) infected leaves I; (c) uredospores U and (d) predator (hyperparasites) P. All parameter values are given in Table 5.1 and correspond to $R = 2.53 > 1$.</figcaption></figure>

### 5.3.3 Asymptotic behaviour for yearly releases

Herein, we investigate the asymptotic behaviour of system (5.22) when predators are released once a year at the beginning of each production period. To do so, we consider the parameter values given by the two dots in the stability regions of Figure 5.3. We point out that the gray dot in subplot (b) of Figure 5.3 is in the instability region of system (5.22) without control which means that the CLR would persist in the coffee plantation if unchecked. We consider the initial condition given in equation (5.21), adding the initial release of predators $P(0^+) = \Lambda_P = 3000$, which yields

$$ (S(0^+), I(0^+), U(0^+), P(0^+), F(0^+), B(0^+)) = (500, 0, 3000, 3000, 0, 0). \quad (5.28) $$

Figure 5.4 compares the temporal dynamics of system (5.1) without control (black curves) and model (5.22) with control (blue curves). From this figure, one can observe that under the actions of biocontrol, the number of healthy leaves increases and becomes stationary (see Figure 5.4(a)), while the number of infectious leaves and uredospores dwindle until extinction (see Figures 5.4(b) and (c)). During the first year, the quantity of predators (hyperparasites) increases through uredospore consumption. From the second year on, the number of uredospores is very close to zero, so that the solution is close to the cPDFS. Biocontrol using hyperparasites reduces drastically CLR in the coffee plantation.
## 5.4 Multiple releases of hyperparasites per year

In this section, we suppose that the yearly quantity $\Lambda_P$ of hyperparasites is released at different times over the year: it is uniformly divided into $m$ releases of size $\frac{\Lambda_P}{m}$, with release interval $\frac{T}{m}$. We thus obtain an impulsive model with a new switching condition for the predator, given by

Rainy season for $t \neq nT$:

$$ 
\left\{
\begin{aligned}
\dot{S} &= \Lambda - \omega\nu U \frac{S}{N} - \mu S; \\
\dot{I} &= \omega\nu U \frac{S}{N} - (\mu + d)I; \\
\dot{U} &= \gamma I - (\nu + \mu_U)U - aP \frac{U}{K+U}; \\
\dot{P} &= eaP \frac{U}{K+U} - \mu_P P, \quad t \neq nT + \frac{j}{m}T; \\
\dot{F} &= \delta_S S + \delta_I I - (\theta + \mu_F)F \\
\dot{B} &= \theta F - \mu_B B.
\end{aligned}
\right.
\qquad (5.29) \\
\text{Predator releases for } j \in \{1, \dots, m-1\}: \\
P \left(nT + \frac{j}{m}T^+\right) = P \left(nT + \frac{j}{m}T\right) + \frac{\Lambda_P}{m};
$$

Dry season:

$$ 
\left\{
\begin{aligned}
S(nT^+) &= \varphi_s S(nT); \\
I(nT^+) &= \varphi_I I(nT); \\
U(nT^+) &= \varphi_U U(nT); \\
P(nT^+) &= \varphi_P P(nT) + \frac{\Lambda_P}{m}; \\
F(nT^+) &= 0 \\
B(nT^+) &= 0.
\end{aligned}
\right.
$$

System (5.29) has the same properties as system (5.6) for the existence, uniqueness, non-negativity, and boundedness of solution. We made the mathematical analysis of system (5.29) without $\dot{F}$ and $\dot{B}$ equations, since B is not present in other equations.

### 5.4.1 Multiple release controlled periodic disease free solution and its stability

To compute the multiple release controlled periodic disease free solution (m-cPDFS), we only consider the first year for readability purposes, which means that $t \in (0, T]$. The m-cPDFS of biocontrol system (5.29) occurs when $I = 0$ and $U = 0$, so healthy leaf dynamics are unchanged compared to the PDFS and still satisfy (5.10); the hyperparasite dynamics are

$$ 
\left\{
\begin{array}{ll}
\dot{P} = -\mu_P P, & t \neq nT; \\
P\left(\frac{j}{m}T^+\right) = P\left(\frac{j}{m}T\right) + \frac{\Lambda_P}{m}, & \text{for } j = 1, \dots, m-1; \\
P(T^+) = \varphi_P P\left(\frac{(m-1)}{m}T\right) + \frac{\Lambda_P}{m}.
\end{array}
\right.
\qquad (5.30)
$$
Solving equation $\dot{P} = -\mu_P P$ for $t \in ]\frac{jT}{m}, \frac{(j+1)T}{m}]$ we obtain

$$P(t) = P\left(\frac{jT^+}{m}\right) e^{-\mu_P(t-\frac{jT}{m})}. \qquad (5.31)$$

Using equation (5.31) and the second equation of (5.30) we have

$$P\left(\frac{(j+1)}{m}T^+\right) = P\left(\frac{jT^+}{m}\right) e^{-\mu_P \frac{T}{m}} + \frac{\Lambda_P}{m}.$$

For $j = 0$

$$P\left(\frac{T^+}{m}\right) = P(0^+) e^{-\mu_P \frac{T}{m}} + \frac{\Lambda_P}{m}.$$

For $j = 1$

$$P\left(\frac{2T^+}{m}\right) = \left(P(0^+) e^{-\mu_P \frac{T}{m}} + \frac{\Lambda_P}{m}\right) e^{-\mu_P \frac{T}{m}} + \frac{\Lambda_P}{m}; \\ = P(0^+) e^{-\mu_P \frac{2T}{m}} + \frac{\Lambda_P}{m} e^{-\mu_P \frac{T}{m}} + \frac{\Lambda_P}{m}.$$

For $j = 2$

$$P\left(\frac{3T^+}{m}\right) = P(0^+) e^{-\mu_P \frac{3T}{m}} + \frac{\Lambda_P}{m} \sum_{i=0}^{2} e^{-\mu_P \frac{iT}{m}}.$$

Generally, for all $j = 0, \dots, m-1$

$$P\left(\frac{jT^+}{m}\right) = P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \sum_{i=0}^{j-1} e^{-\mu_P \frac{iT}{m}}; \\ = P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right),$$

in particular for $j = m - 1$, that is

$$P\left(\frac{(m-1)T^+}{m}\right) = P(0^+) e^{-\mu_P \frac{(m-1)T}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{(m-1)T}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right).$$

Using this latter expression and the third equation of (5.30), we obtain

$$P(T^+) = \varphi_P \left[ P(0^+) e^{-\mu_P T} + \frac{\Lambda_P}{m} \left( \frac{e^{-\mu_P \frac{T}{m}} - e^{-\mu_P T}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right] + \frac{\Lambda_P}{m}.$$

The value for which $P(T^+) = P(0^+)$ is

$$P(0^+) = \frac{\Lambda_P}{m} \left[ \frac{\varphi_P \left( e^{-\mu_P \frac{T}{m}} - e^{-\mu_P T} \right) + 1 - e^{-\mu_P \frac{T}{m}}}{(1 - \varphi_P e^{-\mu_P T}) \left( 1 - e^{-\mu_P \frac{T}{m}} \right)} \right].$$

So we can conclude that, for $t \in ]\frac{jT}{m}, \frac{(j+1)T}{m}]$, the m-cPDFS is given by $Z^T(t) =$
$(S^T(t), 0, 0, P_{m,j}^T(t))$, with $S^T(t)$ given by equation (5.10) and

$$P_{m,j}^T(t) = \left[ P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right] e^{-\mu_P (t - \frac{jT}{m})}, \quad (5.32)$$

$P(0^+)$ being defined above.

Then, the m-cPDFS is $Z^T(t) = (S^T(t), 0, 0, P_{m,j}^T(t))$, for $t \in ]\frac{jT}{m}, \frac{(j+1)T}{m}]$, where $S^T(t)$ is given by (5.10) and $P_{m,j}^T(t)$ by

$$P_{m,j}^T(t) = \left[ P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right] e^{-\mu_P (t - \frac{jT}{m})}, \quad (5.33)$$

with

$$P(0^+) = \frac{\Lambda_P}{m} \left[ \frac{\varphi_P \left( e^{-\mu_P \frac{T}{m}} - e^{-\mu_P T} \right) + 1 - e^{-\mu_P \frac{T}{m}}}{(1 - \varphi_P e^{-\mu_P T}) \left( 1 - e^{-\mu_P \frac{T}{m}} \right)} \right]. \quad (5.34)$$

**Lemma 5.3.** The m-cPDFS $Z^T(t) = (S^T(t), 0, 0, P_{m,j}^T(t))$ of controlled model with multiple releases, is globally asymptotically stable when $\mathcal{R} < 1$

*Proof of Lemma 5.3.* The proof is similar to the proof of Lemma 5.2 which holds for one release per year. □

To present the impact of the multiple release strategy on CLR control, we first compute a proxy of the control intensity over the season, namely the yearly average number of hyperparasites present in the plantation for the m-cPDFS. This quantity depends on $m$, the number of hyperparasite releases per year.

The yearly average number of hyperparasites for the m-cPDFS is computed over the interval $[0, T]$ as follows:

$$g(m) = \int_{0}^{T} P^{T}(t)dt = \sum_{j=0}^{m-1} \int_{\frac{jT}{m}}^{\frac{(j+1)T}{m}} P_{m,j}^{T}(t)dt.$$

Replacing $P_{m,j}^T(t)$ by its expression in equation (5.33), we obtain

$$\begin{align*}
g(m) &= \sum_{j=0}^{m-1} \int_{\frac{jT}{m}}^{\frac{(j+1)T}{m}} \left[ P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right] e^{-\mu_P (t - \frac{jT}{m})} (t) dt; \\
&= \sum_{j=0}^{m-1} \left[ P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right] \int_{\frac{jT}{m}}^{\frac{(j+1)T}{m}} e^{-\mu_P (t - \frac{jT}{m})} (t) dt; \\
&= \sum_{j=0}^{m-1} \left[ P(0^+) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right] \left[ \frac{1 - e^{-\mu_P \frac{T}{m}}}{\mu_P} \right].
\end{align*}$$
Replacing $P(0^+)$ by its expression in equation (5.34), we obtain

$$
\begin{aligned}
g(m) &= \sum_{j=0}^{m-1} \left[ \frac{\Lambda_P}{m} \left( \frac{\varphi_P \left( e^{-\mu_P \frac{T}{m}} - e^{-\mu_P T} \right) + 1 - e^{-\mu_P \frac{T}{m}}}{(1 - \varphi_P e^{-\mu_P T}) \left( 1 - e^{-\mu_P \frac{T}{m}} \right)} \right) e^{-\mu_P \frac{jT}{m}} \right] \\
&\quad + \sum_{j=0}^{m-1} \left[ \frac{\Lambda_P}{m} \left( \frac{1 - e^{-\mu_P \frac{jT}{m}}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \left[ \frac{1 - e^{-\mu_P \frac{T}{m}}}{\mu_P} \right], \\
&= \frac{\Lambda_P}{m \mu_P} \sum_{j=0}^{m-1} \left[ \frac{\varphi_P \left( e^{-\mu_P \frac{T}{m}} - e^{-\mu_P T} \right) + 1 - e^{-\mu_P \frac{T}{m}}}{(1 - \varphi_P e^{-\mu_P T})} e^{-\mu_P \frac{jT}{m}} + 1 - e^{-\mu_P \frac{jT}{m}} \right], \\
&= \frac{\Lambda_P}{m \mu_P} \left[ \left( \frac{\varphi_P \left( e^{-\mu_P \frac{T}{m}} - e^{-\mu_P T} \right) + 1 - e^{-\mu_P \frac{T}{m}}}{(1 - \varphi_P e^{-\mu_P T})} \right) \left( \frac{1 - e^{-\mu_P T}}{1 - e^{-\mu_P \frac{T}{m}}} \right) \right. \\
&\qquad \left. + m - \frac{1 - e^{-\mu_P T}}{1 - e^{-\mu_P \frac{T}{m}}} \right], \\
&= \frac{\Lambda_P}{m \mu_P} \left[ \left( \frac{(\varphi_P - 1) e^{-\mu_P \frac{T}{m}}}{(1 - \varphi_P e^{-\mu_P T})} \right) \left( \frac{1 - e^{-\mu_P T}}{1 - e^{-\mu_P \frac{T}{m}}} \right) + m \right].
\end{aligned}
$$

Simplifying this expression, we finally obtain

$$g(m) = \left( \frac{-\Lambda_P(1 - \varphi_P)(1 - e^{-\mu_P T})}{\mu_P(1 - \varphi_P e^{-\mu_P T})} \right) \frac{1}{m(e^{\mu_P \frac{T}{m}} - 1)} + \frac{\Lambda_P}{\mu_P}.$$

Then,

$$g(m) = \int_{0}^{T} P^{T}(t)dt = \sum_{j=0}^{m-1} \int_{\frac{jT}{m}}^{\frac{(j+1)T}{m}} P_{m,j}^{T}(t)dt = \frac{\Gamma}{m(e^{\mu_P \frac{T}{m}} - 1)} + \frac{\Lambda_P}{\mu_P}, \quad (5.35)$$

with

$$\Gamma = \frac{-\Lambda_P(1 - \varphi_P)(1 - e^{-\mu_P T})}{\mu_P(1 - \varphi_P e^{-\mu_P T})}.$$

Differentiating $g(m)$ with respect to $m$, one has

$$g'(m) = \frac{\Gamma e^{\mu_P \frac{T}{m}} \left( e^{-\mu_P \frac{T}{m}} + \frac{\mu_P T}{m} - 1 \right)}{m^2 (e^{\mu_P \frac{T}{m}} - 1)^2}.$$

Since $e^{-\mu_P \frac{T}{m}}$ is always larger than its first order approximation $1 - \frac{\mu_P T}{m}$ for $\frac{T}{m} > 0$, all factors are positive except $\Gamma$, so that $g'(m) < 0$. This allows to conclude that the yearly average number of hyperparasites decreases with the number of releases per year, without CLR (quantity computed for the m-cPDFS). One might then think that spreading $\Lambda_P$ over several releases is less efficient than releasing everything at the beginning of the season. This is investigated below through simulation.

We use two contrasted predator mortality rates to simulate the solution $P_{m,j}^T(t)$ in the m-cPDFS. Subplots (a) and (b) of Figure 5.5 present the simulation of $P_{m,j}^T(t)$ for various release frequencies and for both mortality rates over three years; they illustrate the important impact these two factors on the m-cPDFS.
<figure><img src="image_26.png" /><figcaption>Figure 5.5: Impact of the predator mortality and its release frequency on its dynamics without CLR. The plots present the predator dynamics of the m-cPDFS (multiple release controlled periodic disease free solution) of system (5.29) for: (a) the default low mortality rate $\mu_P = 0.003/day$; and (b) a high mortality rate $\mu_P = 0.1/day$. Various release frequencies are plotted: $m = 1$ (blue curve), $m = 2$ (brown curve), $m = 5$ (magenta curve) and $m = 50$ (light blue curve) releases per year. Initial condition is $P(0^+)$ given in (5.34), so the solutions are periodic. The yearly released quantity is $\Lambda_P = 3000$ predators. Remaining parameter values are given in Table 5.1.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 5.6: Impact of release frequency on the stability for low predator mortality rate $\mu_P = 0.003$/day. Subplot (a) presents the stable and unstable regions of the m-cPDFS of system (5.29) according to parameter pair $(\Lambda_P, a)$. The regions are separated by $\mathcal{R}_{c,m} = 1$ for $m = 1$ (blue curve) or $m = 5$ (magenta curve) releases/year. Subplot (b) shows the dynamics of infected leaves $I$ for parameter values $(\Lambda_P = 3000$ predators, $a = 0.6$ spores/predator.day) corresponding to the purple dot in subplot (a). Remaining parameter values are given in Table 5.1 and initial conditions by (5.36).</figcaption></figure>

### 5.4.2 Impact of release frequency on CLR control

Herein, we present the impact of the yearly hyperparasite release frequency on CLR control. This is first done using the semi-numerical analysis of the model with multiple releases, to show the $(\Lambda_P, a)$-stability region of the m-cPDFS (as in Figure 5.3), for different values of the release number $m$. Secondly, we present simulations of the dynamic behaviour of system (5.29).

In all simulations, the epidemiological parameters are set to their value in Table 5.1. The only parameters that vary are the ones related to the hyperparasites: $\Lambda_P, a, \mu_P$ and $m$. As in Figure 5.5, we use two contrasted values for the predator mortality $\mu_P$, the low and high values representing respectively a favourable and unfavourable environment for the hyperparasites.

Moreover, we consider that an initial hyperparasite release occurs at time 0 with the usual initial conditions (5.21) for the other variables, so that

$$ (S(0^{+}), I(0^{+}), U(0^{+}), P(0^{+}), F(0^{+}), B(0^{+})) = (500, 0, 3000, \frac{\Lambda_P}{m}, 0, 0) \quad (5.36) $$

#### 5.4.2.1 Impact on stability

Herein, we present the impact of the yearly release frequency on the m-cPDFS stability of controlled system (5.29), for the two contrasted predator mortality rates depicted in Figure 5.5. Only two release frequencies are considered, $m = 1$ and $m = 5$ releases per year, as results obtained for $m = 5$ also hold for higher frequencies.

**Case 1** – Less frequent is more effective (low mortality rate $\mu_P = 0.003$/day) In subplot (a) of Figure 5.6, $\mathcal{R}_{c,1} = 1$ for $m = 1$ (blue curve) and $\mathcal{R}_{c,5} = 1$ for $m = 5$ (magenta curve) separate the instability (below) and stability (above) regions of the m-cPDFS of controlled system (5.29), when parameters $(\Lambda_P, a)$ vary. The stability region is greater for $m = 1$ than for $m = 5$, that is when the yearly release frequency is lower.
<figure><img src="image_28.png" /></figure>

Figure 5.7: Impact of release frequency on the stability for high predator mortality rate $\mu_P = 0.1$/day. Subplot (a) presents the stable and unstable regions of the m-cPDFS of system (5.29) according to parameter pair ($\Lambda_P, a$). The regions separated by $\mathcal{R}_{c,m} = 1$ for $m = 1$ (blue curve) or $m = 5$ (magenta curve) releases/year. Subplot (b) shows the dynamics of infected leaves $I$ for parameter values ($\Lambda_P = 30000$ predators, $a = 1.5$ spores/predator.day) corresponding to the purple dot in subplot (a). Remaining parameter values are given in Table 5.1 and initial conditions by (5.36)

Subplot (b) of Figure 5.6 illustrates the long term dynamics (over 12 years) of the infected leaves of system (5.29), with initial conditions (5.36), for parameter values $\Lambda_P = 3000$ predators and $a = 0.6$ spores/predator.day, which correspond to the purple dot in subplot (a). As expected from subplot (a), CLR goes extinct for one release per year (blue curve), while it persists for five releases per year (magenta curve).

Therefore, if the hyperparasite mortality is very low, the best strategy for eradicating CLR is to release hyperparasites once a year at the beginning of the season, implying that releasing less frequently in larger quantity is the most effective. This is consistent with the control intensity proxy computed in equation (5.35), which is higher for lower frequencies.

**Case 2 – More frequent is more effective (high mortality rate $\mu_P = 0.1$/day)** Figure 5.7 is built as Figure 5.6, but for a high predator mortality. In subplot (a), as opposed to what is observed for low mortality, the stability region is greater for $m = 5$ than for $m = 1$, that is when the yearly release frequency is higher.

Subplot (b) is plotted with parameter values $\Lambda_P = 30000$ predators and $a = 1.5$ spores/predator.day, which correspond to the purple dot in subplot (a). As expected from subplot (a), CLR extinction is achieved for five releases per year but not for one release per year. Infected leaves $I$ behave quite differently. For $m = 1$, $I$ decrease at the beginning of each year and then increase sharply; indeed, the hyperparasite population is large at the beginning of the year, but declines sharply due to the high mortality rate (see Figure 5.5(b)), so that CLR is not kept in check. By contrast, for $m = 5$, hyperparasites do not stay close to 0 for long so that no large increase of $I$ occurs.

For high predator mortality, the best strategy to control CLR is to release hyperparasites five times per year, meaning that more frequent releases are better. This conclusion was not obvious, considering that the control intensity proxy computed in equation (5.35) is lower for more frequent releases.
<figure><img src="image_29.png" /></figure>

Figure 5.8: Impact of release frequency on the transient dynamics with low predator mortality. Subplot (a) presents the stable and unstable regions of the m-cPDFS of system (5.29) according to parameter pair $(\Lambda_P, a)$. The regions are separated by $\mathcal{R}_{c,m} = 1$ for $m = 1$ (blue curve), $m = 2$ (brown curve), $m = 5$ (magenta curve) or $m = 50$ (light blue curve) releases/year. Subplots (b) and (c) show the dynamics of infected leaves I and berries B for parameter values $(\Lambda_P = 3000$ predators, $a = 0.8$ spores/predator.day) corresponding to the purple dot in subplot (a). Remaining parameter values are given in Table 5.1 and initial conditions by (5.36).

### 5.4.2.2 Impact on transient dynamics

Herein, we still look at the impact of the yearly release frequency, for the same contrasted predator mortality rates, but we focus on the first years and assume that the m-cPDFS of controlled system (5.29). Four frequency values are considered, namely $m = 1, 2, 5, 50$ releases per year.

**Case 1 – Less frequent is faster (low mortality rate $\mu_P = 0.003$/day)** Just as Figure 5.6(a), subplot (a) of Figure 5.8 displays the instability (below) and stability (above) regions of the m-cPDFS of controlled system (5.29). These regions are separated by $\mathcal{R}_{c,1} = 1$ for $m = 1$ release/year (blue curve), $\mathcal{R}_{c,2} = 1$ for $m = 2$ releases/year (brown curve), $\mathcal{R}_{c,5} = 1$ for $m = 5$ releases/year (magenta curve) and $\mathcal{R}_{c,50} = 1$ for $m = 50$ releases/year (light blue curve). Subplot (a) confirms that for low hyperparasite mortality, the best strategy for controlling CLR is to release the hyperparasites once a year at the beginning of the season.

Figure 5.8 also shows the temporal evolution of infected leaves (subplot (b)) and berries (subplot (c)), using initial conditions (5.36) and parameter values $\Lambda_P = 3000$ predators and $a = 0.8$ spores/predator.day. These parameter values, which correspond to the purple dot in subplot (a), ensure that the m-CPDFS is stable for all four values of the release frequency $m$. In subplot (b), we can observe that the disease is driven to extinction for $m = 1$ (blue curve). For $m = 5$ (magenta curve) and $m = 50$ (light blue curve), the number of infected leaves after six years is still notably greater than zero (less notably so for $m = 2$, brown curve). However, if we increase the number of years, the extinction will be observed for all values of $m$.

Subplot (c) shows that the number of berries at the end of each year decreases when $m$ increases. This confirms our previous result: the best strategy for low hyperparasite mortality is
<figure><img src="image_30.png" /></figure>

Figure 5.9: Impact of release frequency on the transient dynamics with high predator mortality. Subplot (a) presents the stable and unstable regions of the m-cPDFS of system (5.29) according to parameter pair $(\Lambda_P, a)$. The regions are separated by $\mathcal{R}_{c,m} = 1$ for $m = 1$ (blue curve), $m = 2$ (brown curve), $m = 5$ (magenta curve) or $m = 50$ (light blue curve) releases/year. Subplots (b) and (c) show the dynamics of infected leaves I and berries B for parameter values $(\Lambda_P = 30000$ predators, $a = 2.5$ spores/predator.day) corresponding to the purple dot in subplot (a). Remaining parameter values are given in Table 5.1 and initial conditions by (5.36).

to few releases per year.

**Case 2 – Best strategy depends on the horizon (high mortality rate $\mu_P = 0.1$)** As expected for high hyperparasite mortality, in subplot (a) of Figure 5.9, the $\mathcal{R}_{c,m} = 1$ curves, delimiting the stable and unstable regions of the m-cPDFS, are in reverse order compared to Figure 5.8(a). However, we should note that $\mathcal{R}_{c,5} = 1$ (magenta curve) and $\mathcal{R}_{c,50} = 1$ (light blue curve) are almost indistinguishable, implying that both strategies are almost equivalent in terms of eradication success.

Subplots (b) and (c) of Figure 5.9 are built using parameter values $\Lambda_P = 30000$ predators and $a = 2.5$ spores/predator.day, which correspond to the purple dot in subplot (a)). In subplot (b), we observe CLR extinction in all cases but for $m = 1$, for which extinction takes longer. In the latter case, the high mortality rate drives the hyperparasite population to very small values very quickly so that it is absent for most of the season, which slows down the control process.

To better explain the impact of the high hyperparasite mortality on berries depicted in subplot (c), we consider the number of berries at the end of each season and thus we obtain Table 5.2. The majority of parameter values of the model are certainly assumed, but this qualitative study shows the effectiveness of the biocontrol of CLR using hyperparasite

In terms of berry production, the best strategy from the first year to the second year is to release the hyperparasites once. From year 3 to year 4, the highest berry production occurs with $m = 5$ releases per year. From year 5 and on, the production increases with the release frequency, which is the result expected in the case of high predator mortality. However, except for $m = 1$, the berry production numbers are then almost all identical since CLR is nearly eradicated. Therefore, the best release strategy depends on the time horizon considered.
<table><caption>Table 5.2: Berry production at the end of each season, for different yearly release frequencies and high predator mortality. Production values correspond to subplot (c) of Figure 5.9. The optimal production for each year is enhanced in bold.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Releases per year</th></tr><tr><th>m = 1</th><th>m = 2</th><th>m = 5</th><th>m = 50</th></tr></thead><tbody><tr><td>year 1</td><td><strong>5845</strong></td><td>5554</td><td>5131</td><td>4717</td></tr><tr><td>year 2</td><td><strong>7906</strong></td><td>7756</td><td>7600</td><td>7290</td></tr><tr><td>year 3</td><td>8663</td><td>8866</td><td><strong>8985</strong></td><td>8897</td></tr><tr><td>year 4</td><td>8987</td><td>9399</td><td><strong>9568</strong></td><td>9553</td></tr><tr><td>year 5</td><td>9159</td><td>9657</td><td>9793</td><td><strong>9796</strong></td></tr><tr><td>year 6</td><td>9270</td><td>9785</td><td>9876</td><td><strong>9880</strong></td></tr></tbody></table>

### 5.4.2.3 Summary

From the four cases described above, we can conclude that the best hyperparasite release strategy, in terms of CLR control and berry production, strongly depends on hyperparasite characteristics, namely its mortality rate ($\mu_P$) and its uredospore consumption rate (a). Both parameters determine the yearly released quantity of hyperparasites ($\Lambda_P$) needed to control CLR. Moreover, the mortality rate has an impact on the best hyperparasite release frequency (m). On the long run, more frequent is more effective for high hyperparasite mortality, and vice versa for low mortality. On the shorter run however, less frequent may also be better for high mortality.

## 5.5 Discussion

In this work, we built an original impulsive model of coffee leaf rust dynamics in a coffee plantation, in which the non production (dry) season takes the form of an impulse, it is based in a previous study [38]. There are few CLR models in the literature and they represent different geographical scales, from the individual coffee bush to the country or even the continent. Bebber et al. determined the germination and infection risk depending on the climate in Colombia and neighbouring countries, based upon existing experimental data [15]. In contrast to these static approaches, Vandermeer et al. studied the interaction between the regional and local dynamics of CLR model by representing the evolution of the proportion of infected bushes and farms [169]. Vandermeer et al. also represented the CLR dynamics in a coffee farm in Chiapas using an SI (susceptible–infected) epidemiological model of the host and concluded that the network approach can be a useful way of gaining qualitative insight about disease dynamics in space [165].

The originality of our model, compared to the dynamical approaches above, is the hybrid formalism that is particularly suited to describe the seasonality of coffee and rust dynamics, as well as impulsive biocontrol. Floquet theory was used to compute stability thresholds for the periodic disease free solution. The initial coffee leaf rust model was then coupled with the dynamics of predator-like hyperparasites, as recommended by Zambolin [180] for alternative biocontrol methods; this resulted in a controlled impulsive model. The semi-numerical analysis of this controlled impulsive model allowed to conclude that hyperparasite-based biocontrol can drastically reduce coffee leaf rust. Moreover, we studied how successful biocontrol implementation depends on hyperparasite characteristics, in particular its mortality.

We set the study in a framework where the yearly released quantity of hyperparasites is fixed and only its release frequency varies. We deemed this comparison relevant since, whatever the
release strategy, the same yearly budget may be allocated to biocontrol agents.

Low and high hyperparasite mortality rates were considered to represent different environmental situations and hyperparasite characteristics. High mortality occurs if the hyperparasite is a specialist parasite like *Mycodiplosis* [67], which means that it only consumes CLR uredospores, or if it is a generalist hyperparasite whose alternate food sources are absent; in both cases, the hyperparasite quickly dies out in absence of CLR. Low mortality occurs for generalist hyperparasites in a rich environment. Most hyperparasites of *H. vastatrix* in the literature have alternative preys. For example, the *Lecanicillium lecanii* fungus is a CLR hyperparasite that could be used for biocontrol [168, 143], and it is also a pathogen of plant-parasitic nematodes [57]. Hence, both low and high mortality situations could occur in the field.

For low mortality rates, the best strategy is to release the hyperparasites once a year at the beginning of the production season. This strategy is efficient because it is the one that ensures the largest average hyperparasite density over the year, hyperparasites being always present because of the low mortality rate. However, the cost of such a large release at the beginning of the production season, on top of all the actions that need to be taken for coffee farming at that same time (such as fertiliser and manure application, or trimming of dry or excess leaves), can put too large a financial burden on the coffee producer. This problematic represents a well-known trade-off between yield and affordability in farming practices. Another potential problem is that massively advising farmers to adopt these practices could cause an excessive increase in the local demand of hyperparasites, and consequently a shortage in the market. For both these reasons, a strategy considering two releases could be recommended, even though it is less cost-efficient.

When the hyperparasite mortality rate is high, the best strategy depends on the horizon. For the first year, the most efficient approach is to have a large release at the beginning of the season in order to have an immediate massive impact. For the subsequent years, 5 and then 50 releases per year are advocated in Table 5.2. However, there is little difference in berry production between 5 and 50 releases in years 5 and 6. Moreover, the eradication success is almost identical for both strategies in terms of yearly released quantity of hyperparasites needed to control CLR (see Figure 5.9(a)). Considering the burden and cost associated with having 50 interventions over the production season, five yearly releases should be preferred to 50. Releasing predators twice a year may appear slightly suboptimal in terms of eradication success. It is however fairly close to the berry production achieved when releases occur five times per year and it is less labour consuming. In terms of cost-effectiveness, two releases per year might then be a good compromise. This result is complementary to the results of Henk et al [67], which shows that the larvae of some *Mycodiplosis* species can feed on spores of rust fungi and their frequency play an important role.

In the present work, we have focused our attention on one family of strategies: those that repeat themselves every year and are evenly spread over the production season, with a first release at the beginning of each year. Potential generalisations are numerous, but could not be considered here. One could for example choose to change strategy every year; this would seem to be a good idea if we look at Table 5.2, where we would choose to have a large single release the first year, and potentially five releases per year after that. One could also choose not to have a first release at the beginning of the year, but later in the season, or to have a larger release at the beginning of the year and then some smaller ones during the year. Finally, hyperparasites could be released in reaction to peaks of uredospore densities in the plantation, so as to maximise their impact; the latter approach however requires close monitoring of the disease in the plantation.
In addition to them potentially being more efficient, these approaches have the advantage of not imposing the full burden of acquiring hyperparasites at the beginning of the season.

The continuous part of our model contains most elements that we deem relevant, keeping it reasonably small and tractable. However, uredospore germination and mycelium formation depend on various environmental factors [10], such as temperature [35]. Introducing meteorological data in our model would allow for a better description of the day-to-day dynamics, but constant parameters gives a fairly good representation of the average dynamics for our purpose.

The hybrid model we developed is able to describe the dynamics of other fungal diseases that attack leaves of perennial plants. It appears as a promising tool to explore the efficiency of biocontrol strategies based on the impulsive release of hyperparasites. The impulsive component of the model presented in this work corresponds to plants that are harvested once a year; such is the case in temperate climate or in tropical climate with dry and rainy seasons. However, this model can be perfectly adapted to other countries without such contrasted seasons, where there are more than one harvest per year. In this case, several impulses can occur over the year, corresponding to partial harvests only, and not seasonality. The conclusions that we have drawn for CLR might not translate to other pathogen-plant pairs in other contexts, beyond our analytical results; a full semi-analytical study and numerical simulations should be performed for each case.
# GENERAL CONCLUSION

## Contents

<table><tr><td>6.1 Summary of main results</td><td>119</td></tr><tr><td>6.2 Perspectives</td><td>120</td></tr></table>

## 6.1 Summary of main results

This thesis focused on the mathematical modelling and biocontrol of the dynamics of fungal diseases applied to the particular case of CLR and its impact on coffee production. The main objective was to develop rigorous mathematical models of plant-pathogen interactions, and also to develop an efficient and long-lasting control strategy that reduces the fungal development and preserves the yield. We reviewed literature about the impact of pest and fungal diseases, the biology of *H. vastatrix* and coffee tree in Chapter 1. This chapter also presents existing mathematical models of plant pathogens and particularly of CLR. In Chapter 2, we presented the mathematical tools used in the thesis. The main contributions of the work are found in Chapters 3, 4 and 5. Chapter 3 is based on a manuscript in revision in *Journal of Mathematical Biology*. Chapter 4 is built on a paper published in *Mathematical Modelling of Natural Phenomena* [38]. Chapter 5 is a manuscript in revision in *Mathematical Biosciences*, which is an extension of a paper published in *IFAC-PapersOnLine* [37] and given in Appendix A.

In our first contribution, presented in Chapter 3, we proposed and analysed a mathematical model that describes the propagation of a fungal disease on a perennial plant using an ODE model. The model is based on classical epidemiological models and its originality is the differentiation between young and mature leaves. Our aim was to understand the host-pathogen interactions and study the impact of host development on the aggressiveness of the pathogen. We proved that the model exhibits complex asymptotic properties, that differ from classical epidemiological models: a basic reproduction number smaller than one is insufficient to eradicate the disease as it could persist even though it could not invade, and a basic reproduction number greater than one does not guarantee the persistence of the disease, as large values could result in the destruction of the plantation.

Our second contribution, presented in Chapter 4, consisted in proposing and analysing a PDE model describing the disease propagation in the plantation, differentiating the rainy and dry seasons. The main goal was to understand the impact of seasonality on the propagation
of the fungal disease in order to propose a biocontrol solution, based on a hyperparasite that
inhibits the fungus reproduction. We computed the equilibria of the two subsystems, defined
during the rainy and dry seasons. We proved that the basic reproduction numbers during the two
seasons could determine the dynamics of the full model: when the basic reproduction number is
less than one during the rainy season, then CLR globally decreases until extinction; when it is
greater than one for the dry season, then CLR persists. We showed that a rather high biocontrol
efficiency (75% at least) is necessary in our model to control the disease, but lower efficiencies
still improve coffee production notably.

In our last contribution, presented in Chapter 5, we developed a multi-seasonal model, alter-
nating continuous crop-fungus dynamics during the rainy seasons and discrete events to represent
the simpler dynamics during the dry seasons. The originality of our model, compared to the dy-
namical approaches above, is the hybrid formalism that is particularly suited to describe the
seasonality of coffee and rust dynamics, as well as impulsive biocontrol. Our goal was to develop
a biocontrol strategy based on the impulsive release of hyperparasites in the plantation. We used
semi-numerical analysis to show that the best hyperparasite release strategy strongly depends on
the hyperparasite characteristics, namely its mortality rate and its uredospore consumption rate.
Both parameters determine the yearly released quantity of hyperparasites needed to control the
disease. Moreover, we showed that the mortality rate has an impact on the best hyperparasite
release frequency. On the long run, more frequent is more effective for high hyperparasite mor-
tality, and vice versa for low mortality. On the shorter run however, less frequent may also be
better for high mortality.

## 6.2 Perspectives

We identified two main extensions for this modelling work: the first one based on intra-guild predation and the second one on the deployment of resistant cultivars.

In this thesis, the biological control we proposed to control coffee leaf rust focused on the
deployment of hyperparasites that prey on *H. vastarix* uredospores, such as fungus *L. lecanii*
[168] and *Mycodiplosis* insect [139]. According to Chapter 5, the mortality rate of hyperparasites
determines how often they need to be released in the plantation. The growth of *L. lecanii* is
favoured by the same climatic conditions as *H. vastatrix*. Moreover *L. lecanii* is a generalist [48],
so it can survive in the plantation in the absence of *H. vastarix* uredospores. Its mortality should
hence be relatively low in a coffee plantation affected by coffee leaf rust. Based on our findings,
it could then be released in the plantation at a low frequencies. However, because *L. lecanii* is a
generalist, it may select other preys, such as green coffee scales [166], decreasing the effectiveness
of the control. The *Mycodiplosis* larva is a specialist, that only feeds on *H. vastarix* uredospores
[67]. However, these larvae are eaten by the *Azteca sericeasur* ant [63, 89], so their mortality may
be significant. In that case, according to our findings, they should be released fairly frequently.
Based on field studies in Puerto Rico and Mexico, Hajian-Forooshani at al. suggested that the
joint presence of both *L. lecanii* and *Mycodiplosis* may provide control for *H. vastarix*, with
varying outcomes depending on the community composition. It would be interesting to extend
our work and investigate intra-guild predation, including the dynamics of both predators in the
coffee leaf rust model.

In this thesis, we did not make a difference between the coffee varieties arabica and robusta.
The arabica variety represents 70% of the global production, while robusta accounts for less
than 30%. Indeed arabica has a higher agronomic value, but it is more susceptible to CLR
than robusta. Therefore, developing CLR resistant coffee varieties with a suitable agronomic
value is a current challenge. Most resistant varieties are obtained by combining arabica and
robusta varieties, some are listed in Subsection 1.2.3. Several genes providing resistance to
strains of the fungus *H. vastatrix* were identified in several coffee varieties [128, 129]. However,
the fungus *H. vastatrix*, like many others, is adapting to these resistant varieties, necessitating the
development of new resistant plants or management strategies that reduce or slow the resistance
breakdown. Many modelling studies proposed control methods based on resistant plants, mixing
or alternating resistant and susceptible plants in space and/or time to preserve the resistance
efficiency on the long run [133, 45]. Indeed, the extensive use of resistant plants increases the
selection pressure on virulent pathogens, that can overcome the plant resistance. However,
reducing the use of resistant plants may face acceptability issues, as resistant plants provide
better returns in the short term [114]. Moreover, strategies are better suited for annual crops.
Replacing coffee plants is particularly challenging for farmers because it takes an average of four
years to produce coffee. It is hence particularly interesting to gain some knowledge based on
theoretical studies in the case of coffee leaf rust, to help the deployment of resistant varieties
with mixtures of plants carrying different resistances the most promising approach for perennial
plants..
# IFAC PAPER

## Contents

A.1 Optimal control of the impulsive model of CLR 123

A.1.1 Introduction 124

A.1.2 The model and its basic properties 125

A.1.3 Mathematical analysis 126

A.1.4 Optimal control 132

A.1.5 Conclusion 136

This work is published in IFAC-PapersOnLine [37]

## A.1 Optimal control of the impulsive model of the propagation of the coffee leaf rust

**Abstract:** This paper presents the dynamics and optimal control of the coffee leaf rust (CLR) over the coffee plantation with impulsive effect, which represents the non production period of coffee. Our aim is to estimate the damage cause by the CLR and the optimal control the propagation of the CLR during many years using the biological control. We show the existence of solution and prove using Floquet theory that there exist the threshold calling basic reproduction number, that determine the stability of solutions. Numerical simulations are performed to illustrate the theoretical results. This paper deals with the problem of modelling and optimal control of the coffee leaf rust (CLR) in a coffee plantation. We first formulate and analyze a mathematical model for the propagation of CLR with impulsive effect, which represents the non-production period of coffee. Then, the CLR control is formulated and solved as an optimal control, showing how a biological control strategy can be implemented in a coffee plantation to reduce the damage caused by CLR, therefore increasing the coffee production. Numerical simulations are performed to illustrate and validate the theoretical results.

**Keywords:** Epidemiological modeling; Impulsive system; Floquet theory; Stability; Optimal control, *Hemileia vastatrix*.
A.1.1 Introduction

Coffee is one of the most widely consumed beverages in the world, with their consumption in 2019/20 being estimated at 168 million 60-kilogram bags [72]. Its cultivation is widespread in almost all the countries of the tropic and is an important factor of social stability as it supports no less than twenty-five million small producers and their families worldwide. However, due to the climate of these areas, multiple coffee diseases have persisted for several years. Among these pests, the Coffee Leaf Rust (CLR) is considered the most important leaf disease known to this day for being one of the most devastating pathogens. It is caused by the fungus *Hemileia vastatrix*, first discovered in Sri Lanka in 1869, and nowadays present in the majority of the coffee-growing regions of the world. The pathogenic principle behind CLR is that it causes premature defoliation, which reduces the photosynthetic capacity of coffee leaves, inducing severe yield losses to the sector (up to 70%).

The dynamics of CLR is complex due to the multiple interactions between the plant host and the pathogen in the environment. A better understanding of such CLR dynamics can significantly improve the accuracy of the models, which can be then analyzed through theoretical studies and numerical simulations. Many efforts have been and are still being devoted to the modeling of the propagation of CLR [165]. However, none of them take into account the impulsive effect in the dynamics. Considering discrete events in the model—such as harvesting—could allow for a more comprehensive representation of the plantation process, providing in turn a more precise framework for pest treatment. In this regard, several works have used ordinary impulsive differential equations for representing population dynamics in ecology and epidemiology in the recent years [96, 154, 156, 40].

The capacity to accurately model the time evolution of coffee plantations allows not only to effectively predict trends in the process, but also to act on them. In this context, optimal control theory has proven to be a powerful tool for investigating potential control strategies in pest treatment [1, 178, 18]. Such approaches are commonly based on the well-known PMP (Pontryagin's maximum principle) [125]: for a given cost function to maximize, this theory can provide necessary (and often sufficient) conditions for optimality of control strategies in systems of ordinary differential equations, partial differential equations and hybrid systems with given constraints [39].

The present contribution offers an impulsive perspective to the modeling of the propagation of CLR in a coffee plantation. We first formulate and analyze a non-controlled impulsive model for the propagation of CLR, for which we compute the periodic disease-free solution (PDFS) and a threshold parameter $R$. Then, we conduct an in-depth analysis of the global asymptotic stability of the PDFS when $R < 1$, and we perform numerical simulations to illustrate the theoretical results. Furthermore, we devise an extended model that considers a biological control for pest treatment, intended to maximize the coffee production while minimizing the damages undergone by the CLR. The latter objective is written as an Hybrid OCP (Optimal Control Problem), for which we propose an analogous problem in continuous form, which allows to obtain numerical solutions with Bocop [20]. Later on, we provide a numerical solution of an optimal trajectory of the system for a relevant set of parameters and initial conditions. As a result, we show that a cost-effective inoculative biological control can successfully contain the damages on a coffee plantation.
## A.1.2 The model and its basic properties

### A.1.2.1 Mathematical Modeling of coffee leaf rust

In this section, we formulate a mathematical model for the propagation of the CLR in the coffee plantation. To do so, we consider the dynamics of fungus during the production season—which corresponds to the rainy season for some countries—through ordinary differential equations; and during the non-production period—which is the dry season for some countries—represented by the impulsive effect. More precisely, we consider that in a coffee plantation we can find: susceptible branches S, infected branches I and urediniospores U. During the production period with length $T > 0$, the recruitment of healthy branches occurs at rate $\Lambda$. Urediniospores are deposited on leaves of all branches at rate $\nu$, and a fraction $S/N$ lands on susceptible branches. When the conditions are favorable, urediniospores germinate with efficiency rate $\omega$ and the susceptible branches becomes infected branches. All branches undergo natural mortality with baseline rate $\mu$ and the infected branches have an additional mortality rate $d$ due to the disease. Urediniospores are produced by infected branches at rate $\gamma$ and lose their ability to infect coffee branches at constant rate $\mu_U$. Susceptible and infected branches produce the berries with constant rate $\delta_S$ and $\delta_I$, respectively. At the end of the production period, harvest occurs instantaneously and we consider that the dry season or non-production period corresponds to the time when the number of branches at the end of the production period switches, this switch due to the fact that harvesting reduces the number of susceptible branches, infected branches and urediniospores with rates $\phi_S$, $\phi_I$ and $\phi_U$ respectively. Using the fact that urediniospores lose their ability to infect fastly, we assume that $\phi_U$ is very close to 0. The model is given by the following hybrid system:

$$ 
\left\{
\begin{aligned}
\dot{S} &= \Lambda - \frac{\omega\nu U}{N}S - \mu S, \quad t \neq nT; \\
\dot{I} &= \frac{\omega\nu U}{N}S - (\mu + d)I, \quad t \neq nT; \\
\dot{U} &= \gamma I - (\nu + \mu_U)U, \quad t \neq nT; \\
\dot{B} &= \delta_S S + \delta_I I - \mu_B B, \quad t \neq nT; \\
S(nT^+) &= \phi_S S(nT); \\
I(nT^+) &= \phi_I I(nT); \\
U(nT^+) &= \phi_U U(nT); \\
B(nT^+) &= 0.
\end{aligned}
\right.
\quad (A.1) $$

where $S(t)$, $I(t)$, $U(t)$ and $N(t) = S(t) + I(t)$ represent the number of susceptible branches, infected branches, urediniospores and total branches, respectively at time $t$. With $0 < \phi_I \le \phi_S < 1$ and $\phi_U > 0$ but close to 0.

Table A.1 summarizes the biological meaning and the values of parameter of model (A.1).

### A.1.2.2 Basic properties of the model

Let $\mathbb{R}_+^4 = \{X \in \mathbb{R}^4, X \ge 0\}$. Denote by $f = (f_1, f_2, f_3, f_4)$ the map given by the right-hand side of system (A.1), and consider the initial conditions

$$ (S(0), I(0), U(0), B(0)) \in \mathbb{R}_+^4, S(0) > 0. \quad (A.2) $$
<table><caption>Table A.1: Description and values of parameters for system (A.1)</caption><thead><tr><th></th><th>Biological meaning</th><th>Value</th></tr></thead><tbody><tr><td>Λ</td><td>Recruitment of S</td><td>50/day</td></tr><tr><td>ω</td><td>Inoculum effectiveness</td><td>6.5% [131]</td></tr><tr><td>μ</td><td>Mortality rate of branches</td><td>0.0034/day</td></tr><tr><td>γ</td><td>Production rate of spores</td><td>∈ [0, 20]/day</td></tr><tr><td>d</td><td>Mortality rate due to CLR</td><td>0.056/day</td></tr><tr><td>ν</td><td>Deposition rate</td><td>0.09/day [19]</td></tr><tr><td>μ<sub>U</sub></td><td>Mortality rate of spores</td><td>0.015/day</td></tr><tr><td>φ<sub>S</sub></td><td>Reduction rate of S</td><td>0.7/day</td></tr><tr><td>φ<sub>I</sub></td><td>Reduction rate of I</td><td>0.4/day</td></tr><tr><td>φ<sub>U</sub></td><td>Reduction rate of spores</td><td>0.001/day</td></tr></tbody></table>

The solutions of system (A.1) are positives. Indeed, suppose that one of the variables, that we will denote $x$, is equal to 0 at some instant, with all the others being non-negative. A quick analysis shows that $\dot{x} \ge 0$ cannot become negative. This implies that $x(t) \ge 0$. Then, the solutions are positive. The solutions of system (A.1) are piecewise continuous functions $X: \mathbb{R}_+ \to \mathbb{R}_+^4$, $X$ is continuous on $(nT, (n+1)T]$, $n \in \mathbb{N}$ and $X(nT^+) = \lim_{t \to nT^+} X(t)$. Clearly, the smooth properties of $f$ guarantee the existence and uniqueness of positive solutions of system (A.1). By adding the first and second equations in system (A.1), we see that the dynamics of total branches satisfies

$$ \begin{cases} \dot{N} = \Lambda - \mu N - dI, & t \neq nT \\ N(nT^+) \le \varphi_S N(nT). \end{cases} \quad (A.3) $$

From the equation (A.3), and using the positivity of the variables of the model, we have $\dot{N} \le \Lambda - \mu N$ and $N(nT^+) \le \varphi_S N(nT)$. It follows that $0 \le \lim_{t \to \infty} N(t) \le \frac{\Lambda}{\mu}$. Then, the following region

$$ G = \left\{ (N, U) \in \mathbb{R}_+^2 \mid N(t) \le \frac{\Lambda}{\mu}, U(t) \le \frac{\gamma\Lambda}{\mu(\nu + \mu_U)} \right\}. \quad (A.4) $$

is positively invariant by system (A.1).

### A.1.3 Mathematical analysis

Herein, we present the mathematical analysis of system (A.1). Since the state variable of berries $B$ is not present in the other equations of system (A.1), we do not consider it in the mathematical analysis.

#### A.1.3.1 Periodic disease free solution and its stability

The periodic disease free solution (PDFS) occurs when $I = 0$ and $U = 0$. Replacing these values in the impulsive differential equation (IDE) (A.1), we obtain

$$ \dot{S} = \Lambda - \mu S, \quad t \neq nT; \quad S(nT^+) = \varphi_S S(nT). \quad (A.5) $$

The resolution of the first equation of system (A.5) gives:
For $t \in (nT, (n+1)T]$

$$S(t) = \frac{\Lambda}{\mu} + \left( S(nT^{+}) - \frac{\Lambda}{\mu} \right) e^{-\mu(t-nT)}. \qquad (A.6)$$

This implies that

$$S((n+1)T) = \frac{\Lambda}{\mu} + \left( S(nT^{+}) - \frac{\Lambda}{\mu} \right) e^{-\mu T}. \qquad (A.7)$$

Using the impulsive condition $S((n+1)T^+) = \varphi_S S((n+1)T)$ yields

$$S((n+1)T^+) = \varphi_S \left( \frac{\Lambda}{\mu} + (S(nT^+) - \frac{\Lambda}{\mu})e^{-\mu T} \right). \qquad (A.8)$$

The fixed point of the equation (A.8) is given by

$$S(nT^+) = \frac{\Lambda\varphi_S(1 - e^{-\mu T})}{\mu(1 - \varphi_S e^{-\mu T})} > 0. \qquad (A.9)$$

Substituting the value of $S(nT)$ into the equation (A.6), for $t \in (nT, (n+1)T]$ one has

$$S^P(t) = \frac{\Lambda}{\mu} \left[ 1 - \frac{(1 - \varphi_S) e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right]. \qquad (A.10)$$

Finally, the PDFS is $X^P(t)=(S^P(t), 0, 0)$, where $S^P(t)$ is defined as in equation (A.10).

### A.1.3.2 Stability of the PDFS

We first study the local stability of $X^P(t)$ using small amplitude pertubation methods. Let us denote $\tilde{X}(t) = X(t) - X^P(t)$ where $X(t) = (S(t), I(t), U(t))^T$ and $\tilde{X}(t)$ is understood to be small amplitude pertubations. Substituting the expression of $\tilde{X}(t)$ in the system (A.1) gives

$$\left\{
\begin{aligned}
\dot{\tilde{S}} &= -\frac{\omega\nu\tilde{U}}{\tilde{N} + S^P(t)}(\tilde{S} + S^P(t)) - \mu\tilde{S}, \quad t \neq nT; \\
\dot{\tilde{I}} &= \frac{\omega\nu\tilde{U}}{\tilde{N} + S^P(t)}(\tilde{S} + S^P(t)) - (\mu + d)\tilde{I}, \quad t \neq nT; \\
\dot{\tilde{U}} &= \gamma\tilde{I} - (\nu + \mu_U)\tilde{U}, \quad t \neq nT; \\
\tilde{S}(nT^+) &= \varphi_S\tilde{S}(nT); \\
\tilde{I}(nT^+) &= \varphi_I\tilde{I}(nT); \\
\tilde{U}(nT^+) &= \varphi_U\tilde{U}(nT).
\end{aligned}
\right.
\qquad (A.11)$$

where $\tilde{N}(t) = \tilde{S}(t) + \tilde{I}(t)$ and $N^P(t) = S^P(t)$. The linearization of system (A.11) gives

$$\left\{
\begin{aligned}
\dot{\tilde{X}}(t) &= A\tilde{X}(t), \quad t \neq nT; \\
\tilde{X}(nT^+) &= \mathrm{diag}(\varphi_S, \varphi_I, \varphi_U)\tilde{X}(nT);
\end{aligned}
\right.
\qquad (A.12)$$

where

$$A = \begin{pmatrix}
-\mu & 0 & -\omega\nu \\
0 & -(\mu + d) & \omega\nu \\
0 & \gamma & -(\nu + \mu_U)
\end{pmatrix}.$$
Solving the system (A.12) yields $\tilde{X}(t) = \Phi(t)\tilde{X}(0)$, where $\Phi(0) = I$ with $I$ being the identity matrix of dimension 3 and $\Phi(t)$ the fundamental matrix that satisfies $\frac{d\Phi(t)}{dt} = A\Phi(t)$. A calculation gives

$$ \Phi(t) = \begin{pmatrix} e^{-\mu t} & \psi_{12}(t) & \psi_{13}(t) \\ 0 & \psi_{22}(t) & \psi_{23}(t) \\ 0 & \psi_{32}(t) & \psi_{33}(t) \end{pmatrix}, \qquad (A.13) $$

where

$$ \left\{ \begin{aligned} \psi_{22}(t) &= \frac{1}{2\beta} \left[ (\beta + k_1 - k_2)e^{\lambda_1 t} + (\beta - k_1 + k_2)e^{\lambda_2 t} \right], \\ \psi_{23}(t) &= \frac{\nu\omega}{\beta} \left[ -e^{\lambda_1 t} + e^{\lambda_2 t} \right], \psi_{32}(t) = \frac{\gamma}{\beta} \left[ -e^{\lambda_1 t} + e^{\lambda_2 t} \right], \\ \psi_{33}(t) &= \frac{1}{2\beta} \left[ (\beta - k_1 + k_2)e^{\lambda_1 t} + (\beta + k_1 - k_2)e^{\lambda_2 t} \right]. \end{aligned} \right. \qquad (A.14) $$

with

$$ \left\{ \begin{aligned} k_1 &= \mu + d, & k_2 &= \nu + \mu_U, \\ \alpha &= k_1 + k_2, & \beta &= \sqrt{(k_1 - k_2)^2 + 4\gamma\omega\nu}, \\ \lambda_1 &= -\frac{\alpha}{2} - \frac{\beta}{2} & \text{and} & \lambda_2 &= -\frac{\alpha}{2} + \frac{\beta}{2}. \end{aligned} \right. \qquad (A.15) $$

Using the impulsive conditions, one obtains the impulsive matrix $\tilde{X}(nT^+) = diag(\varphi_S, \varphi_S, \varphi_U)\tilde{X}(nT)$. Finally, the solution of system (A.12) becomes

$$ \tilde{X}(t) = diag(\varphi_S, \varphi_I, \varphi_U)\Phi(t)\tilde{X}(0). \qquad (A.16) $$

From the above equation, the monodromy matrix is $M = diag(\varphi_S, \varphi_I, \varphi_U)\Phi(T)$, i.e

$$ M = \begin{pmatrix} \varphi_S e^{-\mu T} & \varphi_S \xi_1(T) & \varphi_S \xi_2(T) \\ 0 & \varphi_I \psi_{22}(T) & \varphi_I \psi_{23}(T) \\ 0 & \varphi_U \psi_{32}(T) & \varphi_U \psi_{33}(T) \end{pmatrix}. \qquad (A.17) $$

The Floquet multipliers of the matrix $M$ is given by $\chi_1 = \varphi_S e^{-\mu T} < 1$ and the multipliers of the submatrix

$$ M_1 = \begin{pmatrix} \varphi_I \psi_{22}(T) & \varphi_I \psi_{23}(T) \\ \varphi_U \psi_{32}(T) & \varphi_U \psi_{33}(T) \end{pmatrix}. \qquad (A.18) $$

Then, the PDFS $X^P(t)$ is locally asymptotically stable if the eigenvalues of matrix $M_1$ stay in the unit cycle, which need to satisfy verify the following Jury conditions [173]:

$$ \left\{ \begin{array}{l} -\mathrm{tr}(M_1) - \det(M_1) < 1, \\ \det(M_1) < 1, \\ \mathrm{tr}(M_1) - \det(M_1) < 1, \end{array} \right. \qquad (A.19) $$
where

$$ \begin{cases} tr(M_1) = \varphi_I \psi_{22}(T) + \varphi_U \psi_{33}(T), \\ \det(M_1) = \varphi_I \varphi_U (\psi_{22}(T)\psi_{33}(T) - \psi_{23}(T)\psi_{32}(T)), \end{cases} \quad (A.20) $$

The conditions (A.19) hold if $\varphi_U \to 0$ and $\varphi_I \psi_{22}(T) < 1$. Finally, one has the following lemmata.

**Lemma A.1.** The PDFS $X^P(t)$ is locally asymptotically stable provided that $\mathcal{R} < 1$ where

$$ \mathcal{R} = \frac{\varphi_I}{2\beta} \left[ (\beta + k_1 - k_2)e^{\lambda_1 T} + (\beta - k_1 + k_2)e^{\lambda_2 T} \right]. \quad (A.21) $$

We have the following result about the global stability of the PDFS of the system (A.1).

**Theorem A.1.** The PDFS $X^P(t)$ of the system (A.1) is globally asymptotically stable in region $G$ provided that $\mathcal{R} < 1$.

*Proof of Theorem A.1.* We will prove the global attractivity of the PDFS $X^P(t)$. Consider the subsystem following:

$$ \left\{ \begin{aligned} \dot{I} &= \frac{\omega \nu U}{N} S - (\mu + d)I, \quad t \neq nT; \\ \dot{U} &= \gamma I - (\nu + \mu_U)U, \quad t \neq nT; \\ I(nT^+) &= \varphi_I I(nT), \quad U(nT^+) = \varphi_U U(nT). \end{aligned} \right. \quad (A.22) $$

The trivial solution of system (A.22) is $(0,0)$ and the solution $(S(t), I(t), U(t)) \in G$. The system (A.22) is cooperative because $(\frac{\partial I}{\partial U} = \frac{\omega \nu S}{N} > 0, \frac{\partial U}{\partial I} = \gamma > 0)$ and using the fact that $\frac{S(t)}{N(t)} < 1$, one has the following cooperative system:

$$ \left\{ \begin{aligned} \dot{I}_1 &= \omega \nu U_1 - (\mu + d) I_1, \quad t \neq nT; \\ \dot{U}_1 &= \gamma I_1 - (\nu + \mu_U) U_1, \quad t \neq nT; \\ I_1(nT^+) &= \varphi_I I_1(nT), \quad U_1(nT^+) = \varphi_U U_1(nT). \end{aligned} \right. \quad (A.23) $$

Applying the Kamke's theorem [33], one has that $I(t) \le I_1(t), U(t) \le U_1(t)$. The asymptotically, system (A.23) behaves like the system (A.22) [[64]], when $\mathcal{R}_0 < 1$. Then, one has

$$ I_1(t), U_1(t) \to 0^+ \Rightarrow I(t), U(t) \to 0^+. \quad (A.24) $$

Since $U \to 0^+$, for any arbitrary positive $\varepsilon_U$ there exists a $t_0 > 0$ such that for $t \ge t_0$, $U(t) \le \varepsilon_U$. Using the positivity of the solutions and the fact that $S(t)/N(t) \le 1$ into the first equation of system (A.1), one has:

$$ \Lambda - \omega \nu \varepsilon_U - \mu S(t) \le \dot{S}(t) \le \Lambda - \mu S(t) \quad (A.25) $$

Applying the comparison theorem on the above differential inequalities, one has $S_1(t) \le S(t) \le S_2(t)$, where for $t \in [nT; (n+1)T]$,

$$ \left\{ \begin{aligned} S_1(t) &= \frac{\Lambda - \omega \nu \varepsilon_U}{\mu} \left[ 1 - \frac{(1 - \varphi_S) e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right], \\ S_2(t) &= \frac{\Lambda}{\mu} \left[ 1 - \frac{(1 - \varphi_S) e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right], \end{aligned} \right. \quad (A.26) $$
<figure><img src="image_31.png" /><figcaption>Figure A.1: Trajectories of system (A.1) when $\gamma = 1.6$ (so that $\mathcal{R}_0 = 0.36$). All other parameter values are given in Table A.1.</figcaption></figure>

with $S_1$ and $S_2$ the solutions of the following impulsive differential equations:

$$ \dot{S}_1 = \Lambda - \omega\nu\varepsilon_U - \mu S_1, \quad t \neq nT; \quad S_1(nT^+) = \varphi_S S_1(nT), \qquad (A.27) $$

$$ \dot{S}_2 = \Lambda - \mu S_2, \quad t \neq nT; \quad S_2(nT^+) = \varphi_S S_2(nT), \qquad (A.28) $$

One can observe that $S_2(t) \to S^P(t)$ if $\varepsilon_U \to 0$ and $S_2(t) \to S^P(t)$, which implies that $S(t) \to S^P(t)$ asymptotically. In the other words, independently from the initial value $(S_0, I_0, U_0) \in G$, one has that

$$ (S(t), I(t), U(t)) \to (S^P(t), 0, 0) \qquad (A.29) $$

This concludes the proof. □

Numerical results for the global stability of the PDFS and persistence are present in Figure A.1 and A.2 respectively. To do so, we consider the initial conditions $(S_0, U_0) \in [(500, 200), (1000, 400), (1500, 600), (2000, 800)]$, $I_0 = 0$ and $B_0 = 0$. Figure A.1 shows that, for different initial conditions, the number of infected branches and urediniospores converge towards zero, which is the state free disease. Additionally, we can observe that the trajectories of susceptible branches and berries becomes stationary after three seasons, as proved in Theorem A.1. Although the persistence has not been proved analytically, numerical results in Figure A.2 tend to show that the persistence of the CLR appears for system (A.1) when $\mathcal{R}_0 = 5.7 > 1$. From Figures A.1 and A.2 one can deduce that, during the disease free period (when $\mathcal{R}_0 < 1$), the
<figure><img src="image_32.png" /><figcaption>Figure A.2: Trajectories of system (A.1) for different initial condition, when $\gamma = 2.1$ (so that $\mathcal{R}_0 = 5.7 > 1$). All other parameter values are given in Table A.1.</figcaption></figure>
number of berries is 378910, during the persistence of disease state (when $R_0 > 1$) the number of berries is 117698. Then, one can conclude that the yield lost due to the CLR is estimate at 68,93%. Thus, one can ask the following question, how to reduce the damage caused by the CLR in coffee plantation.

### A.1.4 Optimal control

#### A.1.4.1 Optimal problem

Many different methods are commonly used in CLR control, with varying levels of environmental impact. Herein, we choose the biological control, because it is a potentially powerful tool for managing coffee leaf rust that allows organic certification. This biological control uses natural enemies, e.g. mycoparasites such as *Lecanicillium lecanii* [26, 168]. The goal is to introduce the mycoparasite in the plantation so that it will consume the urediniospores. This method represents the inoculative biological control and is denoted by $v$. The new system with control $v(t)$ is the first system(A.1) where equation U changes to

$$ \dot{U} = \gamma I - (\nu + \mu_U)U - v, \quad t \neq nT; \qquad (A.30) $$

where $v(t) \in [0, v_{max}]$, with $v_{max}$ the maximum value of the biological control. Given $k$ seasons of duration $T$, we are interested in maximizing the net revenue, which is directly proportional to the amount of grains of coffee at the end of each season

$$ \sum_{j=1}^{k} B(jT), \qquad (A.31) $$

and inversely proportional to the usage of the control treatment throughout the whole period

$$ \int_{0}^{kT} v \, dt. \qquad (A.32) $$

Thus, we define the cost function to maximize representing the profit for a time period $[0, kT]$, given by

$$ J(v) = \alpha_B \sum_{j=1}^{k} B(jT) - \alpha_v \int_{0}^{kT} v \, dt, \qquad (A.33) $$

where $\alpha_B$ and $\alpha_v$ are the market prices per coffee grain and per spore, respectively. By using the relation

$$ B(jT) = \int_{(j-1)T}^{jT} (\delta_S S + \delta_I I - \mu_B B) \, dt, \qquad (A.34) $$

which takes into account that $B(jT^+) = 0$, we can rewrite the cost function as

$$ J(v) = \int_{0}^{kT} \left( \alpha_B (\delta_S S + \delta_I I - \mu_B B) - \alpha_v v \right) dt. \qquad (A.35) $$
Then, we define the Hybrid Optimal Control Problem

$$ \left\{ \begin{array}{ll} \text{maximize} & J(v) \\ \text{subject to} & \text{dynamics of system (A.30),} \\ & \text{initial conditions (A.2),} \\ & v(\cdot) \in \mathcal{V}, \end{array} \right. \qquad (\text{HOCP}) $$

with $\mathcal{V}$ being the set of admissible controllers, which are Lebesgue measurable real-valued functions defined on the interval $[0, kT]$ and satisfying the constraint $v(t) \in [0, v_{max}]$. We notice that $v(t)$ appears linearly in the dynamics (A.30), and, since $S = 0$ is repulsive, $S + I > 0$ for any trajectory, which means that the dynamics cannot blow up in finite time (and so all trajectories remain in a compact set). An analytical computation of an hybrid optimal control can be derived through the Hybrid Maximum Principle [39], a generalization of the classical PMP (Pontryagin's Maximum Principle). For (HOCP), we define the state $x = (S, I, U, B)$ and adjoint state $p = (p_S, p_I, p_U, p_B)$, and so the Hamiltonian is written as

$$ H(x, p, p_0, v) = \langle \lambda, f(x, v) \rangle - p_0 f^0(x, v) \qquad (\text{A.36}) $$

where $f^0(x, v) = \alpha_B \dot{B} - \alpha_v v$, which corresponds to the integrand of the cost function $J(v)$, $f$ is the right-hand side of (A.30), and $p_0 \le 0$. According to the Maximum Principle, the optimal control $v(t)$ should maximize the Hamiltonian along the optimal trajectory. Since $H$ is linear in the control, the solution of (HOCP) becomes

$$ v_{opt}(t) = \begin{cases} 0 & \text{if } \phi(x, p) < 0; \\ v_{max} & \text{if } \phi(x, p) > 0; \\ v_{sing}(t) & \text{if } \phi(x, p) = 0. \end{cases} \qquad (\text{A.37}) $$

where $\phi(x, p)$ is the switching function defined as

$$ \phi(x, p) \doteq \frac{\partial H}{\partial v}, \qquad (\text{A.38}) $$

and $v_{sing}(t)$ is called a singular control. According to (A.37), any solution of (HOCP) is a concatenation of bangs and singular arcs. A singular arc occurs when the switching function $\phi$ vanishes over some subinterval of time, and can be obtained by successively differentiating $\phi$ until the control $v$ can be computed explicitly in terms of the state and adjoint state. However, a more detailed description of the control law is often hard to obtain, due to the presence of the adjoint state $p$ in the singular arc $v_{sing}$ and in the switching function $\phi(x, p)$. One way to help understand the structure of the optimal control is to simulate numerical trajectories using direct methods, by approximating the problem by a finite dimensional optimization problem. Nevertheless, the use of optimal control solvers is usually restricted to continuous systems. Thus, in order to be able to apply such methods, we rewrite (HOCP) in a continuous form.
### A.1.4.2 Problem in continuous form

(HOCP) can be rewritten in the continuous form by defining $k$ set of states $x_n = (S_n, I_n, U_n, B_n)$ and $k$ controls $v_n(t)$, with dynamics defined for $t \in [0, T]$ as

$$ 
\left\{
\begin{aligned}
\dot{S}_n &= \Lambda - \omega\nu U_n \frac{S_n}{S_n + I_n} - \mu S_n, \\
\dot{I}_n &= \omega\nu U_n \frac{S_n}{S_n + I_n} - (\mu + d)I_n, \\
\dot{U}_n &= \gamma I_n - (\nu + \mu_U)U_n - v_n, \\
\dot{B}_n &= \delta_S S_n + \delta_I I_n - \mu_B B_n,
\end{aligned}
\right.
\qquad (A.39) $$

for $n = 1, \dots, k$, with boundary conditions

$$ 
\left\{
\begin{aligned}
S_1(0) &= S_0 > 0, \quad I_1(0) = I_0 \ge 0, \\
U_1(0) &= U_0 \ge 0, \quad B_1(0) = B_0 > 0 \\
S_{n+1}(0) &= \varphi_s S_n(T), \quad I_{n+1}(0) = \varphi_I I_n(T), \\
U_{n+1}(0) &= \varphi_U U_n(T), \quad B_{n+1}(0) = 0.
\end{aligned}
\right.
\qquad (A.40) $$

for $n = 1, \dots, k-1$. Then, the cost function becomes

$$ J_c(v_1, \dots, v_k) = \alpha_B \sum_{n=0}^{k} B_n(T) - \alpha_v \sum_{n=0}^{k} \int_{0}^{T} v_n dt, \qquad (A.41) $$

and the resulting continuous optimal control problem is

$$ 
\left\{
\begin{array}{ll}
\text{maximize} & J_c(v_1, \dots, v_k) \\
\text{subject to} & \text{dynamics of system (A.39),} \\
& \text{boundary conditions (A.40),} \\
& v_n(\cdot) \in \mathcal{V}, \text{ for } n = 1, \dots, k.
\end{array}
\right.
\qquad (\text{COCP}) $$

While this new OCP includes the boundary conditions relating the consecutive set of states $x_n$, there are no terminal conditions on (HOCP) which means that the set of admissible controllers is not empty. Thus, existence of a solution for (COCP)—and therefore for (HOCP)—holds by Filippov's theorem [2].

### A.1.4.3 Numerical results

Figure A.3 presents the results of a numerical simulation performed with Bocop of the optimal control and the optimal trajectory for problem (HOCP). In this Figure, we compare an open-loop strategy ($v \equiv 0$) with $\mathcal{R}_0 > 1$, and the same initial conditions with the optimal control. For the latter case, the biological treatment successfully eradicates the pest. The optimal strategy is characterized by a bang control $v \equiv v_{max}$ during an initial period of time in the first season, followed by a singular arc in which the use of the pest treatment is gradually reduced until no further control is applied. Indeed, the optimal control exploits the fact that both $U$ and $I$ are greatly reduced at the end of each season due to harvesting, to economize on pest treatment. After that, no further control is required. The latter pattern has also been observed for different set of parameters and initial conditions defined in the previous section. In general, the treatment
<figure><img src="image_33.png" /><figcaption>Figure A.3: Trajectories of HCOP (A.30) when $\gamma = 2.1$ (so that $R_0 = 5.7$ without control). The maximum control $v_{max} = 6.5$. All other parameter values are given in Table A.1. Initial conditions are set to $S(0) = 500$, $U(0) = 200$, $I(0) = B(0) = 0$.</figcaption></figure>
is applied rather for an initial interval of the season, followed by a decrease that takes the control
to $v \equiv 0$. As a conclusion for the analyzed particular case, coffee producers should maximize
the quantity of mycoparasitic (between given boundaries) roughly only during the beginning of
the first year, while discontinuing the treatment in posterior seasons. Naturally, the instant at
which the treatment should be turned off is strongly linked to the initial conditions, and so each
case requires a specific analysis.

A.1.5 Conclusion

This paper presented a comprehensive, impulsive model for the transmission dynamics of the CLR within a coffee plantation. The model has been analyzed to gain insight into its qualitative dynamics. Using the Floquet theory, we have found that there exists a threshold parameter that determines the extinction and persistence of CLR in a coffee plantation. Later on, we presented a controlled model based on an inoculative biological control, for which we defined an OCP targeting the maximization of the net revenue of a coffee plantation. For such problem, we performed a numerical analysis for a relevant set of parameters and initial conditions, intended to provide guidance in how to control the damages caused by CLR. In a follow-up of this work, we will perform a more detailed pest treatment analysis by proposing different control strategies including the dynamics of the mycoparasite in the model.
List of Publications and talks

**Paper 1** C. Djuikem, F. Grognard, R. T. Wafo, S. Bowong, S. Touzeau. Modelling coffee leaf rust dynamics to control its spread. *Mathematical Modelling of Natural Phenomena*, 26(16):25, 2021. https://doi.org/10.1051/mmnp/2021018.

* **Talk 1:** International Conference on Mathematical Methods and Models in Biosciences (BIOMATH 2019) | Jun 16 – Jun 22 | Będlewo | Poland | Conference grant awarded | Oral presentation: Modelling and controlling fungus *Hemileia vastatrix*, a coffee pest. https://hal.inria.fr/hal-02309454

**Paper 2** C. Djuikem, F. Grognard, S. Touzeau. Impulsive modelling of rust dynamics and hyperparasite releases for biocontrol. In revision for *Mathematical Biosciences*, September 2022.

* **Talk 2:** Conference on Mathematical Population Dynamics, Ecology and Evolution (MPDEE 2021) | Apr 26 – Apr 30 | CIRM | Marseille | France | Virtual oral presentation: Impulsive model and hyperparasite-based biocontrol of coffee leaf rust propagation. https://hal.inria.fr/hal-03456344.

* **Talk 3:** 7th IFAC Conference on Analysis and Design of Hybrid Systems (ADHS 2021) | July 7 – July 9 | Brussels | Belgium | Virtual oral presentation: Mathematical modelling and optimal control of the seasonal coffee leaf rust propagation. https://hal.inria.fr/hal-03274865.

* **Conference paper:** C. Djuikem, A. G. Yabo, F. Grognard, S. Touzeau. Mathematical modelling and optimal control of the seasonal coffee leaf rust propagation. *IFAC-PapersOnLine*, 54(5):193-198, 2021. https://doi.org/10.1016/j.ifacol.2021.08.497.

**Paper 3** C. Djuikem, F. Grognard, S. Touzeau. Bifurcation analysis in an epidemiological model of fungal crop disease. In revision for *Journal of Mathematical Biology*, August 2022.

* **Talk 4:** 12th European Conference on Mathematical and Theoretical Biology (ECMTB 2022) | Sept 19 – Sept 23 | Heidelberg | Germany | Oral presentation: Bifurcation analysis in a coffee leaf rust epidemiological model. https://ecmtb2022.org/wp-content/uploads/sites/18/2022/07/list_of_talks.pdf
Bibliography

[1] Z. Abbasi, I. Zamani, A. H. A. Mehra, M. Shafieirad, and A. Ibeas. Optimal control design of impulsive sqeiar epidemic models with application to covid-19. *Chaos, Solitons and Fractals*, 139:110054, 2020. doi:10.1016/j.chaos.2020.110054. (Cited on page 124.)

[2] A. A. Agrachev and Y. Sachkov. *Control theory from the geometric viewpoint*, volume 87. Springer Science & Business Media, 2013. (Cited on page 134.)

[3] G. N. Agrios. *Plant Pathology*. Academic Press, San Diego, USA, 4 edition, 1997. URL: https://books.google.fr/books?id=CnzbgZgby60C&lpg=PP1&ots=FrFqseYEpg&lr&hl=fr&pg=PP1#v=onepage&q&f=false. (Cited on pages 1, 75 and 91.)

[4] L. F. Aristizábal and M. A. Johnson. Monitoring coffee leaf rust (hemileia vastatrix) on commercial coffee farms in hawaii: Early insights from the first year of disease incursion. *Agronomy*, 12(5):1134, 2022. doi:10.3390/agronomy12051134. (Cited on page 7.)

[5] P. A. Arneson. Coffee rust, 2000. The Plant Health Instructor (updated 2011). URL: http://www.apsnet.org/edcenter/intropp/lessons/fungi/Basidiomycetes/Pages/CoffeeRust.aspx, doi:10.1094/PHI-I-2000-0718-02. (Cited on pages 5, 9, 61, 62 and 75.)

[6] J. Arroyo-Esquivel, F. Sanchez, and L. A. Barboza. Infection model for analyzing biological control of coffee rust using bacterial anti-fungal compounds. *Mathematical Biosciences*, 307:13-24, 2019. doi:10.1016/j.mbs.2018.10.009. (Cited on pages 17, 61 and 91.)

[7] J. Avelino, S. Cabut, B. Barboza, M. Barquero, R. Alfaro, C. Esquivel, J.-F. Durand, and C. Cilas. Topography and crop management are key factors for the development of american leaf spot epidemics on coffee in costa rica. *Phytopathology*, 97(12):1532-1542, 2007. doi:10.1094/PHYTO-97-12-1532. (Cited on page 4.)

[8] J. Avelino, M. Cristancho, S. Georgiou, P. Imbach, L. Aguilar, G. Bornemann, P. Läderach, F. Anzueto, A. J. Hruska, and C. Morales. The coffee rust crises in colombia and central america (2008-2013): impacts, plausible causes and proposed solutions. *Food Security*, 7(2):303-321, 2015. doi:10.1007/s12571-015-0446-9. (Cited on pages 7, 9 and 92.)

[9] J. Avelino, S. Gagliardi, I. Perfecto, M. Isaac, T. Liebig, J. Vandermeer, I. Merle, Z. Hajian-Forooshani, and N. Motisi. Tree effects on coffee leaf rust at field and landscape scales. *Plant Disease*, Jan. 2022. doi:10.1094/PDIS-08-21-1804-FE. (Cited on page 17.)

[10] J. Avelino, L. Willocquet, and S. Savary. Effects of crop management patterns on coffee rust epidemics. *Plant pathology*, 53(5):541-547, 2004. doi:10.1111/j.1365-3059.2004.01067.x. (Cited on pages 17, 61 and 118.)
[11] J. Avelino, H. Zelaya, A. Merlo, A. Pineda, M. Ordóñez, and S. Savary. The intensity of a coffee rust epidemic is dependent on production situations. *Ecological modelling*, 197(3-4):431-447, 2006. doi:10.1016/j.ecolmodel.2006.03.013. (Cited on pages 16, 33 and 61.)

[12] D. Bainov and P. Simeonov. *Impulsive differential equations: periodic solutions and applications*, volume 66. CRC Press, 1993. (Cited on pages 23, 24 and 25.)

[13] BBC News: Guatemala's coffee rust 'emergency' devastates crops, 2013. [accessed 21/10/2022]. URL: https://www.bbc.com/news/world-latin-america-21392257. (Cited on pages 3 and 18.)

[14] E. M. Beasley, N. Aristizábal, E. M. Bueno, and E. R. White. Spatially explicit models predict coffee rust spread in fragmented landscapes. *Landscape Ecology*, 37(8):2165-2178, 2022. doi:10.1007/s10980-022-01473-1. (Cited on page 17.)

[15] D. P. Bebber, Á. D. Castillo, and S. J. Gurr. Modelling coffee leaf rust risk in Colombia with climate reanalysis data. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 371(1709):20150458, 2016. doi:10.1098/rstb.2015.0458. (Cited on pages 17, 34, 61 and 116.)

[16] J. A. M. Bedimo, B. P. Dufour, C. Cilas, and J. Avelino. Effets des arbres d'ombrage sur les bioagresseurs de Coffea arabica. *Cahiers Agricultures*, 21(2-3):89-97, 2012. doi: 10.1684/agr.2012.0550. (Cited on page 62.)

[17] J. M. Bedimo, I. Njiayoum, D. Bieysse, M. N. Nkeng, C. Cilas, and J.-L. Notteghem. Effect of shade on arabica coffee berry disease development: Toward an agroforestry system to reduce disease impact. *Phytopathology*, 98(12):1320-1325, 2008. doi:10.1094/PHYTO-98-12-1320. (Cited on page 3.)

[18] S. Belbas and W. Schmidt. Optimal control of impulsive volterra equations with variable impulse times. *Applied mathematics and computation*, 214(2):353-369, 2009. doi:10.1016/j.amc.2009.04.066. (Cited on page 124.)

[19] K. R. Bock. Dispersal of uredospores of Hemileia vastatrix under field conditions. *Transactions of the British Mycological Society*, 45(1):63-74, 1962. doi:10.1016/S0007-1536(62)80035-7. (Cited on pages 36, 64, 100 and 126.)

[20] BOCOP: an open source toolbox for optimal control, 2017. URL: http://bocop.org. (Cited on page 124.)

[21] P. M. Boonekamp. Are plant diseases too much ignored in the climate change debate? *European Journal of Plant Pathology*, 133(1):291-294, 2012. doi:10.1007/s10658-011-9934-8. (Cited on page 33.)

[22] J. B. Burie, A. Calonnec, and A. Ducrot. Singular perturbation analysis of travelling waves for a model in phytopathology. *Mathematical Modelling of Natural Phenomena*, 1(1):49-62, 2006. doi:10.1051/mmnp:2006003. (Cited on page 12.)
[23] J.-B. Burie, A. Calonnec, and M. Langlais. Modeling of the invasion of a fungal disease over a vineyard. In *Mathematical Modeling of Biological Systems*, volume 2, pages 11-21. Springer, 2008. doi:10.1007/978-0-8176-4556-4_2. (Cited on pages 12, 14, 33, 61, 62 and 64.)

[24] J.-B. Burie and A. Ducrot. A field scale model for the spread of fungal diseases in crops: the example of a powdery mildew epidemic over a large vineyard. *Mathematical Methods in the Applied Sciences*, 38(17):3720-3737, 2015. doi:10.1002/mma.3312. (Cited on page 12.)

[25] A. Capucho, L. Zambolim, U. Lopes, and N. Milagres. Chemical control of coffee leaf rust in coffee canephora cv. conilon. *Australasian Plant Pathology*, 42(6):667-673, 2013. doi:10.1007/s13313-013-0242-y. (Cited on pages 9, 61 and 75.)

[26] G. Carrion and V. Rico-Gray. Mycoparasites on the coffee rust in Mexico. *Fungal Diversity*, 11(5):49-60, 10 2002. URL: https://www.fungaldiversity.org/fdp/sfdp/FD_11_49-60.pdf. (Cited on pages 75, 80 and 132.)

[27] C. Castillo-Chavez and B. Song. Dynamical models of tuberculosis and their applications. *Mathematical Biosciences & Engineering*, 1(2):361, 2004. doi:10.3934/mbe.2004.1.361. (Cited on pages 23 and 83.)

[28] G. Champéroux. *Manuel du planteur de café laotien*. CIRAD-IRCC, IRCC, Montpellier, 1991. URL: http://agritrop.cirad.fr/345484/. (Cited on pages 4, 5, 62 and 100.)

[29] A. Charrier and J. Berthaud. Botanical classification of coffee. In *Coffee*, pages 13-47. Springer, 1985. doi:10.1007/978-1-4615-6657-1_2. (Cited on pages 36 and 64.)

[30] A. Charrier and A. B. Eskes. *Botany and Genetics of Coffee*, chapter 2, pages 25-56. John Wiley & Sons, Ltd, 2004. doi:/10.1002/9783527619627.ch2. (Cited on page 4.)

[31] M. N. Clifford and K. C. Willsons, editors. *Coffee: botany, biochemistry and production of beans and beverage*. Springer, New York, 1985. doi:10.1007/978-1-4615-6657-1. (Cited on pages 5, 11, 99 and 100.)

[32] J. S. Coleman. Leaf development and leaf stress: increased susceptibility associated with sink-source transition. *Tree Physiology*, 2(1-2-3):289-299, 12 1986. arXiv:https://academic.oup.com/treephys/article-pdf2/1-2-3/289/4594414/2-1-2-3-289.pdf, doi:10.1093/treephys/2.1-2-3.289. (Cited on page 33.)

[33] W. A. Coppel. *Stability and asymptotic behavior of differential equations*. Heath, 1965. (Cited on pages 98 and 129.)

[34] S. Daivasikamani et al. Biological control of coffee leaf rust pathogen, hemileia vastatrix berkeley and broome using bacillus subtilis and pseudomonas fluorescens. *Journal of Biopesticides*, 2(1):94-98, 2009. URL: http://www.jbiopest.com. (Cited on pages 10 and 92.)

[35] E. De Jong, A. Eskes, J. Hoogstraten, and J. Zadoks. Temperature requirements for germination, germ tube growth and appressorium formation of urediospores of hemileia vastatrix. *Netherlands Journal of Plant Pathology*, 93(2):61-71, 1987. doi:10.1007/BF01998091. (Cited on page 118.)
[36] M.-L. Desprez-Loustau, F. M. Hamelin, and B. Marçais. The ecological and evolutionary trajectory of oak powdery mildew in Europe. In K. Wilson, A. Fenton, and D. Tompkins, editors, *Wildlife Disease Ecology: Linking Theory to Data and Application*, Ecological Reviews, chapter 15, page 429-457. Cambridge University Press, 2019. doi:10.1017/9781316479964.015. (Cited on pages 12, 33 and 91.)

[37] C. Djuikem, A. Gabriel Yabo, F. Grognard, and S. Touzeau. Mathematical modelling and optimal control of the seasonal coffee leaf rust propagation. *IFAC-PapersOnLine*, 54(5):193-198, 2021. 7th IFAC Conference on Analysis and Design of Hybrid Systems ADHS 2021. doi:10.1016/j.ifacol.2021.08.497. (Cited on pages 89, 92, 119 and 123.)

[38] C. Djuikem, F. Grognard, R. Tagne Wafo, S. Touzeau, and S. Bowong. Modelling coffee leaf rust dynamics to control its spread. *Mathematical Modelling of Natural Phenomena*, 16:26, 2021. doi:10.1051/mmnp/2021018. (Cited on pages 34, 59, 116 and 119.)

[39] A. V. Dmitruk and A. M. Kaganovich. The hybrid maximum principle is a consequence of pontryagin maximum principle. *Systems & Control Letters*, 57(11):964-970, 2008. doi:10.1016/j.sysconle.2008.05.006. (Cited on pages 124 and 133.)

[40] A. d'Onofrio. Stability properties of pulse vaccination strategy in seir epidemic model. *Mathematical biosciences*, 179(1):57-72, 2002. doi:10.1016/S0025-5564(02)00095-0. (Cited on page 124.)

[41] e-Krishi Shiksha. Diseases of coffee. [accessed 15/9/2022]. URL: http://ecoursesonline.iasri.res.in/mod/page/view.php?id=9675. (Cited on pages 2, 3 and 4.)

[42] M. El Jarroudi, H. Karjoun, L. Kouadio, and M. El Jarroudi. Mathematical modelling of non-local spore dispersion of wind-borne pathogens causing fungal diseases. *Applied Mathematics and Computation*, 376:125107, 2020. doi:10.1016/j.amc.2020.125107. (Cited on page 12.)

[43] A. Eskes. *Incomplete resistance to coffee leaf rust (Hemileia vastatrix)*. Wageningen University and Research, 1983. (Cited on page 11.)

[44] A. Eskes, M. Mendes, and C. Robbs. Laboratory and field studies on parasitism of hemileia vastatrix with verticillium lecanii and v. leptobactrum. *Café, Cacao, Thé (Francia) v. 35 (4) p. 275-282*, 1991. (Cited on page 10.)

[45] F. Fabre, J.-B. Burie, A. Ducrot, S. Lion, Q. Richard, and R. Djidjou-Demasse. An evolutionary model for predicting the adaptation of spore-producing pathogens to quantitative resistance in heterogeneous environments. *Evolutionary applications*, 15(1):95-110, 2022. doi:10.1111/eva.13328. (Cited on page 121.)

[46] Climate change fans spread of pests and threatens plants and crops, new fao study. [accessed 15/9/2022]. URL: https://www.fao.org/news/story/en/item/1402920/icode/. (Cited on page 1.)

[47] FAOSTAT: Cameroon, agricultural land, 2016. [Online; accessed August 3, 2022]. URL: https://www.fao.org/faostat/en/#country/32. (Cited on page 2.)
[48] M. Fazeli-Dinan, R. Talaei-Hassanloui, and M. Goettel. Virulence of the entomopathogenic fungus lecanicillium longisporum against the greenhouse whitefly, trialeurodes vaporarium and its parasitoid encarsia formosa. *International Journal of Pest Management*, 62(3):251-260, 2016. doi:10.1080/09670874.2016.1182228. (Cited on page 120.)

[49] R. A. Fleming. The potential for control of cereal rust by natural enemies. *Theoretical Population Biology*, 18(3):374-395, 1980. doi:10.1016/0040-5809(80)90060-X. (Cited on pages 12, 13, 33 and 91.)

[50] J. Flood. *Coffee wilt disease*. CABI, 2010. (Cited on page 3.)

[51] Fortune Business Insights. Gree coffee market size, share & covid-19 impact analysis. [accessed 15/9/2022]. URL: https://www.fortunebusinessinsights.com/green-coffee-market-106635. (Cited on page 2.)

[52] Y. Fotso Fotso. *Modelling, analysis and control of coffee berry borer*. Theses, Université de Dschang (Cameroun), Jan. 2022. URL: https://hal.archives-ouvertes.fr/tel-03608463. (Cited on page 2.)

[53] D. Fróna, J. Szenderák, and M. Harangi-Rákos. The challenge of feeding the world. *Sustainability*, 11(20):5816, 2019. doi:10.3390/su11205816. (Cited on page 1.)

[54] L. Galbusera, M. P. E. Marciandi, P. Bolzern, and G. Ferrari-Trecate. Control schemes based on the wave equation for consensus in multi-agent systems with double-integrator dynamics. In *2007 46th IEEE Conference on Decision and Control*, pages 1498-1503. IEEE, 2007. doi:10.1109/CDC.2007.4434358. (Cited on page 72.)

[55] W. Gams and A. Van Zaayen. Contribution to the taxonomy and pathogenicity of fungi-colous verticillium species. i. taxonomy. *Netherlands Journal of Plant Pathology*, 88(2):57-78, 1982. doi:10.1007/BF01977339. (Cited on page 92.)

[56] B. Ghosh, F. Grognard, and L. Mailleret. Natural enemies deployment in patchy environments for augmentative biological control. *Applied Mathematics and Computation*, 266:982-999, 2015. doi:10.1016/j.amc.2015.06.021. (Cited on pages 26 and 104.)

[57] M. S. Goettel, M. Koike, J. J. Kim, D. Aiuchi, R. Shinya, and J. Brodeur. Potential of Lecanicillium spp. for management of insects, nematodes and plant diseases. *Journal of Invertebrate Pathology*, 98(3):256-261, 2008. Special Issue for SIP 2008. doi:10.1016/j.jip.2008.01.009. (Cited on page 117.)

[58] G. Grée. Epidemiology of coffee leaf rust in the Eastern Highlands. *Coffee Research Institute Newsletter*, 2:16-20, 1993. URL: https://agritrop.cirad.fr/399280/. (Cited on pages 61 and 92.)

[59] F.-M. Gumpert, H. Geiger, and U. Stähle. A mathematical model of the epidemics in homogeneous and heterogeneous host stands/ein mathematisches modell für epidemien in homogenen und heterogenen wirtsbeständen. *Zeitschrift für Pflanzenkrankheiten und Pflanzenschutz/Journal of Plant Diseases and Protection*, pages 206-215, 1987. URL: https://www.jstor.org/stable/43383228. (Cited on page 12.)
[60] S. Hacquard, B. Petre, P. Frey, A. Hecker, N. Rouhier, and S. Duplessis. The poplar-poplar rust interaction: insights from genomics and transcriptomics. *Journal of pathogens*, 2011, 2011. doi:10.4061/2011/716041. (Cited on page 33.)

[61] F. Haddad, L. A. Maffia, E. S. G. Mizubuti, and H. Teixeira. Biological control of coffee rust by antagonistic bacteria under field conditions in Brazil. *Biological Control*, 49(2):114-119, 2009. doi:10.1016/j.biocontrol.2009.02.004. (Cited on page 75.)

[62] F. Haddad, R. M. Saraiva, E. S. G. Mizubuti, R. S. Romeiro, and L. A. Maffia. Antifungal compounds as a mechanism to control Hemileia vastatrix by antagonistic bacteria. *Tropical Plant Pathology*, 38(5):398-405, 2013. doi:10.1590/S1982-56762013000500004. (Cited on pages 10 and 75.)

[63] Z. Hajian-Forooshani, I. S. Rivera Salinas, E. Jiménez-Soto, I. Perfecto, and J. Vandermeer. Impact of regionally distinct agroecosystem communities on the potential for autonomous control of the coffee leaf rust. *Environmental Entomology*, 45(6):1521-1526, 10 2016. arXiv:https://academic.oup.com/ee/article-pdf/45/6/1521/8660201/nvw125.pdf, doi:10.1093/ee/nvw125. (Cited on pages 10, 92 and 120.)

[64] J. Hale. *Ordinary differential equations, 1969*. Wiley-Interscience, New York, 1969. (Cited on page 129.)

[65] W. H. Hamer. *Epidemic disease in England: the evidence of variability and of persistency of type*. Bedford Press, 1906. (Cited on page 11.)

[66] Hawaii Department of Agriculture Divisions. Board of agriculture expands coffee quarantine to o'ahu and lana'i. [accessed 20/10/2022]. URL: https://hdoa.hawaii.gov/blog/main/nr21-08coffeequarantineexpanded/. (Cited on pages 3, 4 and 7.)

[67] D. Henk, D. Farr, and M. Aime. Mycodiplosis (diptera) infestation of rust fungi is frequent, wide spread and possibly host specific. *Fungal Ecology*, 4(4):284-289, 2011. doi:10.1016/j.funeco.2011.03.006. (Cited on pages 10, 18, 92, 117 and 120.)

[68] A. Heydari, M. Pessarakli, et al. A review on biological control of fungal plant pathogens using microbial antagonists. *Journal of biological sciences*, 10(4):273-290, 2010. doi: 10.3923/jbs.2010.273.290. (Cited on page 92.)

[69] F. D. Hinnah, P. C. Sentelhas, F. R. Alves Patrício, R. N. Paiva, and M. V. Parenti. Performance of a weather-based forecast system for chemical control of coffee leaf rust. *Crop Protection*, 137:105225, 2020. doi:10.1016/j.cropro.2020.105225. (Cited on page 17.)

[70] P. Hrelia, C. Fimognari, F. Maffei, F. Vigagni, R. Mesirca, L. Pozzetti, M. Paolini, and G. C. Forti. The genetic and non-genetic toxicity of the fungicide vinclozolin. *Mutagenesis*, 11(5):445-453, 1996. doi:10.1093/mutage/11.5.445. (Cited on page 91.)

[71] International Coffee Organization. Coffee production by exporting countries. [accessed 23/2/202]. URL: http://www.ico.org/prices/po-production.pdf. (Cited on pages 61 and 64.)

[72] International Coffee Organization. World coffee consumption, 2020. [accessed 12/11/2020]. URL: http://www.ico.org/prices/new-consumption-table.pdf. (Cited on page 124.)
[73] D. Jackson, J. Skillman, and J. Vandermeer. Indirect biological control of the coffee leaf rust, hemileia vastatrix, by the entomogenous fungus lecanicillium lecanii in a complex coffee agroecosystem. *Biological Control*, 61(1):89-97, 2012. doi:10.1016/j.biocontrol.2012.01.004. (Cited on page 10.)

[74] J. Jaramillo, E. Muchugu, F. E. Vega, A. Davis, C. Borgemeister, and A. Chabi-Olaye. Some like it hot: the influence and implications of climate change on coffee berry borer (hypothenemus hampei) and coffee production in east Africa. *PloS one*, 6(9):e24528, 2011. doi:10.1371/journal.pone.0024528. (Cited on page 2.)

[75] M. J. Jeger, P. Jeffries, Y. Elad, and X.-M. Xu. A generic theoretical model for biological control of foliar plant diseases. *Journal of Theoretical Biology*, 256(2):201-214, 2009. doi: 10.1016/j.jtbi.2008.09.036. (Cited on page 91.)

[76] M. J. Jeger, P. J. Wijngaarden, and R. F. Hoekstra. Adaptation to the cost of resistance in a haploid clonally reproducing organism: The role of mutation, migration and selection. *Journal of Theoretical Biology*, 252(4):621-632, 2008. doi:10.1016/j.jtbi.2008.02.023. (Cited on page 91.)

[77] John Megahan. Life cycle of the coffee leaf rust. [accessed 8/9/2022]. URL: https://feralatlas.supdigital.org/poster/coffee-rust-spreads-together-with-coffee-plantations. (Cited on page 8.)

[78] J.-C. Kamgang. *Contribution à la stabilisation des systèmes mécaniques: contribution à l'étude de la stabilité des modèles épidémiologiques*. PhD thesis, Université Paul Verlaine-Metz, 2003. URL: https://hal.univ-lorraine.fr/tel-01749936. (Cited on page 40.)

[79] J. C. Kamgang and G. Sallet. Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (dfe). *Mathematical biosciences*, 213(1):1-12, 2008. doi:10.1016/j.mbs.2008.02.005. (Cited on page 42.)

[80] W. O. Kermack and A. G. McKendrick. A contribution to the mathematical theory of epidemics. *Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character*, 115(772):700-721, 1927. doi:10.1098/rspa.1927.0118. (Cited on page 11.)

[81] H. Kielhöfer. Stability and semilinear evolution equations in Hilbert space. *Archive for Rational Mechanics and Analysis*, 57(2):150-165, 1974. doi:10.1007/BF00248417. (Cited on pages 80 and 81.)

[82] M. Komárek, E. Čadková, V. Chrastný, F. Bordas, and J.-C. Bollinger. Contamination of vineyard soils with fungicides: a review of environmental and toxicological aspects. *Environment international*, 36(1):138-151, 2010. doi:10.1016/j.envint.2009.10.005. (Cited on page 91.)

[83] A. Koutouleas, H. J. L. Jørgensen, B. Jensen, J.-P. B. Lillesø, A. Junge, and A. Ræbild. On the hunt for the alternate host of Hemileia vastatrix. *Ecology and Evolution*, 9(23):13619-13631, 2019. doi:10.1002/ece3.5755. (Cited on page 33.)
[84] A. Kushalappa, M. Akutsu, S. Oseguera, G. Chaves, C. Melles, J. Miranda, and G. Bartolo. Equations for predicting the rate of coffee rust development based on net survival ratio for monocyclic process of hemileia vastatrix *Coffea arabica*. *Fitopatologia Brasileira* (Brazil), 1984. URL: https://www.apsnet.org/publications/phytopathology/backissues/Documents/1983Articles/Phyto73n01_96.PDF. (Cited on page 16.)

[85] A. Kushalappa, G. Chaves, et al. An analysis of the development of coffee rust in the field. *Fitopatologia Brasileira*, 5(1):95-103, 1980. (Cited on page 16.)

[86] A. C. Kushalappa. Linear models applied to variation in the rate of coffee rust development. *Journal of Phytopathology*, 101(1):22-30, 1981. doi:https://doi.org/10.1111/j.1439-0434.1981.tb03317.x. (Cited on page 16.)

[87] A. C. Kushalappa and B. E. Albertus. *Coffee rust: Epidemiology, resistance and management*. CRC Press, Florida, USA, 1989. (Cited on pages 3 and 9.)

[88] J. P. LaSalle. *The stability of dynamical systems*, volume 25. Siam, 1976. (Cited on page 71.)

[89] K. Li, J. H. Vandermeer, and I. Perfecto. Disentangling endogenous versus exogenous pattern formation in spatial ecology: a case study of the ant azteca sericeasur in southern Mexico. *Royal Society open science*, 3(5):160073, 2016. doi:10.1098/rsos.160073. (Cited on page 120.)

[90] T. Liebig, F. Ribeyre, P. Läderach, H.-M. Poehling, P. van Asten, and J. Avelino. Interactive effects of altitude, microclimate and shading system on coffee leaf rust. *Journal of Plant Interactions*, 14(1):407-415, 2019. doi:10.1080/17429145.2019.1643934. (Cited on page 10.)

[91] T.-K. Lim and W. Z. W. Nik. Mycoparasitism of the coffee rust pathogen, hemileia vastatrix, by verticillium psalliotae in Malaysia. *Pertanika.*, 6(2):23-25, 1983. doi:http://agris.upm.edu.my:8080/dspace/handle/0/13206. (Cited on page 10.)

[92] U. Lopes, L. Zambolim, P. Neto, H. Duarte, J. Ribeiro, A. Souza, and F. Rodrigues. Silicon and Triadimenol for the management of coffee leaf rust. *Journal of Phytopathology*, 162, 02 2014. doi:doi:10.1111/jph.12166. (Cited on page 64.)

[93] J. A. Lucas, N. J. Hawkins, and B. A. Fraaije. The evolution of fungicide resistance. *Advances in applied microbiology*, 90:29-92, 2015. doi:10.1016/bs.aambs.2014.09.001. (Cited on page 9.)

[94] L. V. Madden and F. Van Den Bosch. A population-dynamics approach to assess the threat of plant pathogens as biological weapons against annual crops: Using a coupled differential-equation model, we show the conditions necessary for long-term persistence of a plant disease after a pathogenic microorganism is introduced into a susceptible annual crop. *BioScience*, 52(1):65-74, 2002. doi:10.1641/0006-3568(2002)052[0065:APDATA]2.0.CO;2. (Cited on page 15.)
[95] M. Mahfud, M. A. ZA, S. Meon, and J. Kadir. In vitro and in vivo tests for parasitism of verticillium psalliotae treschow on hemileia vastatrix berk. and br. *Malaysian Journal of Microbiology*, 2(1):46-50, 2006. doi:10.21161/MJM.210608. (Cited on page 10.)

[96] L. Mailleret and F. Grognard. Global stability and optimisation of a general impulsive biological control model. *Mathematical Biosciences*, 221(2):91-100, 2009. doi:10.1016/j.mbs.2009.07.002. (Cited on pages 80 and 124.)

[97] Y. Mammeri, J. Burie, M. Langlais, and A. Calonnec. How changes in the dynamic of crop susceptibility and cultural practices can be used to better control the spread of a fungal pathogen at the plot scale? *Ecological Modelling*, 290:178-191, 2014. Special Issue of the 4th International Symposium on Plant Growth Modeling, Simulation, Visualization and Applications (PMA'12 ). doi:10.1016/j.ecolmodel.2014.02.017. (Cited on pages 12, 33 and 91.)

[98] L. Manuel, P. Talhinhas, V. Varzea, and J. Neves-Martins. Characterization of col-letotrichum kahawae isolates causing coffee berry disease in angola. *Journal of phytopathol-ogy*, 158(4):310-313, 2010. doi:10.1111/j.1439-0434.2009.01613.x. (Cited on page 3.)

[99] D. Masaba and J. M. Waller. Coffee berry disease: the current status. *Coffee berry disease: the current status.*, pages 237-249, 1992. (Cited on page 3.)

[100] A. Maublanc and L. Roger. Une nouvelle rouille du caféier au Cameroun. *Bulletin de la Société Mycologique de France*, 50:193-202, 1934. (Cited on page 4.)

[101] A. Maupetit, R. Larbat, M. Pernaci, A. Andrieux, C. Guinet, A.-L. Boutigny, B. Fabre, P. Frey, and F. Halkett. Defense compounds rather than nutrient availability shape ag-gressiveness trait variation along a leaf maturity gradient in a biotrophic plant pathogen. *Frontiers in Plant Science*, 9, 2018. doi:10.3389/fpls.2018.01396. (Cited on pages 12 and 33.)

[102] S. McCook and J. Vandermeer. The big rust and the red queen: Long-term perspec-tives on coffee rust research. *Phytopathology*, 105(9):1164-1173, 2015. doi:10.1094/PHYTO-04-15-0085-RVW. (Cited on page 9.)

[103] J. McDonald. A preliminary account of a disease of green coffee berries in kenya colony. *Transactions of the British Mycological Society*, 11(1):145-154, 1926. doi:10.1016/S0007-1536(26)80033-6. (Cited on page 3.)

[104] K. Mendgen. Growth of verticillium lecanii in pustules of stripe rust (puccinia striiformis). *Phytopathologische Zeitschrift*, 102(3-4):301-309, 1981. doi:10.1111/j.1439-0434.1981.tb03391.x. (Cited on page 92.)

[105] X. Meng and Z. Li. The dynamics of plant disease models with continuous and impulsive cultural control strategies. *Journal of Theoretical Biology*, 266(1):29-40, 2010. doi:10.1016/j.jtbi.2010.05.033. (Cited on pages 13 and 91.)

[106] I. Merle. *Effets du microclimat sur le développement de l'épidémie de rouille orangée du caféier Arabica (Hemileia vastatrix-Coffea arabica) dans une gamme de situations de production*. PhD thesis, Université de Montpellier, 2019. (Cited on pages 18 and 34.)
[107] I. Merle, P. Tixier, E. de Melo Virginio Filho, C. Cilas, and J. Avelino. Forecast models of coffee leaf rust symptoms and signs based on identified microclimatic combinations in coffee-based agroforestry systems in costa rica. *Crop Protection*, 130:105046, 2020. doi:10.1016/j.cropro.2019.105046. (Cited on page 17.)

[108] R. E. Mickens. Nonstandard finite difference schemes for reaction-diffusion equations. *Numerical Methods for Partial Differential Equations: An International Journal*, 15(2):201-214, 1999. doi:10.1080/1023619021000000807. (Cited on pages 28 and 73.)

[109] L. C. Monaco et al. Consequences of the introduction of coffee rust into brazil. *The genetic basis of epidemics in agriculture*, pages 57-71, 1977. (Cited on page 9.)

[110] L. B. Mondiale. Agriculture et alimentation. [accessed 14/9/2022]. URL: https://www.banquemondiale.org/fr/topic/agriculture/overview. (Cited on page 1.)

[111] D. Moore, G. D. Robson, and A. P. J. Trinci. *21st Century Guidebook to Fungi*. Cambridge University Press, 2000. doi:10.1017/CBO9780511977022. (Cited on page 11.)

[112] R. A. Muller, D. Berry, J. Avelino, and D. Bieysse. *Coffee Diseases*, chapter 4, pages 491-545. John Wiley & Sons, Ltd, 2008. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/9783527619627.ch18, doi:10.1002/9783527619627.ch18. (Cited on page 61.)

[113] C. Nembot, P. Takam Soh, G. M. Ten Hoopen, and Y. Dumont. Modeling the temporal evolution of cocoa black pod rot disease caused by phytophthora megakarya. *Mathematical Methods in the Applied Sciences*, 41(18):8816-8843, 2018. doi:10.1002/mma.5206. (Cited on pages 13 and 91.)

[114] S. Nilusmas, M. Mercat, T. Perrot, C. Djian-Caporalino, P. Castagnone-Sereno, S. Touzeau, V. Calcagno, and L. Mailleret. Multi-seasonal modelling of plant-nematode interactions reveals efficient plant resistance deployment strategies. *Evolutionary applications*, 13(9):2206-2221, 2020. doi:10.1111/eva.12989. (Cited on page 121.)

[115] S. Nundloll, L. Mailleret, and F. Grognard. Two models of interfering predators in impulsive biological control. *Journal of Biological Dynamics*, 4(1):102-114, 2010. doi:10.1080/17513750902968779. (Cited on pages 13 and 91.)

[116] F. J. Nutman, F. M. Roberts, and R. T. Clarke. Studies on the biology of hemileia vastatrix berk. & br. *Transactions of the British Mycological Society*, 46(1):27-44, 1963. doi:10.1016/S0007-1536(63)80005-4. (Cited on pages 33 and 62.)

[117] E.-C. Oerke. Crop losses to pests. *The Journal of Agricultural Science*, 144(1):31-43, 2006. (Cited on page 1.)

[118] T. O. of Economic Complexity (OEC). Historical data, 2017. [Online; accessed August 3, 2022]. URL: https://oec.world/en/profile/country/cmr#yearly-trade. (Cited on page 2.)

[119] C.-V. Pao. *Nonlinear parabolic and elliptic equations*. Springer Science & Business Media, 2012. (Cited on page 26.)
[120] J. Papaïx, K. Adamczyk-Chauvat, A. Bouvier, K. Kiêu, S. Touzeau, C. Lannou, and H. Monod. Pathogen population dynamics in agricultural landscapes: The Ddal modelling framework. *Infection, Genetics and Evolution*, 27:509–520, 2014. (Cited on pages 12 and 61.)

[121] A. Pazy. *Semigroups of linear operators and applications to partial differential equations*, volume 44. Springer Science & Business Media, 2012. (Cited on pages 26, 27, 68, 69 and 70.)

[122] F. Pellegrin, B. Seivert, F. Kohler, C. Van Bercie, and B. Boccas. La rouille orangée du caféier Arabica en Nouvelle Calédonie : historique et épidémiologie. *Café, Cacao, Thé*, 27(1):27-40, 1983. URL: https://www/documentation.ird.fr/hor/fdi:03582. (Cited on pages 18 and 34.)

[123] D. R. Pereira, D. H. Nadaleti, E. C. Rodrigues, A. D. da Silva, M. R. Malta, S. P. de Carvalho, and G. R. Carvalho. Genetic and chemical control of coffee rust (hemileia vastatrix b erk et b r.): impacts on coffee (coffea arabica l.) quality. *Journal of the Science of Food and Agriculture*, 101(7):2836-2845, 2021. doi:10.1002/jsfa.10914. (Cited on page 9.)

[124] S. Pivonia and X. Yang. Relating epidemic progress from a general disease model to seasonal appearance time of rusts in the united states: Implications for soybean rust. *Phytopathology*, 96(4):400-407, 2006. doi:10.1094/PHYTO-96-0400. (Cited on pages 12, 33 and 91.)

[125] L. S. Pontryagin. *Mathematical theory of optimal processes*. Routledge, 2018. (Cited on page 124.)

[126] N. S. Prakash, J. Devasia, Jayarama, and R. K. Aggarwal. Chapter 8 - coffee industry in india: Production to consumption—a sustainable enterprise. In V. R. Preedy, editor, *Coffee in Health and Disease Prevention*, pages 61–70. Academic Press, San Diego, 2015. doi:10.1016/B978-0-12-409517-5.00008-5. (Cited on page 2.)

[127] S. Rahman, A. R. Anik, and J. R. Sarker. Climate, environment and socio-economic drivers of global agricultural productivity growth. *Land*, 11(4):512, 2022. doi:10.3390/land11040512. (Cited on page 1.)

[128] D. Ramiro, A. Jalloul, A.-S. Petitot, M. F. Grossi De Sá, M. P. Maluf, and D. Fernandez. Identification of coffee wrky transcription factor genes and expression profiling in resistance responses to pathogens. *Tree Genetics & Genomes*, 6(5):767-781, 2010. doi:10.1007/s11295-010-0290-1. (Cited on page 121.)

[129] D. A. Ramiro, O. Guerreiro-Filho, and P. Mazzafera. Phenol contents, oxidase activities, and the resistance of coffee to the leaf miner leucoptera coffeella. *Journal of chemical ecology*, 32(9):1977-1988, 2006. doi:10.1007/s10886-006-9122-z. (Cited on page 121.)

[130] V. Ravigné, V. Lemesle, A. Walter, L. Mailleret, and F. M. Hamelin. Mate limitation in fungal plant parasites can lead to cyclic epidemics in perennial host populations. *Bulletin of mathematical biology*, 79(3):430-447, 2017. doi:10.1007/s11538-016-0240-7. (Cited on pages 12, 33 and 91.)