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

HAL
open science

# Mathematical modelling and control of perennial plant phytopathogens

Clotilde Djuikem

► To cite this version:

Clotilde Djuikem. Mathematical modelling and control of perennial plant phytopathogens. Automatic Control Engineering. Université Côte d'Azur, 2023. English. NNT: 2023COAZ4002. tel-04026720

HAL Id: tel-04026720

https://theses.hal.science/tel-04026720

Submitted on 13 Mar 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.
UNIVERSITÉ
CÔTE D'AZUR

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

ÉCOLE DOCTORALE
SCIENCES ET
TECHNOLOGIES DE
L'INFORMATION ET DE
LA COMMUNICATION

## THÈSE DE DOCTORAT

# Modélisation mathématique et contrôle de phytopathogènes de plantes pérennes

**Clotilde DJUIKEM**

BIOCORE, Centre Inria d'Université Côte d'Azur

**Présentée en vue de l'obtention du grade de docteur en**

Automatique, Traitement du Signal et des Images
d'Université Côte d'Azur

**Dirigée par :**

Jean-Luc GOUZÉ

**Co-encadrée par :**

Frédéric GROGNARD
Suzanne TOUZEAU
Samuel BOWONG

**Soutenue le :** 5 janvier 2023

**Devant le jury, composé de :**

Julien ARINO, Professeur, University of Manitoba, Canada

Samuel BOWONG, Professeur, Université de Douala, Cameroun

Arnaud DUCROT, Professeur, Université Le Havre Normandie

Jean-luc GOUZÉ, Directeur de recherche, Inria

Frédéric GROGNARD, Chargé de recherche, Inria

Fabien HALKETT, Chargé de recherche, INRAE

Béatrice LAROCHE, Directrice de recherche, INRAE

Jean Jules TEWA, Professeur, Université de Yaoundé I, Cameroun

Suzanne TOUZEAU, Chargée de recherche, INRAE

Inria
# Modélisation mathématique et contrôle de phytopathogènes de plantes pérennes

## Jury :

Président du jury

Julien ARINO, Professeur, University of Manitoba, Canada

Rapporteurs :

Arnaud DUCROT, Professeur, Université Le Havre Normandie

Jean Jules TEWA, Professeur, Université de Yaoundé I, Cameroun

Examineurs/Examinatrice :

Samuel BOWONG, Professeur, Université de Douala, Cameroun

Jean-Luc GOUZÉ, Directeur de recherche, Inria

Fabien HALKETT, Chargé de recherche, INRAE

Béatrice LAROCHE, Directrice de recherche, INRAE

Suzanne TOUZEAU, Chargée de recherche, INRAE

Invité :

Frédéric GROGNARD, Chargé de recherche, Inria
Dedication

To my mother NGOULLA Thérèse.
Acknowledgments

This work like others could not have been realised without the numerous consultation of scholarly
works and collaboration from different people. I was always thought not to say «thank you» but
to know how to show gratefulness to those who gave you support. I would therefore show my
thankfulness to all those who supported me throughout this work.

First of all, I thank the Almighty God for giving me the strength to finish my three-year
thesis despite the surgeries I went through.

I would like to thank Jean-Luc Gouzé for accepting to direct my thesis. Thank you for your
help over the years.

I would like to express my particular gratitude to supervisors Suzanne Touzeau and Frédéric Grognard for years of incredible collaboration with them in rigorous scientific research. Moreover, I want to express my gratitude to them for acting as my parents over the years by assisting me during my surgeries and driving me to and from the research facilities and to my hospital appointments. Again, I thank them for the wonderful environment they helped to develop during my journey at Inria in Sophia Antipolis.

I am grateful to Samuel Bowong for his advices and support during all my academic and
research period. Thank you for believing in me and encouraging me in this way. On a daily basis,
I much value your advises.

I would want to express my gratitude to Ludovic Mailleret and Fabien Hallket for agreeing to
lead my Individual Survival Committee. Your advises was really helpful in completing my thesis.

I would like to thank the jury's members composed of Julien Arino as President, Jean Jules Tewa, Arnaud Ducrot the Rapporteur and Beatrice and Fabien Hallket as members of this jury, for their high contribution in the improvement of the quality of this manuscript.

I would like to thank all the members of my team (Biocore), as well as all those who exchanged a smile, a sweet word, in the hallway or at lunch. A special thank goes to the Team Leader Olivier Bernard.

Johan Monagnat deserves a special thank you for being there all these years and welcoming me into his house as his own daughter. I am grateful to you and Suzanne for welcoming me back after these operations. The completion of my thesis would not have been feasible without your help throughout the years.

I appreciate and will never forget the people I made at Inria, Inrae and others place who kept the three years pleasant and lighthearted : Bruno, Yan, Ali, Ignacio, Juan-Carlos, Lucie, Alesia, Carlos, Nicolas, Elena, Lamberto, Agústin, Walid, Israël, Nicolás, Yan, Ali, Ignacio, Juan-Carlos, Adel, Odile, Joseph, Fransheca, Jineth, David, Ludovic, Josué, Pauline, Méline, Hannah, Adel, David, Maxime, Marine, Beryl, Tehreem . Thank you for our time spending together at the volleyball beach and coffee break.

I would also like to warmly thank Marielle for her friendship, advice, presence at my several
hospital appointments, and her role as Stella's aunty during this time.

A special thank you to Masimba for his friendship and advice on how to improve my English
during these three years.

A special thank you to Yves for being there from the beginning of this experience with his
assistance concerning my health issues and his scientific advice during my internship.

A special thanks to my mother Thérèse Ngoulla for all the sacrifices and efforts made during all these years to encourage her children to go as far as possible in their achievement. Your character as a hard-working woman helped me to better understand life and its difficulties.
I would like to thank Stella, my lovely daughter for the strength and the support she gives me everyday by her presence here in my life.

I am indebted to my brothers and sisters Cyrille, Delphine, Gilbertine, Aubin, Rodrigue, Carole, Chanceline who are dearest in my heart. Thanks for their support and permanent assistance. It is a blessing to have a family that supports you in what you do and values what you do. My brothers and sisters, this journey today is made possible in part by you, since your adventures motivate me every day, and you are my first role models of success.

A special thanks to Yvon for being there, for encouraging me, and for reminding me when things get tough that I can accomplish anything.

Special thanks go to Fanou for his advices, support and encouragements in all my undertakings.

A particular thanks to Marc, my physiotherapist; the several operations I had throughout my thesis years required a patient and attentive individual to enhance my health, and you were such person.

A particular thank you to my best friend Gaïtan, who I affectionately call Gaï; you are one of those people who has been there since the beginning of this study experience, since Babété, and who is always there to encourage me and provide a helping hand when required.

I cannot my friends Daniel, Rodin, Naomie, Foutse, Charles, Vanessa, Irène, Auriole, Ralph, Mitterand, Anita, Myriam, Damarique, Crepin, Roméo, Guerlaine, Bienvenu, Cedric, Camelle, Arantes, Pablo, Jerry, Léon, Harold, Franklin, Arnold, Alphonse, Pierre who also stand by me in one way and the other during these years.

To everyone whose is not mentioned here, your prayers, encouragements, patience and tolerance will always be remembered.
Modélisation mathématique et contrôle de
phytopathogènes de plantes pérennes

**Résumé** Les maladies fongiques causent d'importants dégâts dans les cultures et menacent la sécurité alimentaire mondiale. Les fongicides chimiques sont couramment utilisés, mais ils sont nocifs pour l'environnement, la santé humaine et peuvent entraîner des résistances des pathogènes. Il est donc nécessaire de développer des méthodes de contrôle plus durables. Pour aborder cette question, nous avons développé des modèles mathématiques représentant l'interaction plante-champignon, afin de proposer des méthodes alternatives. Nous avons obtenu des résultats qualitatifs génériques, qui ont été appliqués à la rouille orangée du caféier causée par *Hemileia vastatrix*, une importante maladie du café. Cette culture de rente faisant vivre de nombreux petits producteurs et leurs familles, le contrôle de cette maladie est un enjeu socio-économique majeur.

Tout d'abord, nous avons développé un modèle représentant l'interaction plante-champignon en utilisant des équations différentielles ordinaires. Comme dans les modèles épidémiologiques classiques, la population hôte a été subdivisée par état de santé, l'individu étant la feuille de la plante. L'infection étant transmise par les spores libérées par le champignon, leur dynamique a également été incluse. De plus, les feuilles jeunes et matures ont été différenciées, pour tenir compte de variations d'agressivité du champignon selon l'âge de l'hôte. Nous avons calculé le taux de reproduction de base $R_0$, qui détermine classiquement la stabilité de l'équilibre sans maladie. Ce modèle a exhibé des propriétés asymptotiques complexes, différentes des modèles classiques : $R_0 < 1$ ne suffisait pas à obtenir la stabilité de l'équilibre sans maladie, un équilibre endémique stable pouvant exister ; $R_0 > 1$ ne garantissait pas l'existence et la stabilité d'un équilibre endémique, car de fortes valeurs de $R_0$ pouvaient conduire à la destruction de la plantation.

Puis, les spores de champignons tels que *H. vastatrix* étant dispersées par le vent, nous avons considéré un modèle spatio-temporel décrivant la propagation de la maladie dans une plantation pendant la saison des pluies et la saison sèche, à base d'équations aux derivées partielles. Nous avons calculé deux seuils, les taux de reproduction de base en saison des pluies et saison sèche, qui caractérisent la stabilité des équilibres des sous-systèmes saisonniers. Pour illustrer ces résultats théoriques, des simulations numériques ont été réalisées, en utilisant une méthode non-standard de différences finies pour intégrer le modèle. Nous avons également étudié numériquement l'impact d'un agent de biocontrôle qui réduit la reproduction du champignon. Nous avons déterminé son seuil d'efficacité afin d'assurer l'éradication de la maladie.

Enfin, nous avons développé un modèle multi-saison, avec des dynamiques continues pour l'hôte et le champignon en saisons des pluies et des événements discrets pour les dynamiques plus simples des saisons sèches. En outre, nous avons implémenté une stratégie de biocontrôle fondée sur des hyperparasites se nourrissant de spores, dont la dynamique a été explicitement représentée. Les hyperparasites étaient lâchés une ou plusieurs fois pendant les saisons des pluies. Des études analytiques et semi-numériques ont été réalisées afin de déterminer la quantité et la fréquence des lâchers nécessaires pour contrôler efficacement la maladie. Nous avons montré que la meilleure stratégie dépendait de la mortalité des hyperparasites : les parasites à faible mortalité ne devaient être introduits qu'une fois par an, tandis que les parasites avec une mortalité élevée devaient l'être plus fréquemment afin d'assurer leur persistance dans la plantation.

Grâce à la modélisation mathématique, ce travail fournit des bases qualitatives et quantitatives pour la compréhension des interactions plante-champignon, ainsi que pour la mise en place d'alternatives aux fongicides chimiques afin de lutter contre les maladies fongiques de cultures pérennes.

**Mots-clés:** modélisation mathématique, épidémiologie végétale, culture pérenne, champignon phytopathogène, rouille orangée du caféier, méthodes de lutte durables
Mathematical modelling and control of
perennial plant phytopathogens

**Abstract** Fungal diseases cause serious damages in crops, which threatens food security world-
wide. In response, chemical fungicides are commonly used, yet they are harmful to the envi-
ronment and human health and they may induce fungus resistance. There is a definite need
for more sustainable control methods. To tackle this issue, our approach consisted in developing
mathematical models representing crop-fungus interactions, in order to design alternative control
methods. We obtained generic qualitative results, which we applied to coffee leaf rust caused by
*Hemileia vastatrix*, a major coffee disease. Coffee being an important cash crop which provides a
living for numerous small producers and their family, controlling coffee leaf rust is hence a major
socio-economic issue.

First, we developed a crop-fungus interaction model using ordinary differential equations. As in classical epidemiological models, the host population was subdivided by health status, considering the crop leaf as an individual. The infection being mediated by fungus spores released in the plantation, their dynamics were also included. Moreover, this model differentiated between young and mature leaves, to take into account variations in fungus aggressiveness according to host development. We computed the basic reproduction number $R_0$, which classically determines the stability of the disease free equilibrium. This model exhibited complex asymptotic properties, that differed from classical epidemiological models: $R_0 < 1$ was not sufficient to obtain the stability of the disease free equilibrium, as a stable endemic equilibrium could exist; $R_0 > 1$ did not guarantee the existence and stability of an endemic equilibrium, because high $R_0$ values could lead to the destruction of the plantation.

Then, based on the knowledge that spores of fungi such as *H. vastatrix* are dispersed by wind, we considered a spatio-temporal model describing the disease propagation in a plantation during the rainy and dry seasons, using partial differential equations. We computed two threshold parameters, the rainy and dry period basic reproduction numbers, that characterised the stability of the equilibria for seasonal subsystems. To illustrate these theoretical results, numerical simulations were performed, using a non-standard finite method to integrate the pest model. We also numerically investigated the impact of a biocontrol agent which reduces the reproduction of the fungus. We determined its efficiency threshold in order to ensure disease eradication.

Finally, we developed a multi-seasonal model, alternating continuous crop-fungus dynamics
during the rainy seasons and discrete events to represent the simpler dynamics during the dry
seasons. Moreover, we implemented a biocontrol strategy based on hyperparasites, whose dy-
namics were explicitly represented. Hyperparasites were introduced through one or more discrete
events over the rainy seasons. Analytical and semi-numerical studies were performed to deter-
mine how much and how frequently hyperparasites had to be introduced to efficiently control
the disease. We showed that the best strategy depended on the hyperparasite mortality: low
mortality parasites needed be released only once a year, while high mortality parasites had to
be released more frequently to ensure their persistence in the plantation.

Through mathematical modelling, this work hence provides qualitative and quantitative bases
for the understanding of interactions between a crop and a phytopathogenic fungus, as well as
for the implementation of alternatives to chemical fungicides for the control of fungal diseases in
perennial crops.

**Keywords:** mathematical modelling, plant epidemiology, perennial crop, phytopathogenic fungus, coffee leaf rust, sustainable control methods
# Modélisation mathématique et contrôle de phytopathogènes de plantes pérennes

**Résumé grand public** Les maladies fongiques causent d'importants dégâts dans les cultures et menacent ainsi la sécurité alimentaire mondiale. Les fongicides sont largement utilisés, malgré leur impact négatif sur l'environnement et la santé humaine. En outre, ils induisent des résistances qui limitent leur utilisation à long terme. Des méthodes de lutte alternatives et durables doivent donc être développées. Cette thèse contribue à cet enjeu par une approche de modélisation mathématique. Les modèles représentent les dynamiques saisonnières des interactions plantechampignon et intègrent une méthode de lutte biologique par des lâchers d'hyperparasites. Les résultats obtenus ont permis de déterminer la quantité et la fréquence de lâchers nécessaires pour contrôler la maladie. Ils ont été appliqués à la rouille orangée du caféier, une importante maladie du café. Cette culture de rente faisant vivre de nombreux petits producteurs, le contrôle de cette maladie est un enjeu socio-économique majeur.

# Mathematical modelling and control of perennial plant phytopathogens

**Lay summary** Fungal diseases cause significant damages in crops and thus threaten global food security. Fungicides are widely used, despite their negative impact on the environment and human health. In addition, they induce resistance that limits their long-term use. Alternative and sustainable control methods therefore need to be developed. This thesis contributes to this issue through a mathematical modelling approach. The models represent the seasonal dynamics of plant-fungus interactions and integrate a biological control method using hyperparasite releases. The results obtained allowed to determine the quantity and frequency of releases necessary to control the disease. These results were applied to coffee leaf rust, a major coffee disease. Coffee being an important cash crop which provides a living for numerous small producers, controlling coffee leaf rust is a major socio-economic issue.
# Contents

<table><tr><td>Dedication</td><td>i</td></tr><tr><td>Acknowledgements</td><td>iii</td></tr><tr><td>Résumé</td><td>v</td></tr><tr><td>Abstract</td><td>vi</td></tr><tr><td>Abstract et résumé grand public</td><td>vii</td></tr><tr><td>Contents</td><td>xi</td></tr><tr><td>List of figures</td><td>xiii</td></tr><tr><td>List of tables</td><td>xv</td></tr><tr><td>1 Introduction</td><td>1</td></tr><tr><td>1.1 Context</td><td>1</td></tr><tr><td>1.1.1 Socio-economic impact of crop pests and diseases</td><td>1</td></tr><tr><td>1.1.2 Coffee pests and diseases</td><td>2</td></tr><tr><td>1.2 Biology of coffee leaf rust (CLR)</td><td>4</td></tr><tr><td>1.2.1 Coffee development</td><td>4</td></tr><tr><td>1.2.2 Fungus Hemileia vastatrix</td><td>5</td></tr><tr><td>1.2.3 CLR control methods</td><td>9</td></tr><tr><td>1.3 Mathematical models for fungal disease propagation</td><td>11</td></tr><tr><td>1.3.1 Models for fungal crop diseases</td><td>11</td></tr><tr><td>1.3.2 CLR models</td><td>16</td></tr><tr><td>1.4 Objectives</td><td>18</td></tr><tr><td>1.5 Thesis structure</td><td>18</td></tr><tr><td>2 Mathematical preliminaries</td><td>21</td></tr><tr><td>2.1 Stability of epidemiological models</td><td>21</td></tr><tr><td>2.1.1 Computation of the basic reproduction number</td><td>21</td></tr><tr><td>2.1.2 Bifurcation analysis</td><td>23</td></tr><tr><td>2.2 Stability of semi-discrete models</td><td>23</td></tr><tr><td>2.2.1 Floquet theory</td><td>24</td></tr><tr><td>2.2.2 Jury conditions</td><td>25</td></tr><tr><td>2.2.3 Semi-numerical analysis</td><td>26</td></tr><tr><td>2.3 Analysis and simulation of parabolic PDE</td><td>26</td></tr><tr><td>2.3.1 Existence of solution</td><td>26</td></tr><tr><td>2.3.2 Nonstandard finite difference schemes</td><td>28</td></tr></table>
<table><tbody><tr><td><b>3</b></td><td><b>Bifurcation analysis in an epidemiological model</b></td><td><b>31</b></td></tr><tr><td>3.1</td><td>Introduction</td><td>33</td></tr><tr><td>3.2</td><td>Model formulation and basic properties</td><td>34</td></tr><tr><td>3.2.1</td><td>Formulation of the CLR model</td><td>34</td></tr><tr><td>3.2.2</td><td>Basic properties</td><td>36</td></tr><tr><td>3.3</td><td>The disease free equilibrium and its stability</td><td>39</td></tr><tr><td>3.3.1</td><td>The basic reproduction number and the local stability of the DFE</td><td>40</td></tr><tr><td>3.3.2</td><td>Global stability of DFE</td><td>40</td></tr><tr><td>3.4</td><td>Endemic equilibria and their stability</td><td>43</td></tr><tr><td>3.4.1</td><td>Bifurcation analysis</td><td>44</td></tr><tr><td>3.4.2</td><td>Forward bifurcation and stability</td><td>47</td></tr><tr><td>3.4.3</td><td>Backward bifurcation and stability</td><td>49</td></tr><tr><td>3.5</td><td>Other types of bifurcation</td><td>56</td></tr><tr><td>3.6</td><td>Conclusion</td><td>56</td></tr><tr><td><b>4</b></td><td><b>Modelling CLR dynamics to control its spread</b></td><td><b>59</b></td></tr><tr><td>4.1</td><td>Introduction</td><td>61</td></tr><tr><td>4.2</td><td>Coffee leaf rust dynamics</td><td>62</td></tr><tr><td>4.2.1</td><td>Biological background</td><td>62</td></tr><tr><td>4.2.2</td><td>Model formulation</td><td>62</td></tr><tr><td>4.3</td><td>Mathematical analysis</td><td>65</td></tr><tr><td>4.3.1</td><td>Basic properties</td><td>65</td></tr><tr><td>4.3.2</td><td>Equilibria and their stability</td><td>70</td></tr><tr><td>4.3.3</td><td>Comparison of the dynamics in the dry and rainy seasons</td><td>72</td></tr><tr><td>4.4</td><td>Numerical simulations</td><td>73</td></tr><tr><td>4.4.1</td><td>CLR dynamics</td><td>74</td></tr><tr><td>4.4.2</td><td>Control of CLR</td><td>75</td></tr><tr><td>4.5</td><td>Conclusion</td><td>80</td></tr><tr><td>4.6</td><td>Appendix: Mathematical analysis</td><td>80</td></tr><tr><td>4.7</td><td>Appendix: Numerical scheme</td><td>86</td></tr><tr><td>4.7.1</td><td>Discretization of model system</td><td>86</td></tr><tr><td>4.7.2</td><td>Discretization of boundary conditions</td><td>88</td></tr><tr><td>4.7.3</td><td>Consistence of scheme</td><td>88</td></tr><tr><td><b>5</b></td><td><b>Impulsive modelling of rust and biocontrol</b></td><td><b>89</b></td></tr><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><tr><td><b>6</b></td><td><b>General conclusion</b></td><td><b>119</b></td></tr><tr><td>6.1</td><td>Summary of main results</td><td>119</td></tr><tr><td>6.2</td><td>Perspectives</td><td>120</td></tr></tbody></table>
<table><tr><td><b>A Ifac paper</b></td><td><b>123</b></td></tr><tr><td>A.1 Optimal control of the impulsive model of CLR</td><td>123</td></tr><tr><td>A.1.1 Introduction</td><td>124</td></tr><tr><td>A.1.2 The model and its basic properties</td><td>125</td></tr><tr><td>A.1.3 Mathematical analysis</td><td>126</td></tr><tr><td>A.1.4 Optimal control</td><td>132</td></tr><tr><td>A.1.5 Conclusion</td><td>136</td></tr><tr><td><b>List of Publications and Talks</b></td><td><b>137</b></td></tr><tr><td><b>Bibliography</b></td><td><b>139</b></td></tr></table>
List of Figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Coffee maturation, from flower to ripe fruit</td>
      <td>5</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>World distribution of coffee leaf rust</td>
      <td>6</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Timeline of coffee leaf rust, from 1865 to 2020</td>
      <td>7</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Life cycle of fungus H. vastatrix</td>
      <td>8</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Diagram of the CLR propagation in the coffee plantation: ODE model</td>
      <td>35</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Bifurcation diagram and dynamics of susceptible young leaves for the forward bifurcation case</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Persistence of CLR in the coffee plantation</td>
      <td>50</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>CLR impact according to host development</td>
      <td>51</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Bifurcation diagram and dynamics of susceptible young leaves for the backward bifurcation case</td>
      <td>53</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Extinction and persistence of CLR in the coffee plantation</td>
      <td>54</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Destruction of the plantation by CLR</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>The (R<sub>0</sub>, J<sub>S</sub>)-bifurcation diagram using parameter values in Table 3.4 and in Table 3.1.</td>
      <td>56</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Diagram of the CLR PDE model</td>
      <td>63</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Spatio-temporal simulation of system (4.1,4.4,4.30) when &gamma;<sub>D</sub> = &gamma;<sub>R</sub> = 1.5, which leads to R<sub>0</sub><sup>(D)</sup> = 0.52 and R<sub>0</sub><sup>(R)</sup> = 0.95.</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Spatio-temporal simulation of system (4.1,4.4,4.30) when &gamma;<sub>D</sub> = &gamma;<sub>R</sub> = 8, which leads to R<sub>0</sub><sup>(R)</sup> = 5.09 and R<sub>0</sub><sup>(D)</sup> = 2.77.</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Temporal simulation, for three values of x, of system (4.1,4.4,4.30) when &gamma;<sub>D</sub> = &gamma;<sub>R</sub> = 8, which leads to R<sub>0</sub><sup>(R)</sup> = 5.09 and R<sub>0</sub><sup>(D)</sup> = 2.77.</td>
      <td>77</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Temporal simulations, for various efficiencies q of the biocontrol mycoparasite, of controlled system (4.31,4.4,4.30).</td>
      <td>79</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Diagram of the CLR multi-seasonal model</td>
      <td>93</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Impact of CLR on the production of coffee berries</td>
      <td>101</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Stability regions of periodic disease free solutions as functions of model parameters</td>
      <td>105</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Impact of biocontrol on the dynamics of CLR</td>
      <td>106</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Impact of the predator mortality and its release frequency on its dynamics without CLR</td>
      <td>111</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Impact of release frequency on the stability for low predator mortality rate &mu;<sub>P</sub> = 0.003/day</td>
      <td>112</td>
    </tr>
    <tr>
      <td>5.7</td>
      <td>Impact of release frequency on the stability for high predator mortality rate &mu;<sub>P</sub> = 0.1/day</td>
      <td>113</td>
    </tr>
    <tr>
      <td>5.8</td>
      <td>Impact of release frequency on the transient dynamics with low predator mortality</td>
      <td>114</td>
    </tr>
    <tr>
      <td>5.9</td>
      <td>Impact of release frequency on the transient dynamics with high predator mortality</td>
      <td>115</td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>Trajectories of system (A.1) when &gamma; = 1.6 (so that R<sub>0</sub> = 0.36).</td>
      <td>130</td>
    </tr>
    <tr>
      <td>A.2</td>
      <td>Trajectories of system (A.1) for different initial condition, when &gamma; = 2.1 (so that R<sub>0</sub> = 5.7 &gt; 1).</td>
      <td>131</td>
    </tr>
    <tr>
      <td>A.3</td>
      <td>Trajectories of HCOP (A.30) when &gamma; = 2.1 (so that R<sub>0</sub> = 5.7 without control).</td>
      <td>135</td>
    </tr>
  </tbody>
