<figure><img src="image_1.png" /></figure>

HAL
open science

# Impulsive modelling of rust dynamics and predator releases for biocontrol

Clotilde Djuikem, Frédéric Grognard, Suzanne Touzeau

► To cite this version:

Clotilde Djuikem, Frédéric Grognard, Suzanne Touzeau. Impulsive modelling of rust dynamics and predator releases for biocontrol. Mathematical Biosciences, 2023, 356, pp.108968. 10.1016/j.mbs.2023.108968. hal-03952381

HAL Id: hal-03952381

https://inria.hal.science/hal-03952381

Submitted on 23 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Graphical Abstract

Impulsive modelling of rust dynamics and predator releases for
biocontrol

Clotilde Djuikem, Frédéric Grognard, Suzanne Touzeau

<figure><img src="image_3.png" /></figure>
Highlights

Impulsive modelling of rust dynamics and predator releases for
biocontrol

Clotilde Djuikem, Frédéric Grognard, Suzanne Touzeau

* Original impulsive model of multi-seasonal fungal disease dynamics

• Explicit dynamics of predators, acting as predators of spores

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

* predator release frequency modulates the control efficiency
Impulsive modelling of rust dynamics and predator
releases for biocontrol

Clotilde Djuikem<sup>a,*</sup>, Frédéric Grognard<sup>a</sup>, Suzanne Touzeau<sup>a,b</sup>

<sup>a</sup>Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, BIOCORE, France

<sup>b</sup>Université Côte d'Azur, INRAE, CNRS, ISA, France

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 dy-
namics during the dry season. Biological control using predators through
one or more discrete introduction events over the year is then added. Ana-
lytical and semi-numerical studies are performed to identify how much and
how frequently predators 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 predator mortality:
low mortality parasites need be released only once a year, while high mortal-
ity 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 predator-based biocontrol, a promising alternative to
fungicides for rust control.

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

<footnote>*Corresponding author</footnote>

<footnote>Email address: clotilde.djuikem@inria.fr (Clotilde Djuikem)</footnote>
1

# 1. Introduction

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

Many of the most serious crop diseases are caused by fungi such as rusts [1]. 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 [2]. 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 [3]. 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 [4]. Mammeri et al. studied the impact of spatial heterogeneities on the spread and control of grapevine powdery mildew [5], 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 be observed for banana Sigatoka diseases [6]. 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 [7].

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 [8]. 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
38

39

40

41

42

43

44

45

46

47

the minimal predator rate required to eradicate the pest [9]. 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 [10]. 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 [11]. Semi-discrete models are hence particularly appropriate for seasonal plant and pathogen dynamics, as well as impulsive control strategies.

48

49

50

51

52

53

54

55

56

57

58

59

60

Control methods used to fight fungal diseases include chemical fungicide application, cultural practices, and the use of resistant cultivars [12]. These methods induce significant labour and/or material costs. Moreover, chemical fungicides are harmful to the environment [13], potentially also to farmers and consumers [14]; in addition, they may affect non-target organisms and induce pest resistance [15]. Therefore, researchers are currently investigating alternative control methods such as biocontrol [16, 17]. 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 [18, 19, 20, 21]. Other investigations on biocontrol focused on hyperparasites and predators. For example, several insects such as *Mycodiplosis* (Diptera) [22] and fungi such as *Lecanicillium lecanii* (previously called *Verticillium lecanii*) [23, 24] feed on rust spores.

61

62

63

64

65

66

67

68

69

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 [25, 26]. Indirect impacts include increased costs to control the disease.

70

71

72

73

74

75

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 spread of CLR (up to 68%) under in vitro conditions [21]. The fungus *Cladosporium hemileiae*, inhibits the evolution of *H. vastatrix* mycelium [27].
76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

The fungus *Lecanicillium lecanii* [28] and the insect *Mycodiplosis* [29, 30] can feed on spores. In all these biocontrol studies, the dose and time of application of the bacteria, hyperparasites or predators affected their effec- tiveness in controlling coffee leaf rust development. In this work, we chose to focus on a *Mycodiplosis*-like predator, 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 predator releases, using a mathematical modelling approach. To achieve this goal, we extend an original impulsive CLR model that we previously developed [31]. In particular, we explicitly represent the predator dynamics, the latter acting as predators of spores. Moreover, we consider in this work discrete releases of the predator instead of a continuous control.

This paper is organised as follows. Section 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 math- ematical analysis of the model using Floquet theory [32] and simulations of this impulsive model to illustrate the theoretical results. In Section 3, we introduce the biocontrol in the model, based on predator yearly releases. We semi-analytically study the stability of the controlled model to obtain the sta- bility regions as functions of parameter values. Section 4 presents the impact of multiple yearly releases of predators, with contrasting results according to the value of the predator mortality rate. The last Section 5 concludes the paper with a discussion of the main results and possible perspectives for future work.

102

## 2. Coffee leaf rust dynamics

103

### 2.1. *Modeling of CLR*

104

105

106

107

108

109

110

111

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-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 all leaves at rate $\nu$. They germinate with efficiency $\omega$ ([33, 34]) 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 on infected leaves by the fungus at rate $\gamma$ and lose their ability to infect coffee leaves at constant rate $\mu_U$. The production of flowers and berries depends on the photosynthesis, which in turn depends on the leaves, infected leaves synthesising less carbon due to rust lesions. Moreover, coffee flowers are carried at the leaf nodes on the branches and are hence also related to leaves. We translate this in the model by flower production at constant rates $\delta_S$ and $\delta_I$ "by" susceptible and infected leaves respectively, with $\delta_S > \delta_I$. The flowers become berries at rate $\theta$. Flowers and berries have a natural mortality rate $\mu_F$ and $\mu_B$, respectively. 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 proportions $\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 1, with the continuous dynamics during the rainy season and switching event for the dry season.

Using flowchart in Figure 1, we can write the following impulsive differential system:
<figure><img src="image_4.png" /><figcaption>Figure 1: Diagram of the CLR multi-seasonal model. Model equations are given in system (1). State variables are: healthy leaves (S), infected leaves (I), uredospores (U), flowers (F) and berries (B). T corresponds to the length of rainy season and n the year number.</figcaption></figure>

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, \\
\dot{F} &= \delta_S S + \delta_I I - (\theta + \mu_F)F, \\
\dot{B} &= \theta F - \mu_B B;
\end{aligned}
\right.
\qquad (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.
$$

<footnote n="141">141 where $S(t)$, $I(t)$, $U(t)$, $F(t)$, $B(t)$ and $N(t) = S(t) + I(t)$ represent the</footnote>

<footnote n="142">142 number of healthy leaves, infected leaves, uredospores, flowers, berries and</footnote>
total number of leaves, respectively, at time $t$; with $0 < \varphi_S, \varphi_I, \varphi_U < 1$ and $\varphi_U$ close to 0.

## 2.2. Mathematical analysis

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

### 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 (1) when $t \ne nT$, and consider the initial conditions $(S(0), I(0), U(0), F(0), B(0)) \in \mathbb{R}_+^4 \times \mathbb{R}_+^4$.

The solutions of system (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 (1), we obtain the dynamics of total leaves which satisfies:

$$ 
\left\{ 
\begin{array}{l} 
\dot{N} = \Lambda - \mu N - dI, \quad t \ne nT; \\ 
N(nT^+) = \varphi_S S(nT) + \varphi_I I(nT). 
\end{array} 
\right. 
\qquad (2) 
$$

From equation (2), and using the non-negativity of variable $I$, we obtain

$$ 
\left\{ 
\begin{array}{l} 
\Lambda - (\mu + d)N \le \dot{N} \le \Lambda - \mu N, \quad t \ne nT; \\ 
0 \le N(nT^+) \le N(nT). 
\end{array} 
\right. 
$$

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

$$ \frac{\Lambda}{\mu + d} + \left( N(nT^+) - \frac{\Lambda}{\mu + d} \right) e^{-(\mu+d)(t-nT)} \le N(t) \le \frac{\Lambda}{\mu} + \left( N(nT^+) - \frac{\Lambda}{\mu} \right) e^{-\mu(t-nT)} $$

It follows that $\min(\Lambda/(\mu + d), N(nT^+)) \le N(t) \le \max(\Lambda/\mu, N(nT^+))$ for $t \in (nT, (n+1)T]$. The impulse with $0 < \varphi_S, \varphi_I < 1$ then imposes:

$$ \min(\varphi_S, \varphi_I) N((n+1)T) \le N((n+1)T^+) \le N((n+1)T). $$

and hence

$$ \min(\varphi_S, \varphi_I) \min(\Lambda/(\mu + d), N(nT^+)) \le N(t) \le \max(\Lambda/\mu, N(nT^+)) \quad (3) $$
for $t \in (nT, (n+1)T^+]$. Applying the R.H.S. recursively, we get

$$N(t) \leq \max(\Lambda/\mu, N(0)) = \Gamma_N$$

162 The recursive application of the L.H.S. of (3) does not yield a single lower bound, but ensures that there is a positive lower bound to $N$ over each $T$-interval.

165 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:

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

$$\begin{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)$$

167 Identifying the values where the upper bounds of (4) and (5) 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

$$\begin{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 (6)$$

170 Identifying the values where the upper bounds of (6) 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\{ \begin{array}{l} (S, I, U, F, B) \in \mathbb{R}_+^5 \mid S(t) + I(t) \leq \frac{\Lambda}{\mu}, U(t) \leq \frac{\gamma\Lambda}{\mu(\nu + \mu_U)}, \\ F(t) \leq \frac{\Lambda(\delta_S + \delta_I)}{\mu(\theta + \mu_F)}, B(t) \leq \frac{\Lambda\theta(\delta_S + \delta_I)}{\mu(\theta + \mu_F)\mu_B} \end{array} \right\}$$

173 is positively invariant for impulsive system (1).

174 We have hence shown the non-negativity and boundedness of the solutions. Moreover, the smooth properties of the right side of system (1) guarantee the existence and uniqueness of the solutions of this system.
Since the flower and berry state variables F and B are not present in the other equations of system (1), we do not 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 (7) $$

### 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 (7), 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 (8) $$

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 (9) $$

The fixed point of equation (9) is given by

$$ S(nT^+) = \frac{\Lambda\varphi_S(1 - e^{-\mu T})}{\mu(1 - \varphi_S e^{-\mu T})} > 0. \quad (10) $$
187 Substituting the value of $S(nT^+)$ into equation (8), for $t \in (nT, (n+1)T]$,
188 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 (11)$$

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

190 For stability analysis, we use Floquet theory and its implementation on
191 imolusive systems as developed in [32]; it transforms the analysis of the im-
192 pulsive system into that of a discrete mapping over the period. We first
193 study the local stability of the PDFS $X^T(t)$ using small amplitude per-
194 turbation methods. For that, let us denote $\tilde{X}(t) = X(t) - X^T(t)$, where
195 $X(t) = (S(t), I(t), U(t))^T$ and $\tilde{X}(t)$ is understood to be small amplitude
196 perturbations. Substituting the expression of $\tilde{X}(t)$ in system (7) 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}, && 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.
\quad (12) $$

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

198 The linearization of system (12) in the neighbourhood of $0_{\mathbb{R}^3}$ is

$$ 
\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.
\quad (13) $$

199 where

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

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

$$\tilde{X}(t) = \Phi_A(t)\tilde{X}(0^+), \quad (14)$$
201 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} $$

202 and

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

203 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. $$

204 Observe that $\psi_{22}(t) > 0$, indeed $\beta^2 > (k_1 - k_2)^2$, which implies $\beta - k_1 + k_2 > 0$ and $\beta + k_1 - k_2 > 0$.

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

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

208 Using [32], 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}. $$

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

211 According to [32], 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}. \quad (16) $$

214 The PDFS $X^T(t)$ is locally asymptotically stable if the Floquet multipliers of the monodromy matrix $M$ stay inside the unit circle [32], which
216 is equivalent to the Floquet multipliers of sub-matrix $M_1$ staying inside the
217 unit circle. To obtain this result, the submonodromy matrix $M_1$ needs to
218 satisfy the following Jury conditions [35]:

$$ egin{cases} -\operatorname{tr}(M_1) - \det(M_1) < 1, \\ \det(M_1) < 1, \\ \operatorname{tr}(M_1) - \det(M_1) < 1, \end{cases} \tag{17} $$

219 where

$$ egin{cases} \operatorname{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} $$

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

223 **Lemma 1.** The PDFS $X^T(t)$ is locally asymptotically stable provided that
224 $\mathcal{R} < 1$ and unstable otherwise, 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]. \tag{18} $$

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