</table>
List of Tables

<table><tbody><tr><td>3.1</td><td>Description and values of parameters for system (3.1)</td><td>36</td></tr><tr><td>3.2</td><td>Descartes rule for λ-polynomial (3.19).</td><td>44</td></tr><tr><td>3.3</td><td>Force of infection λ and susceptible young leaves J<sub>S</sub> at equilibrium using parameter values in Table 3.1.</td><td>49</td></tr><tr><td>3.4</td><td>Parameter values of system (3.1) chosen to obtain a backward bifurcation.</td><td>50</td></tr><tr><td>3.5</td><td>Force of infection and susceptible young leaves at equilibrium</td><td>52</td></tr><tr><td>4.1</td><td>Description of parameters for system (4.1) – numerical values are used in section 4.4.</td><td>64</td></tr><tr><td>4.2</td><td>Impact of biocontrol efficiency q on reproduction numbers, defined in (4.32), and on berry production during the 8th and last year.</td><td>79</td></tr><tr><td>5.1</td><td>Description and values of parameters for system (5.22).</td><td>100</td></tr><tr><td>5.2</td><td>Berry production at the end of each season, for different yearly release frequencies and high predator mortality.</td><td>116</td></tr><tr><td>A.1</td><td>Description and values of parameters for system (A.1)</td><td>126</td></tr></tbody></table>
# INTRODUCTION

## Contents

<table><tr><td>1.1 Context</td><td>1</td></tr><tr><td>1.1.1 Socio-economic impact of crop pests and diseases</td><td>1</td></tr><tr><td>1.1.2 Coffee pests and diseases</td><td>2</td></tr><tr><td>1.2 Biology of coffee leaf rust (CLR)</td><td>4</td></tr><tr><td>1.2.1 Coffee development</td><td>4</td></tr><tr><td>1.2.2 Fungus Hemileia vastatrix</td><td>5</td></tr><tr><td>1.2.3 CLR control methods</td><td>9</td></tr><tr><td>1.3 Mathematical models for fungal disease propagation</td><td>11</td></tr><tr><td>1.3.1 Models for fungal crop diseases</td><td>11</td></tr><tr><td>1.3.2 CLR models</td><td>16</td></tr><tr><td>1.4 Objectives</td><td>18</td></tr><tr><td>1.5 Thesis structure</td><td>18</td></tr></table>

## 1.1 Context

### 1.1.1 Socio-economic impact of crop pests and diseases

Agriculture plays a key role in reducing poverty, raising incomes and improving food security [110]. The rapid growth of global population raises the need for agricultural production while simultaneously increasing the need for agricultural land. As agriculture struggles to feed this growing population, climate change, pest and crop diseases reduce quantity and quality of crop production [110, 53, 3]. According to FAO estimates, up to 40% of the world crop production is lost each year due to pests, and plant diseases cost the global economy more than $220 billion yearly [46]. Global yield losses for crops that dominate the world in production, such as rice, wheat, barley, maize, potatoes, soybeans, cotton, and coffee in various regions and countries range from 20% to 40% [117]. Due to the intricate interactions between social and economic issues, the primary problem facing agriculture is ensuring food security while increasing production and reducing environmental cost [127]. In order to boost productivity, it is necessary to manage pests and plant diseases, but due to the resilience of them, these efforts are not only expensive but also incredibly difficult. Consequently, it is necessary to continuously develop the control methods.
Control of crop diseases and pests is a major issue, particularly in African countries such as Cameroon, where agriculture is a very important sector in terms of income and employment. More precisely, Cameroon's tropical climate is very favourable to agriculture. With 21% of agricultural land, this sector represents 30% of the country export revenues and 14% of its GDP [47, 118]. It employs almost 62% of the working population. The main export crops are cocoa (4th largest producing country in 2016, 15% of exports), bananas and plantains (7.2%), cotton (4.2%) and coffee (1.5%). Food crops include bananas and plantains, cassava, yam, maize, sorghum, millet, etc. Amount these crops, coffee is an economically important crop in the tropics and subtropics. Globally, more than 50 countries across Latin America, Africa, Asia and Oceania grow coffee on a commercial scale [126]. Coffee exports are a significant source of income for developing countries, and the worldwide green coffee market is expected to increase from $35.40 billion in 2021 to $47.22 billion in 2028 [51]. However, the production of coffee was attacked by many pests and pathogens.

### 1.1.2 Coffee pests and diseases

Coffee is grown in practically every tropical country [177]. Optimal coffee-growing conditions include cold to warm tropical climates, rich soils, and theses optimal conditions favour pests or diseases. The world coffee belt spans the globe along the equator, with cultivation in North, Central, and South America; the Caribbean; Africa; the Middle East; and Asia. Brazil is now the world largest coffee-producing country, accounting for one-third of all the coffee produced worldwide [171]. The climate change increased intense heat waves, drought, and excessive rainfall in some coffee cultivation areas which make more difficult to grow coffee, and increase the pathogens and pests that affect it. The list provided here does not include all of the diseases and pests that are damaging coffee farms around the world, but it does represent a significant portion of them. We classify two category, coffee pests and diseases [170].

*Coffee pests*

* The coffee berry borer (CBB), whose scientific name is *Hypothenemus Hampei*, is an insect of the *Scolytidae* family that feeds and grows at the expense of the berries. According to historical records, this insect pest was originally identified in Gabon in Central Africa in 1901. But this pest currently presents a risk to every coffee plantation in the world, starting in Africa [74]. Accidental introductions of CBB in coffee plantations have led to its rapid geographical expansion worldwide. This insect can adapt to a variety of environmental circumstances and defy a variety of eradication programmes implemented by coffee growers since it spends the most of its life inside a berry [52, 138].

* Root knot nematodes (*Meloidogyne spp.*) develop inside the roots of plants, their impact reduce the ability of the coffee tree to absorb water and nutrients. This is in turn reflected in a decline in the general health and vigour of the tree and hence yield and bean quality [41, 138].

* Nursery disease, caused by *Rhizoctonia solani*, occurs on seedlings of 1-3 months age in the nursery if the conditions are favourable for the fungus growth [41].

* *Monochamus leuconotus*, commonly named white coffee stem borer, is a pest whose larvae, by feeding on the bark and tunnelling into the wood, cause serious damage to the tree.
This pest is endemic to Africa and has been reported in many coffee-growing countries
[138].

* *Coccus viridis*, named green scale, is an important widespread coffee scale present in most coffee-producing countries. It affects plants during any stage of their growth, whether vegetative, flowering, fruiting, or post-harvest [138].

* *Pseudococcus коffеae*, named coffee root mealybug, is a tiny insect that consumes plants by draining the sap from their roots. Mealybugs can multiply rapidly, especially in dry weather. Root mealybugs have been known to cause severe damage to the crop, especially in Kenya [138].

Coffee diseases

* Coffee berry disease (CBD) is a disease caused by fungus *Colletotrichum kahawae* that produces areas of black necrosis and premature dropping of the green coffee berries. High humidity, relatively warm temperatures, and high altitude are ideal for disease formation [98]. The first report of CBD dates back to 1922 in western Kenya, where it led to the destruction of coffee plantations [103]. It has only been found to date in African production zones, as in the west of Cameroon [17], and it can lead to 60% harvest losses [99].

* Coffee leaf rust (CLR) is a disease caused by an obligate fungus, *Hemileia vastatrix*, which attacks leaves and reduces the photosynthesis rate by either defoliating the plant or covering a portion of the leaf surface area. It has an impact on yield because berry production is correlated with the amount of foliage [87]. It was first reported by an English explorer on wild *Coffea* species in the Lake Victoria region of East Africa in 1861 [171]. Since this initial report, CLR has continued to cause significant economic losses. As of 1990, CLR has become endemic in all major coffee-producing countries [171]. In 2012, ten nations in Latin America and the Caribbean experienced a severe CLR outbreak. The disease spread rapidly, and the consequent crop losses caused the supply to drop significantly. Between 2012 and 2014, CLR damaged over $1 billion and impacted over 2 million people in Latin America [160]. During this period in Honduras, 80,000 hectares of coffee farms were infected [175]. Coffee crops in Guatemala were ruined and a state of emergency was declared in February 2013 [13]. Sanitary emergency was also declared in Peru by the government. In the Hawaiian Islands, coffee leaf rust was reported in January 2021 [66].

* There are four different root diseases that affect coffee plants: brown root disease, red root disease (by *Poria hypolateritia*), black root disease (by *Rosellinia bunodes*) and Santavery root disease (by *Fusarium oxysporum*). These diseases damage the tree roots: they stop the trunk supply of nutrients and water, which kills the tree [41].

* The soil-borne fungus *Fusarium stilbioides*, which causes coffee bark disease, only seems to affect *Coffea arabica*. It was initially identified in Tanzania in 1932, then later in Malawi, where it significantly damaged the economy of the country [138].

* Coffee wilt disease, also known as “fusarium wilt”, is a vascular disease caused by fungus *Fusarium xylarioides*. The fungus blocks the xylem system, causing host reactions that progressively lead to plant death [50]. Coffee wilt disease was first observed in the Central
African Republic in 1928 where it developed slowly and went on to cause two epidemics between the 1930s and the 1960s [138].

* American leaf spot of coffee is a disease caused by gemmiferous fungus *Mycena citricolor* [7].

* The coffee brown-eye spot, a disease caused by the fungus *Cercospora coffeicola*, has been reported in all coffee-growing countries. The disease is usually a problem when coffee plants are not growing well, because of poor nutrition or too little shade. In nurseries, it can cause leaf fall of seedlings and, in severe cases, stem dieback. Generally, it is less important on mature plants, although epidemics may occur on well-maintained trees [41].

* Powdery rust of coffee is a disease caused by fungus *Hemileia coffeicola*. It was first recorded on *Coffea arabica* in Cameroon in 1932. Infected leaves eventually turn yellow and are desiccated [100].

In contrast to other harmful coffee diseases like coffee berry disease, which is more common in Africa, CLR affects all coffee plantations worldwide. Moreover, new epidemics still happen nowadays, as was the case in Hawaii in 2021 [66]. A major epidemic like the one that struck America in 2012 could occur at any time and bring about the collapse of the coffee industry, in addition to the losses that the disease continues to cause today. Estimates of yield loss due to CLR vary by country and might be as high as 80% [171]. The costs of controlling the disease are also high. They are largely due to the reliance on fungicides, which are estimated, globally, between US$1 billion and US$3 billion per year [138]. This is why we chose to concentrate on CLR.

## 1.2 Biology of coffee leaf rust (CLR)

### 1.2.1 Coffee development

The coffee tree is a tropical perennial plant of the *Rubiaceae* genus. There are around 100 species of this plant. The two main species of coffee tree cultivated worldwide are *Coffea arabica* (arabica) and *Coffea canephora* (robusta) [30]. Less cultivated species include *Coffea liberica* and *Coffea excelsa*, which are mainly restricted to West Africa and Asia, and account for only 1-2% of the global production [177]. The taxonomic classification has become increasingly complex due to the discovery of several new species during the twentieth century in West Africa, Central Africa, Madagascar, and East Africa [30]. The growth of a coffee plant from seed to first flowering and fruit production takes around five years. The fruit of the coffee tree is called a cherry or berry. The beans that grow within the berry are used to make roast and grind coffee, soluble coffee powders, coffee liquor, etc. The economic lifespan of a coffee plantation is rarely more than 30 years, although a well-managed coffee tree can be productive for 80 years or more, but fruit production decreases over time [177]. There are two main periods in the life cycle of a coffee tree, described below.

* Growth period: from the germination of the seed until adulthood. The coffee bean may be seeded fast and can germinate right away after harvest. After germination it takes around three to seven years to become an adult plant [28]. The coffee tree often has an upright
<figure><img src="image_3.png" /><figcaption>Figure 1.1: Coffee maturation, from flower to ripe fruit - from [159].</figcaption></figure>

main shoot (trunk) with lateral branches that are classified as primary, secondary, and tertiary. On the sides of the main stem and branches of the coffee tree, leaves develop on petioles in opposing pairs [31].

* Productivity period: from the time the coffee tree initially flowers until it dies. The flowers begin to form in the axils of the leaves at the age of three. There are 15 to 50 blossoms at each node, which consists of two leaves [28, 177]. After pollination, 65% of flowers will produce fruit. However, this proportion varies among species and trees, ranging from 3 to 91% [31, 177]. Figure 1.1 shows the fruiting period, which ranges from flowering to fruit maturity, for a single season of coffee production.

Usually, the number of fruits is higher for arabica than robusta, due to the fact that the former is self-fertile, while the latter is self-sterile [28]. Varieties of coffee differ in the size and shape of the coffee bean (seed) but, on average, beans are approximately 10 mm long and 6 mm wide. The weight of a seed is 0.45–0.50 g for arabica and 0.37–0.40 g for robusta [177].

### 1.2.2 Fungus *Hemileia vastatrix*

Coffee leaf rust (CLR) is caused by the fungus known scientifically as *Hemileia vastatrix*. This fungus is an obligate parasite of coffee, meaning that it can only complete its life cycle on plants of the *Coffea* genus.

#### 1.2.2.1 History of *Hemileia vastatrix*

After the first report in 1861, CLR was detected in 1867 on the Island of Ceylon (now Sri Lanka), where it decimated all coffee trees [5]. In Figure 1.2, which depicts the progression of CLR worldwide since its first report until 2014, we can see that practically all of the nations that produce coffee have been affected.

Cameroon went through an epidemic from 1930 to 1935, Brazil in 1970, Costa Rica in 1989, Nicaragua in 1995, El Salvador from 2002 to 2003, Colombia from 2008 to 2013, Peru from 2008 to 2011, and Ecuador in 2013. Coffee production declined in many nations in 2012 as a result of the serious CLR outbreak in Central America. El Salvador, the nation most affected, reported a 70% decrease in coffee yields, and Honduras, whose production had reached a record level in 2012, saw a 23% decrease between 2012 and 2013. In Guatemala, production dropped by 18% between 2012 and 2014, while Nicaragua, the country least affected, met an 11% decline between
<figure><img src="image_4.png" /><figcaption>Figure 1.2: World distribution of coffee leaf rust. Source: Bionovelus, 2014.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 1.3: Timeline of coffee leaf rust, from 1865 to 2020 – from [176].</figcaption></figure>

2012 and 2014 [8]. The first CLR reports in Hawaii appeared in 2020 [66, 4]. Figure 1.3 presents the timeline of CLR since its discovery.

1.2.2.2 Life cycle of *Hemileia vastatrix*

The life cycle of fungus *H. vastatrix* involves several steps, presented in Figure 1.4.

**Step 1:** Uredospore dispersal is generally caused by wind, rain, human impact, or a combination of these factors [153].

**Step 2:** Under favourable conditions uredospores germinate. Although temperature and moisture are key factors for infection, dispersal, and colonisation, plant resistance is also important in determining whether *H. vastatrix* will penetrate the leaves [153].

**Step 3:** The germinal tubes of the uredospores penetrate the open stomata, which are present only on the underside of the leaves. In order to create a biotrophic connection, the fungus colonises the leaf through an intercellular mycelium<fnref n="1" /> growth hypha<fnref n="2" /> that forms haustoria<fnref n="3" /> within the host cells [153].

**Step 4:** Within 24–48 hours, infection is completed. A successful infection will cause yellow

<footnote n="1">¹The mycelium is a root-like structure of a fungus</footnote>

<footnote n="2">²The hypha is a long, branching, filamentous structure of a fungus</footnote>

<footnote n="3">³The haustoria is a rootlike structure that grows into another structure to absorb water or nutrients</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 1.4: Life cycle of fungus *H. vastatrix*: (1) Dissemination, (2) germination, (3) intercellular mycelium, (4) yellow spots and (5) sporulation – adapted from [77].</figcaption></figure>
lesions to appear on the leaf within 2-3 weeks, and these lesions will colonise and produce spores.
Each lesion is made up of multiple pustules that will continue to develop for several months [87].

**Step 5:** Fungus *H. vastatrix* produces uredospores asexually and teliospores sexually. However, because teliospores do not infect coffee leaves and have no known hosts, there is no evidence for an active sexual reproductive cycle [136].

1.2.2.3 Impact of coffee leaf rust

Coffee leaf rust has direct and indirect economic impacts on coffee production. Direct impacts include decreased quantity and quality of yield produced by the diseased plants. Indeed, CLR cause:

* premature defoliation of leaves, which reduces the plant ability to derive energy through photosynthesis [87];

* abortion of flowers and fruits at the premature stage, which decreases their quantity and quality, and which ultimately reduces the beverage quality [123];

* desiccation of young shoots and branches, which reduces the production for the following years [8].

Indirect impacts include increased costs to control the disease [102].

**1.2.3 CLR control methods**

Farmers use integrated pest management to reduce the spread of rust. Integrated pest management is a combination of biological, biotechnological, chemical, physical, cultural methods, as well as the use of resistant cultivars [145]. CLR control methods can be classified into five categories.

**Quarantine** During more than a century, strict quarantine measures prevented coffee rust from invading the Americas. In Central America, new infections were removed by spraying diseased coffee plants as well as asymptomatic plants within 30 metres with a herbicide mixed with diesel fuel [5]. The downside is that coffee plant take several years to produce. Thus, removing infected plants implies to wait at least three years before the replacement plants start producing coffee.

**Chemical control** Among the chemical fungicides available, the ones usually used against CLR are copper fungicides and dithiocarbamate, which must be present on the leaves before the infection occurs [25]. One disadvantage of using fungicides, beside the cost, is that the copper accumulates in the soil, particularly in the organic matter, and can reach levels toxic to plants and to other organisms in the environment [5]. Furthermore, the application of fungicides leads to resistance of fungi over time, which considerably limits their effectiveness and disqualifies them as sustainable control methods [93, 109].

**Cultural practices** Although cultural practices are crucial, there are no clear-cut guidelines to follow. Nutrition is an important aspect in CLR management, since the susceptibility to rust is associated with the nutritional status of the plant. Fertilisation with nitrogen and phosphorus tends to reduce susceptibility to rust, but excessive potassium increases susceptibility [5]. One of
the key cultural management decisions is to produce coffee either in full sun or with some degree of shade [90]. This decision is often more sociopolitical than agronomic. Some say that rust is easier to control on properly spaced plants in full sun, since they dry faster and, therefore, have shorter periods of leaf wetness. Others argue that coffee in the shade has less rust because the closed canopy prevents the formation of dew on coffee leaves and thus reduces the infection. Sun may hinder the development of the pathogen, but it increases the plant productivity. The plant then needs an extra input of nutrients, particularly nitrogen. Furthermore, a higher production may provide more scope to the spread of CLR [177].

**Biological control** It is an approach that employs natural pest enemies. Several biological agents have been tested to control CLR, among which various bacteria. Antagonist bacteria *Bacillus subtilis* isolated from the rhizosphere of coffee crops significantly decreased uredospore development (up to 68%) in vitro [34]. In greenhouse and field conditions, antagonist bacteria *Pseudomonas putida* (P286) and *Bacillus thuringiensis* (B157) completely inhibited uredospore germination; electron microscopy revealed that the B157 isolate prevented the production of germination tubes. An investigation done in Ervália in 2007 revealed that isolate B157 decreased rust intensity as efficiently as copper hydroxide. As a result, it is being investigated as a possible biocontrol agent for coffee rust in organic agricultural systems in Brazil [62]. Forty endophytic bacteria isolates were tested for their capacity to inhibit *H. vastatrix* uredospore germination, among which twenty three reduced more than 40% of the germination; in general, the endophytes were more effective when applied before the inoculation of *H. vastatrix* uredospores [146].

Other investigations on biocontrol focused on hyperparasites, which are parasites whose host, often an insect, is also a parasite. Several potential hyperparasites emerged to be natural enemies of fungus *H. vastatrix*, among which *Verticillium psalliotae*, shown to reduce the germination of uredospores [95, 91], and fungus *Cladosporium hemileiae*, which inhibits the evolution of *H. vastatrix* mycelium [149]. Fungus *Lecanicillium lecanii* (previously called *Verticillium lecanii*) affects the viability of *H. vastatrix* uredospores and was reported as a promising biocontrol agent in several studies [168, 73, 44]. In southern Ethiopia, investigations conducted during the dry and wet seasons for three consecutive years at 60 showed that *L. lecanii* inhibits the rust growth when the seasons change from wet to dry [181]. Finally, the larvae of some *Mycodiplosis* insect species can feed on spores of rust fungus *H. vastatrix* and can therefore be considered as rust spore predators [139, 63, 67].

**Resistant cultivars** Recent biotechnology advances made it possible to detect and choose the best cultivars for coffee farming that should be resistant to fungus *H. vastatrix* [145] and to other fungi [141]. Combining CLR resistance with desirable agronomic traits and high calibre coffee is a problem for breeders. The next issue is to use these resistance genes such that new strains of *H. vastatrix* do not quickly overcome them. The Major resistant cultivars were developed from hybridisation between *C. arabica* and *C. canephora* or their descendants [177]. There are several types of resistant coffee, with varying degrees of CLR resistance (high, moderate, or low) [142]. Some resistant cultivars are cited below.

* Hibrido de Timor is the result of natural hybridisation between *C. arabica* and *C. canephora* [135] and grown widely on the island of Timor. It is the most important progenitor for resistance to coffee leaf rust.
* S288, S333, S795 and BA selections of the Balehonnur Coffee Research Station in India are resistant to many rust races and are believed to have originated from a natural hybrid between *C. arabica* and *C. canephora* [31].

* Icatú is a backcross<fnref n="4" /> between *C. canephora* and Bourbon Vermelho. Its has a resistance genes to many rust races of *H. vastatrix* [43].

* Sarchimor is a cross between the Timor Hybrid and Costa Rican Villa Sarchi that is commonly cultivated in both Costa Rica and India [142].

* Catimor is a cross between the Caturra and the Timor Hybrid. The ideal combination of high yields, excellent disease resistance, and tiny plant size was discovered by scientists in Portugal in 1959. In Brazil, Catimor was initially introduced in 1970 [31].

It is important to remember that the disease has a probability of getting past even these defences. For long-lasting resistance, it is advised to combine several resistance genes in a single cultivar. There are signs that *H. vastatrix* is sometimes unable to combine certain virulence components; when the rust is still able to overcome resistance based on multiple genes, the resulting virulence is anticipated to be modest [43, 31].

## 1.3 Mathematical models for fungal disease propagation

### 1.3.1 Models for fungal crop diseases

Due to the increased spread of fungal diseases, researchers have developed several mathematical models to understand the dynamics of these pathogens, but also to propose control strategies against these diseases. Plant pathogens can be generally divided into two main groups: the necrotrophs which kill their host and feed on their content, and the biotrophs which establish a long term feeding relationship with the living cells of the host [111]. During the infection stage, the penetrating necrotrophic fungi release a group of fungal enzymes which disrupt the cell integrity of the host tissue, causing the cell death which thereby constitutes the food supply. In contrast biotrophic fungi form specialised physiological and morphological adaptations to the living host to ensure the supply of nutrients [150]. Biotrophic fungi also require a living host to complete their life cycle. The existing mathematical models are based on the life cycle of pathogens and their interactions with their host (cf. Section 1.2.2.2). Many of these models are rooted in the epidemic SIR models developed by pioneers Hamer [65], Ross [137], and Kermack and McKendrick [80] as the disease is transmitted by spores produced on infected plants and carried by wind, water, insects, human and impact currents and deposited on new susceptible host plants.

The first models describing the temporal progress of a plant disease epidemic was developed by Van der Plank in the 1960s [163, 164]. Since then, mathematical models of fungal diseases have received a lot of attention from researchers. For instance, Van den Bosch et al. developed a framework for modelling the continental spread of focal plant disease. Therefore, they restricted their first attempt at modelling pandemic spread to simple parameter-sparse submodels. The authors conclude that the simple analytically tractable models, are useful tools to gain qualitative