227 **Theorem 2.** The PDFS $X^T(t)$ of system (7) is globally asymptotically stable
228 if $\mathcal{R} < 1$ and unstable otherwise.

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

$$ egin{cases} \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; \\ I(nT^+) = \varphi_I I(nT); \\ U(nT^+) = \varphi_U U(nT). \end{cases} \tag{19} $$

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

$$ egin{cases} 
\dot{I}_1 = \omega\nu U_1 - (\mu + d)I_1, & t \neq nT; \\
\dot{U}_1 = \gamma I_1 - (\nu + \mu_U)U_1, & t \neq nT; \\
I_1(nT^+) = \varphi_I I_1(nT); \\
U_1(nT^+) = \varphi_U U_1(nT). 
\end{cases} 
\quad (20) $$

234 Applying Kamke's theorem [36], 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^+). $$

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

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

240 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:

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

244 and

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

245 One can observe that $S_1(t) \to S_{\varepsilon_U}^T(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]$ and $S_2(t) \to S^T(t)$ asymptotically, where $S^T(t)$ is given in equation (11). This implies that

$$ S_{\varepsilon_U}^T(t) \le S(t) \le S^T(t). \quad (21) $$

248 Since we can take $\varepsilon_U$ as small as we want since $U(t) \to 0$, we obtain $S_{\varepsilon_U}^T(t) \to S^T(t)$ when $\varepsilon_U \to 0$. Finally, using (21), one has that $S(t) \to$
250 $S^T(t)$. This 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)) \rightarrow (S^T(t), 0, 0) $$

252 This concludes the proof of global asymptotic stability. □

## 2.3. Numerical simulations of CLR dynamics

254 Herein, we present the results of numerical simulations of system (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 [37] 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). \quad (22) $$

264 Table 1 summarises the parameter definitions and values used in this paper.

265 In Theorem 2, 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 2, drawn for $\mathcal{R} = 2.53 > 1$. 266

267 The dashed curve of Figure 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. 270 271 272 273 274 275 276 277 278 279
<table><caption>Table 1: Description and values of parameters for system (23).</caption><thead><tr><th>Symbol</th><th>Biological meaning</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td><i>T</i></td><td>Rainy season duration</td><td>250 days</td><td>[37]</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 efficiency</td><td>0.065 leaf/spore</td><td>Assumed</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>[38]</td></tr><tr><td><i>d</i></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>[38]</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>Survival proportion of S</td><td>0.7</td><td>Assumed</td></tr><tr><td>φ<sub>I</sub></td><td>Survival proportion of I</td><td>0.4</td><td>Assumed</td></tr><tr><td>φ<sub>U</sub></td><td>Survival proportion 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>[39]</td></tr><tr><td>δ<sub>I</sub></td><td>Flower production rate by I</td><td>0.04 flowers/leaf.day</td><td>[39]</td></tr><tr><td>μ<sub>F</sub></td><td>Mortality rate of F</td><td>2, 4 × 10<sup>-3</sup>/day</td><td>[40]</td></tr><tr><td>θ</td><td>Maturation rate of F</td><td>0.0055/day</td><td>[41]</td></tr><tr><td>μ<sub>B</sub></td><td>Mortality rate of B</td><td>2, 4 × 10<sup>-3</sup>/day</td><td>[40]</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><i>a</i></td><td>Consumption rate</td><td>variable spores/predator.day</td><td>Assumed</td></tr><tr><td><i>e</i></td><td>Biomass transformation rate</td><td>0.7 predators/spore</td><td>Assumed</td></tr><tr><td><i>K</i></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>Survival proportion of P</td><td>0.3</td><td>Assumed</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 2: Impact of CLR on the production of coffee berries. The plots present the trajectories of impulsive differential system (7) 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) infected leaves I, (c) uredospores U and (d) berries B. Remaining parameter values are given in Table 1 and initial conditions by (22).</figcaption></figure>
## 3. Biocontrol of CLR using predators

To limit the impact of CLR on coffee berry production, we now introduce a biocontrol agent in the model, more specifically a predator such as *Mycodiplosis*, which consumes uredospores of *H. vastatrix*. The new model is then an extension of system (7), to which we add an equation for the predators (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 predators in order to control CLR. In this section, we assume that a quantity *Λ*<sub>*P*</sub> of predators is released once a year, at the start of the production season. This translates into a *Λ*<sub>*P*</sub> jump in the predator population at each switching moment. The model with biocontrol is then given by

$$ 
\begin{array}{l}
\text{Rainy season for } t \neq nT: \\
\left\{
\begin{array}{l}
\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} = e a P \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{array}
\right.
\end{array}
\qquad (23)
$$

Dry season:

$$ 
\left\{
\begin{array}{l}
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{array}
\right.
$$

System (23) has the same properties as system (7) for the existence, uniqueness, non-negativity and boundness of solutions. We pursue the mathematical analysis and simulations of system (23) without the $\dot{F}$ and $\dot{B}$ equations, since berry variable *B* is not present in the other equations.
### 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 (11) while the predator dynamics are

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

Solving equation (24) 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 predator at the cPDFS yields

$$ P^T(t) = \frac{\Lambda_P}{1 - \varphi_P e^{-\mu_P t}} e^{-\mu_P(t-nT)}. \qquad (25) $$

Thus, the cPDFS for system (23) is $Y^T(t) = (S^T(t), 0, 0, P^T(t))$, where $S^T(t)$ and $P^T(t)$ are given in equation (11) and (25), respectively. We have the following result for the global stability of cPDFS.

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

*Proof.* 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
311 expression of $\tilde{Y}(t)$ in the system (7) 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.
\tag{26}
$$

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

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

$$
\left\{
\begin{array}{ll}
\dot{\tilde{Y}}(t) = G(t)\tilde{Y}(t), & t \neq nT; \\
\tilde{Y}(nT^{+}) = \text{diag}(\varphi_S, \varphi_I, \varphi_U, \varphi_P)\tilde{Y}(nT);
\end{array}
\right.
\tag{27}
$$

314 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}.
$$

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

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

317 where $\Phi_G(0) = I$. Re-ordering the variables as follows, $(\tilde{S}, \tilde{P}, \tilde{I}, \tilde{U})$, we note
318 that the system is block-triangular. The stability is then determined by
319 considering separately $(\tilde{S}, \tilde{P})$ and the pair $(\tilde{I}, \tilde{U})$. The Floquet multipliers
320 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 (27) reduces to the stability of the following sub-system

$$ \begin{cases} \dot{\tilde{Y}}_1(t) = G_1(t)\tilde{Y}_1(t), & t \neq nT; \\ \tilde{Y}(nT^+) = \text{diag}(\varphi_I, \varphi_U)\tilde{Y}(nT); \end{cases} \qquad (28) $$

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 (16), which is stable if $\mathcal{R} < 1$, where $\mathcal{R}$ is given in equation (18). Using the fact that system (28) 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. $\square$

## 3.2. Semi-numerical analysis of the controlled model

In order to analyse system (23) 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 (28) 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 (23). We term this method “semi-numerical” as it requires numerical computations, but not extensive simulations of the system [42].

In Figure 3, we plot the threshold level of the spectral radius $\mathcal{R}_c = 1$ of controlled system (23) for parameter pairs $(\Lambda_P, a)$ (subplot (a)) and $(\nu, \gamma)$
<figure><img src="image_6.png" /><figcaption>Figure 3: Stability regions of periodic disease-free solutions as functions of model pa- rameters. Blue curves represent the threshold $\mathcal{R}_C = 1$ separating the stable and unstable regions of the model with predator (23), for different values of parameter pair $(\Lambda_P, a)$ in subplot (a) and $(\nu, \gamma)$ in subplot (b); moreover, the black curve $\mathcal{R} = 1$ in subplot (b) separates the stable and unstable regions of the model without predator (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 1.</figcaption></figure>

(subplot (b)), as well as the threshold stability level $\mathcal{R} = 1$ of uncontrolled system (1) for parameter pair $(\nu, \gamma)$. In subplot (a), the regions below and above the blue curve are, respectively, the unstable and stable regions of the cPDFS of controlled system (23). 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 extinc- tion 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

In subplot (b), as in subplot (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.
### 3.3. Asymptotic behaviour for yearly releases

Herein, we investigate the asymptotic behaviour of system (23) 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 3. We point out that the gray dot in subplot (b) of Figure 3 is in the instability region of system (23) without control which means that the CLR would persist in the coffee plantation if unchecked. We consider the initial condition given in equation (22), 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 (29) $$

Figure 4 compares the temporal dynamics of system (1) without control (black curves) and model (23) 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 4(a)), while the number of infected leaves and uredospores dwindles until extinction (see Figures 4(b) and (c)). During the first year, the quantity of predators (predators) 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 predators reduces drastically CLR in the coffee plantation.

## 4. Multiple releases of predators per year

In this section, we suppose that the yearly quantity $\Lambda_P$ of predators 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
<figure><img src="image_7.png" /><figcaption>Figure 4: Impact of biocontrol on the system dynamics. 6-year simulations of model (1) without predator (black curve) and model (23) with predator (blue curve), using initial conditions (29): (a) healthy leaves S; (b) infected leaves I; (c) uredospores U and (d) predators P. All parameter values are given in Table 1 and correspond to $R = 2.53 > 1$.</figcaption></figure>

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. 
\quad (30) 
$$

Predator releases for $j \in \{1, ..., 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); \quad 23 \\
P(nT^+) &= \varphi_P P(nT) + \frac{\Lambda_P}{m}; \\
F(nT^+) &= 0 \\
B(nT^+) &= 0.
\end{aligned}
\right. 
$$
System (30) has the same properties as system (7) for the existence, uniqueness, non-negativity, and boundedness of solution. We made the mathematical analysis of system (30) without $\dot{F}$ and $\dot{B}$ equations, since B is not present in other equations.

### 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 (30) occurs when $I = 0$ and $U = 0$, so healthy leaf dynamics are unchanged compared to the PDFS and still satisfy (11); the predator dynamics are

$$ egin{cases} 
\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{cases} 
\quad (31) $$

The expression of the m-cPDFS is established in Appendix B.1. We obtain $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 (11) 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 \left( t - \frac{jT}{m} \right)}, \quad (32) $$

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 (33) $$

**Lemma 4.** 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.* The proof is similar to the proof of Lemma 3 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 predators present in the plantation for the m-cPDFS. This quantity depends on $m$, the number of predator releases per year. Its expression, established in Appendix B.2, is

$$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 (34)$$

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 predators 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 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.

## 4.2. Impact of release frequency on CLR control

Herein, we present the impact of the yearly predator 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 3), for different values of the release number $m$. Secondly, we present simulations of the dynamic behaviour of system (30).

In all simulations, the epidemiological parameters are set to their value in Table 1. The only parameters that vary are the ones related to the predators: $\Lambda_p, a, \mu_p$ and $m$. As in Figure 5, we use two contrasted values for the predator
<figure><img src="image_8.png" /><figcaption>Figure 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 (30) 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 (33), so the solutions are periodic. The yearly released quantity is $\Lambda_P = 3000$ predators. Remaining parameter values are given in Table 1.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 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 (30) 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 1 and initial conditions by (35).</figcaption></figure>

431 mortality $\mu_p$, the low and high values representing respectively a favourable and unfavourable environment for the predators.

432 Moreover, we consider that an initial predator release occurs at time 0 with the usual initial conditions (22) 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 (35) $$

### 4.2.1. Impact on stability

433 Herein, we present the impact of the yearly release frequency on the m-cPDFS stability of controlled system (30), for the two contrasted predator mortality rates depicted in Figure 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.

434 Case 1 – Less frequent is more effective (low mortality rate $\mu_P = 0.003$/day). In subplot (a) of Figure 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 (30), 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.

440 Subplot (b) of Figure 6 illustrates the long term dynamics (over 12 years) of the infected leaves of system (30), with initial conditions (35), for pa-
<figure><img src="image_10.png" /></figure>

Figure 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 (30) 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 1 and initial conditions by (35)

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 predator mortality is very low, the best strategy for eradicating CLR is to release predators 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 (34), which is higher for lower frequencies.

*Case 2 – More frequent is more effective (high mortality rate $\mu_P = 0.1$/day).* Figure 7 is built as Figure 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 sub-plot (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 predator population is large at the
### 4.2.2. Impact on transient dynamics

<footnote n="468">468 beginning of the year, but declines sharply due to the high mortality rate (see Figure 5(b)), so that CLR is not kept in check. By contrast, for $m = 5$, 469 predators do not stay close to 0 for long so that no large increase of $I$ occurs. 470</footnote>

<footnote n="471">471 For high predator mortality, the best strategy to control CLR is to release 472 predators five times per year, meaning that more frequent releases are better. 473 This conclusion was not obvious, considering that the control intensity proxy 474 computed in equation (34) is lower for more frequent releases. 475</footnote>

<footnote n="476">476 Herein, we still look at the impact of the yearly release frequency, for 477 the same contrasted predator mortality rates, but we focus on the first six 478 years. Four frequency values are considered, namely $m = 1, 2, 5, 50$ releases 479 per year. We choose the predator parameter values corresponding to a stable 480 m-cPDFS for all release frequencies considered. 481</footnote>

<footnote n="481">481 Case 1 – Less frequent is faster (low mortality rate $\mu_P = 0.003$/day). Just 482 as Figure 6(a), subplot (a) of Figure 8 displays the instability (below) and 483 stability (above) regions of the m-cPDFS of controlled system (30). These 484 regions are separated by $R_{c,1} = 1$ for $m = 1$ release/year (blue curve), $R_{c,2} = 485$ 486 1 for $m = 2$ releases/year (brown curve), $R_{c,5} = 1$ for $m = 5$ releases/year 487 (magenta curve) and $R_{c,50} = 1$ for $m = 50$ releases/year (light blue curve). 488 Subplot (a) confirms that for low predator mortality, the best strategy for 489 controlling CLR is to release the predators once a year at the beginning of 490 the season. 491</footnote>

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

<footnote n="503">503 Subplot (c) shows that the number of berries at the end of each year 504 decreases when $m$ increases. This confirms our previous result: the best 505 strategy for low predator mortality is to few releases per year. 506</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 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 (30) 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 1 and initial conditions by (35).</figcaption></figure>
<figure><img src="image_12.png" /></figure>

Figure 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 (30) 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 1 and initial conditions by (35).

*Case 2 – Best strategy depends on the horizon (high mortality rate $\mu_P = 0.1$). As expected for high predator mortality, in subplot (a) of Figure 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 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 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 predator 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 predator mortality on berries depicted in subplot (c), we consider the number of berries at the end of each season and thus we obtain Table 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 predator

<table><caption>Table 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 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>

In terms of berry production, the best strategy from the first year to the second year is to release the predators 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.

### 4.2.3. Summary

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

## 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 [43]. 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 [44]. 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 [45]. 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 [46].

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 predators, as recommended by Zambolin [33] for alternative biocontrol methods; this resulted in a controlled impulsive model. The semi-numerical analysis of this controlled impulsive model allowed to conclude that predator-based biocontrol can drastically reduce coffee leaf rust. Moreover, we studied how successful biocontrol implementation depends on predator characteristics, in particular its mortality.

We set the study in a framework where the yearly released quantity of predators 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 predator mortality rates were considered to represent different environmental situations and predator characteristics. High mortality occurs if the predator is a specialist parasite like *Mycodiplosis* [22], which means that it only consumes CLR uredospores, or if it is a generalist predator whose alternate food sources are absent; in both cases, the predator quickly dies out in absence of CLR. Low mortality occurs for generalist predators in a rich environment. Most predators 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 [28, 47], and it is also a pathogen of plant-parasitic nematodes [48]. Hence, both low and high mortality situations could occur in the field.
For low mortality rates, the best strategy is to release the predators once a year at the beginning of the production season. This strategy is efficient because it is the one that ensures the largest average predator density over the year, predators 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 predators, 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 predator 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 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 predators needed to control CLR (see Figure 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 [22], 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 2, where we would choose to have a large single release the first year, and potentially five releases per year after
619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

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, predators 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 predators 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 [49], such as temperature [50]. 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 predators. 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.

To conclude, we were able to obtain analytical results on the qualitative behaviour of our model without control. Compared to purely numerical analyses, our semi-numerical analyses hold for large parameter regions. This is particularly true for one of our main results, regarding the best predator release frequency: for a low predator mortality, one release per season is better than five, and vice versa for high predator mortality. However, the minimum yearly quantity of predators that need to be released to control the disease depends on the model parameter values. Therefore, even if this research provides valuable insights on predator-based biocontrol, further research, most notably around parameter calibration, would be needed to complete the step from theory to practice. Furthermore, 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 sim-
ulations should be performed for each case.

## Acknowledgments

This work is supported by EPITAG, an Inria Associate team part of the LIRIMA (https://team.inria.fr/epitag/)

## Appendix A. Parameter values

We present in this appendix how parameter values presented in Table 1 were set.

• The maturation of flowers takes 180 days, which represents stage I + stage II + stage III presented by Torres Castillo et al. [41], so the maturation rate is $θ = 1/180 \approx 0.0055$/ day.

• In the coffee tree, a node produces 30 flowers during the flowering period [40, page 59]. Each node has 2 leaves, so a leaf corresponds to 15 flowers. The production of flowers in the model is hence $δ_S = 15/180 \approx 0.08$ flowers/leaf.day. We suppose that $δ_I = 0.04 < \text{flowers/leaf.day}$ due to the disease.

• Out of 100 flowers, 65 set fruit after pollination [40, page 28]. This implies that $e^{-180μ_F} = 0.65$, so $μ_F = 2.4 \times 10^{-3}$ per day. We suppose that the berries have the same mortality rate, so $μ_B = 2.4 \times 10^{-3}$ per day.

• Leaf mortality during the rainy season is fairly low. The rainy season lasts 250 days. We assumed that leaves survive on average 194 days. Then $μ = 1/194$ day = $0.0034$ day$^{-1}$.

• Leaf mortality during the dry season is naturally occurring, caused by either the dry season or by man pruning the trees. We set the leaf survival to the dry season to $φ_S = 0.7$ for susceptible leaves and $φ_I = 0.4$ for infected leaves, branches with CLR-affected leaves being preferentially targeted for pruning.

• Lesions form on the leaves 2 to 3 weeks after the spores germinate. These lesions will directly cause the leaf to die; it may take longer, but we have assumed that the average time is $(14+21)/2=18$ days. Then $d = 1/18$ days = $0.056$ day$^{-1}$.
688 • The values of parameters $\Lambda$, $\mu_U$, $\phi_U$, $\phi_P$, and $K$ were assumed.

689 • The parameter values $\mu_P$ are chosen to have high and low mortality of

690 predator

691 Appendix B. Properties of the m-cPDFS

692 Appendix B.1. Expression of the m-cPDFS

693 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\left(t - \frac{jT}{m}\right)}. \qquad (B.1)$$

694 Using equation (B.1) and the second equation of (31) 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}.$$

695 For $j = 0$

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

696 For $j = 1$

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

697 For $j = 2$

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

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

$$\begin{aligned} P\left(\frac{jT^+}{m}\right) &= P\left(0^+\right) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \sum_{i=0}^{j-1} e^{-\mu_P \frac{iT}{m}}; \\ &= P\left(0^+\right) e^{-\mu_P \frac{jT}{m}} + \frac{\Lambda_P}{m} \left(1 - e^{-\mu_P \frac{jT}{m}}\right) \end{aligned}$$
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 (31), 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) =$

$P(S^T(t), 0, 0, P_{m,j}^T(t))$, with $S^T(t)$ given by equation (11) 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} \left( t - \frac{jT}{m} \right)}, \quad (B.2)$$

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