<footnote n="4">⁴Backcrossing aims at transferring one or two major genes.</footnote>
insight. But they have also been able to conclude that a parameter-rich models is difficult to study due to parameter dependence of the velocity of pandemic spread [161]. 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. He also prove that control at low rust density can delay epidemic development and thus substantially reduce yield losses. However, the model only describes the dynamics of pathogen and it is not possible to estimate the impact on the cereal production [49]. Pivonia et al. studied the seasonal appearance of rusts in the United States, in particular soybean rust. The authors used the mathematical models developed by Gumpert et al. [59] and Van der Plank's [163] to proved the existence of the potential lag period which favours a development of different rusts over a time [124]. Ravigné et al. looked at the impact of sexual and asexual reproduction on the epidemiological dynamics of fungal plant parasites and showed that they could induce cyclic persistence of the disease, which can occur with Sigatoka diseases of banana [130].

Among the models that have focused on airborne spread of these fungal diseases, one can cite the DDAL framework developed by Papaïx et al. that focuses on the deployment of susceptible and resistant crop hosts in an agricultural landscape. The authors proposed a modelling approach for landscape epidemiology that takes advantage of theoretical results developed in the metapopulation context while considering realistic landscape structures [120]. Burie et al. developed two models coupling host and pathogen dynamics for mildew infestation in a vineyard. First, they investigated the structure of travelling waves of the PDE [22]. In a second instance, the authors proved the existence of solutions of their model, the stability of the disease free equilibrium and showed numerically that the disease establishes [23]. Based on this last work, Mammeri et al. studied the impact of spatial heterogeneities on the spread and control of grapevine powdery mildew, but only during a cropping season [97]. Burie et al. proposed a field scale model for the spread of a powdery mildew epidemic over a large vineyard and showed the well-posedness of the asymptotic model and obtained a convergence result confirmed by numerical simulations [24]. Sapoukhina et al. developed a model to simulate the propagation of a fungal disease in a 2D field, including a reaction-diffusion model for short-distance disease dispersal, and a stochastic model for long-distance dispersal. The authors demonstrated that the spatial arrangement of the mixture components in both agro-ecosystems and forestry has to accord with the modes of pathogen dispersal if the efficacy of crop diversification is to be maximized [140]. Mustapha et al. proposed and analysed a mathematical model for describing the dispersal by wind of fungal pathogens in plant populations. The authors modelled the dispersal of pathogen spores using a non-local diffusion equation which took into account variations in wind velocity components [42]

Plant growth and disease spread may be affected by seasonal patterns, which then need to be included in epidemiological models. Among these models, one can cite Rimbaud et al. who 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]. 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]. Maupetit et al. tested which of these factors have a dominant effect on the pathogens development and applied to the particular case of poplar rust [101]. For the crops cultivated in temperate climates or tropical regions that alternate dry and rainy seasons, the pathogens dynamics can substantially differ between the seasons with possibly rapid transitions. Impulsive or semi-discrete models have
been developed for several plant fungal diseases. Among these models, one can also cite Nembot et al. who built 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]. Blériot Tchienkou-Tchiengang et al. build two multi-seasonal mathematical models to show the impact of Maize Stalk Borer (*Busseola fusca*) on maize tissues [157]. Tankam-Chedjou et al. developed 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 economic than continuous removing [105].

Among these models presented in this subsection, we have detailed three cases, which represent the modelling approaches that we have considered in the different chapters of the thesis.

### 1.3.1.1 A cereal rust model

Fleming proposed a simple model that describes the dynamics of the population densities of the cereal rust (in pustules/area) [49]. The model considers low-density rust dynamics in an aerobiological zone, then the population density of the cereal rust N (in pustules/area) is described as

$$ \frac{dN}{dT} = G(N) - L(N) $$

where G(N) and L(N) are the rates per unit area at which the rust population gains and loses pustules, respectively. Time is represented by T.

The rate at which the rust population gains is defined by

$$ G(N) = B \cdot f(N) \cdot T_L + IT_I, $$

where B, T_L and T_I are, respectively, the rate of spore production per pustule, the transmission factors for local and immigrant spores. Author suppose that f(N) = N is the effective number of pustules per unit area. I is the density-independent rate of spore immigration per unit area.

The density-dependent loss rate is defined by

$$ L(N) = f_D(N) + f_E(N) + f_P(N), $$

where f_D(N) = DN is the sum of losses due to natural host deaths and D is the intrinsic rate of natural rust mortality. f_B(N) = EN is the sum of losses due to the essentially passive and wind-mediated process of emigration and E is the instantaneous rate of emigration. The function f_P(N) = AN(Z+N)^{-1} is a sum of losses of the rust fungi to the complex of background predators and parasites, A is the maximum attack rate (pustules/time) and Z is the rust density at which the rate of losses to the complex of background predators and parasites is A/2.
Finally the model is given

$$ \frac{dN}{dT} = (BT_L - D - E)N + IT_I - AN(Z + N)^{-1} \quad (1.1) $$

Fleming used this mathematical model to prove that a complex of polyphagous non-synchronised predators and parasites is likely to control only low-density cereal rust populations. However, the model only describes the dynamics of the pathogen and it is not possible to estimate the impact on the cereal production. This model can be adapted to the dynamics of the uredospores of fungus *H. vastatrix* but we still need to couple it with the dynamics of coffee tree development.

### 1.3.1.2 A model of mildew in a vineyard

Burie et al. proposed a model of the dynamics of mildew in a vineyard [23]. Their modelling hypotheses are that, when spores fall upon the plant tissue, they may create a new colony which will produce spores after some latency period and during some sporulating period. They assumed for simplicity that the time variation of the surface of a colony can be neglected.

The authors proposed a SEIR model in which the total density N of sites than can host a colony of fungus at location x and time t is divided into four categories: healthy H, latent L, sporulating I, and removed (postinfectious) R. The total spore density is separated into two categories according to the range of dispersal: short range $S_1(x, t)$ and longer range $S_2(x, t)$. A sporulating colony produces them at a rate of $r_p > 0$, and they can spread over short distances with a probability of $F \in [0, 1]$ and over longer distances with a probability of $(1 - F)$. The authors assumed that the spores disperse according to a diffusion process with Fickian diffusion coefficient $D_1 > 0$ (short range) or $D_2 > D_1 > 0$ (long range). Spores fall upon the vineyard with some deposition rate $\delta_1 > 0$ or $\delta_2 > 0$ ($\delta_1 = \delta_2$ in numerical simulations). The authors thus obtain the equations that describe the production of spores by the colonies and their dispersal:

$$ \left\{ \begin{array}{l} \frac{\partial S_1}{\partial t}(x, t) = \nabla \cdot (D_1 \nabla S_1(x, t)) - \delta_1 S_1(x, t) + r_p F I(x, t) \\ \frac{\partial S_2}{\partial t}(x, t) = \nabla \cdot (D_2 \nabla S_2(x, t)) - \delta_2 S_2(x, t) + r_p (1 - F) I(x, t) \end{array} \right. $$

for $x \in \Omega$, a regular spatial domain and $t > 0$.

Moreover, they assumed that no spores come from outside the vineyard, so that they imposed Dirichlet conditions on the boundary:

$$ S_1(x, t) = S_2(x, t) = 0 \text{ for } x \in \partial\Omega \text{ and } t > 0 $$

with non-negative initial conditions:

$$ S_1(x, 0) = S_1^0(x) \geq 0, S_2(x, 0) = S_2^0(x) \geq 0 \text{ for } x \in \Omega $$

The set $\Omega_r \subset \Omega$ denotes the area covered by the vine rows. They devised their model in such a way that for all $t > 0$ and $x \in \Omega$, $N(x, t)$ equals 0 if $x \notin \Omega_r$.

According to the authors, the host development is unaffected by the powdery mildew outbreak. New places become accessible for colonisation as a result of this increase. They examined the epidemic over the period of a single season and made the assumption that the logistic rule
governs the variation in time of the total number of colony sites inside the rows.

$$ \frac{\partial N}{\partial t}(x,t) = rN(x,t) \left(1 - \frac{N(x,t)}{K}\right), \quad \text{for } x \in \Omega_r $$

where $r > 0$ is the growth rate and $K > 0$ the carrying capacity.

The local evolution of the disease at some point $x \in \Omega_r$ (inside a row) then obeys the classical SEIR model, whereas we set $N(x,t) = L(x,t) = I(x,t) = R(x,t) = 0$ for $t \ge 0$ if $x \notin \Omega_r$. Let $p$ and $i$ denote the mean duration of the latency and infectious period respectively. Let $E$ be the inoculum effectiveness (probability for the spores succeed in creating a new colony upon a site). The second set of equations of the model for $x \in \Omega_r$ is then

$$ \left\{ \begin{aligned} \frac{\partial H}{\partial t}(x,t) &= -E(\delta_1 S_1(x,t) + \delta_2 S_2(x,t)) \frac{H(x,t)}{N(x,t)} + rN(x,t) \left(1 - \frac{N(x,t)}{K}\right) \\ \frac{\partial L}{\partial t}(x,t) &= +E(\delta_1 S_1(x,t) + \delta_2 S_2(x,t)) \frac{H(x,t)}{N(x,t)} - \frac{1}{p}L(x,t) \\ \frac{\partial I}{\partial t}(x,t) &= +\frac{1}{p}L(x,t) - \frac{1}{i}I(x,t) \\ \frac{\partial R}{\partial t}(x,t) &= +\frac{1}{i}I(x,t) \end{aligned} \right. $$

supplemented with non-negative initial conditions

$$ \begin{aligned} H(x,0) &= H^0(x) \ge 0, L(x,0) = L^0(x) \ge 0 \\ I(x,0) &= I^0(x) \ge 0, R(x,0) = R^0(x) \ge 0 \quad \text{for } x \in \Omega_r \end{aligned} $$

Burie et al. proved the existence of solutions of the full model, computed the basic reproduction number and concluded about the stability of the disease free equilibrium. The model would probably be a good starting point for CLR dynamics modelling, but the authors only showed numerically that the epidemic reach a maximal rate of infection. Their model has to be reduced for the global mathematical analysis because it is difficult, especially with the Fickian diffusion.

### 1.3.1.3 A multi-seasonal model

Madden and Van den Bosch developed a coupled differential equation model for the multi-seasonal dynamics of plant disease for annual crops in the early 2000s [94]. Their model is a semi-discrete SEIR model of a plant disease introduced in an annual cropping system. They looked at the possibility of using plant infections as biological weapons and described the circumstances in which they may endure from season to season. Naming $S$ the disease free plant individuals, $E$ the latent infected individuals, $I$ the infectious individuals, such that $Y = E + I + R$ is the total diseased individual, and naming $P$ the abundance of pathogen inoculum, the model of Madden and Van den Bosch is given by the following equation :

$$ \left\{ \begin{aligned} \frac{dS}{dt} &= -\nu SP - \beta SI \\ \frac{dE}{dt} &= \nu SP + \beta SI - \sigma E \\ \frac{dI}{dt} &= \sigma E - \gamma I \\ \frac{dR}{dt} &= \gamma I \\ \frac{dP}{dt} &= -(\mu_g' + \nu') P \end{aligned} \right. \qquad (1.2) $$

where $P(0) = P_0, S(0) = S_0$. The constant $\nu$ is the primary infection rate parameter, such that $\nu SP$ is the rate of occurrence of new infections, and $\nu S$ is the mean number of new infected plant individuals produced per unit of inoculum per time. Infectious plants are first in latent
state and become infectious at rate $\sigma$. Infectious plants become post-infectious or remove from the epidemic at the rate $\gamma$. New infection occurs from infectious plants at a rate $\beta SI$, with $\beta$ being the secondary infection rate parameter. Thus, $\beta S$ is the mean number of new infected individuals per infected individual per time.

The inoculum decays exponentially at a natural mortality rate $\mu_g P$ and a depletion rate $\nu' P$. In the latter expression, $\mu' = \mu$ because each unit of inoculum that produces a plant infecting is one less unit of inoculum available for infecting another plant.

In this model, an epidemic ends at $t = T_g$, the time of crop harvest. All abundances are set to 0 immediately after $T_g$: $E(T_g^+) = I(T_g^+) = R(T_g^+) = S(T_g^+) = 0$; where the “+” and the “-” superscript represent the instant after the end of the season and the instant before the end of the season, respectively.

Diseased individuals $E, I$ and $R$ produce each an amount of pathogen inoculum at the end of the growing season. If $\theta_E, \theta_I, \theta_R$ represent the amount of pathogen inoculum produced per diseased individuals, then $k_E = \theta_E E(T_g^-)$, $k_I = \theta_I I(T_g^-)$, $k_R = \theta_R R(T_g^-)$ represent the total amount of inoculum at $t = T_g^+$ from latent, infectious and removed diseased individuals at time $t = T_g^-$. Hence, $P(T_g^+) = k_P + k_E + k_I + k_R$; $k_p = P_0 \exp(-(\nu' + \mu_g)T_g^-)$, where $k_p$ is the inoculum from the beginning of the season left at the end of the season.

Between growing seasons, $P$ decays exponentially with the rate $\mu_b(\varepsilon)$. This correspond to $dP/dt = -\mu_b(\varepsilon)P$. As a result, the single-season model (1.2) is extended to account for several seasons, with initial susceptible individuals set to $S_0$ and new inoculum of the growing season provided by the inoculum that survived from the previous season.

Madden and Van den Bosch determined the basic reproduction number and the infection persistence thresholds using this model. This model, that covers multiple seasons, can be adapted to CLR but would not be sufficient as it does not contain a spore compartment.

### 1.3.2 CLR models

Coffee leaf rust epidemics are based on a seemingly simple infection cycle, but develop poly-cyclic epidemics in a season and polyetic epidemics over successive seasons. To study these outbreaks, several authors have published results of their investigations about the germination of the uredospores of *H. vastatrix*, the propagation of CLR in the whole plantation and many other approaches in order to understand the dynamics of the pathogen and to control it.

Modelling coffee rust epidemics necessitates a variety of approaches, thus CLR models in the literature represent a range of scales, from the individual coffee bush to the country or even the continent. Kushalappa developed a method to quantify the rate of CLR development using a linear model. Stepwise regression identified combination of biological and meteorological variables that explained significant variation in the rate of CLR development [86]. Kushalappa et al. used equations for predicting the rate of coffee rust development based on net survival ratio for monocyclic process of *Hemileia vastatrix* [84]. Later, they analysed the development of CLR in the field and using the proportion of infected leaves, pustules/leaf and percentage of infected area to determine the intensity of *H. vastatrix* [85]. Avelino et al. used multivariate analysis to show that CLR development is linked to three sets of factors: the environment, plant growth and development, and grower's practices. The authors used a statistical segmentation technique to define some recommendation domains to control CLR in Honduras, based on these three sets of factors [11]. Avelino et al. also sowed how other trees affect CLR at the field and landscape
scales. The authors identified several canopy characteristics and leaf traits that can help manage
coffee leaf rust at the plot scale: namely, thin canopies with high openness, short base height,
horizontal branching, and small, dentate leaves [9]. Bebber et al. developed a model of germina-
tion and infection risk, and drove this model using estimates of leaf wetness duration and canopy
temperature. They modelled germination and infection as Weibull functions with different tem-
perature optima, based upon existing experimental data. They found that the germination and
infection risks do not depend on the climate in Colombia and neighbouring countries. They
concluded that biological modelling need to improve using spatial resolution and accuracy of
climate reanalyses [15]. Beasley et al. developed a probabilistic model to simulate within-patch
and between-patch transmission using a random walk with spore movement inhibited by canopy
cover. The authors concluded that increasing the spacing between coffee farms and reducing for-
est fragmentation in coffee-growing regions can benefit biodiversity conservation and reduce the
economic impacts of coffee rust [14]. Merle et al. statistically identified the complex combinations
of microclimatic variables responsible for changes in lesion status to construct three models pre-
dicting lesion emergence probability, lesion sporulation probability and growth of its infectious
area [107]. 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. The
model suggested that it could be the larger ecological structure of the agro-ecosystem, both local
and regional aspects, that needs to be considered [169]. Vandermeer et al. also represented the
CLR dynamics in a coffee farm in Chiapas using an SI epidemiological model of the host. They
showed that the network approach can be a useful way of gaining qualitative insight on spatial
disease dynamics [165].

Several researchers were intrigued by the incorporation of integrated crop management meth-
ods into coffee rust models. A survey conducted in Honduras illustrated how crop management,
different combinations of shade, coffee tree density, fertilisation and pruning may strongly influ-
ence coffee rust epidemics [10]. Hannah et al. used Gompertz growth model to describe CLR
epidemics. Their forecast system showed to be promising for recommending sprays for coffee
plantations in Brazil; but authors argue that more field trials are essential to evaluate its perfor-
mance in other coffee regions and under different fungicide management situations [69]. To study
biological control of CLR using hyperpasites, Vandermeer et King used a basic Lotka–Volterra
model, with a consumer depending on two resources, one fixed, another renewable [167]. They
applied the model to the particular case of the fungus *H. vastatrix*, coffee green scales and a my-
coparasite *Lecanicillium lecanii* which predates the former two [166]. Arroyo et al. developed a
stochastic model based on [167] that considers the interactions between bacteria *B. thuringiensis*
B157 and fungus *H. vastatrix*. They studied the effect of the spatial distribution of the disease
and its spread on the overall dynamics of the disease. The authors found that at equilibrium,
the proportion of tree leaves infected with coffee rust is directly proportional to how much the
coffee rust both grows and spreads inside the coffee plantation, and inversely proportional to how
much bacteria population is present in each tree [6]. Zewdie et al. studied CLR dynamics in the
presence of hyperparasites during the dry and wet seasons for three consecutive years at 60 sites
across the southwestern Ethiopia. Their findings highlight the potential of the hyperparasite to
suppress the rust growth rate from the wet season to dry season transition when the rust severity
could otherwise be at its peak. However, they argue that more detailed knowledge is needed on
the interaction of these species to assess its importance for reducing rust induced yield losses or
the risk of rust outbreaks [181].
All these existing CLR models do not explicitly incorporate coffee plant development. In particular, the differences between young and mature leaves are not considered while they clearly have an impact on the severity of infection. In fact, the mature leaves become infected more quickly than the young ones and it is only the latter that produce the coffee berries [122, 106].

## 1.4 Objectives

World agriculture is still affected by serious fungal diseases. Even though several studies have been conducted and a variety of management strategies have been developed, crop losses caused by these diseases are still quite severe. Sometimes the pathogen dynamics are poorly understood and they also become more resistant over time to the chemical fungicides that are still widely used, which makes the disease control more challenging. As a result, the focus of this thesis is to employ mathematical techniques to understand and manage fungal diseases.

The initial goal of the thesis is to develop rigorous mathematical model of plant-pathogen interactions. The models range from ordinary differential equations, which represent the global dynamics of plants and pathogens, to partial differential equations which study aerial dispersal and to hybrid models with discrete jumps, which represent (almost) instantaneous events like harvesting, the introduction of biocontrol agents or seasonality. The qualitative analysis of these models helps to understand their asymptotic behaviour. Additionally, numerical simulations of the models illustrate these theoretical findings and help to better understand the disease transient phases.

The thesis ultimate goal is to develop efficient and long-lasting control strategies that reduce the fungal development and preserve the yield. We favour biological control methods, which are less harmful for the environment and human health than chemical fungicides.

In this thesis, we examine the specific case of coffe leaf rust. Indeed, due to the economic importance of coffee, particularly for developing countries like Cameroon, managing this rust is still a major problem. Latin America and the Caribbean suffered from a major outbreak in 2012 despite the many management measures implemented. The emergency rust summit meeting held in Guatemala in April 2013 compiled a comprehensive list of shortcomings, among which a lack of resources to combat CLR and inefficient fungicide application methods [13]. In this thesis, we will explore alternative control approaches such as the release of hyperparasites. Different hyperparasites have been identified for CLR: for instance fungus *Lecanicillium lecanii* [181] and insect *Mycodiplosis* [67], which feed on CLR spores.

## 1.5 Thesis structure

The thesis is divided into 5 chapters:

* Introduction 1 reviews the impact of fungal diseases and pests in agriculture, and presents different existing fungal disease models. It also describes the biological background on coffee leaf rust, paying special attention to recent developments and CLR modelling.

* Chapter 2 presents some major mathematical tools that were used in the thesis.

* In Chapter 3, we built our baseline model, an ODE model in which the host population is subdivided by health status, considering the crop leaf as an individual. The infection being
mediated by fungus spores released in the plantation, their dynamics were also included.
Moreover, this model differentiated between young and mature leaves, to take into account
variations in fungus aggressiveness according to host development. We then studied an-
alytically the existence and the stability of the equilibria. We also explored numerically
the disappearance or persistence of CLR in the plantation, or the whole destruction of the
coffee plantation by CLR.

* Chapter 4 is based on the ODE model of Chapter 3, which was simplified by combining young and mature leaves into branches. We introduced spore dispersal and we took into account the seasonality which is important factor of coffee production by distinguishing between dry and rainy seasons. The well-posedness of the PDE model with seasonality obtained was addressed. We studied analytically the existence and the stability of the equilibria. We investigated the CLR control by varying a parameter which reduces spore production.

* In Chapter 5, we presented a seasonal model which does not distinguish between young and mature leaves. The dynamics are continuous during the rainy seasons and represented by discrete events during the rainy seasons. For the biological control of the disease, we included the dynamics of hyperparasites which feed on uredospores. We used Floquet theory to analyse the stability of the periodic disease free solution. We used semi-numerical analysis to find the frequency and quantity of predators to release to control the disease.

* Chapter 6 presents the main results and the perspectives resulting from this work.
# MATHEMATICAL PRELIMINARIES

## Contents

<table><tr><td><strong>2.1 Stability of epidemiological models</strong></td><td><strong>21</strong></td></tr><tr><td>2.1.1 Computation of the basic reproduction number</td><td>21</td></tr><tr><td>2.1.2 Bifurcation analysis</td><td>23</td></tr><tr><td><strong>2.2 Stability of semi-discrete models</strong></td><td><strong>23</strong></td></tr><tr><td>2.2.1 Floquet theory</td><td>24</td></tr><tr><td>2.2.2 Jury conditions</td><td>25</td></tr><tr><td>2.2.3 Semi-numerical analysis</td><td>26</td></tr><tr><td><strong>2.3 Analysis and simulation of parabolic PDE</strong></td><td><strong>26</strong></td></tr><tr><td>2.3.1 Existence of solution</td><td>26</td></tr><tr><td>2.3.2 Nonstandard finite difference schemes</td><td>28</td></tr></table>

In this chapter, we present the main mathematical tools that are used in this thesis. As evidenced by the table of content, they are at the crossing point between mathematical epidemiology, dynamical systems and partial differential equations.

## 2.1 Stability of epidemiological models

The transmission of diseases in animal (underlay humans) or plants, such as fungi, can be described using epidemiological models. Once the model is established, we can investigate the behaviour of the disease using mathematical analysis.

### 2.1.1 Computation of the basic reproduction number

The basic reproduction number is computed in classical epidemiological models to determine whether the disease will disappear or establish itself. Here we present a method of Van den Driessche et Watmough [162] to compute that number.

Let us note $x = (x_1, ..., x_n)^t$, with each $x_i \ge 0$, be the number of individuals in each compartment. Let us define the domain

$$X_s = \{x \ge 0 | x_i = 0, i = 1, ..., m\}$$
where the disease is not present.

Let us consider the following epidemiological model with a non-negative initial conditions

$$ \dot{x}_i = f_i(x) = \mathcal{F}_i(x) - \mathcal{V}_i(x), \quad i = 1, \dots, n, \qquad (2.1) $$

where $\mathcal{F}_i(x)$ is the rate of appearance of new infections in compartment $i$, and $\mathcal{V}_i = \mathcal{V}_i^+ - \mathcal{V}_i^-$ with $\mathcal{V}_i^+(x)$ is the rate of transfer of individuals into compartment $i$ by all other means, and $\mathcal{V}_i^-(x)$ is the rate of transfer of individuals out of compartment $i$. It is assumed that each function is continuously differentiable at least twice in each variable.

Let us suppose the following assumptions $(A_1)-(A_5)$

(A_1) if $x \ge 0$, then $\mathcal{F}_i(x), \mathcal{V}_i^+, \mathcal{V}_i^- \ge 0$ for $i = 1, \dots, n$.

(A_2) if $x_i = 0$ then $\mathcal{V}_i^- = 0$. In particular, if $x \in X_s$ then $\mathcal{V}_i^- = 0$ for $i = 1, \dots, m$.

(A_3) $\mathcal{F}_i = 0$ if $i > m$.

(A_4) if $x \in X_s$ then $\mathcal{V}_i^+ = 0$ for $i = 1, \dots, m$.

(A_5) If $\mathcal{F}(x)$ is set to zero, then all eigenvalues of $Df(x_0)$ have negative real parts. where $Df(x_0)$ is the Jacobian matrix $(\frac{\partial f_i}{\partial x_j})$ evaluated at a disease free equilibrium (DFE) $x_0$.

**Lemma 2.1.** [162]

If $x_0$ is a DFE of equation (2.1) and $f_i(x)$ satisfies assumptions $(A_1)-(A_5)$, then the Jacobian matrices $D\mathcal{F}(x_0)$ and $D\mathcal{V}(x_0)$ are partitioned as

$$ D\mathcal{F}(x_0) = \begin{pmatrix} F & 0 \\ 0 & 0 \end{pmatrix} \quad \text{and} \quad D\mathcal{V}(x_0) = \begin{pmatrix} V & 0 \\ J_3 & J_4 \end{pmatrix} $$

where $F$ and $V$ are the $m \times m$ matrices defined by

$$ F = \left( \frac{\partial \mathcal{F}_i}{\partial x_j}(x_0) \right) \quad \text{and} \quad V = \left( \frac{\partial \mathcal{V}_i}{\partial x_j}(x_0) \right), \quad \text{with } 1 \le i, j \le m. $$

Further, $F$ is non-negative, $V$ is a non-singular M-matrix and all eigenvalues of $J_4$ have a positive real part.

From these elements, the basic reproduction number is defined:

**Definition 2.1.** The basic reproduction number, denoted $\mathcal{R}_0$, is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual. $\mathcal{R}_0$ is computed as the spectral radius of matrix $FV^{-1}$:

$$ \mathcal{R}_0 = \rho(FV^{-1}) \qquad (2.2) $$

where $FV^{-1}$ is called the next generation matrix for model (2.1).

We have the following theorem about the local stability of the disease free equilibrium

**Theorem 2.1** (Local stability of a DFE). Consider the disease transmission model given by (2.1) with $f(x)$ satisfying conditions $(A_1)-(A_5)$. If $x_0$ is a DFE of the model, then $x_0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$, and unstable if $\mathcal{R}_0 > 1$.
### 2.1.2 Bifurcation analysis

The Castillo-Chavez and Song theorem [27] allows for the analysis of the different types of bifurcations that the epidemiological model exhibits, giving then the stability of the disease free and endemic equilibria.

**Theorem 2.2.** 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 (2.3) $$

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

A$_1$: $A = D_x f(0, 0) = (\frac{\partial f_i}{\partial x_j}(0, 0))$ is the linearisation matrix of system (2.3) 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 (2.3) around 0 are totally determined by a and b.

1. $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;

2. $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;

3. $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;

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

## 2.2 Stability of semi-discrete models

In the present manuscript, semi-discrete models are developed to model the phenomena related to seasonality or the introduction of hyperparasites to control fungal diseases in crop. For example, plant development and plant-pathogen interactions during the wet season may be modelled using continuous dynamics, whereas simpler dynamics during the dry season allow for their modelling as a discrete event. In this work, we consider semi-discrete models in which impulses occur at fixed times, so that they can be described as follows [12]:
$$ egin{cases} \frac{dx(t)}{dt} = f(t, x), \quad t \neq t_k, \\ \Delta x = I_k(x(t)), \quad t = t_k, \quad \text{for } k \in \mathbb{Z}. \\ x(t_0^+) = x_0, \quad t_0 \ge 0. \end{cases} \qquad (2.4) $$

where $\Delta x = x(t_k^+) - x(t_k)$ with $t_k^+$ denoting the instant just after $t = t_k$ so that $x(t_k^+) = \lim_{h \to 0^+} x(t_k + h)$. $f$ is locally Lipschitz continuous on $\mathbb{R} \times \Omega$, with $\Omega$ an open set of $\mathbb{R}^n$, $I_k : \Omega \to \Omega$.

In the particular case where $t_k = kT$ define regularly spaced discrete events, we define a periodic solution for such a system (2.4) as a map $x : \mathbb{R}^+ \to \mathbb{R}^n$ is said to be a periodic solution of equation (2.4) with period $T$ if it satisfies $x(t+T) = x(t)$ for $t \neq t_k$ and $x(t_k+T^+) = x(t_k^+)$, for $k \in \mathbb{Z}$.

### 2.2.1 Floquet theory

We apply Floquet's theory to study semi-discrete models. Floquet theory, in general, is a branch of the theory of ordinary differential equations that deals with solutions to periodic linear differential equations. These can arise from semi-discrete systems like (2.4) when $f$ is $T$-periodic in $t$ and the $t_{k+q} = t_k + T$ for a fixed $q$. Let us consider the following linear impulsive differential system associated with non-linear model (2.4):

$$ egin{cases} \frac{dx(t)}{dt} = A(t)x(t), \quad t \neq t_k, \\ x(t_k^+) = x(t_k) + B_k x(t_k), \quad t = t_k, \quad t_k < t_{k+1}, \quad k \in \mathbb{Z}, \\ x(t_0^+) = x_0, \quad t_0 \ge 0. \end{cases} \qquad (2.5) $$

where $A(t)$ and $B$ are $n \times n$ matrices. These equations are obtained by linearising system (2.4) around one of its periodic solutions, in order to consider the stability of the periodic solution (the $x$ of (2.5) is then the difference between the $x$ of (2.4) and the periodic solution).

Let us consider the following hypotheses :

$H_1$: $A(.) \in PC(\mathbb{R}, \mathcal{M}_n)$ and $A(t+T) = A(t)$, where $\mathcal{M}_n$ is the set of $(n \times n)$-matrices

$H_2$: $B_k \in \mathcal{M}_n$ and $\det(I + B_k) \neq 0, k \in \mathbb{Z}$

$H_3$: there exists $q \in \mathbb{N}$ such that $B_{k+q} = B_k$ and $t_{k+q} = t_k + T, k \in \mathbb{Z}$.

where $PC(\mathbb{R}, \mathcal{M}_n)$ denotes the set of functions $h : \mathbb{R} \to \mathcal{M}_n$ that are continuous $\forall t \in \mathbb{R}, t \neq t_k$, and admitting jump discontinuities at $t_k$.

**Theorem 2.3** (Fundamental matrix [12]). Suppose that condition $(H_1) - (H_3)$ hold. Then each fundamental matrix of equation (2.5) can be represented in the form $X(t) = \phi(t)e^{\Lambda t}, (t \in \mathbb{R})$ for a non-singular, $T$-periodic matrix $\phi(.) \in PC^1(\mathbb{R}, \mathcal{M}_n(\mathbb{C}))$ and a constant matrix $\Lambda \in \mathcal{M}_n(\mathbb{C})$.

**Definition 2.2** (Monodromy matrix [12]). Let $X(t)$ be a fundamental matrix of equation (2.5). Then the matrix $X(t+T)$ is also fundamental, and there corresponds a unique non-singular matrix $M \in \mathcal{M}_n(\mathbb{C})$ such that

$$ X(t+T) = X(t)M, $$
for all $t \in \mathbb{R}$. Here, M is called monodromy matrix of (2.5) corresponding to the fundamental matrix $X(t)$. All monodromy matrices of equation (2.5) are similar and have the same eigenvalues. Their eigenvalues $\lambda_1, \dots, \lambda_n$ are called Floquet multipliers.

In the order to calculate the Floquet multipliers of (2.5), we have to choose an arbitrary fundamental matrix $X(t)$ and calculate the eigenvalues of the matrix $M = X(t_0+T)X^{-1}(t_0)$ where $t_0 \in \mathbb{R}$ is fixed.

If $X(0) = I$ (or $X(0^+) = I$), then we can choose $M = X(T)$ (or $M = X(T^+)$) as the monodromy matrix. It follows that for any monodromy matrix M of (2.5), we have by [12]

$$ \det M = \lambda_1 \cdot \lambda_2 \dots \lambda_n = \prod_{k=1}^{q} \det (I + B_k) \exp \left( \int_{0}^{T} \mathrm{Tr} A(s) ds \right). $$

Under these hypotheses, the following stability result holds.

**Theorem 2.4** (Stability [12]). Suppose that the conditions $(H_1) - (H_3)$ hold. Then the linear $T$-periodic impulsive equation (2.5) is:

(1) stable if and only if all Floquet multipliers $\lambda_j$, ($j = 1, \dots, n$) of equation (2.5) satisfy $|\lambda_j| \le 1$ and moreover, to those $\lambda_j$ for which $|\lambda_k| = 1$, there corresponds a simple elementary divisor.

(2) asymptotically stable if and only if all Floquet multipliers $\lambda_j$, ($j = 1, \dots, n$) of the equation (2.5) satisfy the inequality $|\lambda_j| < 1$.

(3) unstable if there is a Floquet multiplier $\lambda_j$, ($j = 1, \dots, n$) such that $|\lambda_j| > 1$.

### 2.2.2 Jury conditions

Through Floquet theory and the computation of the monodromy matrix, the stability study often requires the analysis of a discrete-time system, corresponding to the map from time $kT^+$ to time $(k+1)T^+$. To achieve that in dimension 2, we will need a simple tool that circumvent the computation of the eigenvalues, that is the Jury condition.

Let $A$ a $2 \times 2$-monodromy matrix for which we are not able to compute explicitly the Floquet multiplications. The conditions for these multipliers to stay in the unit circle can be established using the Jury criteria.

Defining $A$ as

$$ A = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} $$

the Jury conditions [173] guarantees that the eigenvalues of matrix $A$ stay inside the unit circle if and only if

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

where

$$ \left\{ \begin{array}{l} \mathrm{Tr}(A) = a_{11} + a_{22}, \\ \det(A) = a_{11}a_{22} - a_{12}a_{21}. \end{array} \right. $$
### 2.2.3 Semi-numerical analysis

In order to analyse system (2.5) when it is not possible to compute analytically the Floquet multipliers. For that, we consider the simplified situation where $B_k = B$ and $t_k = kT$. The monodromy matrix associated with system (2.5) is evaluated in two steps: we first compute the state transition matrix of system $\dot{x} = A(t)x$ by numerically solving the linear system (2.5) using the initial conditions $x(0) = e_i$ with $e_i$ the elements of the orthonormal basis. The $n$ solutions evaluated at time $T$, which is a period, are then put together to obtain the fundamental matrix $\Phi(T)$. Then we use the impulsive condition to obtain the monodromy matrix. We name this method "semi-numerical" as it requires numerical computations, but not extensive simulations of the system [56].

The algorithm is given by the following steps:

* Matrix $A(t)$ is a $n \times n$-matrix

* The fundamental matrix verifies $\frac{d\Phi(t)}{dt} = A(t)\Phi(t)$

* We use initial conditions $\Phi(0) = I$ to compute numerically the fundamental matrix $\Phi(T)$

* The monodromy matrix associated with system (2.5) is then

$$x((k+1)T^+) = (I + B)\Phi(T)x(kT^+)$$

* We compute $\mathcal{R}$, the spectral radius of this monodromy matrix

* $\mathcal{R} < 1$ implies stability of the zero solution

The monodromy matrix associated with $x$ is

$$M = (I + B)\Phi(T).$$

## 2.3 Some results for the analysis and simulation of reaction-diffusion partial differential equations

For the positivity of PDE systemen, we use the following Lemma defined by [119] and given as follow:

**Lemma 2.2.** Suppose $w \in C([0, T] \times \bar{\Omega}) \cap C^{1,2}((0, T] \times \Omega)$ that satisfies:

$$\left\{ \begin{array}{ll} w_t - D\Delta w \geq -c(t,x)w(t,x), & x \in \Omega, \quad 0 < t \leq T, \\ \frac{\partial w}{\partial \eta} \geq 0, & \text{on } \partial\Omega, \\ w(0,x) \geq 0, & x \in \Omega, \end{array} \right.$$

where $c(t, x)$ is any bounded function in $[0, T] \times \Omega$. Then, $w(t, x) \geq 0$ on $[0, T] \times \bar{\Omega}$. Moreover, $w(t, x) > 0$ or $w \equiv 0$ in $(0, T] \times \Omega$.

### 2.3.1 Existence of solution

In this section we defined the initial value problem and the existence of solution for the abstract problem following the approach described by Pazy [121].
Let $1 < p < \infty$ and let $\Omega$ be a bounded domain with smooth boundary $\partial\Omega$ in $\mathbb{R}^n$. Consider the initial value problem:

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial u}{\partial t} + A(t, x, D) = f(t, x), & \text{in } \Omega \times [0, T], \\
\\
\displaystyle D^{\alpha}u(t, x) = 0, & \text{|\alpha| < m, on } \partial\Omega \times [0, T], \\
\\
u(0, x) = u_0(x) & \text{in } \Omega
\end{array}
\right.
\qquad (2.7) $$

where

$$ A(t, x, D) = \sum_{|\alpha| \le 2m} a_{\alpha}(t, x) D^{\alpha} \qquad (2.8) $$

Let us define the following assumptions:

(H$_1$) The operator $A(t, x, D)$, $t \ge 0$ are uniformly strongly elliptic in $\Omega$.

(H$_2$) The coefficients $a_\alpha(t, x)$ are smooth functions of the variables $x \in \bar{\Omega}$ for every $0 \le t \le T$ and satisfy for some constants $C_1 > 0$ and $0 < \beta \le 1$

$$ |a_{\alpha}(t, x) - a_{\alpha}(s, x)| \le C_1 |t - s|^{\beta} $$

for $x \in \bar{\Omega}$, $0 \le s, t \le T$ and $|\alpha| \le 2m$

With the family $A(t, x, D)$, $t \in [0, T]$, of strongly elliptic operators, we associate a family of linear operators $A_p(t)$, $t \in [0, T]$, in $\mathcal{L}^p(\Omega)$, $1 < p < \infty$. This is done as follows:

$$ D(A_p(t)) = D = W^{2m,p}(\Omega) \cap W_0^{m,p}(\Omega) $$

and

$$ A_p(t)u = A(t, x, D)u \text{ for } u \in D. $$

If $u_0 \in \mathcal{L}^p(\Omega)$ and $f(t, x) \in \mathcal{L}^p(\Omega)$ for every $0 \le t \le T$ then a classical solution $u$ of the following abstract initial value problem

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{du}{dt} + A_p(t)u = f, \\
\\
u(0) = u_0
\end{array}
\right.
\qquad (2.9) $$

in $\mathcal{L}^p(\Omega)$ is defined to be a generalized solution of the initial value problem (2.7).

**Theorem 2.5** (Existence of solution [121]). Let the family $A(t, x, D)$, $0 \le t \le T$, satisfy the conditions (H$_1$) and (H$_2$) and let $f(t, x) \in \mathcal{L}^p(\Omega)$ for $0 \le t \le T$ satisfy

$$ \left( \int_{\Omega} |f(t, x) - f(s, x)|^p dx \right)^{\frac{1}{p}} \le C |t - s|^{\gamma} \qquad (2.10) $$

for some constants $C > 0$ and $0 \le \gamma < 1$. Then for every $u_0(x) \in \mathcal{L}^p(\Omega)$ the evolution equation (2.7) possesses a unique generalized solution $u \in D$.
### 2.3.2 Nonstandard finite difference scheme for parabolic equations

Here, we present a finite difference scheme for a particular class of reaction-diffusion partial differential equations (PDEs) in dimension one [108]. Consider the following equation

$$ \frac{\partial u_i}{\partial t} = D_i \frac{\partial^2 u_i}{\partial x^2} - f_i(u_1, u_2, \dots, u_n)u_i + g_i(u_1, u_2, \dots, u_n)u_i + \Lambda_i, \quad i = 1, 2, \dots, n, \quad (2.11) $$

where $D_i$ are nonnegative diffusion coefficients, and the $f_i$ and $g_i$ are polynomial or rational functions of $(u_1, u_2, \dots, u_n)$ with positive coefficients. With Neumann boundary condition

$$ \frac{\partial u_i}{\partial x} = 0 \qquad (2.12) $$

We suppose that,

$$ u_i(x, 0) \geq 0 \implies u_i(x, t) \geq 0. $$

Let us introduce the notation:

$$ t \to t_k = \psi(\Delta t)k, \quad x \to x_m = \phi(\Delta x)m, $$

where $k$ is a non-negative integer, $m$ is an integer; $\psi(\Delta t)$ and $\phi(\Delta x)$ are, respectively, the time and space step-sizes, with $\psi$ and $\phi$ satisfying

$$ \begin{cases} \psi(h) = h + O(h^2) \\ \phi(h) = h^2 + O(h^4). \end{cases} \qquad (2.13) $$

The discrete approximations to the dependent variables are

$$ u_i(x, t) = [u_i]_m^k, \quad i = 1, 2, \dots, n. $$

A positivity preserving finite difference scheme is one for which

$$ [u_i]_m^k \geq 0 \Rightarrow [u_i]_{m-1}^{k+1} \geq 0, \quad i = 1, 2, \dots, n. $$

One consequence of this requirement is that it can hold only if a functional relation exists between the time and space step-sizes. In particular, for the PDE system given by equation (2.11), the following relation holds:

$$ \psi(\Delta t) \leq \alpha(\phi(\Delta x))^2, \qquad (2.14) $$

where $\alpha$ is a constant determined by the diffusion coefficients ($D_1, D_2, \dots, D_n$).

The discrete forms for the first-order time derivative and the second-order space derivative are given, respectively, by the usual forward Euler and central difference is given by:

$$ \begin{aligned} \frac{\partial u_i}{\partial t} &= \frac{[u_i]_m^{k+1} - [u_i]_m^k}{\psi(\Delta t)} + O(\psi(\Delta t)) \\ \frac{\partial^2 u_i}{\partial x^2} &= \frac{[u_i]_{m+1}^k - 2[u_i]_m^k + [u_i]_{m-1}^k}{(\phi(\Delta x))^2} + O((\phi(\Delta x))^2). \end{aligned} $$

A non-standard finite difference scheme for the PDE system of equation (2.7) is obtained by making the substitutions of the above expressions. Doing this gives ($i = 1, 2, \dots, n$)
$$ \frac{[u_i]_m^{k+1} - [u_i]_m^k}{\psi(\Delta t)} = D_i \frac{[u_i]_{m+1}^k - 2[u_i]_m^k + [u_i]_{m-1}^k}{(\phi(\Delta x))^2} - f_i([u]_m^k) [u_i]_{m+1}^k + g_i([u]_m^k) [u_i]_{m-1}^k + \Lambda_i, \quad (2.15) $$

where $[u]_m^k = ([u_1]_m^k, [u_2]_m^k, \dots, [u_n]_m^k)$ represent u evaluated at $t = t_k$ and $x = x_m$,
If $R_i$ is defined as

$$ R_i \equiv \frac{D_i \psi(\Delta t)}{(\phi(\Delta x))^2} $$

then equation (2.15) can be written as

$$ \begin{aligned} \left[1 + \psi(\Delta t) f_i([u]_m^k)\right] [u]_m^{k+1} = R_i \left([u]_{m+1}^k + [u]_{m-1}^k\right) + (1 - 2R_i) [u]_m^k + \psi(\Delta t) \Lambda_i \\ + \psi(\Delta t) g_i([u]_m^k) [u]_m^k. \end{aligned} \quad (2.16) $$

Inspection of the second term in equation (2.16) shows that the positivity condition

$$ [u]_m^k \ge 0 \Rightarrow [u]_{m-1}^{k+1} \ge 0 $$

is satisfied, provided that

$$ \frac{D\psi(\Delta t)}{(\phi(\Delta x))^2} \le \frac{1}{2}, \quad D = \max(D_1, D_2, \dots, D_n) \quad (2.17) $$

Then the time-step size should be selected such that the condition given in equation (2.17) holds and the non-standard scheme is (conditionally) stable. However, in general, the solutions of this condition do not satisfy this principle. The imposing of this min-max principle on the solutions to equation (2.17) gives the conditional stability requirement.

Consider Neumann condition in equation (2.12), the discrete form is given by:

$$ \left\{ \begin{aligned} \left. \frac{\partial u_i}{\partial x} \right|_{m=0} &= \frac{[u]_{m+1}^k - [u]_m^k}{(\phi(\Delta x))^2} = 0 &\implies & [u]_1^k &= [u]_0^k \\ \left. \frac{\partial u_i}{\partial x} \right|_{m=N_x} &= \frac{[u]_{m+1}^k - [u]_m^k}{(\phi(\Delta x))^2} = 0 &\implies & [u]_{N_x+1}^k &= [u]_{N_x}^k \end{aligned} \right. $$

where $N_x$ is the number of space point. Replacing the above Neumann condition in equation (2.16), we obtain

$$ \left\{ \begin{aligned} (1 + \psi(\Delta t) f_i([u]_0^k)) [u]_0^{k+1} &= [u]_0^k + \psi(\Delta t) (\Lambda_i + g_i([u]_0^k)) [u]_0^k, \\ (1 + \psi(\Delta t) f_i([u]_{N_x}^k)) [u]_{N_x}^{k+1} &= [u]_{N_x}^k + \psi(\Delta t) (\Lambda_i + g_i([u]_{N_x}^k)) [u]_{N_x}^k. \end{aligned} \right. \quad (2.18) $$

Combining equations (2.16) and (2.18), we obtain a non-standard schemes for a reaction-diffusion partial differential equations with Neumann boundary condition.
BIFURCATION ANALYSIS IN AN
EPIDEMIOLOGICAL MODEL OF FUNGAL
CROP DISEASE

In this chapter, we develop an epidemiological host-pathogen model that describes the propaga-
tion of a fungal disease, applied to CLR, 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 (susceptible leaves $S$, latent leaves $L$, infectious leaves $I$ and dry leaves $R$), considering the crop leaf as an individual. This model differentiates between young and mature leaves, to take into account variations in fungus aggressiveness according to host development. The model includes crop production (berries $B$) by mature leaves. The infection being mediated by fungus spores (uredospores $U$) released in the plantation, their dynamics are also included.

$H_2$: The model takes into account the non-constant recruitment of young susceptible leaves, which are generated by mature leaves. Due to competition, the dynamics of young leaves are negatively impacted by mature leaves.

$H_3$: Although the germination rate (which determines the infection success) differs between young and mature leaves, both face the same force of infection. The latter is proportional to the number of uredospores and inversely proportional to the total number of leaves. Thus, the force of infection is density-dependent.

$H_4$: The model does not take seasonality into account.

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

$H_6$: The model uses an ODE formalism.

$H_7$: No control methods are implemented.

Based on these hypotheses we obtain a complex ODE model with nine compartments. We use it to understand the dynamics of fungal diseases applied to the particular case of CLR. We study the model behaviour for different values of parameters and initial states of the pathogen.
Contents

<table><tr><td>3.1</td><td>Introduction</td><td>33</td></tr><tr><td>3.2</td><td>Model formulation and basic properties</td><td>34</td></tr><tr><td>3.2.1</td><td>Formulation of the CLR model</td><td>34</td></tr><tr><td>3.2.2</td><td>Basic properties</td><td>36</td></tr><tr><td>3.3</td><td>The disease free equilibrium and its stability</td><td>39</td></tr><tr><td>3.3.1</td><td>The basic reproduction number and the local stability of the DFE</td><td>40</td></tr><tr><td>3.3.2</td><td>Global stability of DFE</td><td>40</td></tr><tr><td>3.4</td><td>Endemic equilibria and their stability</td><td>43</td></tr><tr><td>3.4.1</td><td>Bifurcation analysis</td><td>44</td></tr><tr><td>3.4.2</td><td>Forward bifurcation and stability</td><td>47</td></tr><tr><td>3.4.3</td><td>Backward bifurcation and stability</td><td>49</td></tr><tr><td>3.5</td><td>Other types of bifurcation</td><td>56</td></tr><tr><td>3.6</td><td>Conclusion</td><td>56</td></tr></table>

**Abstract:** Most crops are subject to many pests and pathogens, which reduce their yield. In this paper, we develop an SEIR-U epidemiological model for a fungal crop disease, where U stands for spores. The originality of this model is a differentiation between young and mature leaves, that takes into account the impact of host development on the aggressiveness of the pathogen. This is motivated by the particular case of coffee leaf rust (CLR), which can drastically reduce coffee berry production. We compute the basic reproduction number $R_0$, which classically determines the stability of the disease-free equilibrium (DFE), and identify two bifurcation cases at its threshold value 1. The first classical case corresponds to a forward bifurcation. Above the threshold, the DFE becomes unstable and a stable endemic equilibrium (EE) appears. In the second case, there is backward bifurcation. Below the threshold, there are two EE, one stable and one unstable, and a stable DFE. Above the threshold, the DFE becomes unstable, one EE disappears and the other EE remains stable. In both cases, when $R_0$ is notably larger than 1, either the EE persists and remains stable or it disappears and the trivial equilibrium becomes stable. This model hence exhibits complex asymptotic properties: $R_0 < 1$ is not sufficient to eradicate the disease, as a stable EE may exist; moreover, high $R_0$ values may lead to the destruction of the coffee plantation.

**Keywords:** Epidemiological model, Host age classes, Coffee leaf rust, Basic reproduction number, Stability, Backward bifurcation

**MSC Classification:** 92D30, 37N25, 34D20, 34C23, 34C60
**3.1 Introduction**

Plant pest and pathogens account for up to 40% of yield loss worldwide [21]. In the current context of ever growing population and climate change, this problem needs to be tackled at scale. Hence, the interactions between pest/pathogens and their hosts need to be better understood. Mathematical models of fungal diseases have received a lot of attention from researchers. For instance, [49] used a mathematical model to describe the transition of a cereal rust population from endemic to epidemic densities and prove that natural enemies are likely to control only low-density cereal rust populations. [124] studied the seasonal appearance of rusts in the United States, in particular soybean rust, thanks to a general disease model. [134] investigated how the spatial deployment of resistant cultivars affected the resistance efficiency and durability, using a demogenetic model, for a seasonal crop infected by a fungal-like pathogen.

There are fewer models of fungal diseases that target perennial hosts. For instance, [23] explored the dynamical behaviour of mildew in a vineyard, while [97] studied the impact of spatial heterogeneities on its spread and control during a cropping season. [140] studied susceptible and resistant crop mixtures for a fungal disease propagated by airborne spores in a field. [36] 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. [130] looked at the impact of sexual and asexual reproduction on the epidemiological dynamics of fungal plant parasites and showed that they could induce cyclic persistence of the disease, which can occur with Sigatoka diseases of banana.

In this paper, we develop an epidemiological model to describe the spread of a fungal disease on a perennial plant in an agriculture context. The originality of this model is that it is structured to represent host developmental stages, as they often impact the aggressiveness of the pathogen, stemming from nutrient availability [32] or presence of different defense compounds [101] in the host.

The present work is motivated by the particular case of coffee leaf rust (CLR) on the coffee tree, which is a perennial plant that is exploited commercially for about 30 years [177]. CLR has a major impact on coffee production [180]. However, our study can be applied to other fungal diseases that target perennial hosts since most fungal diseases present the same type of life cycle, going through spore dispersion for asexual and/or sexual reproduction. CLR is caused by a basidiomycete fungus, *Hemileia vastatrix*, which infects coffee leaves of all stages, and dries them up before they fall down. The life cycle of *H. vastatrix* starts with the dispersion of uredospores, which land on the leaves. They germinate and penetrate the leaf, the infection process requires 48 hours to 10 days [116] depending of the age of the leaf and the climatic conditions [132]. Two to three weeks after the infection, yellow spots appear on the leaves. They grow and sporulate for two weeks to several months to produce two types of spores: uredospores for asexual reproduction and teliospores for sexual reproduction [172]. As the teliospore part of the cycle is not well known, except that these spores do not infect coffee leaves and might even be a dead end [83], it is not included in our model. For other rusts such as poplar rust, this part of the cycle is crucial for overwintering [60].

Mathematical models have been developed to study the epidemiology of CLR. [11] investigated the factors (coffee tree characteristics, crop management patterns, environment) that affect CLR intensity in several plots in Honduras. They showed that coffee rust development is linked
to three sets of factors: the environment, plant growth and development, and grower practices. [15] determined the germination and infection risk depending on the climate in Colombia and neighbouring countries, based upon existing experimental data. [165] represented the CLR dynamics in a coffee farm in Chiapas using an SI epidemiological model of the host. They showed that the network approach can be a useful way of gaining qualitative insight on spatial disease dynamics. They also showed in their study that coffee rust pattern suggests something other than random dispersal of spores. [38] studied the dynamics of CLR in a coffee plantation. However, existing models do not explicitly incorporate coffee plant development. In particular, in these models the differences between young and mature leaves are not considered while they clearly have an impact on the severity of infection. In fact, the mature leaves become infected more quickly than the young ones and it is only the latter that produce the coffee berries [122, 106].

Our aim is hence first to build a mathematical model describing the interactions between a fungal pathogen and a perennial plant, taking into account the development stages of the host and the epidemiological status of the leaves. It is based on the SEIR-U epidemiological model, where U stands for uredospores, presented in [38]. Compared to this previous work, this new model does not include the spatially-explicit spore dispersal, but it differentiates young and mature leaves and the young leaves grow logistically. We explicitly include a berry compartment to study the impact of CLR on crop production. Secondly, we analyse the equilibria of this complex ODE model and their stability.

This paper is organized as follows. Section 3.2 is devoted to the formulation of the model and derive its basic properties. In Section 3.3, we compute the disease free equilibrium and the basic reproduction number which determines its stability. In Section 3.4, we prove the existence of endemic equilibria and study their stability. We also presents different types of bifurcations that occurs in the model. Finally, we conclude the paper and propose several perspectives for future work.

## 3.2 Model formulation and basic properties

### 3.2.1 Formulation of the CLR model

Herein, we formulate a mathematical model for the CLR in the coffee plantation. To do so, we propose an ODE coffee-CLR interaction model that take into account a logistic growth of young leaves, the different stages of development of a leaf and the production of uredospores by all types of leaves. More precisely, we consider that in a coffee plantation we can find: susceptible young and mature leaves ($J_S$ and $M_S$), that have not (yet) been affected, latent young and mature leaves ($J_L$ and $M_L$), that are infected but not (yet) infectious, infectious young and mature leaves ($J_I$ and $M_I$), dry leaves ($M_R$), uredospores ($U$), which the fungus *H. vastatrix* uses for its asexual reproduction, and berries ($B$).

Figure 3.1 gives the compartmental representation of the corresponding ODE model.

The process begins with a logistic-like production of susceptible young leaves at rate $r_0 M_1 - r_1 M J_S$. In the first term, $M_1 = M_S + M_L + \varepsilon(M_I + M_R)$ and $\varepsilon < 1$, since all mature leaves can produce young ones, with the infectious and dry leaves producing less due to infection. The second term represents the competition between susceptible young leaves and all mature leaves with $M = M_S + M_L + M_I + M_R$. All young leaves become mature at rate $\beta$. Uredospores land on leaves of all leaves at deposition rate $\nu$, with fractions $\frac{J_S}{J+M}$ and $\frac{M_S}{J+M}$ on susceptible
<figure><img src="image_7.png" /><figcaption>Figure 3.1: Diagram of the CLR propagation in the coffee plantation corresponding to system (3.1). State variables are: susceptible young and mature leaves ($J_S$ and $M_S$), latent young and mature leaves ($J_L$ and $M_L$), infectious young and mature leaves ($J_I$ and $M_I$), dry leaves ($M_R$), uredospores ($U$) and berries ($B$).</figcaption></figure>

young and mature leaves, respectively, with $J = J_S + J_L + J_I$. Uredospores germinate at rates $\omega_1$ and $\omega_2$, respectively. Young and mature latent leaves become infectious at rate $\theta$, where $1/\theta$ corresponds to the average latency period. Infectious leaves sporulate over on average period $1/\alpha$; young leaves then die, while mature ones become dry. All young and mature leaves undergo natural mortality with baseline rate $\mu_J$ and $\mu_M$, respectively. Dry leaves have an additional mortality rate $d$. Uredospores are produced by infectious young and mature leaves at rate $\gamma_1$ and $\gamma_2$, respectively, and lose their infection ability at rate $\mu_U$. Berries are produced by all types of mature leaves at rates $\delta_S \ge \delta_L \ge \delta_I \ge \delta_R$ and die at rate $\mu_B$.

Table 3.1 summarises the parameter definitions and values of system (3.1).

From the compartmental model in Figure 3.1, we have the following system of ordinary differential equations:

$$ 
\left\{ 
\begin{aligned}
\dot{J}_S &= r_0 M_1 - r_1 M J_S - \frac{\omega_1 \nu U}{J+M} J_S - (\beta + \mu_J) J_S, \\
\dot{J}_L &= \frac{\omega_1 \nu U}{J+M} J_S - (\beta + \theta + \mu_J) J_L, \\
\dot{J}_I &= \theta J_L - (\beta + \alpha + \mu_J) J_I, \\
\dot{M}_S &= \beta J_S - \frac{\omega_2 \nu U}{J+M} M_S - \mu_M M_S, \\
\dot{M}_L &= \beta J_L + \frac{\omega_2 \nu U}{J+M} M_S - (\theta + \mu_M) M_L, \\
\dot{M}_I &= \beta J_I + \theta M_L - (\alpha + \mu_M) M_I, \\
\dot{M}_R &= \alpha M_I - (\mu_M + d) M_R, \\
\dot{U} &= \gamma_1 J_I + \gamma_2 M_I - (\nu + \mu_U) U, \\
\dot{B} &= \delta_1 M_S + \delta_2 M_L + \delta_3 M_I + \delta_4 M_R - \mu_B B,
\end{aligned}
\right.
\quad (3.1) $$

Since $B$ is not present in the other equations of system (3.1), we do not consider this state variable in the mathematical analysis. With this in mind, system (3.1) reduces to
<table><caption>Table 3.1: Description and values of parameters for system (3.1)</caption><thead><tr><th>Symbol</th><th>Biological meaning</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>r<sub>0</sub></td><td>Production rate of J<sub>S</sub></td><td>2 /day</td><td>Assumed</td></tr><tr><td>r<sub>1</sub></td><td>Competition rate</td><td>0.02 /leaf.day</td><td>Assumed</td></tr><tr><td>ε</td><td>Infection-related production ratio</td><td>0.06</td><td>Assumed</td></tr><tr><td>ω<sub>1</sub></td><td>Germination rate on J<sub>S</sub></td><td>0.055 leaves/spore</td><td>ref<sup>1</sup></td></tr><tr><td>ω<sub>2</sub></td><td>Germination rate on M<sub>S</sub></td><td>0.065 leaves/spore</td><td>ref<sup>1</sup></td></tr><tr><td>γ<sub>1</sub></td><td>Sporulation rate by J<sub>I</sub></td><td>5 spores/leaf.day</td><td>ref<sup>2</sup></td></tr><tr><td>γ<sub>2</sub></td><td>Sporulation rate by M<sub>I</sub></td><td>7 spores/leaf.day</td><td>ref<sup>2</sup></td></tr><tr><td>ν</td><td>Deposition rate</td><td>0.09 /day</td><td>ref<sup>2</sup></td></tr><tr><td>β</td><td>Maturation rate</td><td>0.05 /day</td><td>Assumed</td></tr><tr><td>θ</td><td>1 / duration of latency period</td><td>0.033 /day</td><td>ref<sup>3</sup></td></tr><tr><td>α</td><td>1 / duration of sporulation period</td><td>0.066 /day</td><td>Assumed</td></tr><tr><td>μ<sub>J</sub></td><td>Mortality rate of J</td><td>0.0054 /day</td><td>Assumed</td></tr><tr><td>μ<sub>M</sub></td><td>Mortality rate of M</td><td>0.0034 /day</td><td>Assumed</td></tr><tr><td>μ<sub>U</sub></td><td>Mortality rate of U</td><td>0.035 /day</td><td>ref<sup>4</sup></td></tr><tr><td>d</td><td>Extra mortality rate of M<sub>R</sub></td><td>0.056 /day</td><td>Assumed</td></tr><tr><td>δ<sub>S</sub></td><td>Berry production rate by M<sub>S</sub></td><td>0.7 berries/leaf.day</td><td>ref<sup>5</sup></td></tr><tr><td>δ<sub>L</sub></td><td>Berry production rate by M<sub>L</sub></td><td>0.5 berries/leaf.day</td><td>ref<sup>5</sup></td></tr><tr><td>δ<sub>I</sub></td><td>Berry production rate by M<sub>I</sub></td><td>0.3 berries/leaf.day</td><td>ref<sup>5</sup></td></tr><tr><td>δ<sub>R</sub></td><td>Berry production rate by M<sub>R</sub></td><td>0.05 berries/leaf.day</td><td>ref<sup>5</sup></td></tr><tr><td>μ<sub>B</sub></td><td>Mortality rate of B</td><td>0.0021 /day</td><td>Assumed</td></tr></tbody></table>

ref¹: [131], ref²: [19], ref³: [172],

ref⁴: [180] and ref⁵: [29].

$$ 
\left\{
\begin{aligned}
\dot{J}_S &= r_0 M_1 - r_1 M J_S - \frac{\omega_1 \nu U}{J+M} J_S - (\beta + \mu_J) J_S, \\
\dot{J}_L &= \frac{\omega_1 \nu U}{J+M} J_S - k_1 J_L, \\
\dot{J}_I &= \theta J_L - k_2 J_I, \\
\dot{M}_S &= \beta J_S - \frac{\omega_2 \nu U}{J+M} M_S - \mu_M M_S, \\
\dot{M}_L &= \beta J_L + \frac{\omega_2 \nu U}{J+M} M_S - k_3 M_L, \\
\dot{M}_I &= \beta J_I + \theta M_L - k_4 M_I, \\
\dot{M}_R &= \alpha M_I - k_5 M_R, \\
\dot{U} &= \gamma_1 J_I + \gamma_2 M_I - k_6 U.
\end{aligned}
\right.
\qquad (3.2) $$

where

$$ k_1 = \beta + \theta + \mu_J, \quad k_2 = \beta + \alpha + \mu_J, \quad k_3 = \theta + \mu_M, $$

$$ k_4 = \alpha + \mu_M, \quad k_5 = \mu_M + d, \quad k_6 = \nu + \mu_U, $$

System (3.2) is defined in

$$ \Gamma = \mathbb{R}_+^8 \setminus \{(0, 0, 0, 0, 0, 0, 0, U) \text{ with } U \ge 0\} $$

where $\mathbb{R}_+$ is the set of non-negative real numbers.

### 3.2.2 Basic properties

We obtain the following results for the positivity of solutions of system (3.2).

**Lemma 3.1.** For any initial condition in $\Gamma$, system (3.2) has a unique solution in $\Gamma$.
Proof of Lemma 3.1. Let $x = (J_S, J_L, J_I, M_S, M_L, M_I, M_R, U)$ and $x(0)$ be an initial condition in $\Gamma$ of system (3.2).

Suppose that one of the variables, that we will denote $x_i$ is equal to 0 at some time, with all others being non-negative. A quick analysis shows that $\dot{x}_i \ge 0$ so that $x_i$ cannot become negative. To ensure that any solution remains non-negative, we still need to prove that $J + M$ cannot become negative in finite time.

Using system (3.2), we obtain

$$ \begin{aligned} \dot{J} + \dot{M} &= r_0 M_1 - r_1 M J_S - \mu_J J - \mu_M M - \alpha J_I - d M_R, \\ &\ge -r_1 M J_S - \mu_J J - \mu_M M - \alpha J_I - d M_R. \end{aligned} $$

Using the fact that $J_S < J < J + M$, $M < J + M$, $J_I < J + M$ and $M_R < J + M$, one has

$$ \dot{\overline{J}} + \overline{M} \ge -r_1(J + M)^2 - \Upsilon(J + M) \quad (3.3) $$

where $\Upsilon = \mu_J + \mu_M + \alpha + d$. Using the comparison theorem to solve equation (3.3) yields

$$ J(t) + M(t) \ge \frac{\Upsilon(J(0) + M(0))}{\Upsilon e^{\Upsilon t} + (J(0) + M(0))(e^{\Upsilon t} - 1)}. $$

The above expression proves that $J(t) + M(t)$ cannot reach 0 in finite time. Therefore, the right hand-side of system (3.2) is Lipschitz continuous along all solutions, which implies that this system admits a unique solution that stays in $\Gamma$ at all time.

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

We also obtain the following result about the boundedness of solutions of system (3.2).

**Lemma 3.2.** All solutions of system (3.2) are bounded.

*Proof of Lemma 3.2.*

(1) Firstly we will show that the compact region

$$ \Omega = \left\{ \begin{array}{l} (J_S, J_L, J_I, M_S, M_L, M_I, M_R, U) \in \mathbb{R}_+^8, J_S \le \overline{J_S}, J_L \le \overline{K_{J_L}} \sqrt{U_{\max}}, \\ J_I \le \overline{K_{J_I}} \sqrt{U_{\max}}, M_S \le \overline{M_S}, M_L \le \overline{K_{M_L}} \sqrt{U_{\max}}, M_I \le \overline{K_{M_I}} \sqrt{U_{\max}}, \\ M_R \le \overline{K_{M_R}} \sqrt{U_{\max}} \quad \text{and} \quad U \le U_{\max}, \text{ with } U_{\max} \ge \overline{K_U} \end{array} \right\}, $$

is positively invariant for system (3.2), where

$$ \begin{aligned} \overline{J_S} &= \max(J_S(0), \frac{r_0}{r_1}), \quad \overline{M_S} = \max\left(M_S(0), \frac{\beta}{\mu_M} \overline{J_S}\right), \quad \overline{K_{J_L}} = \sqrt{\frac{\omega_1 \nu \overline{J_S}}{k_1}}, \\ \overline{K_{J_I}} &= \frac{\theta \overline{K_{J_L}}}{k_2}, \quad \overline{T_1} = \sqrt{\beta^2 \overline{K_{J_L}^2} + 4k_3 \omega_2 \nu \overline{M_S}}, \quad \overline{K_{M_L}} = \frac{\beta \overline{K_{J_L}} + \overline{T_1}}{2k_3}, \\ \overline{K_{M_I}} &= \frac{\beta \overline{K_{J_L}} + \theta \overline{K_{M_L}}}{k_4}, \quad \overline{K_{M_R}} = \frac{\alpha \overline{K_{M_I}}}{k_5}, \quad \overline{T_2} = \gamma_1 \overline{K_{J_I}} + \gamma_2 \overline{K_{M_I}}, \quad \text{and} \quad \overline{K_U} = \frac{\overline{T_2}^2}{k_6^2}. \end{aligned} $$

Consider the first equation of system (3.2):

$$ \dot{J}_S = r_0 M_1 - r_1 M J_S - \frac{\omega_1 \nu U}{J + M} J_S - (\beta + \mu_J) J_S. $$
Using the fact that $\frac{\omega_1 \nu U}{J+M} J_S \ge 0$ and $M_1 \le M$, we obtain the following equation

$$ \dot{J}_S \le r_0 M \left( 1 - \frac{J_S}{\frac{r_0}{r_1}} \right) - (\beta + \mu_J) J_S. $$

If $J_S \ge \frac{r_0}{r_1}$, we obtain $\dot{J}_S \le 0$, this implies that $\forall t \ge 0$, $J_S(t) \le \overline{J_S}$.

Similarly, let us consider the fourth equation of system (3.2).

$$ \dot{M}_S = \beta J_S - \frac{\omega_2 \nu U}{J+M} M_S - \mu_M M_S, $$

Using the fact that $\frac{\omega_2 \nu U}{J+M} M_S \ge 0$ and $J_S \le \overline{J_S}$, we obtain

$$ \dot{M}_S \le \beta \overline{J_S} - \mu_M M_S. $$

If $M_S \ge \frac{\beta}{\mu_M} \overline{J_S}$, one has $\dot{M}_S \le 0$, this allows to conclude for all $t \ge 0$ and $M_S(t) \le \overline{M_S}$.

Now for the boundness of the other variables, we use the hypothesis that $\forall t \ge 0$, $U(t) \le U_{\max}$ whose validity will be confirmed later. Let us consider the second equation of system (3.2):

$$ \dot{J}_L = \frac{\omega_1 \nu U}{J+M} J_S - k_1 J_L $$

Using the fact that $\frac{J_S}{J+M} \le \frac{\overline{J_S}}{\overline{J_L}}$ and $U \le U_{\max}$, we obtain

$$ \dot{J}_L \le \frac{\omega_1 \nu U_{\max} \overline{J_S}}{J_L} - k_1 J_L. \qquad (3.4) $$

If $J_L \ge \mathcal{K}_{J_L} \sqrt{U_{\max}}$ in the equation (3.4) we obtain $\dot{J}_L \le 0$. This implies that, if $J_L(0) \le \mathcal{K}_{J_L} \sqrt{U_{\max}}$, then $J_L(t) \le \mathcal{K}_{J_L} \sqrt{U_{\max}}$, for all $t \ge 0$.

Now, replacing the upper-bound value of $J_L$ into the third equation of system (3.2), gives

$$ \dot{J}_I \le \theta \mathcal{K}_{J_L} \sqrt{U_{\max}} - k_2 J_I. \qquad (3.5) $$

In the equation (3.5), one can observe that if $J_I \ge \mathcal{K}_{J_I} \sqrt{U_{\max}}$, we have $\dot{J}_I \le 0$ so that, when $J_I(0) \le \mathcal{K}_{J_I} \sqrt{U_{\max}} = \frac{\theta \mathcal{K}_{J_L}}{k_2} \sqrt{U_{\max}}$, we obtain $J_I(t) \le \mathcal{K}_{J_I} \sqrt{U_{\max}}$, for all $t \ge 0$.

Replacing the upper-bound of $J_L$ and $U$ into the fifth equation of system (3.2) and using the fact $\frac{M_S}{J+M} \le \frac{\overline{M_S}}{\overline{M_L}}$, yields

$$ \begin{aligned} \dot{M}_L &\le \beta \mathcal{K}_{J_L} \sqrt{U_{\max}} + \frac{\omega_2 \nu U_{\max} \overline{M_S}}{M_L} - (\theta + \mu_M) M_L \\ &= \frac{-k_3 M_L^2 + \beta \mathcal{K}_{J_L} \sqrt{U_{\max}} M_L + \omega_2 \nu U_{\max} \overline{M_S}}{M_L}. \end{aligned} \qquad (3.6) $$

The sign of the right-hand side of equation (3.6) is given by the sign of the second order equation of the numerator. Hence, it is negative when $M_L$ is larger than the largest of its two roots. This largest root is given by $\mathcal{K}_{M_L} \sqrt{U_{\max}} = \frac{\beta \mathcal{K}_{J_L} + T_1}{2k_3} \sqrt{U_{\max}}$. Then, one can conclude that when $M_L(0) \le \mathcal{K}_{M_L} \sqrt{U_{\max}}$, one has $M_L(t) \le \mathcal{K}_{M_L} \sqrt{U_{\max}}$ for all $t \ge 0$.
Replacing the upper-bound of $J_I$ and $M_L$ into sixth equation of system (3.2), gives

$$ \dot{M}_I \le \beta K_{J_L} \sqrt{U_{\max}} + \theta K_{M_L} \sqrt{U_{\max}} - k_4 M_I. \quad (3.7) $$

In the equation (3.7), one can observe that if $M_I \ge K_{M_I} \sqrt{U_{\max}} = \frac{\beta K_{J_L} + \theta K_{M_L}}{k_4} \sqrt{U_{\max}}$, we have $\dot{M}_I \le 0$ so that, when $M_I(0) \le K_{M_I} \sqrt{U_{\max}}$, we obtain $M_I(t) \le K_{M_I} \sqrt{U_{\max}}$, for all $t \ge 0$.

Replacing the upper-bound value of $M_I$ into the seventh equation of system (3.2), yields

$$ \dot{M}_R \le \alpha K_{M_I} \sqrt{U_{\max}} - k_5 M_R. \quad (3.8) $$

Thus, if $M_R \ge K_{M_R} \sqrt{U_{\max}} = \frac{\alpha K_{M_I}}{k_5} \sqrt{U_{\max}}$ in the equation (3.8) we obtain $\dot{M}_R \le 0$, this implies, if $M_R(0) \le K_{M_R} \sqrt{U_{\max}}$, then $M_R(t) \le K_{M_R} \sqrt{U_{\max}}$, for all $t \ge 0$.

Replacing the upper-bound value of $J_I$ and $M_I$ into the eighth equation of system (3.2), gives

$$ \dot{U} \le \gamma_1 K_{J_I} \sqrt{U_{\max}} + \gamma_2 K_{M_I} \sqrt{U_{\max}} - k_6 U. \quad (3.9) $$

Let $T_3 = \gamma_1 K_{J_I} + \gamma_2 K_{M_I}$, from the equation (3.9), one can observe that if $U \ge \frac{T_3}{k_6} \sqrt{U_{\max}}$, we obtain $\dot{U} \le 0$. Thus, if $U_{\max} \ge \frac{T_3}{k_6^2}$, we obtain $\dot{U} \le 0$ for $U = U_{\max}$. Hence, the positive invariance of $\Omega$ is shown.

(2) Secondly, for any given initial conditions, one can choose $U_{\max}$ large enough i.e

$$ U_{\max} \ge \max \left\{ \frac{J_L(0)^2}{K_{J_L}^2}, \frac{J_I(0)^2}{K_{J_I}^2}, \frac{M_L(0)^2}{K_{M_L}^2}, \frac{M_I(0)^2}{K_{M_I}^2}, \frac{M_R(0)^2}{K_{M_R}^2}, K_U \right\}, $$

such that the initial condition belongs to the corresponding positively invariant compact region $\Omega$. Hence the solution initiated there is bounded.

Finally, using these two items one can conclude that all solutions are bounded. This achieves the proof.

□

**Remark 3.1.** System (3.2) is not defined for $(J_S, J_L, J_I, M_S, M_L, M_I, M_R) = 0_{\mathbb{R}^7}$ but is extended to the trivial equilibrium $0_{\mathbb{R}^8}$ which makes sense biologically. However, the stability of this trivial equilibrium cannot be analysed as the Jacobian matrix of system (3.2) is not defined at this point.

## 3.3 The disease free equilibrium and its stability

The disease free equilibrium (DFE) of system (3.2) occurs in the absence of CLR, so the DFE is

$$ x^1 = (J_S^1, 0, 0, M_S^1, 0, 0, 0, 0) $$

where

$$ J_S^1 = \frac{r_0 \beta - \mu_M (\beta + \mu_J)}{\beta r_1} \quad \text{and} \quad M_S^1 = \frac{\beta J_S^1}{\mu_M}. \quad (3.10) $$
The above expression shows that the DFE exists under the following condition:

$$r_0\beta - \mu_M(\beta + \mu_J) > 0. \qquad (3.11)$$

If condition (3.11) does not hold, the plantation of coffee cannot survive.

### 3.3.1 The basic reproduction number and the local stability of the DFE

To compute the basic reproduction number and the study local stability of the DFE of system (3.2), we use the standard method of the next generation matrix developed by [162]. In system (3.2), we consider only the variables in which the infection is in progression and uredospores, i.e. $J_L, J_I, M_L, M_I, M_R$ and $U$. The corresponding equations can be rewritten in the following way

$$\frac{dx_I}{dt} = \mathcal{F}(x) - \mathcal{V}(x),$$

where $x_I = (J_L, J_I, M_L, M_I, M_R, U)$, $\mathcal{F}_i(x)$ represents the rate of appearance of the new infections in compartment $i$ and $\mathcal{V}(x)$ the transfer terms between compartments or the outside. We have

$$\mathcal{F}(x) = \begin{pmatrix} \frac{\omega_1 \nu U}{J+M} J_S \\ 0 \\ \frac{\omega_2 \nu U}{J+M} M_S \\ 0 \\ 0 \\ 0 \end{pmatrix} \quad \text{and} \quad \mathcal{V}(x) = \begin{pmatrix} k_1 J_L \\ -\theta J_L + k_2 J_I \\ -\beta J_L + k_3 M_L \\ -\beta J_I - \theta M_L + k_4 M_I \\ -\alpha M_I + k_5 M_R \\ -\gamma_1 J_I - \gamma_2 M_I + k_6 U \end{pmatrix}.$$

The Jacobian matrices $F$ of $\mathcal{F}$ and $V$ of $\mathcal{V}$ evaluated at the DFE $x^1$ are respectively

$$F = \begin{pmatrix} 0 & 0 & 0 & 0 & \frac{\omega_1 \nu J_S^1}{J_S^1+M_S^1} \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{\omega_2 \nu M_S^1}{J_S^1+M_S^1} \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \quad \text{and} \quad V = \begin{pmatrix} k_1 & 0 & 0 & 0 & 0 & 0 \\ -\theta & k_2 & 0 & 0 & 0 & 0 \\ -\beta & 0 & k_3 & 0 & 0 & 0 \\ 0 & -\beta & -\theta & k_4 & 0 & 0 \\ 0 & 0 & 0 & -\alpha & k_5 & 0 \\ 0 & -\gamma_1 & 0 & -\gamma_2 & 0 & k_6 \end{pmatrix}.$$

The basic reproduction number $R_0$ is obtained as the spectral radius of the positive matrix $FV^{-1}$, that is

$$R_0 = \frac{\nu\theta[\omega_1\mu_M[\beta\gamma_2(k_2+k_3) + \gamma_1k_3k_4] + \beta\omega_2\gamma_2k_1k_2]}{k_1k_2k_3k_4k_6(\beta + \mu_M)} \qquad (3.12)$$

Moreover, we have the following result about the local stability of the DFE.

**Lemma 3.3.** The DFE $x^1$ of system (3.2) is locally asymptotically stable (LAS) in $\Gamma$ when $R_0 < 1$, and unstable if $R_0 > 1$.

### 3.3.2 Global stability of DFE

Using a result obtained by Kamgang and Sallet [78], we have the following result about the global stability of the DFE.
**Theorem 3.1.** If $\mathcal{R}_0 \le \xi < 1$, where

$$ \xi = \frac{\mu_M}{\beta + \mu_M} \left( 1 + \frac{\nu\theta\omega_2\gamma_2}{k_3k_4k_6} \left( \frac{\beta}{\mu_M} - 1 \right) \right), $$

then the DFE $x^1$ of system (3.2) is globally asymptotically stable (GAS) in $\Gamma$.

*Proof of Theorem 3.1.* For the proof of this global stability, using the property of DFE, we rewrite system (3.2) in the following compact form

$$ \left\{ \begin{aligned} \dot{x}_S &= G_{11}(x_S, x_I)(x_S - x_S^1) - G_{12}(x_S, x_I)x_I \\ \dot{x}_I &= G_{22}(x_S, x_I)x_I \end{aligned} \right. \qquad (3.13) $$

where $x_S$ is the vector representing the state of different compartments of non-transmitting leaves (e.g. susceptible leaves) and the vector $x_I$ represents the state of compartments of different transmitting leaves (e.g. infected leaves). Here, we have $x_S = (J_S, M_S)^T$, $x_I = (J_L, J_I, M_L, M_I, M_R, U)^T$, $x_S^1 = (J_S^1, M_S^1)^T$, $x = (x_S, x_I)$,

$$ G_{11}(x_S, x_I) = \begin{pmatrix} -r_1M - (\beta + \mu_J) & r_0 - r_1J_S^1 \\ \beta & -\mu_M \end{pmatrix}, $$

$$ G_{12}(x_S, x_I) = \begin{pmatrix} 0 & 0 & r_0 - r_1 J_S^1 & \varepsilon r_0 - r_1 J_S^1 & \varepsilon r_0 - r_1 J_S^1 & -\frac{\omega_1 \nu J_S}{J+M} \\ 0 & 0 & 0 & 0 & 0 & -\frac{\omega_2 \nu M_S}{J+M} \end{pmatrix} $$

and

$$ G_{22}(x_S, x_I) = \begin{pmatrix} -k_1 & 0 & 0 & 0 & 0 & \frac{\omega_1 \nu J_S}{J+M} \\ \theta & -k_2 & 0 & 0 & 0 & 0 \\ \beta & 0 & -k_3 & 0 & 0 & \frac{\omega_2 \nu M_S}{J+M} \\ 0 & \beta & \theta & -k_4 & 0 & 0 \\ 0 & 0 & 0 & \alpha & -k_5 & 0 \\ 0 & \gamma_1 & 0 & \gamma_2 & 0 & -k_6 \end{pmatrix} $$

Let $\Omega \subset \mathcal{U} = \mathbb{R}_+^2 \times \mathbb{R}_+^6$. The right-hand side of system (3.13) is of class $C^1$ on the open set $\mathcal{U}$. Let us verify the following five hypotheses for the application of the Kamgang and Sallet's theorem.

$H_1$: System (3.13) is defined on a positively invariant set $\Omega$ of the nonnegative orthant. System (3.13) is dissipative on $\Omega$

$H_2$: Subsystem $\dot{x}_I = G_{11}(x_S, 0)(x_S - x_S^1)$ is globally asymptotically stable at the equilibrium $x_S^1 = (J_S^1, M_S^1)$ on the canonical projection $\mathcal{D} = \Omega \cap (\mathbb{R}_+^2 \times \{0_{\mathbb{R}^6}\})$ of $\Omega$ on $\mathbb{R}_+^2$

Indeed, this subsystem can be written as

$$ \left\{ \begin{aligned} \frac{dJ_S(t)}{dt} &= r_0 M_S - r_1 J_S M_S - (\beta + \mu_J) J_S, \\ \frac{dM_S(t)}{dt} &= \beta J_S - \mu_M M_S \end{aligned} \right. \qquad (3.14) $$

System (3.14) has two equilibria: a trivial which is unstable and the DFE $x_S^1 = (J_S^1, M_S^1)$, (where the values of $J_S^1$ and $M_S^1$ are giving above and exist under the condition (3.11)) which is globally asymptotically stable.
$H_3$: Matrix $G_{22}(x_S, x_I)$ is Metzler and irreducible for any given $x \in \Omega$.

$H_4$: There exists an upper-bound matrix $\hat{G}_{22}$ for $\mathfrak{M} = \{G_{22}(x) \in \mathcal{M}_6(\mathbb{R})/x \in \Omega\}$ with the property that either $\hat{G}_{22} \notin \mathfrak{M}$ or if $\hat{G}_{22} \in \mathfrak{M}$, (i.e., $\hat{G}_{22} = \max_{\Omega}\mathfrak{M}$), then for any $\bar{x} \in \Omega$ such that $\hat{G}_{22} = G_{22}(\bar{x})$, $\bar{x} \in \mathbb{R}_+^2 \times \{0\}$.

Indeed, consider the subsystem $\dot{x}_I = G_{22}(x_S, x_I)x_I$, on the matrix $G_{22}(x)$, we have $\frac{J_S}{J+M} \le 1$ and $\frac{M_S}{J+M} \le 1$, this implies that $\dot{x}_I \le \hat{G}_{22}x_I$, where $\hat{G}_{22}$ can be written as $\hat{G}_{22} = \hat{F} - V$, with $V$ given above and

$$ \hat{F} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & \omega_1\nu \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \omega_2\nu \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}. $$

$H_5$: The spectral radius of the matrix $\hat{G}_{22}$ verify $\sigma(\hat{G}_{22}) \le 0$

Indeed, consider the following auxiliary problem:

$$ \frac{d\hat{x}_I}{dt} = (\hat{F} - V)\hat{x}_I \qquad (3.15) $$

The matrix $\hat{F}$ and $V$ verified the condition of Theorem 2 in [162], then we have $\sigma(\hat{G}_{22}) = \sigma(\hat{F} - V) < 0$, which implies that $\mathcal{R}_C = \sigma(\hat{F}V^{-1}) \le 1$, where $\sigma(\hat{F} - V)$ represents the spectral radius of the matrix $\hat{F} - V$ (i.e. the largest real part of all the eigenvalues of this matrix) and $\hat{F}V^{-1}$ is the next generation matrix. The computation gives

$$ \mathcal{R}_C = \frac{\nu\theta[\omega_1[\beta\gamma_2(k_2 + k_3) + \gamma_1k_3k_4] + \omega_2\gamma_2k_1k_2]}{k_1k_2k_3k_4k_6} $$

Thus $\mathcal{R}_C \le 1$ ensures that all eigenvalues of $\hat{G}_{22}$ have a negative real parts. We have computed $\mathcal{R}_C$ and we have seen that the hypotheses $H_1, H_2, H_3, H_4$ and $H_5$ are satisfied. Then, by [79], we have proven that the DFE $x^1$ for the system (3.2) is GAS if $\mathcal{R}_C \le 1$.

One can observe that

$$ \mathcal{R}_0 = \frac{\mu_M}{\beta + \mu_M} \left( \mathcal{R}_C + \frac{\nu\theta\omega_2\gamma_2}{k_3k_4k_6} \left( \frac{\beta}{\mu_M} - 1 \right) \right), $$

and using the fact that $\mathcal{R}_C \le 1$, we obtain

$$ \mathcal{R}_0 \le \xi, $$

where

$$ \xi = \frac{\mu_M}{\beta + \mu_M} \left( 1 + \frac{\nu\theta\omega_2\gamma_2}{k_3k_4k_6} \left( \frac{\beta}{\mu_M} - 1 \right) \right) < 1. $$

Then, Lemma 3.3 and hypothesis $H_1 - H_5$ allow to conclude that the DFE is GAS if $\mathcal{R}_0 \le \xi < 1$, where $\xi$ is given above. This completes the proof. $\square$
## 3.4 Endemic equilibria and their stability

Herein we prove the existence of endemic equilibria (EE) of system (3.2) and analysed their stability.

We denote by $\lambda$ the force of infection given by

$$ \lambda = \frac{\nu U}{J + M}. \qquad (3.16) $$

Let $x^* = (J_S^*, J_L^*, J_I^*, M_S^*, M_S^*, M_L^*, M_I^*, M_R^*, U^*)$ be an endemic equilibrium of system (3.2). Setting the right hand side of this system to 0 yields

$$ \left\{ \begin{aligned} 0 &= r_0 M_1^* - r_1 M^* J_S^* - \omega_1 \lambda^* J_S^* - (\beta + \mu_J) J_S^* \\ J_L^* &= \lambda^* \frac{\omega_1}{k_1} J_S^*, \\ J_I^* &= \lambda^* \frac{\theta \omega_1}{k_1 k_2} J_S^*, \\ M_S^* &= \frac{\beta J_S^*}{(\mu_M + \lambda^* \omega_2)}, \\ M_L^* &= \lambda^* \left( \frac{\beta \omega_1}{k_1 k_3} J_S^* + \frac{\omega_2}{k_3} M_S^* \right), \\ M_I^* &= \lambda^* \left( \frac{\beta \theta \omega_1 (k_2 + k_3)}{k_1 k_2 k_3 k_4} J_S^* + \frac{\theta \omega_2}{k_3 k_4} M_S^* \right), \\ M_R^* &= \lambda^* \left( \frac{\beta \theta \alpha \omega_1 (k_2 + k_3)}{k_1 k_2 k_3 k_4 k_5} J_S^* + \frac{\theta \alpha \omega_2}{k_3 k_4 k_5} M_S^* \right), \\ U^* &= \lambda^* \left( \frac{\theta \omega_1 (\beta \gamma_2 (k_2 + k_3) + \gamma_1 k_3 k_4)}{k_1 k_2 k_3 k_4 k_6} J_S^* + \frac{\theta \gamma_2 \omega_2}{k_3 k_4 k_6} M_S^* \right). \end{aligned} \right. \qquad (3.17) $$

One can observe that there exists an endemic equilibrium $x^*$ if and only if $\lambda^*$ and $J_S^*$ are positive.

Substituting the expressions of $J_{\{L,I\}}^*$ and $M_{\{S,L,I,R\}}^*$ into the first equation of system (3.17) and solving this equation for $J_S^*$, we obtain either $J_S^* = 0$ which corresponds to the trivial equilibrium or

$$ J_S^* = \frac{b_2(\lambda^*)^2 + b_1\lambda^* + b_0}{\beta r_1(c_2(\lambda^*)^2 + c_1\lambda^* + c_0)}, \qquad (3.18) $$

where

$$ \left\{ \begin{aligned} b_2 &= \omega_1 \omega_2 (\beta r_0 (\theta \varepsilon (\alpha + k_5) (k_2 + k_3) + k_2 k_4 k_5) - k_1 k_2 k_3 k_4 k_5), \\ b_1 &= \beta r_0 k_2 [(k_5 (\varepsilon \theta + k_4) + \varepsilon \theta \alpha) (k_1 \omega_2 + \mu_M \omega_1)] + \beta r_0 \mu_M \theta \varepsilon \omega_1 k_3 (\alpha + k_5), \\ & \quad - k_1 k_2 k_3 k_4 k_5 (\omega_2 (\beta + \mu_J) + \mu_M \omega_1), \\ b_0 &= (\beta r_0 - \mu_M (\beta + \mu_J)) k_1 k_2 k_3 k_4 k_5, \\ c_2 &= \omega_1 \omega_2 [\theta (k_2 + k_3) (\alpha + k_5) + k_2 k_4 k_5], \\ c_1 &= k_2 (k_1 \omega_2 + \mu_M \omega_1) [k_5 (\theta + k_4) + \alpha \theta] + \mu_M \theta \omega_1 k_3 (\alpha + k_5), \\ c_0 &= k_1 k_2 k_3 k_4 k_5. \end{aligned} \right. $$

Note that, in the equation (3.18) when $\lambda^*$ is positive, the sign of $J_S^*$ is determined by its numerator in (3.18), as all c. are positive.

Firstly, let us study the positivity of $\lambda^*$. From equation (3.16) one has
$$ \lambda^* = \frac{\nu U^*}{J^* + M^*} \Rightarrow (J^* + M^*) = \frac{\nu U^*}{\lambda^*}, $$

Replacing the expression of $U^*$ given in system (3.17), we obtain

$$ J^* + M^* = \nu\theta \left[ \frac{\omega_1[\beta\gamma_2(k_2 + k_3) + \gamma_1k_3k_4]}{k_1k_2k_3k_4k_6} + \frac{\beta\omega_2\gamma_2}{k_3k_4k_6(\mu_M + \lambda^*\omega_2)} \right] J_S^*. $$

As $J^* + M^* = J_S^* + J_L^* + J_I^* + M_S^* + M_L^* + M_I^* + M_R^*$, replacing the expressions given in equation (3.17) into the left-hand side of the above equation, we obtain either $J_S^* = 0$ or the following second degree polynomial governed by $\lambda^*$:

$$ P(\lambda^*) = m_2(\lambda^*)^2 + m_1\lambda^* + m_0 \quad (3.19) $$

where

$$ \left\{ \begin{array}{l} m_2 = \omega_1\omega_2\left[\frac{\theta+k_2}{k_1k_2} + \frac{\beta}{k_1k_3} + \frac{\beta\theta(k_2+k_3)(\alpha+k_5)}{k_1k_2k_3k_4k_5}\right], \\ m_1 = \mu_M\frac{a_2}{\omega_2} + \omega_2\left[1 + \frac{\beta}{k_3} + \frac{\beta\theta(\alpha+k_5)}{k_3k_4k_5}\right] - \frac{\theta\omega_1\omega_2\nu[\beta\gamma_2(k_2+k_3)+\gamma_1k_3k_4]}{k_1k_2k_3k_4k_6}, \\ m_0 = (\beta + \mu_M)(1 - \mathcal{R}_0). \end{array} \right. $$

Now, using the Descartes sign rules for $\lambda$-polynomial (3.19), the number of positive solutions is given in Table 3.2.

<table><caption>Table 3.2: Descartes rule for $\lambda$-polynomial (3.19).</caption><thead><tr><th>m<sub>2</sub></th><th>m<sub>1</sub></th><th>m<sub>0</sub></th><th>R<sub>0</sub></th><th>Number of positive solutions</th></tr></thead><tbody><tr><td>+</td><td>+</td><td>+</td><td>R<sub>0</sub> &lt; 1</td><td>No solution</td></tr><tr><td>+</td><td>-</td><td>+</td><td>R<sub>0</sub> &lt; 1</td><td>No or two solutions</td></tr><tr><td>+</td><td>+</td><td>-</td><td>R<sub>0</sub> &gt; 1</td><td>One solution</td></tr><tr><td>+</td><td>-</td><td>-</td><td>R<sub>0</sub> &gt; 1</td><td>One solution</td></tr></tbody></table>

A positive value for $\lambda^*$ is not sufficient to obtain an endemic equilibrium of system (3.2). We also need to verify if the corresponding value of $J_S^*$, determined from equation (3.18), is positive. If both $\lambda^*$ and $J_S^*$ are positive, according to system (3.17), all state variables are positive, so it corresponds to an endemic equilibrium.

From the above proof we have the following Lemma

**Lemma 3.4.** If $\mathcal{R}_0 > 1$ and $J^* > 0$, the system (3.2) has one endemic equilibrium

### 3.4.1 Bifurcation analysis

To study the stability of the endemic equilibria around $\mathcal{R}_0$, we use Theorem 2.2 given in Subsection 2.1.2. We first introduce a positive parameter $\eta$ and set $\gamma_2 = \eta\gamma_1$. We then use $\gamma_1$ as bifurcation parameter, as its variations cover all $\mathcal{R}_0$ values. Indeed,

$$ \frac{\mathcal{R}_0}{\gamma_1} = \frac{\nu\theta(\omega_1\mu_M(\beta\eta(k_2+k_3)+k_3k_4)+\beta\omega_2\eta k_1k_2)}{k_1k_2k_3k_4k_6(\beta+\mu_M)} := \frac{1}{\zeta} > 0. \quad (3.20) $$
The linearisation of system (3.2) in the neighbourhood of DFE $x^1$ and around $\mathcal{R}_0 = 1$, using expression in equation (3.20) setting $\gamma_1$ to $\zeta$, yields

$$ \frac{dx(t)}{dt} = J(x^1)x(t), \qquad (3.21) $$

where

$$ J(x^1) = \begin{pmatrix} -\frac{r_0\beta}{\mu_M} & 0 & 0 & \frac{\mu_M(\beta+\mu_J)}{\beta} & \frac{\mu_M(\beta+\mu_J)}{\beta} & Q & Q & -\frac{\omega_1\nu\mu_M}{\beta+\mu_M} \\ 0 & -k_1 & 0 & 0 & 0 & 0 & 0 & \frac{\omega_1\nu\mu_M}{\beta+\mu_M} \\ 0 & \theta & -k_2 & 0 & 0 & 0 & 0 & 0 \\ \beta & 0 & 0 & -\mu_M & 0 & 0 & 0 & -\frac{\omega_2\nu\beta}{\beta+\mu_M} \\ 0 & \beta & 0 & 0 & -k_3 & 0 & 0 & \frac{\omega_2\nu\beta}{\beta+\mu_M} \\ 0 & 0 & \beta & 0 & \theta & -k_4 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \alpha & -k_5 & 0 \\ 0 & 0 & \zeta & 0 & 0 & \eta\zeta & 0 & -k_6 \end{pmatrix}, $$

$$ \text{where } Q = \frac{\mu_M(\beta+\mu_J)+r_0(\varepsilon-1)\beta}{\beta}. $$

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

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

Solving the above equation, yields

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

Similarly, the right eigenvector of matrix $J(x^1)$, denoted by $u = (u_1, u_2, u_3, u_4, u_5, u_6, u_7, u_8)^T$, satisfies $J(x^1)u = 0$,
Solving the above equation, we obtain

$$
\left\{
\begin{aligned}
u_1 &= \frac{u_8 \nu \mu_M \chi}{(\beta + \mu_M)((\beta + \mu_M)\mu_M - r_0 \beta)k_1 k_2 k_3 k_4 k_5}, \\
u_2 &= \frac{\omega_1 \nu \mu_M}{(\beta + \mu_M)k_1} u_8, \\
u_3 &= \frac{\theta \omega_1 \nu \mu_M}{(\beta + \mu_M)k_1 k_2} u_8, \\
u_4 &= \frac{\beta(-\nu \omega_2 u_8 + (\beta + \mu_M)u_1)}{(\beta + \mu_M)\mu_M}, \\
u_5 &= \frac{\beta \nu (\mu_M \omega_1 + \omega_2 k_1)}{(\beta + \mu_M)k_1 k_3} u_8, \\
u_6 &= \frac{\beta \nu \theta (\mu_M \omega_1 (k_2 + k_3) + \omega_2 k_1 k_2)}{(\beta + \mu_M)k_1 k_2 k_3 k_4} u_8, \\
u_7 &= \frac{\nu \theta \alpha \beta (\mu_M \omega_1 (k_2 + k_3) + \omega_2 k_1 k_2)}{(\beta + \mu_M)k_1 k_2 k_3 k_4 k_5} u_8, \\
u_8 &> 0.
\end{aligned}
\right.
$$

where

$$
\left\{
\begin{aligned}
\chi &= \chi_1 + \chi_2 + \chi_3 + \chi_4, \\
\chi_1 &= (\beta + \mu_J)\omega_1((k_2 + k_3)(\alpha + k_5)\theta + k_2k_4k_5)\mu_M^2, \\
\chi_2 &= ((r_0(\varepsilon - 1)\omega_1 + k_1\omega_2)k_2 + \omega_1k_3r_0(\varepsilon - 1))\beta + \mu_J\omega_2k_1k_2, \\
\chi_3 &= ((\alpha + k_5)\theta + k_1k_2k_4k_5(\beta\omega_2 - k_3\omega_1 + \mu_J\omega_2))\mu_M, \\
\chi_4 &= k_1\omega_2k_2(\beta r_0(\varepsilon - 1)(\alpha + k_5)\theta - k_3k_4k_5(\beta + \mu_J)).
\end{aligned}
\right.
$$

The expression b defined in Assumption A₂ of Theorem 2.2 for system (3.2) is rewritten by

$$
b = \sum_{k,i=1}^{8} v_k u_i \frac{\partial^2 F_k}{\partial x_i \partial \gamma_1}(x^1, \zeta),
$$

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

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

$$
v_8 \frac{\partial^2 F_8}{\partial x_3 \partial \gamma_1}(x^1, \zeta) = v_8 \quad \text{and} \quad v_8 \frac{\partial^2 F_8}{\partial x_6 \partial \gamma_1}(x^1, \zeta) = v_8 \eta.
$$

Substituting the above expression into the expression of b, yields

$$
b = v_8 u_3 \frac{\partial^2 F_8}{\partial x_3 \partial \gamma_1}(x^1, \zeta) + v_8 u_6 \frac{\partial^2 F_8}{\partial x_6 \partial \gamma_1}(x^1, \zeta) = v_8(u_3 + u_6\eta),
$$

Then

$$
b = v_8(u_3 + u_6\eta) > 0 \qquad (3.22)
$$

Let us compute the expression of a, defined in Assumption A₂ of Theorem 2.2 for system (3.2),
which is

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

The only terms that are non-null correspond to

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

Substituting these above expressions into the expression of a, yields

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

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

$$ \begin{aligned} a = & \frac{2\zeta\mu_M r_1\omega_1\nu\beta\theta u_8 v_8 (\mu_M \sum_{i=2}^7 u_i - \beta u_1) (\gamma_{1r}k_4k_3 + \gamma_{2r}\beta(k_2+k_3))}{k_1k_2k_3k_4(\beta + \mu_M)^2((\beta + \mu_J)\mu_M - \beta r_0)} \\ & + \frac{2r_1\theta\mu_M\zeta\beta\nu u_8 v_8 \gamma_{2r} (\beta(u_1 + u_2 + u_3 + u_5 + u_6 + u_7) - \mu_M u_4)}{(\beta + \mu_M)^2((\beta + \mu_J)\mu_M - \beta r_0)k_3k_4}. \end{aligned} \quad (3.23) $$

The expression b given in equation (3.22) is always positive, then the direction of bifurcation depends on the sign of a given by (3.23). This means that we are in the cases **1** and **4** of Theorem 2.2. Then, the direction of the bifurcation is determined by the sign of parameter a defined in equation (3.23). Using parameter values in Table 3.1, with $\eta = 7/5$, lead to $a < 0$ which corresponds to a forward bifurcation.

### 3.4.2 Forward bifurcation and stability

Using Maple to calculate the solutions of the $\lambda$-polynomial in equation (3.19), we obtain Table 3.3.

From this table, one can conclude that system (3.2) has one endemic equilibrium when $\mathcal{R}_0 > 1$ and that the DFE always exists.

According to the results in Section 3.4.1, when the bifurcation parameter changes from $\gamma_1 < \zeta$ to $\gamma_1 > \zeta$, i.e. the basic reproduction number changes from $\mathcal{R}_0 < 1$ to $\mathcal{R}_0 > 1$, the DFE changes from LAS to unstable. Moreover, the endemic equilibrium $J_{S,1}^*$ in Table 3.3 appears and is locally asymptotically stable.

To complement this analytical result, Figure 3.2 shows the existence of equilibria and their
<figure><img src="image_9.png" /><figcaption>Figure 3.2: Bifurcation diagram and dynamics of susceptible young leaves for the forward bifurcation case. Subplot (a) presents the ($R_0, J_S$)-bifurcation diagram with stable (blue curve) and unstable (red curve) equilibria of system (3.1). Subplot (b) shows the dynamics of susceptible young leaves $J_S$ for parameter values $R_0 = 0.8$ (black curve) and $R_0 = 4$ (orange curve) corresponding to the dotted lines in subplot (a). Bifurcation parameter is $\gamma_1$, $R_0 = \gamma_1/\zeta$ and $\gamma_2 = \eta\gamma_1$, with $\zeta$ defined in equation (3.20) and $\eta = 7/5$. All remaining parameters are given in Table 3.1, leading to $\zeta = 1.185$ spores/leaf.day and initial conditions defined by (3.24) with $U_0 = 1000$ spores.</figcaption></figure>
<table><caption>Table 3.3: Force of infection $\lambda$ and susceptible young leaves $J_S$ at equilibrium using parameter values in Table 3.1.</caption><thead><tr><th rowspan="2">R<sub>0</sub></th><th colspan="2">Equilibrium values for:</th></tr><tr><th>&lambda;</th><th>J<sub>S</sub></th></tr></thead><tbody><tr><td rowspan="4">R<sub>0</sub> = 0.8</td><td>&lambda;<sub>0</sub> = 0</td><td>J<sub>S</sub><sup>0</sup> = 0</td></tr><tr><td>&lambda;<sub>1</sub><sup>*</sup> &lt; 0</td><td>J<sub>S</sub><sup>1</sup> = 99.81</td></tr><tr><td>&lambda;<sub>2</sub><sup>*</sup> &lt; 0</td><td>J<sub>S,1</sub><sup>*</sup> &lt; 0</td></tr><tr><td></td><td>J<sub>S,2</sub><sup>*</sup> &lt; 0</td></tr><tr><td rowspan="4">R<sub>0</sub> = 4</td><td>&lambda;<sub>0</sub> = 0</td><td>J<sub>S</sub><sup>0</sup> = 0</td></tr><tr><td>&lambda;<sub>1</sub><sup>*</sup> = 1.27</td><td>J<sub>S</sub><sup>1</sup> = 99.81</td></tr><tr><td></td><td>J<sub>S,1</sub><sup>*</sup> = 81.94</td></tr><tr><td>&lambda;<sub>2</sub><sup>*</sup> &lt; 0</td><td>J<sub>S,2</sub><sup>*</sup> &lt; 0</td></tr></tbody></table>

stability. The initial conditions are chosen to be

$$x(0) = \begin{cases} J_S(0) = J_S^1, \quad J_L(0) = J_I(0) = 0, \quad M_S(0) = M_S^1, \\ M_L(0) = M_I(0) = M_R(0) = 0, \quad U(0) = U_0 > 0, \quad B(0) = 0 \end{cases}. \quad (3.24)$$

They correspond to the DFE except for the initial presence of uredospores and absence of berries.

For $R_0 = 0.8 < 1$ (black curves), subplot (a) shows that system (3.1) has two equilibria, the DFE which is stable and the trivial equilibrium which we showed to be unstable numerically; subplot (b) shows the convergence of susceptible young leaves to their DFE value $J_S^1$.

For $R_0 = 4 > 1$ (orange curves), subplot (a) shows that system (3.1) has three equilibria: the trivial equilibrium which is unstable, the DFE which is unstable and the endemic equilibrium which is stable; subplot (b) shows the convergence of susceptible young leaves to their endemic equilibrium value $J_{S,1}^*$.

Figure 3.3 shows the convergence of system (3.1), initially at the disease free equilibrium, towards the endemic equilibrium for different initial uredospore values. Susceptible leaves decrease over time. Infected leaves and uredospores overshoot their endemic steady state values. There are globally less leaves at the endemic equilibrium, compared to the initial disease free state. Compared to the DFE value ($B^1 = 48929$ berries), the number of berry at endemic equilibrium ($B^* = 8212$ berries) is drastically reduced by 83%.

Figure 3.4 differentiates CLR impact between young and mature leaves. Young leaves are less susceptible to the disease and hence less affected. There are globally less young leaves, but more young susceptible leaves than mature susceptible leaves at the endemic equilibrium. Indeed, as there are less mature leaves, young leaves are less impacted by the competition with mature leaves.

### 3.4.3 Backward bifurcation and stability

In this subsection we study the behavior of system (3.1) when there is a backward bifurcation, corresponding to $a > 0$ in equation (3.23). We choose ad hoc parameters given in Table 3.4.

Firstly, using Maple to calculate the equilibria of system (3.2), we obtain Table 3.5. In this table, negative or complex values of $\lambda^*$ and/or $J_S^*$ are included (text in gray), but they do not correspond to biological equilibria.
<figure><img src="image_10.png" /><figcaption>Figure 3.3: Persistence of CLR in the coffee plantation. Simulations of system (3.1) for $\mathcal{R}_0 = 4.218$, using initial conditions (3.24) with $U_0 = 10$ (blue curves), $U_0 = 100$ (red curves), $U_0 = 200$ (yellow curves) and $U(0) = 300$ spores (purple curves): (a) susceptible leaves $J_S + M_S$, (b) infected leaves $J_L + J_I + M_L + M_I + M_R$, (c) uredospores $U$ and (d) berries $B$. All parameter values are given in Table 3.1.</figcaption></figure>

<table><caption>Table 3.4: Parameter values of system (3.1) chosen to obtain a backward bifurcation.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>r<sub>0</sub></td><td>0.7094 /day</td><td>θ</td><td>0.2336 /day</td></tr><tr><td>r<sub>1</sub></td><td>0.0065 /day</td><td>α</td><td>0.8261 /day</td></tr><tr><td>ω<sub>1</sub></td><td>0.5092 leaves/spore</td><td>μ<sub>J</sub></td><td>0.0114 /day</td></tr><tr><td>ω<sub>2</sub></td><td>0.4709 leaves/spore</td><td>μ<sub>M</sub></td><td>0.0236 /day</td></tr><tr><td>η</td><td>0.0396</td><td>μ<sub>U</sub></td><td>0.0931 /day</td></tr><tr><td>ν</td><td>0.6269 /day</td><td>d</td><td>0.8961 /day</td></tr><tr><td>β</td><td>0.0555 /day</td><td>ε</td><td>0.9953 /day</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 3.4: CLR impact according to host development. Simulations of system (3.1) for $R_0 = 4.218$, using initial conditions (3.24) with $U_0 = 300$: (a) susceptible young $J_S$ and mature $M_S$ leaves, (b) latent young $J_L$ and mature $M_L$ leaves, (c) infectious young $J_I$, infectious $M_I$ and dry $M_R$ mature leaves, and (d) total young $J = J_S + J_L + J_I$ and mature $M = M_S + M_L + M_I + M_R$ leaves. Leaves are represented in a log scale that preserves 0 ($\log(1 + J_S)$ instead of $J_S$). All parameter values are given in Table 3.1.</figcaption></figure>
<table><caption>Table 3.5: Force of infection $\lambda$ and susceptible young leaves $J_S$ at equilibrium using parameter values in Table 3.4.</caption><thead><tr><th>R<sub>0</sub></th><th colspan="2">Equilibrium values for:</th></tr><tr><th></th><th>λ</th><th>J<sub>S</sub></th></tr></thead><tbody><tr><td rowspan="3">R<sub>0</sub> = 0.3</td><td>λ<sub>0</sub> = 0</td><td>J<sub>S</sub><sup>0</sup> = 0</td></tr><tr><td>λ<sub>1</sub><sup>*</sup> complex</td><td>J<sub>S,1</sub><sup>*</sup> complex</td></tr><tr><td>λ<sub>2</sub><sup>*</sup> complex</td><td>J<sub>S,2</sub><sup>*</sup> complex</td></tr><tr><td rowspan="3">R<sub>0</sub> = 0.8</td><td>λ<sub>0</sub> = 0</td><td>J<sub>S</sub><sup>0</sup> = 0</td></tr><tr><td>λ<sub>1</sub><sup>*</sup> = 0.31</td><td>J<sub>S</sub><sup>1</sup> = 104.76</td></tr><tr><td>λ<sub>2</sub><sup>*</sup> = 0.039</td><td>J<sub>S,1</sub><sup>*</sup> = 65.75</td></tr><tr><td rowspan="3">R<sub>0</sub> = 2</td><td>λ<sub>0</sub><sup>*</sup> = 0</td><td>J<sub>S,2</sub><sup>*</sup> = 100.11</td></tr><tr><td>λ<sub>1</sub><sup>*</sup> = 1.70</td><td>J<sub>S</sub><sup>0</sup> = 0</td></tr><tr><td>λ<sub>2</sub><sup>*</sup> < 0</td><td>J<sub>S</sub><sup>1</sup> = 104.76</td></tr><tr><td rowspan="3">R<sub>0</sub> = 4</td><td>λ<sub>0</sub> = 0</td><td>J<sub>S</sub><sup>0</sup> = 0</td></tr><tr><td>λ<sub>1</sub><sup>*</sup> = 3.92</td><td>J<sub>S</sub><sup>1</sup> = 104.76</td></tr><tr><td>λ<sub>2</sub><sup>*</sup> < 0</td><td>J<sub>S,1</sub><sup>*</sup> = -6.46</td></tr><tr><td></td><td></td><td>J<sub>S,2</sub><sup>*</sup> < 0</td></tr></tbody></table>
<figure><img src="image_12.png" /></figure>

Figure 3.5: Bifurcation diagram and dynamics of susceptible young leaves for the backward bifurcation case. Subplot (a) presents the ($\mathcal{R}_0, J_S$)-bifurcation diagram with stable (blue curve) and unstable (red curve) equilibria of system (3.1). Subplot (b) shows the dynamics of susceptible young leaves $J_S$ for parameter values $\mathcal{R}_0 = 0.3$ (blue curves), $\mathcal{R}_0 = 0.8$ (black curves), $\mathcal{R}_0 = 2$ (purple curves) and $\mathcal{R}_0 = 4$ (orange curves) corresponding to the dotted lines in subplot (a). The bifurcation parameter is $\gamma_1$, $\mathcal{R}_0 = \gamma_1/\zeta$ and $\gamma_2 = \eta\gamma_1$, with $\zeta$ defined in equation (3.20). All remaining parameters are given in Table 3.4, leading to $\zeta = 7.771$ spores/leaf.day and initial conditions defined by (3.24) with $U_0 = 100$ (dashed curves) and $U_0 = 1000$ spores (solid curves).

From this table, we conclude that system (3.2) has two endemic equilibria when $\mathcal{R}_0 < 1$ and one endemic equilibrium when $\mathcal{R}_0 > 1$.

When $\mathcal{R}_0 < 1$, the DFE is stable. Moreover, according to the results in Section 3.4.1, when $\mathcal{R}_0$ is close to one, there exist two endemic equilibria among which one is stable and the other unstable. When $\mathcal{R}_0 > 1$, the DFE is unstable. Moreover, when $\mathcal{R}_0$ is close to one, there is a unique endemic equilibrium, which is stable. For larger $\mathcal{R}_0$ values, the endemic equilibrium disappears.

Figure 3.5 illustrates the existence of the endemic equilibria and their stability. The diagram in subplot (a) presents a sequence of three bifurcations: a saddle-node bifurcation in $\mathcal{R}_0 \approx 0.7$, a backward transcritical bifurcation in $\mathcal{R}_0 = 1$ and an other transcritical bifurcation in $\mathcal{R}_0 \approx 2.7$. Subplot (b) presents the dynamics of susceptible young leaves $J_S$ for representative parameter values.

For $\mathcal{R}_0 = 0.3$ (blue curves), subplot (a) shows that system (3.1) has two equilibria: the trivial equilibrium which we showed to be unstable numerically and the DFE which is stable; subplot (b) shows the convergence of susceptible young leaves to their DFE value $J_S^1$.

For $\mathcal{R}_0 = 0.8$ (black curves), subplot (a) shows that system (3.1) has four equilibria: the trivial equilibrium which is unstable, the DFE which is stable, the endemic equilibrium $J_{S,1}^*$ which is stable and the endemic equilibrium $J_{S,2}^*$ which is unstable; subplot (b) shows the convergence of susceptible young leaves to their DFE value $J_S^1$ (for $U_0 = 100$ spores) and to their endemic equilibrium value $J_{S,1}^*$ (for $U_0 = 1000$ spores).

For $\mathcal{R}_0 = 2$ (orange curves), subplot (a) shows that system (3.1) has three equilibria: the trivial equilibrium which is unstable, the DFE which is unstable and the endemic equilibrium which is stable; subplot (b) shows the convergence of susceptible young leaves to their endemic equilibrium value $J_{S,1}^*$.

For $\mathcal{R}_0 = 4$ (orange curves), subplot (a) shows that system (3.1) has three equilibria: the trivial equilibrium which is stable and the DFE is unstable; subplot (b) shows the convergence
<figure><img src="image_13.png" /><figcaption>Figure 3.6: Extinction and persistence of CLR in the coffee plantation. Simulations of model (3.1) for $R_0 = 0.8$, using initial conditions (3.24) with $U_0 = 10$ (blue curves), $U_0 = 100$ (red curves), $U_0 = 200$ (yellow curves) and $U(0) = 300$ spores (purple curves): (a) susceptible leaves $J_S + M_S$, (b) infected leaves $J_L + J_I + M_L + M_I + M_R$, (c) uredospores $U$ and (d) berries $B$. All parameter values are given in Table 3.4 except $\gamma_1 = 6.2171$ spores/leaf.day.</figcaption></figure>

of susceptible young leaves to their value at trivial equilibrium.

For $R_0 = 0.8$, Figure 3.6 shows the convergence of system (3.1) towards the DFE and towards the stable endemic equilibrium for different initial uredospore values. For $U_0 = 10$ (blue curves) and $U_0 = 100$ spores (red curves), susceptible leaves decrease and then increase towards their DFE values. Simultaneously, the infected leaves and uredospores decrease towards zero and the berries increase. This convergence corresponds to the local stability of the DFE. For $U_0 = 200$ (yellow curves) and $U(0) = 300$ spores (purple curves), susceptible leaves decrease towards their endemic values. Infected leaves and uredospores first decrease, as for smaller initial uredospore values, and then increase towards their endemic values. Berries increase and then decrease towards their endemic values. This convergence corresponds to the local stability of one of the two endemic equilibria.

For $R_0 = 4$, Figure 3.7 shows the convergence of system (3.1) towards the trivial equilibrium for different initial uredospore values. Susceptible leaves decrease over time towards zero. Infected leaves, uredospores and berries first increase and then decrease towards zero. This shows the stability of the trivial equilibrium. This case corresponds to a very virulent pathogen that totally destroys the coffee plantation. Note that the dynamics does not depend on the initial conditions, except for the berries whose dynamics are amplified by the small differences
<figure><img src="image_14.png" /><figcaption>Figure 3.7: Destruction of the plantation by CLR. Simulations of model (3.1) for $\mathcal{R}_0 = 4$, using initial conditions (3.24) with $U_0 = 10$ (blue curves), $U_0 = 100$ (red curves), $U_0 = 200$ (yellow curves) and $U(0) = 300$ spores (purple curves): (a) susceptible leaves $J_S + M_S$, (b) infected leaves $J_L + J_I + M_L + M_I + M_R$, (c) uredospores U and (d) berries B. All parameter values are given in Table 3.4 except $\gamma_1 = 31.0856$ spores/leaf.day.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 3.8: Subplots present the ($\mathcal{R}_0, J_S$)-bifurcation diagram with stable (blue curve) and unstable (red curve) equilibria of system (3.1). The bifurcation parameter is $\gamma_1$, $\mathcal{R}_0 = \zeta\gamma_1$ and $\gamma_2 = \eta\gamma_1$, with $\zeta$ defined in equation (3.20). Subplot (a) shows a forward bifurcation using parameter values in Table 3.4 except $\mu_M = 0.1236$/day leading to $\zeta = 3.676$ spores/leaf.day. Subplot (b) shows the backward bifurcation using parameter values in Table 3.1 except $d = 0.956$/day and $\eta = 0.3/20$ spores/leaf.day leading to $\zeta = 59.306$ spores/leaf.day.</figcaption></figure>

in susceptible leaves.

## 3.5 Other types of bifurcation

Beside the forward and backward bifurcation diagrams presented in section 3.4. Figure 3.8 shows two different bifurcation diagrams. Subplot (a) presents a forward bifurcation, as in Figure 3.2, but where the endemic equilibrium disappears for large $\mathcal{R}_0$ values while the trivial equilibrium becomes stable. Subplot (b) presents a backward bifurcation, as in Figure 3.5, but where the endemic equilibrium persists for large $\mathcal{R}_0$ values.

## 3.6 Conclusion

In this paper, we proposed and analysed a mathematical model that describes the propagation of a fungal disease on a perennial plant. The originality of the model is the differentiation between young and mature leaves, to take into account the impact of host development on the aggressiveness of the pathogen. The model is applied to the particular case of coffee leaf rust in a coffee plantation. It is fairly generic and can be adapted to other fungal diseases for which the host age structure plays an important role.

We computed the disease-free and endemic equilibria and we studied their stability as a function of the basic reproduction number. We identified different bifurcation diagrams. They showed two behaviours for the model when the basic reproduction number is close to one. The first case is a forward bifurcation which corresponds to a stable disease free equilibrium when the basic reproduction number is less than one, and a stable endemic equilibrium when it is greater than one. The second case is a backward bifurcation in which, on top of the stable disease free equilibrium, there is a stable endemic equilibrium when the basic reproduction number is less than one. This means that, disease may persist even in these unfavourable conditions. Moreover, when the basic reproduction number is notably larger than one, be it in the forward or backward
bifurcation case, the endemic equilibrium may disappear and the trivial equilibrium then becomes
stable, which corresponds to the destruction of the plantation.

All bifurcation cases allow to conclude that, when the basic reproduction number is greater
that one, the disease can destroy the plantation and even when the basic reproduction number
is less than one, the disease can persist. These results provide useful insights on the disease
dynamics and its control. For instance, in order to eradicate a disease that spreads with basic
reproduction number greater than one, it might not be sufficient to bring it below one, as an
endemic state still persists. The basic reproduction number should be brought below the level
of the saddle-node bifurcation to achieve effective disease control.
# MODELLING COFFEE LEAF RUST
DYNAMICS TO CONTROL ITS SPREAD

In this chapter, we simplify the model in Chapter 3 by combining young and mature leaves into branches. Moreover, we distinguish between dry and rainy season and introduce spore dispersal. 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 branches $S$, latent branches $L$, infectious branches $I$ and leafless branches $J$), considering the crop branch as an individual. The model includes crop production (berries $B$) by branches during the rainy 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 branches 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 branches.

$H_4$: The model alternates dry and rainy seasons. The same model is used during both seasons, but with different parameter values: crop development and disease progression are reduced during the dry seasons.

$H_5$: Space is taken into account in this model. The branches are homogeneously distributed. The 1D dispersal of uredospores is considered.

$H_6$: The model uses a PDE formalism.

$H_7$: A parameter which reduces spore production is added to control the disease.

We then obtain a PDE model with six compartments. The objective of this second study is to understand the impact of spore diffusion and seasonality on the dynamics of CLR.

This work is published in *Mathematical Modelling of Natural Phenomena* [38]
Contents

<table><tr><td>4.1</td><td>Introduction</td><td>61</td></tr><tr><td>4.2</td><td>Coffee leaf rust dynamics</td><td>62</td></tr><tr><td>4.2.1</td><td>Biological background</td><td>62</td></tr><tr><td>4.2.2</td><td>Model formulation</td><td>62</td></tr><tr><td>4.3</td><td>Mathematical analysis</td><td>65</td></tr><tr><td>4.3.1</td><td>Basic properties</td><td>65</td></tr><tr><td>4.3.2</td><td>Equilibria and their stability</td><td>70</td></tr><tr><td>4.3.3</td><td>Comparison of the dynamics in the dry and rainy seasons</td><td>72</td></tr><tr><td>4.4</td><td>Numerical simulations</td><td>73</td></tr><tr><td>4.4.1</td><td>CLR dynamics</td><td>74</td></tr><tr><td>4.4.2</td><td>Control of CLR</td><td>75</td></tr><tr><td>4.5</td><td>Conclusion</td><td>80</td></tr><tr><td>4.6</td><td>Appendix: Mathematical analysis</td><td>80</td></tr><tr><td>4.7</td><td>Appendix: Numerical scheme</td><td>86</td></tr><tr><td>4.7.1</td><td>Discretization of model system</td><td>86</td></tr><tr><td>4.7.2</td><td>Discretization of boundary conditions</td><td>88</td></tr><tr><td>4.7.3</td><td>Consistence of scheme</td><td>88</td></tr></table>

**Abstract:** Coffee leaf rust (CLR) is one of the main diseases that affect coffee plantations worldwide. It is caused by the fungus *Hemileia vastatrix*. Damages induce severe yield losses (up to 70%). Its control mainly relies on cultural practices and fungicides, the latter having harmful ecological impact and important cost. Our goal is to understand the propagation of this fungus in order to propose a biocontrol solution, based on a mycoparasite that inhibits *H. vastatrix* reproduction. We develop and explore a spatio-temporal model that describes CLR propagation in a coffee plantation during the rainy and dry seasons. We show the existence of a solution and prove that there exists two threshold parameters, the dry and rainy basic reproduction numbers, that determine the stability of the equilibria for the dry and rainy season subsystems. To illustrate these theoretical results, numerical simulations are performed, using a non-standard finite method to integrate the pest model. We also numerically investigate the biocontrol impact. We determine its efficiency threshold in order to ensure CLR eradication.

**Keywords:** spatio-temporal model, coffee leaf rust, basic reproduction number, stability

**MSC Classification:** 35K57, 93D05, 65M06, 92D30
## 4.1 Introduction

Coffee is one of the most widely consumed beverages in the world; its trade satisfies the regular consumption of more than two billion people and exceeds 10$ billion worldwide [158]. Its cultivation is an important factor of social stability as it sustains the living of not less than twenty-five million small producers and their families worldwide by ICO [71]. The total production of all exporting countries in 2018 was more than 172 million 60-kilogram bags [71].

Coffee leaf rust (CLR) is a disease caused by a basidiomycete fungus, *Hemileia vastatrix*, that affects coffee trees. The fungus *H. vastatrix* is a compulsory parasite. In other words, it is a fungus that lives and develops only on coffee leaves. 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. CLR is the most serious leaf disease of the coffee tree known to date [112]. 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]. Indirect impacts include increased costs to control the disease. Control methods include fungicide application [25], cultural practices [10] including stumping of diseased plants, the use of resistant cultivars [5], and biocontrol by fungal and bacterial parasites [166, 6]. These methods induce significant labour and material costs and, in the case of stumping, a year-long decline in production. To better control this fungus, it is necessary to understand the disease propagation and what conditions favour its development and dispersal.

The spread of crop diseases, in particular airborne pathogens such as fungi, has received a lot of attention from researchers. Among models that represent the pathogen spatial dispersal, one can cite the DDAL framework that focuses on the deployment of susceptible and resistant crop hosts in an agricultural landscape [120]. Fewer models represent the pathogen spread by a diffusion term in partial differential equations (PDE). For instance, Sapoukhina et al. also study susceptible and resistant crop mixtures for a fungal disease propagated by airborne spores in a field [140], while Burie et al. explore the dynamical behaviour of mildew in a vineyard [23]. These disease dynamics are relevant for CLR modelling: they include a latency period, a sporulation period, spore dispersal and germination.

CLR models in the literature represent different scales, from the individual coffee bush to the country or even the continent. Avelino et al. investigate the factors (coffee tree characteristics, crop management patterns, environment) that affect CLR intensity in several plots in Honduras [11]. Bebber et al. determine the germination and infection risk depending on the climate in Colombia and neighbouring countries [15]. In contrast to these static approaches, Vandermeer et al. study 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 represent the CLR dynamics in a coffee farm in Chiapas using an SI epidemiological model of the host [165]. In these two latter studies, the fungus life cycle is not represented. Some other models investigate CLR control. Vandermeer et al. look at the interaction between *H. vastatrix* and a mycoparasite *Lecanicillium lecanii* [166], while Arroyo et al. consider interactions with antifungal bacteria [6]. However, no existing CLR model considers *H. vastatrix* dynamics together with its interaction with the coffee host. In particular none considered the impact of CLR on berry production, which is the variable of agronomic interest.

Our aim is to understand the CLR propagation in a field and to propose a biocontrol solution
based on the *Lecanicillium lecanii* mycoparasite, by a mathematical modelling approach. Bio-
control methods are still in development but could provide ecologically friendly alternatives to
fungicides, which can induce pest resistance. Moreover, biocontrol can be coupled with cultural
practices and does not require host replacement as for resistant plant deployment. To achieve
our goal, an original CLR model is needed, based on existing crop-fungus interaction models
with diffusion [23, 140] and including berry production. As coffee growth is climate-dependent
and harvest is seasonal, we have to develop a hybrid spatio-temporal PDE model.

This paper is organized as follows. Section 4.2 is devoted to the biological background and the formulation of the spatio-temporal CLR propagation model with two seasons per year, the dry and rainy seasons. Section 4.3 presents the stability analysis of the subsystems during the two seasons. In Section 4.4, numerical simulations are performed to illustrate and validate theoretical results. We also present numerical results when mycoparasites are used to control CLR spread. Finally, we conclude the paper and propose several possible perspectives for future work.

4.2 Coffee leaf rust dynamics

4.2.1 Biological background

The coffee tree is a perennial plant, belonging to the Rubiaceae family, which has persistent
leaves [16]. The productivity stage (flowers and fruits) begins after 3 years and lasts for 20
years on average. Flowers appear at the end of branches and require a lot of humidity for their
development [28]. The annual productivity period lasts 7 to 11 months depending on coffee
species and cultivars, as well as climate. In the case of Cameroon, with a single rainy period per
year from April to November, there is a single harvest between November and December [151].

The life cycle of *H. vastatrix* in favourable humidity and temperature conditions lasts 5 weeks on average [172]. Uredospores are dispersed by rain and wind. They germinate and penetrate through stomates on the underside of the leaf. This infection process requires 24 to 48 hours [116]. The first symptom is a pale yellow lesion that appears 1 to 3 weeks after infection. Sporulation, i.e. the production of Uredospores and teliospores, occurs 2 weeks to several months after infection [172]. A single lesion produces 4 to 6 crops of spores, releasing 300,000 to 400,000 Uredospores over a period of 3 to 5 months [116].

As most rust fungi, *H. vastatrix* produces Uredospores for asexual reproduction and teliospores for sexual reproduction. Teliospores do not infect coffee leaves and have no known host [136], so it is not established that there is an effective sexual reproduction cycle.

During the dry season, *H. vastatrix* survives primarily as mycelium in the living tissues of the coffee leaves. As infected leaves drop prematurely during the dry months, a large amount of potential inoculum for the next rainy season disappears. However, a few green leaves always persist. Moreover, uredospores can survive about 6 weeks, so there is always some viable inoculum to infect the newly formed leaves at the start of the next rainy season [5].

4.2.2 Model formulation

Based on this biological knowledge, we propose a spatio-temporal coffee-CLR interaction model within a coffee plantation. We represent space and uredospore diffusion inside a bounded domain. At each point in space, we describe the evolution of the number of coffee branches, according to their epidemiological state, based on the CLR life cycle: healthy, latent with spores germinating
<figure><img src="image_16.png" /><figcaption>Figure 4.1: Diagram of the CLR propagation model in the coffee plantation corresponding to system (4.1). State variables are: healthy branches (S), latent branches (L), infectious branches (I), leafless branches (J), Uredospores (U) and berries (B). $\varepsilon\Delta U$ corresponds to uredospore diffusion.</figcaption></figure>

or non sporulating lesions, infectious with sporulating lesions, and leafless at the end of the infection process. The infection of a branch corresponds to the infection of its leaves. We consider two periods: the rainy season during which berries are produced, and the dry season. This implies that most parameters are time-dependent, as they take two values according to the season.

We obtain the following model, that is illustrated in Figure 4.1:

$$ 
\left\{
\begin{aligned}
\partial_t S(t, x) &= \Lambda(t) - \frac{\omega(t)\nu U(t, x)}{N(t, x)} S(t, x) - \mu(t)S(t, x), \\
\partial_t L(t, x) &= \frac{\omega(t)\nu U(t, x)}{N(t, x)} S(t, x) - (\theta(t) + \mu(t))L(t, x), \\
\partial_t I(t, x) &= \theta(t)L(t, x) - (\alpha(t) + \mu(t) + d)I(t, x), \\
\partial_t J(t, x) &= \alpha(t)I(t, x) - \mu(t)J(t, x), \\
\partial_t U(t, x) &= \varepsilon\Delta U(t, x) + \gamma(t)I(t, x) - (\nu + \mu_U)U(t, x), \\
\partial_t B(t, x) &= \delta_S(t)S(t, x) + \delta_L(t)L(t, x) + \delta_I(t)I(t, x) + \delta_J(t)J(t, x) - \mu_B B(t, x),
\end{aligned}
\right.
\qquad (4.1) $$

where $S(t, x)$, $L(t, x)$, $I(t, x)$, $J(t, x)$, $U(t, x)$ and $B(t, x)$ represent the densities of healthy branches, latent branches, infectious branches, leafless branches, Uredospores and berries respectively, at time $t$ and location $x$ defined in $\Omega$, a simply connected bounded domain (of $\mathbb{R}$ or $\mathbb{R}^2$) with smooth boundary $\partial\Omega$. The densities are $/m$ or $/m^2$ depending on the choice of $\Omega$. The total density of branches is $N(t, x) = S(t, x) + L(t, x) + I(t, x) + J(t, x)$.

The recruitment of healthy branches occurs at rate $\Lambda(t)$. Uredospores are deposited on leaves of all branches at rate $\nu$ and a fraction $S/N$ lands on healthy branches. These spore covered healthy branches become latent branches with rate $\omega(t)$, which represents the germination ef-
ficacy, that is the number of healthy branches which become latent branches, per deposited spore. The latent branches become infectious at rate $\theta(t)$, where $1/\theta(t)$ corresponds to the latency period, and in turn become leafless branches at rate $\alpha(t)$, where $1/\alpha(t)$ corresponds to the sporulation period. All branches undergo natural mortality with baseline rate $\mu(t)$ and the infectious branches have an additional constant mortality rate $d$ due to the disease. Uredospores are produced by infectious branches at rate $\gamma(t)$ and lose their ability to infect coffee branches at constant rate $\mu_U$. Parameter $\varepsilon$ is the uredospore diffusion coefficient. Berries are produced by all types of branches at different rates: $\delta_S(t) \ge \delta_L(t) \ge \delta_I(t) \ge \delta_J(t)$. They have a constant mortality rate $\mu_B$. All parameters are assumed to be non negative.

Having two seasons, let us consider $T$ the yearly period and $\tau$ the duration of the dry season. We set the initial time $t = 0$ at the beginning of a dry season. Then, for $n$ years, we assume that time-dependent parameters $m(t)$, where $m(t) \in \{\Lambda(t), \omega(t), \theta(t), \alpha(t), \gamma(t), \mu(t)\}$, and $\delta_i(t)$, with $i \in \{S, L, I, J\}$, are constant during each season and can be written as follows:

$$m(t) = \begin{cases} m_D & \text{for } t \in [nT, nT + \tau[, \\ m_R & \text{else,} \end{cases} \quad \text{and} \quad \delta_i(t) = \begin{cases} 0 & \text{for } t \in [nT, nT + \tau[, \\ \bar{\delta}_i & \text{else,} \end{cases} \quad (4.2)$$

where $D$ represents the dry season and $R$ the rainy season.

Table 4.1 summarizes the biological meaning of parameter values for system (4.1).

<table><caption>Table 4.1: Description of parameters for system (4.1) – numerical values are used in section 4.4.</caption><thead><tr><th rowspan="2">Symbol</th><th rowspan="2">Biological meaning</th><th colspan="2">Value</th><th rowspan="2">Source</th></tr><tr><th>Dry season (D)</th><th>Rainy season (R)</th></tr></thead><tbody><tr><td>Λ<sub>p</sub></td><td>Recruitment rate</td><td>6 day<sup>-1</sup></td><td>8 day<sup>-1</sup></td><td>Assumed</td></tr><tr><td>ω<sub>p</sub></td><td>Inoculum effectiveness</td><td>4.5%</td><td>5.5%</td><td>[131]</td></tr><tr><td>δ&#773;<sub>S</sub></td><td>Production rate of berries by S</td><td>0</td><td>0.7 day<sup>-1</sup></td><td>[71, 29]</td></tr><tr><td>δ&#773;<sub>L</sub></td><td>Production rate of berries by L</td><td>0</td><td>0.5 day<sup>-1</sup></td><td>[71, 29]</td></tr><tr><td>δ&#773;<sub>I</sub></td><td>Production rate of berries by I</td><td>0</td><td>0.3 day<sup>-1</sup></td><td>[71, 29]</td></tr><tr><td>δ&#773;<sub>J</sub></td><td>Production rate of berries by J</td><td>0</td><td>0.05 day<sup>-1</sup></td><td>[71, 29]</td></tr><tr><td>1/θ<sub>p</sub></td><td>Latency period duration</td><td>30 days</td><td>21 days</td><td>[172]</td></tr><tr><td>μ<sub>p</sub></td><td>Mortality of mature branches</td><td>0.0134 day<sup>-1</sup></td><td>0.0034 day<sup>-1</sup></td><td>Assumed</td></tr><tr><td>1/α<sub>p</sub></td><td>Sporulation period duration</td><td>150 days</td><td></td><td>[172]</td></tr><tr><td>γ<sub>p</sub></td><td>Production rate of Uredospores</td><td>∈ [0, 20] day<sup>-1</sup></td><td></td><td>[19]</td></tr><tr><td>μ<sub>B</sub></td><td>Mortality rate of berries</td><td>0.0021 day<sup>-1</sup></td><td></td><td>Assumed</td></tr><tr><td>d</td><td>Mortality rate due to the disease</td><td>0.056 day<sup>-1</sup></td><td></td><td>Assumed</td></tr><tr><td>ε</td><td>Diffusion coefficient</td><td>5000 m<sup>2</sup> day<sup>-1</sup></td><td></td><td>[23]</td></tr><tr><td>ν</td><td>Deposition rate</td><td>0.09 day<sup>-1</sup></td><td></td><td>[19]</td></tr><tr><td>μ<sub>U</sub></td><td>Mortality rate of Uredospores</td><td>0.015 day<sup>-1</sup></td><td></td><td>[92]</td></tr></tbody></table>

Initial conditions are $(\phi_S(x), \phi_L(x), \phi_I(x), \phi_J(x), \phi_U(x)) \in (L^2(\Omega))^5$ where

$$\begin{cases} S(0, x) = \phi_S(x) > 0, & L(0, x) = \phi_L(x) \ge 0, & I(0, x) = \phi_I(x) \ge 0, \\ J(0, x) = \phi_J(x) \ge 0, & U(0, x) = \phi_U(x) \ge 0, & \forall x \in \Omega. \end{cases} \quad (4.3)$$

with $B(0, x) = 0 \forall x \in \Omega$. We assume that the harvest occurs instantaneously at the end of each rainy season. Thus $B(t, x)$ is brought back to 0 at $t = nT$, with other variables unchanged.

We suppose that the domain $\Omega$ is isolated, which implies that spores do not get in or out and that there is no interaction with the outside environment. To ensure this property, we assume
that there is a large enough buffer zone without coffee plants around the coffee plantation (hence around the domain $\Omega$), so that no Uredospores are introduced from outside the plantation. Moreover, we neglect the spores that could potentially get out of the plantation. This absence of transfer is classically represented through the homogeneous Neumann boundary:

$$ \frac{\partial U}{\partial \eta} = 0, \quad \text{on} \quad \partial \Omega \qquad (4.4) $$

where $\eta$ denotes the unit outward normal on $\partial \Omega$.

**Remark 4.1.** When considered in one dimension, the space domain simply is an interval $\Omega = [x_{\min}, x_{\max}]$, whose boundary as considered in (4.4) is the pair of points $\partial\Omega = \{x_{\min}, x_{\max}\}$, and condition (4.4) becomes $\frac{\partial U}{\partial x} = 0$ in $x = x_{\min}$ and $x = x_{\max}$.

## 4.3 Mathematical analysis

Herein, we present the basic properties of the subsystems of system (4.1), defined over each season. For any parameter $m$, we denote by $m_p$ the seasonal constant value of the parameter $m$ where $p = R$ for the rainy season and $p = D$ for the dry season. Replacing the value of every parameter $m(t)$ of system (4.1), with the corresponding constant value $m_p$ with respect to the season, and removing the $\partial_t B$ equation since $B$ is not present in the other equations of system (4.1), we obtain:

$$ \left\{ \begin{array}{l} \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 = \gamma_p I - (\nu + \mu_U)U, \\ \frac{\partial U}{\partial \eta} = 0, \quad \text{on} \quad \partial \Omega, \\ S(0, x) = \phi_S, L(0, x) = \phi_L, I(0, x) = \phi_I, J(0, x) = \phi_J, U(0, x) = \phi_U, \forall x \in \Omega. \end{array} \right. \qquad (4.5) $$

### 4.3.1 Basic properties

In this section, we will first establish positivity of the solutions if they exist in Lemma 4.1, their boundedness in Lemma 4.2, and finally their existence and uniqueness in Lemma 4.4. Through these, we obtain positivity, boundedness and global existence of solutions of subsystem (4.5) in Theorem 4.1.

In order to show positivity, we will be using the Lemma 2.2 defined in section 2.3.

Since the initial values are non-negative and the growth functions on the right-hand side of system (4.5) are assumed to be sufficiently smooth in $\mathbb{R}_+^5$, we have the following result over the positivity of a solution of subsystems (4.5), if it exists:

**Lemma 4.1.** Any solution $(S(t, x), L(t, x), I(t, x), J(t, x), U(t, x)) \in [C([0, T_{max}) \times \bar{\Omega}) \cap C^{1,2}((0, T_{max}) \times \Omega)]^5$ of subsystems (4.5) is positive over $[0, T_{max})$, the largest interval over which the solution of subsystems (4.5) exists.
Proof of Lemma 4.1. For the positivity of S, L, I, J, let us define the supremum $t_1 = \sup \{t > 0, S(t,x) > 0, L(t,x) > 0, I(t,x) > 0, J(t,x) > 0, U(t,x) > 0 \forall x \in \Omega\}$, and suppose $t_1 < T_{max}$.

Using the definition of supremum, there does not exist $t \in [0, t_1)$ and $x \in \Omega$ such that any variable is equal to zero.

Among S, L, I, J, U, let us first consider S to be the variable such that there exists x with $S(t_1, x) = 0$, and let us define:

$$ \lambda(t, x) = \exp \left( \int_{0}^{t} \frac{\omega_{p} \nu U(a, x)}{N(a, x)} da + \mu_{p} t \right). $$

One has:

$$ \begin{aligned} \partial_t(\lambda(t,x)S(t,x)) &= \lambda(t,x) \left[ \partial_t S(t,x) + \frac{\omega_p \nu U}{N} S(t,x) + \mu_p S(t,x) \right] \\ &= \Lambda_p \lambda(t,x) \\ &> 0 \end{aligned} $$

where $\Lambda_p$ is obtained from the expression of $\partial_t S(t,x)$ in subsystems (4.5). By integration of this expression over $[0, t_1]$, we obtain:

$$ S(t_1, x)\lambda(t_1, x) - \phi_S(x)\lambda(0, x) > 0. $$

Since $\lambda(t_1, x) > 0$ and $\phi_S(x)\lambda(0, x) > 0$, we have $S(t_1, x) > 0$, which contradicts the hypothesis that $S(t_1, x) = 0$. Hence $S(t, x) > 0$ for all $(t, x) \in [0, t_1] \times \Omega$.

Also, the $\partial_t L$, $\partial_t I$ and $\partial_t J$ equations are linear scalar equations with forcing terms $\frac{\omega_p \nu U}{N} S$, $\theta_p L$ and $\alpha_p I$ that are positive over the interval $[0, t_1)$; hence $L(t_1, x)$, $I(t_1, x)$ and $J(t_1, x)$ are positive.

Hence, if $t_1 < T_{max}$ as we supposed, we must have, for some $x$, $U(t_1, x) = 0$.

As $I(t, x) > 0$ over $[0, t_1]$, the last equation of subsystems (4.5) then gives rise to the following inequality system:

$$ \left\{ \begin{array}{l} \partial_t U - \varepsilon \Delta U \ge -(\nu + \mu_U)U, \quad (t,x) \in (0, t_1] \times \Omega \\ \frac{\partial U}{\partial \eta} \ge 0 \quad \text{on } \partial \Omega, \\ U(0, x) \ge 0, \quad \forall x \in \Omega. \end{array} \right. \qquad (4.6) $$

Since $\nu + \mu_U$ is bounded in $[0, t_1] \times \Omega$, it follows from Lemma 2.2 that $U(t, x) > 0$ in $[0, t_1] \times \Omega$. This contradicts $U(t_1, x) = 0$, and hence the assumption $t_1 < T_{max}$. Therefore all variables are positive for all times smaller than $T_{max}$.

This achieves the proof. □

**Lemma 4.2.** Let $(S, L, I, J, U) \in C^0([0, T_{max}) \times \bar{\Omega}) \cap C^{1,2}([0, T_{max}) \times \Omega)^5$ be the solution of subsystems (4.5) with bounded initial conditions. Then this solution is bounded, $T_{max} = +\infty$, $0 < N(t, x) = S + L + I + J \le N_m$ and $0 < U(t, x) \le U_m$, where:

$$ N_m = \max \left\{ \frac{\Lambda_p}{\mu_p}, \|\phi\|_{\infty} \right\} \quad \text{and} \quad U_m = \max \left\{ \max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\} \frac{\gamma_p}{(\nu + \mu_U)}, \|\phi_U\|_{\infty} \right\}. $$
Proof of Lemma 4.2. Let $\phi(x) = \phi_S(x) + \phi_L(x) + \phi_I(x) + \phi_J(x)$. The N-dynamics then satisfy:

$$ \left\{ \begin{array}{l} \partial_t N \le \Lambda_p - \mu_p N, \quad (t, x) \in [0, \mathcal{T}_{\max}) \times \Omega \\ N(0, x) = \phi(x) \le \|\phi\|_{\infty} = \max_{x \in \Omega} \phi(x). \end{array} \right. \qquad (4.7) $$

upon which we build the upper-bounding system

$$ \left\{ \begin{array}{l} \frac{dW_1}{dt} = \Lambda_p - \mu_p W_1, \\ W_1(0) = \|\phi\|_{\infty}. \end{array} \right. \qquad (4.8) $$

The comparison principle then yields:

$$ N(t, x) \le W_1(t) = \frac{\Lambda_p}{\mu_p} + \left( \phi(x) - \frac{\Lambda_p}{\mu_p} \right) e^{-\mu_p t} $$

for all $(t, x) \in [0, \mathcal{T}_{\max}) \times \Omega$. Since $W_1(t) \le \max\{\frac{\Lambda_p}{\mu_p}, \|\phi\|_{\infty}\}$ for $t \in [0, \infty)$, one has that $N(t, x)$ stays bounded and is hence defined for all times:

$$ 0 < N(t, x) \le \max \left\{ \frac{\Lambda_p}{\mu_p}, \|\phi\|_{\infty} \right\} = N_m, \quad \forall (t, x) \in [0, \infty) \times \Omega. $$

This upper-bound also holds for $S$, $L$, $I$ and $J$ which are integral parts of $N$.

From the upper-bound on $I$ and the last equation of subsystems (4.5), we can deduce that:

$$ \left\{ \begin{array}{l} \partial_t U - \varepsilon \Delta U \le \max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\}\gamma_p - (\nu + \mu_U)U, \quad (t, x) \in [0, \mathcal{T}_{\max}) \times \Omega \\ U(0, x) = \phi_U(x) \le \|\phi_U\|_{\infty} = \max_{x \in \Omega} \phi_U(x), \\ \frac{\partial U}{\partial \eta} = 0, \quad x \in \partial\Omega, \end{array} \right. \qquad (4.9) $$

and the upper-bounding system:

$$ \left\{ \begin{array}{l} \frac{dW_2}{dt} = \max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\}\gamma_p - (\nu + \mu_U)W_2, \\ W_2(0) = \|\phi_U\|_{\infty}. \end{array} \right. \qquad (4.10) $$

It then follows that for all $(t, x) \in [0, \mathcal{T}_{\max}) \times \Omega$:

$$ U(t, x) \le W_2(t) = \max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\} \frac{\gamma_p}{(\nu + \mu_U)} + \left( \phi_U(x) - \max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\} \frac{\gamma_p}{(\nu + \mu_U)} \right) e^{-(\nu + \mu_U)t}. $$

Hence, as was done for $N$, we have:

$$ U(t, x) \le \max \left\{ \max\{\Lambda_p/\mu_p, \|\phi\|_{\infty}\} \frac{\gamma_p}{(\nu + \mu_U)}, \|\phi_U\|_{\infty} \right\} = U_m, \quad \forall (t, x) \in [0, \infty) \times \Omega. $$

This completes the proof

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

Now we prove existence and uniqueness of a local solution of subsystems (5). The idea is to write the system in the abstract form on a suitable Hilbert space and then to use the Hille-Yosida Theorem to conclude