## Appendix B.2. Expression of the yearly average number of predators

The yearly average number of predators 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 (32), we obtain

$$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} \left( t - \frac{jT}{m} \right)} (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} \left( t - \frac{jT}{m} \right)} (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].$$
709 Replacing $P(0^+)$ by its expression in equation (33), we obtain

$$
\begin{align*}
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{align*}
$$

710 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}.
$$

711 References

712 [1] G. N. Agrios, Plant pathology, Elsevier, 2005.

[2] S. Pivonia, 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) (2006) 400–407.
doi:10.1094/PHYTO-96-0400.

[3] L. Rimbaud, J. Papaïx, L. G. Barrett, J. J. Burdon, P. H. Thrall, Mo-
saics, mixtures, rotations or pyramiding: What is the optimal strategy
to deploy major gene resistance?, Evolutionary Applications 11 (10)
(2018) 1791–1810. doi:10.1111/eva.12681.
721

722

723

[4] R. A. Fleming, The potential for control of cereal rust by natural enemies, Theoretical Population Biology 18 (3) (1980) 374-395. doi: 10.1016/0040-5809(80)90060-X.

724

725

726

727

728

729

[5] Y. Mammeri, J. Burie, M. Langlais, 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 (2014) 178-191, 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.

730

731

732

733

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

734

735

736

737

738

[7] M.-L. Desprez-Loustau, F. M. Hamelin, B. Marçais, The ecological and evolutionary trajectory of oak powdery mildew in Europe, in: K. Wilson, A. Fenton, D. Tompkins (Eds.), Wildlife Disease Ecology: Linking Theory to Data and Application, Ecological Reviews, Cambridge University Press, 2019, Ch. 15, p. 429-457. doi:10.1017/9781316479964.015.

739

740

741

742

[8] I. Tankam-Chedjou, F. Grognard, J. J. Tewa, S. Touzeau, Optimal and sustainable management of a soilborne banana pest, Applied Mathematics and Computation 397 (2021) 125883. doi:10.1016/j.amc.2020.125883.

743

744

745

[9] S. Nundloll, L. Mailleret, F. Grognard, Two models of interfering predators in impulsive biological control, Journal of Biological Dynamics 4 (1) (2010) 102-114. doi:10.1080/17513750902968779.

746

747

748

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

749

750

751

752

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

[12] E. G. Sharvelle, Plant disease control, The AVI Publishing Company, 1979.

755
756
757
758

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

759
760
761
762

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

763
764
765
766

[15] M. J. Jeger, P. J. Wijngaarden, 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) (2008) 621–632. doi:10.1016/j.jtbi.2008.02.023.

767
768
769

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

770
771
772
773

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

774
775
776

[18] J. Shafi, H. Tian, M. Ji, Bacillus species as versatile weapons for plant pathogens: a review, Biotechnology & Biotechnological Equipment 31 (3) (2017) 446–459. doi:10.1080/13102818.2017.1286950.

777
778
779

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

780
781
782

[20] V. Singh, B. Deverall, Bacillus subtilis as a control agent against fungal pathogens of citrus fruit, Transactions of the British Mycological Society 83 (3) (1984) 487–490. doi:10.1016/S0007-1536(84)80045-5.
783

784

785

[21] S. Daivasikamani, et al., Biological control of coffee leaf rust pathogen, hemileia vastatrix berkeley and broome using bacillus subtilis and pseu- domonas fluorescens., Journal of Biopesticides 2 (1) (2009) 94-98.

786

787

788

[22] D. Henk, D. Farr, M. Aime, Mycodiplosis (diptera) infestation of rust fungi is frequent, wide spread and possibly host specific, Fungal Ecology 4 (4) (2011) 284-289. doi:10.1016/j.funeco.2011.03.006.

789

790

791

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

792

793

794

795

[24] W. Gams, A. Van Zaayen, Contribution to the taxonomy and pathogenicity of fungicolous verticillium species. i. taxonomy, Netherlands Journal of Plant Pathology 88 (2) (1982) 57-78. doi:10.1007/BF01977339.

796

797

[25] G. Grée, Epidemiology of coffee leaf rust in the Eastern Highlands, Coffee Research Institute Newsletter 2 (1993) 16-20.

798

799

800

801

802

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

803

804

805

806

[27] R. L. Steyaert, Cladosporium hemileiae n. spec. un parasite de l'hemileia vastatrix berk. et br., Bulletin de la Société Royale de Botanique de Belgique/Bulletin van de Koninklijke Belgische Botanische Vereniging 63 (Fasc. 1) (1930) 46-48.

807

808

809

810

[28] J. Vandermeer, I. Perfecto, H. Liere, Evidence for hyperparasitism of coffee rust (Hemileia vastatrix) by the entomogenous fungus, Lecanicil- lium lecanii, through a complex ecological web, Plant Pathology 58 (4) (2009) 636-641. doi:10.1111/j.1365-3059.2009.02067.x.

811

812

813

814

815

[29] E. Santiago-Elena, E. J. Zamora-Macorra, M. Zamora-Macorra, K. G. Elizalde-Gaytan, Interaction between mycodiplosis and hemileia vasta- trix in three scenarios of coffee crop management (coffea arabica), Revista mexicana de fitopatología 38 (3) (2020) 320-336. doi:10.18781/r.mex.fit.2005-2.
816

817

818

819

820

821

[30] Z. Hajian-Forooshani, I. S. Rivera Salinas, E. Jiménez-Soto, I. Perfecto, J. Vandermeer, Impact of Regionally Distinct Agroecosystem Communities on the Potential for Autonomous Control of the Coffee Leaf Rust, Environmental Entomology 45 (6) (2016) 1521-1526. arXiv:https://academic.oup.com/ee/article-pdf/45/6/1521/8660201/nvw125.pdf, doi:10.1093/ee/nvw125.

822

823

824

825

826

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

827

828

[32] D. Bainov, P. Simeonov, Impulsive differential equations: periodic solutions and applications, Vol. 66, CRC Press, 1993.

829

830

831

[33] L. Zambolim, Current status and management of coffee leaf rust in brazil, Tropical Plant Pathology 41 (1) (2016) 1-8. doi:10.1007/s40858-016-0065-9.

832

833

834

[34] R. W. Rayner, Germination and penetration studies on coffee rust (Hemileia vastatrix B. & Br.), Annals of Applied Biology 49 (3) (1961) 497-505. doi:10.1111/j.1744-7348.1961.tb03641.x.

835

[35] H. Wang, Mathematical modeling I-preliminary, Bookboon, 2012.

836

837

[36] W. A. Coppel, Stability and asymptotic behavior of differential equations, Heath, 1965.

838

839

840

[37] J.-B. Suchel, Quelques remarques à propos de la répartition des pluies au Cameroun durant la période sèche 1969-1973, Hommes et Terres du Nord 3 (1) (1983) 24-28.

841

842

843

[38] K. R. Bock, Dispersal of uredospores of Hemileia vastatrix under field conditions, Transactions of the British Mycological Society 45 (1) (1962) 63-74. doi:10.1016/S0007-1536(62)80035-7.

844

845

846

[39] G. Champéroux, Manuel du planteur de café laotien, CIRAD-IRCC, IRCC, Montpellier, 1991.
URL http://agritrop.cirad.fr/345484/
847

848

849

[40] M. N. Clifford, K. C. Willsons (Eds.), Coffee: botany, biochemistry and production of beans and beverage, no. 1, Springer New York, NY, 1985. doi:https://doi.org/10.1007/978-1-4615-6657-1.

850

851

852

853

854

855

856

857

[41] N. E. Torres Castillo, E. M. Melchor-Martínez, J. S. Ochoa Sierra, R. A. Ramirez-Mendoza, R. Parra-Saldívar, H. M. Iqbal, Impact of climate change and early development of coffee rust – an overview of control strategies to preserve organic cultivars in Mexico, Science of The Total Environment 738 (2020) 140225. doi:https://doi.org/10.1016/j.scitotenv.2020.140225. URL https://www.sciencedirect.com/science/article/pii/S0048969720337463

858

859

860

[42] B. Ghosh, F. Grognard, L. Mailleret, Natural enemies deployment in patchy environments for augmentative biological control, Applied Mathematics and Computation 266 (2015) 982–999.

861

862

863

[43] C. Djuikem, F. Grognard, R. Tagne Wafo, S. Touzeau, S. Bowong, Modelling coffee leaf rust dynamics to control its spread, Mathematical Modelling of Natural Phenomena 16 (2021) 26. doi:10.1051/mmnp/2021018.

864

865

866

867

[44] D. P. Bebber, Á. D. Castillo, 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) (2016) 20150458. doi:10.1098/rstb.2015.0458.

868

869

870

[45] J. Vandermeer, P. Rohani, The interaction of regional and local in the dynamics of the coffee rust disease, arXiv 1407.8247 (2014). doi:10.48550/arXiv.1407.8247.

871

872

873

874

[46] J. Vandermeer, Z. Hajian-Forooshani, I. Perfecto, The dynamics of the coffee rust disease: an epidemiological approach using network theory, European Journal of Plant Pathology 150 (4) (2018) 1001–1010. doi: 10.1007/s10658-017-1339-x.

875

876

877

878

[47] R. Setiawati, A. Widiastuti, A. Wibowo, A. Priyatmojo, Variability of lecanicillium spp. mycoparasite of coffee leaf rust pathogen (hemileia vastatrix) in Indonesia., Pakistan Journal of Biological Sciences: PJBS 24 (5) (2021) 588–598. doi:10.3923/pjbs.2021.588.598.
879

880

881

882

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

883

884

885

[49] J. Avelino, L. Willocquet, S. Savary, Effects of crop management patterns on coffee rust epidemics, Plant pathology 53 (5) (2004) 541–547. doi:10.1111/j.1365-3059.2004.01067.x.

886

887

888

889

[50] E. De Jong, A. Eskes, J. Hoogstraten, J. Zadoks, Temperature requirements for germination, germ tube growth and appressorium formation of urediospores of hemileia vastatrix, Netherlands Journal of Plant Pathology 93 (2) (1987) 61–71. doi:10.1007/BF01998091.