Infection model for analyzing biological control of coffee rust using bacterial anti-fungal compounds
Jorge Arroyo Esquivel Fabio Sanchez Luis Barboza
q-bio.PE

Frustrated Kondo chains and glassy magnetic phases on the pyrochlore lattice
    Gia-Wei Chern
    December 30, 2023
============================================================================




Jorge Arroyo Esquivel[Corresponding author. <jorge.arroyoesquivel@ucr.ac.cr>][Present address: Universidad de Costa Rica. San Pedro de Montes de Oca, San José, Costa Rica, 11501.]

 Universidad de Costa Rica
   Escuela de Matemática
   San Pedro, San José, Costa Rica
 



Fabio Sanchez

 Universidad de Costa Rica
   Escuela de Matemática
   San Pedro, San José, Costa Rica




Luis A. Barboza

 Universidad de Costa Rica
   Escuela de Matemática
   San Pedro, San José, Costa Rica












Coffee rust is one of the main diseases that affect coffee plantations worldwide <cit.>. This causes an important economic impact in the coffee production industry in countries where coffee is an important part of the economy. A common method for combating this disease is using copper hydroxide as a fungicide, which can have damaging effects both on the coffee tree and on human health <cit.>. A novel method for biological control of coffee rust using bacteria has been proven to be an effective alternative to copper hydroxide fungicides as anti-fungal compounds <cit.>. In this paper, we develop and explore a spatial stochastic model for this interaction in a coffee plantation. We analyze equilibria for specific control strategies, as well as compute the basic reproductive number, ℛ_0, of individual coffee trees, conditions for local and global stability under specific conditions, parameter estimation of key parameters, as well as sensitivity analysis, and numerical experiments under local and global control strategies for key scenarios.


Keywords: Infection model, biological control, epidemiology, spatial model, coffee rust, parameter estimation, sensitivity analysis.




§ INTRODUCTION

Throughout the last two centuries, coffee (Coffea arabica) has become one of the most extensively produced commodities worldwide, consequently it has helped several developing countries to enter the global market <cit.>. During that time, coffee production has been affected by a complex network of diseases and pests <cit.>, from which the coffee rust, a leaf disease caused by the infection of the fungi Hemileia vastatrix, has become a major nuisance in the coffee production worldwide, causing serious defoliation in coffee trees, and as a result it causes losses in the crop yield <cit.>.

A typical method used for prevention is using fungicides based on copper compounds <cit.>. The use of these fungicides brings potential problems to the ecosystem through copper enrichment of the soil and health effects on the high consumption of copper <cit.>, as well as the increase in the population of other coffee pests such as the green scale (Coccus viridis) through the destruction of their parasite, the fungi Lecanicillium lecanii <cit.>.

As an alternative method of reducing coffee rust and other fungal pathogens in cash crops, several studies have explored the use of bacteria as a biological control method <cit.>. The strain of Bacillus thuringiensis B157 used as a biological control method has shown to be as effective as copper hydroxide fungicides <cit.>. These bacteria control the rust uredospores by blocking their germination tubes with antibiotic compounds, mainly lipopeptides such as iturine, surfactin, and fengicine <cit.>.

Typically, it is enough to use fungicides based on these compounds to control rust on coffee trees. However, these compounds have been shown to be volatile <cit.>, making them difficult to be effectively used as an independent fungicide. Moreover, considering the interaction between the bacteria and coffee rust in a coffee plantation, it is possible that one can determine the critical size of bacteria population, this in turn will reduce coffee rust on coffee trees.

Several models have been developed to describe the spread of crop diseases in a field <cit.> and some others have been developed to describe the dynamics between the biological control agents and their related pathogens <cit.>. However, these models do not propose an interaction between the pathogens and a control agent, in our case the bacteria, while the latter does not consider the spatial distribution of the interacting populations. In this article we develop a stochastic model that considers the interactions between bacteria and coffee rust populations and consider the effect of the spatial distribution of the disease and its spread on the overall dynamics of the disease.

The rest of the paper is organized as follows. In Section 2, we describe the dynamics between the coffee tree, coffee rust, and bacteria. In Section 3, we develop a stochastic epidemic model that illustrates the spread of coffee rust through a coffee plantation, as well as the interaction between the bacteria and coffee rust. Section 4, shows both analytic results and numerical experiments, including theoretical stability analysis for our model, simulations of the coffee rust and bacteria populations, and parameter estimation with its corresponding uncertainty analysis. Lastly, we discuss model results and their biological significance in Section 5.



§ COFFEE RUST DYNAMICS



 §.§ Coffee


Coffee (Coffea arabica) is a perennial tree of the Rubiaceae family originated in the southeast of the Arabic peninsula <cit.>. Several types and mutations of this tree have been introduced in several regions of Latin America, Southeast Asia, and Africa since the nineteenth century, from which the berry is used as a food product <cit.>.

Typical conditions for coffee to be grown depends mostly on altitude and rainfall. Coffee is a highland crop, usually growing between 1000 and 2000 meters of altitude <cit.>. In addition, it needs a proper amount of rainfall to grow, where more than four months without rain can cause significant damage to crops <cit.>.

The yearly rainfall distribution properly defines how the tree produces berries. If there is a single rainfall period in the year, the coffee will only have a single period of harvest, while if there are two periods, it will have two periods <cit.>. For example, in the case of Costa Rica, where the rainy season goes from May to October, there is a single harvest season of coffee ranging from June to November, the exact moment depending on the altitude and composition of the soil at which the coffee is being grown <cit.> (see Figure <ref>).





 §.§ Coffee rust


Coffee trees are affected by several diseases, which may cause significant yield losses in coffee plantations. One of these is the coffee rust, caused by the fungus Hemileia vastatrix, which is an obligate parasite that uses coffee as its main host. It was first recorded in 1869 in Ceylon (currently Sri Lanka), where an outbreak suddenly reduced the production of coffee in the island by around 90% of the total production, causing the coffee industry to cease and made them cultivate tea instead, which is not as profitable <cit.>. This disease is usually recognizable by yellow spots appearing in the blade of the leaves as shown in Figure <ref>.



The coffee rust disease is caused by the infection of uredospores of H. vastatrix in the leaves. These spores are mainly produced by asexual reproduction of the hyphae <cit.>. For a leaf to be infected, it requires a minimum number of spores to properly perform germination, that number being between 15 and 20 spores <cit.>. The germination of uredospores depends of several environmental and physiological variables. The main factors influencing the germination process are rainfall density, temperature, and age on both the uredospores and leaves <cit.>.

For the environmental factors, in low rainfall densities, the germination of uredospores and growth on their mycelium is slower in comparison to the growth in more humid conditions. In cases of extreme dryness, the mycelium of H. vastatrix is unable to grow <cit.>. The optimal temperature for uredospore germination, mycelium growth, and lesion formation in the leaves is 22^∘C. However, with a previous low-temperature stimulation, followed by an increased temperature, the uredospore germination is faster compared to a constant exposition to 22^∘C. This is the case when the plantation receives rainfall during the night, where the temperature is lower, and then warms up during daylight.

In the case of physiological factors, the germination process of uredospores is both slower and have a lower infection rate in the more mature leaves compared to younger leaves <cit.>. The age of the uredospores plays a role, since the effectiveness of germination of uredospores decrease exponentially.

The uredospores are mostly transmitted inside a single plantation through two main mechanisms, rain splash and wind <cit.>. However, there is a negative correlation between uredospore age and its ease of being transported through the wind. As it becomes easier to transport the uredospores through the wind, the less effective they are for germination. Hence, we can argue that the principal method for local dispersal is rainfall.

In Costa Rica, an epidemic of coffee rust occurred between 2008 and 2013. During this period, coffee production and price decreased by 16% and 55%, respectively. This had an important economic impact on small producers and coffee gatherers <cit.>.

The typical method for controlling coffee rust in coffee plantations is by spraying the plantation with several fungicides. In Costa Rica, the fungicide is normally applied at least twice a year, once at the beginning of the rainy season around May, and during September, in the middle of the rainy season as a precautionary measure <cit.>. For stricter control measures, more applications can be done.



 §.§ Interaction between bacteria and coffee rust


The use of bacteria as a biological control method can lead to the reduction of coffee rust in coffee plantations. It has been shown that biological control through the bacteria Bacillus thuringiensis results in a process as effective as the application of fungicides of copper hydroxide compounds, which are the standard method of control of coffee rust and has the potential to be harmful to humans <cit.>.

This interaction is produced by the effect of several antibiotic compounds produced by the bacteria. In the case of B. thuringiensis, these compounds are mostly iturine, surfactine, and fengicine <cit.>. They interact with the germination pores of the uredospores of H. vastatrix, blocking the germination tubes inside the uredospore, and therefore stopping the growth of the mycelium.

In the following section we develop a stochastic model that describes this process, allowing us to make predictions about more effective control measures of coffee rust in coffee plantations.



§ MODEL

In this section we introduce a nonlinear differential equation model that describes the spread of coffee rust through a coffee plantation and its interaction with the bacteria B. thuringiensis. Previous models have studied a similar phenomenon <cit.>. Our model is moderately based on the two-population model presented in <cit.>. The latter is constructed as follows, let C be the consumer population with carrying capacity K, which consumes the resource population with an effective rate a and has a mortality rate m. Let R be the resource population, which grows at a rate r and is consumed by the consumer population with a conversion rate c. Then the interaction between these two populations is presented by the following system of differential equations:


    dC/dt =    a(K-C/K)CR-mC,
    dR/dt =    rR(1-R)-a(K-C/K)CR.


Within our model we consider the consumer's growth (in our case the bacteria population) is not affected by the population of the source (in our case the coffee rust population), which means that a=1 and R does not affect the logistic growth factor of C. Moreover, the source's growth is similar to the growth of the consumer (this means that the growth of the source is also logistic). We also consider the spatial distribution of the coffee rust (source). 

We consider a coffee plantation of n× m coffee trees, distributed in a rectangular array, where the tree in the i-th row and j-th column occupies the position {i,j}. We denote by H_{i,j} and B_{i,j} the size of the populations of coffee rust spores H. vastatrix and the bacteria B. thuringiensis B157 on tree {i,j}, respectively.

For each i∈{1,…,n},j∈{1,…,m}, both H_{i,j} and B_{i,j} are assumed to grow logistically, where b and h are the growth rates and K_B and K_H are the carrying capacities of the bacteria and coffee rust, respectively. Given that coffee rust appears primarily during rainy season, we will consider the carrying capacities constant. We assume that B_{i,j} has a natural death rate d and the bacteria limits the growth of coffee rust with a constant conversion rate γ.

Because of the evidence presented in Section 2.2, in our model, the dispersal of coffee rust through wind is neglected. Therefore, the principal method of dispersal is through rainfall, which arrives to the neighbor trees at rate α and is removed from their original tree by rainfall at rate β. We suppose that the bacteria are being irrigated into the plantation at rate μ_{i,j}.

The model is then described by the system of nonlinear differential of equations:







    dB_{i,j}/dt =    bB_{i,j}(1-B_{i,j}/K_B)+μ_{i,j}(t)-dB_{i,j},
    dH_{i,j}/dt =    hH_{i,j}(1-H_{i,j}/K_H)+α I_{i,j}-(γ B_{i,j}+β) H_{i,j}.


The I_{i,j} function is defined as a random variable dependent of the coffee tree's neighbors, which correspond to the trees in the cardinal and intermediate directions of the {i,j}-th tree. This function is the number of spores that can go from any neighbor coffee tree to the {i,j}-th tree. I_{i,j} is then defined as:

    I_{i,j} = ∑_k=-1^1∑_l=-1^1 H_{i+k,j+l}(t)·τ_{i+k,j+l}(k,l,t)-H_{i,j}(t)·τ_{i,j}(0,0,t),


where H_{i,j}=0 if i∉{1,…,n}, j∉{1,…,m}. Here, τ_{ij}(k,l,t) is a random variable that represents the proportion of coffee rust spores H_{i,j} that is passed from tree {i,j} to tree {i+k,j+l}.

To define the value of τ_{i,j}, we will assume that the only method of local dispersal of coffee rust is through rain splash. Therefore, we will base the probability of the random variables τ_{i,j} on the dispersal of plant pathogens by the rain splash model provided in <cit.>. In the model, spores splash from a point source and are displaced over a flat ground surface. Each spore is suspended in a water layer, which can be hit by a water splash with probability λ and the spore can be lost in the process and stuck on the ground with probability .

If the initial distribution of the spores is δ(x), then the probability per unit length of finding a spore at position x at time t is:

    P(t,x) = e^-λ tδ(x)+∑_i=1^∞(λ t)^ie^-λ t/i!^iD(x)^⋆ i,


where D(x)^⋆ i is the normal distribution convoluted by itself i times, and the normal distribution is given by:

    D(x) = 1/√(2πσ^2)exp(-x^2/2σ^2).


The model describes how at a point in space x, the probability of the spores that were initially at point x that remain at that point decays in an exponential manner and is compensated by the spores from other points that reach this point through rain splash following a Poisson process.

We make additional assumptions to the model. First, we make a discrete version of the model, considering only the spread from a tree to its neighbors, separated by a distance (k,l)_{i,j}, where ·_{i,j} is a distance function, not necessarily a norm function.


Since the distance between two trees is smaller if they are in the same row, independent of the column, we have that 

    (k,0)_{i,j}≤(0,l)_{i,j}≤(k,l)_{i,j}

for k,l∈{± 1}.


Additionally, if a rust spore is in one tree at time t, we assume that in time t+Δ t there will be at most one splash, where Δ t is a small change in the time scale (for example, if t is scaled in hours, Δ t can correspond to a minute). Furthermore, since we only analyze the spread of coffee rust in time t+Δ t, we do not consider the exponential decay value, which begins at time 0. Hence, we only consider the case i = 1 in Equation <ref>.

In the case of no rain, the dispersal is negligible. To consider the chance that there is no rain, assume there is a probability P_λ(t) at time t of raining. Given these assumptions, our random variables τ_{i,j} satisfy that

    τ_{i,j}(k,l,t)∼Bernoulli(P_{i,j}(k,l,t)),

where:

    P_{i,j}(k,l,t) = ψλ(t)H_{i,j}(t)D((k,l)_{i,j}),


and ψ is a proportionality constant, λ(t)∼Bernoulli(P_λ(t)) is the probability that it is raining at time t, and D is the normal distribution with mean 0 and variance σ^2 as in Equation <ref>. Note that P_λ depends on time, which is important since rainfall amount is dependent on the period of the year, as shown in Figure <ref>.

For the values of μ_{i,j}, each value will define a specific control strategy. First, we define a global control strategy where the whole plantation is irrigated at the same rate, this means, for all {i,j}:

    μ_{i,j}≡μ


for some constant μ. This strategy works as a preemptive measure, preparing the coffee trees to counter coffee rust. Then the local control strategy, which defines a rate that is proportional to the amount of coffee rust in the tree itself and its neighbors, for all {i,j}, is given by:

    μ_{i,j}(t) = ρ H_{i,j}(t)+δ(∑_k=-1^1∑_l=-1^1 H_{i+k,j+l}(t)-H_{i,j}(t)),


for some constants ρ,δ defined as the intensity of spraying on a tree and its neighbors, respectively. This control strategy works as a reactive strategy, where the trees are irrigated depending on its neighbors’ infectious status.



§ RESULTS

In this section we refer to the mathematical and numerical results of our study. First, we show the theoretical equilibria of the model under necessary conditions using the basic reproductive number and prove global stability of the disease-free equilibrium, as well as the local stability of the endemic equilibrium in a specific plantation size. Moreover, we perform several numerical experiments under local and global strategies for the control of coffee rust on coffee plantations to study the temporal and spatial dynamics of our model, as well as parameter estimation of the model under certain conditions to analyze the sensitivity of the parameters in the overall behavior of the system.



 §.§ Bacteria and coffee rust equilibrium

Given the recursive nature of the model, we only study the following case. Consider a global control strategy as shown in Equation <ref>, and take P_λ≡ 1 as the probability for λ in Equation <ref>, which implies the plantation has a constant amount of rainfall. More so, we assume that the difference between the distances of neighbor trees from a specific tree are negligible, this means |D((k_1,l_1)_{i,j})-D((k_2,l_2)_{i,j})|< for small enough  and k_i,l_i∈{±1,0}.

Since μ_{i,j} is constant, we get that:

    dB_{i,j}/dt = dB_1,1/dt


for all i∈{1,…,n},j∈{1,…,m}. Then, for all i,j, if b>d, the equilibrium point of the bacteria population is the following:

    B_{i,j}^* = B = 1/2bK_B(b-d)(1+√(1+4bμ/K_B(b-d)^2)).


We only consider positive equilibria; hence it is clear from Equation <ref> that this equilibrium exists whenever the birth rate of bacteria is bigger than its death rate.

The bacteria equilibrium shown in Equation <ref> exists if and only if b> d. 

For the coffee rust population, since P_λ≡ 1, and the difference in tree distances D((k,l)_{i,j}), are negligible, the probability of dispersal is homogeneous throughout the coffee plantation. Therefore, we can conclude that in equilibrium the value of I_{i,j} behaves in a similar manner to a constant Π for all i,j, where Π corresponds to the proportion of coffee rust spores each tree has if the only dynamic mechanism considered is dispersal (this means, growth and competition between populations is not considered). Hence, the equilibria of each coffee rust population is given by:

    H_{i,j}^* = K_H/γ Bh(h+Πα-β-γ B).


Define F_0 = h+Πα-β as the net growth rate of the coffee rust spore population on each tree. Then, if B>0 (this means, there is bacteria in equilibrium), the basic reproductive number of coffee rust is defined as the average number of infected trees the population of coffee rust in a single tree can produce and is given by:

    R_0 = F_0/γ B,


which turns Equation <ref> into:

    H_{i,j}^* = K_H/h(R_0-1).


From this equation we get the following theorem.

Given R_0 as in Equation <ref>, then the endemic coffee rust equilibrium exists if and only if R_0>1.


Note that, by taking H_{i,j}≡0 for all {i,j}, then it is shown that the disease-free equilibrium always exists.



 §.§ Global stability of disease-free equilibrium

In this section we prove that, for any plantation, the disease-free equilibrium is globally asymptotically stable. To do this, we use the method described on the main result in <cit.>:


Let X∈ℝ^n be the uninfected individuals and Z∈ℝ^m the infected individuals in the system such that the System <ref> is rewritten as:


    dX/dt =    F(X,Z),
    dZ/dt =    G(X,Z).


Then, if the three following conditions are met:


    
  * ℛ_0<1
    
  * For dX/dt = F(X,0), the disease-free equilibrium X^* is globally asymptotically stable.
    
  * G(X,Z) = AZ-Ĝ(X,Z), where A = G_Z(X^*,0) and Ĝ(X,Z)≥ 0 for all (X,Z) where the model makes sense.


Then the disease-free equilibrium is globally asymptotically stable.


In this section we will prove that conditions (2) and (3) are met by our model. To do this, starting from the {1,1}-th tree, relabel the trees from the first column down to the {n,1}-th tree as 1,…,n. Then, starting from the {1,2}-th tree relabel the trees from the second column down to the {n,2}-th tree as n+1,…,2n and repeat this process up to the m-th column. Then we write X=(B_1,…,B_mn) and Z=(H_1,…,H_mn).

To prove (2), note that for all i, by repeatedly derivating we get that


    d^kB_i/dt^k=0
 

in the equilibrium (when B_i=B^* with B^* as in Equation <ref>). Hence, we have that in the equilibrium, B_i becomes constant. Therefore B_i→B^* when t→∞ and the equilibrium X^* is globally asymptotically stable.

To prove (3), given that we consider only the asymptotic behavior, we can consider that I_{i,j} =Π and B_i= B_i^*. Note that for all i, the i-th entry of the vector G equals


    G(X,Z)_i=hH_i(1-H_i/K_H)+αΠ∑_ neighbor_iH_j-(β+B_i)H_i,


where j∈neighbor_i if the following condition is met:


    i-j∈{1,n-1,n+1} if i≢0,1 n
    i-j∈{n-1,n+1} or i-j=-1 if i≡1 n
    i-j∈{n-1,n+1} or i-j=1 if i≡0 n
    .


Then, it follows that:


    (A)_i,j =
    h-(β+B_i^*) if i=j
    αΠ if j∈neighbor_i
    
    0 if else.


Therefore


    (AZ)_i = hH_i+αΠ+∑_ neighbor_iH_j-(β+B_i^*)H_i.


Take


    Ĝ(X,Z)_i = hH_i/K_H


Therefore, we have that the system satisfies condition (3). This lets us conclude the following theorem:


If ℛ_0<1, the disease-free equilibrium is globally asymptotically stable.




 §.§ Local stability of endemic equilibrium in a 1× 2 system

The structure of the model gives an additional difficulty when analyzing the stability of the endemic equilibrium when ℛ_0>1. However, the following result can be derived in a 1×2 system.


If (n,m)∈{(1,2),(2,1)} and ℛ_0>1, then if


    1<2/γ B^*


then the endemic equilibrium is locally asymptotically stable. Otherwise, if


    1>2/γ B^*


then the endemic equilibrium is unstable.


To prove this result, we will consider that B_i=B^* for i∈{(1,2)}. Then, the linearization of the H_i values comes as:


    J=([ h(1-2H^*/K_H)-β-γ B^*                    αΠ;                    αΠ h(1-2H^*/K_H)-β-γ B^* ]).


This matrix has eigenvalues


    h±αΠ-(β+γ B^*+2hH^*/K_H).


Since all parameters are positive, notice that the biggest eigenvalue is found when ±αΠ is positive. Notice that this eigenvalue is negative when


    h+αΠ<β+γ B^*+2hH^*/K_H.


If we order the terms to get ℛ_0 as in Equation <ref> we get that:


    ℛ_0<1+2hH^*/γ K_HB^*.


By making the substitution of H^* as in Equation <ref>, we get that:


    ℛ_0(1-2/γ B^*)<1-2/γ B^*.


Which is equivalent to ℛ_0>1 if 1-2/γ B^*<0. Otherwise if 1-2/γ B^*>0, then we get that ℛ_0<1, which is a contradiction, therefore both eigenvalues are positive, and the endemic equilibrium is unstable.



 §.§ Numerical simulations

Numerical experiments are performed using NetLogo <cit.>. In each simulation, for simplicity of computation and visualization we will consider a 11× 18 plantation, which corresponds to approximately 4% of the recommended amount of coffee trees in a hectare <cit.>. The simulation display in Netlogo is shown in Figure <ref>. Model parameters are summarized in Table <ref>.






For all simulations, we will assume that initially 50% of the plantation is infected, which corresponds to a value like that reported in <cit.>, with different initial values of coffee rust for each {i,j}, as shown in Figure <ref>. We will analyze both local and global control strategies in different settings (high/low growth rate of coffee rust, high/low dispersal of coffee rust, high/low conversion rate of bacteria).

The time scale of the simulations is in hours. We explore the behavior of the system in a small-time interval, hence, the duration of each simulation will be one week (this means t∈ [0,168]). In each simulation we analyze the proportion of leaves infected by coffee rust (with the assumption that on average all trees have the same amount of leaves) and the proportion of trees being infected by coffee rust. An advanced state of this simulation in Netlogo (this means, after 72 hours of execution) is seen in Figure <ref>. Figures <ref>, <ref>, and <ref> show time series of simulations where global control is applied, while Figures <ref>, <ref>,  <ref>, and <ref> show time series of simulations where local control is applied. We also calculate the basic reproductive number ℛ_0 on each simulation using the Maximum Likelihood estimation method presented in the R0 package of R <cit.>. This method provides a workaround to the difficulty of estimating ℛ_0 in an analytical manner.



Figure <ref> shows a scenario where the coffee rust growth rate is small, but its dispersal rate and probability of dispersal are high. These values are reflected by how strongly coffee rust is spread throughout the first days of the simulation, but the inhibition produced by bacteria with a high conversion rate gets into a breaking point, where the proportion of infected trees decrease rapidly. This is also accompanied by a slow decrease of coffee rust through the plantation. This suggests that although the dispersal is effective, the amount of coffee rust spores in the original trees decays, and under some threshold value, the coffee rust population is small enough that it is no longer able to support itself in a particular tree, therefore reducing the proportion of infected trees considerably.


In Figure <ref> we see a different behavior. This is explained by the high coffee rust growth rate, which compensates for the loss obtained by the conversion produced by bacteria, avoiding the coffee rust on each individual tree to pass through a critical point where coffee rust spores can no longer sustain itself on the leaves. Also notice that for small values of ψ (overall probability of dispersal by rain splash) makes the spread of the coffee rust spores behave in stochastic jumps, rather than a continuous function. 


Figure <ref> showcases how the coffee rust decaying rate is affected primarily by the conversion rate. Here, both coffee rust growth rate and conversion rate are small, hence the main source of coffee rust decay is the emigration rate (β). Since the values of ψ and immigration rate are both high, the coffee rust spores propagate through the plantation rapidly. Notice that in the end of the simulation the proportion of infected trees starts to go down, possibly as a delayed decay that happened after coffee rust reaches its threshold.


In Figure <ref>, the coffee rust growth rate is in the middle of the considered range, with a low emigration rate. This appears to be enough to compensate the high conversion rate, which makes the coffee rust decay. Furthermore, the high immigration rate (α) and overall probability of dispersal by rain splash (ψ) help the coffee rust spores reach most of the plantation.


Figure <ref> shows how the trend of the previous simulations does not appear when there are rainless periods. In this scenario the proportion of infected trees remains relatively constant throughout several days and suddenly decays given the small amount of coffee rust spores. This happens even if the value of ψ is high, which is related to an increased dispersal of coffee rust spores. It is important to note that P_λ<1 is a more realistic scenario for ranges of time like weeks, where the rainfall is not constant throughout the day.


Figure <ref> shows a scenario where the coffee rust growth rate is high, which avoids the compensation of other values that inhibit the coffee rust's growth, such as emigration rate and conversion rate. Note that in this scenario, there was more rain than in Figure <ref>, which can be seen as having more opportunities of dispersal (hence more jumps in the graph), even though the value of ψ was smaller. In this scenario coffee rust was successful and managed to get through the control of the bacteria, suggesting that with a high growth rate, coffee rust could not be controlled by this method.



In the case of Figure <ref>, both the coffee rust growth rate and the conversion rate of bacteria are high. In this scenario coffee rust eventually takes over the plantation, and reaches its carrying capacity in two days. This suggests that a high growth rate of coffee rust is not feasible, since no conversion rate from the bacteria seems to be able to compensate the rate at which coffee rust grows on each tree.

Our simulations show that when the plantation is at coffee rust-free equilibrium, that is, H_{i,j}=0 for all {i,j}, it is indicative that a coffee rust threshold exists where the coffee rust incidence on each tree can sustain itself. This suggests there is a minimum viable population of coffee rust in the coffee plantation.





We also performed experiments for high values of h and γ under values of the control parameters above our defined framework and with different initial conditions of coffee rust in the plantation. To do this, we let h = 0.15,K_H = 0.284,α = 0.9,β = 0.1,γ = 0.5, ψ = 20 for all simulations. For global control, we take μ = 2.5 and for local control take ρ = 7,δ = 7. The results for 50% of the plantation initially infected are shown in Figure <ref>, while the results for 10% of the plantation initially infected are shown in Figure <ref>. These results show that under these scenarios the model is not sensitive to the coffee rust initial conditions, and both strategies are effective to control coffee rust. However, in the case of global control, for both sets of initial conditions presented, the coffee rust went to zero, which suggests that in these scenarios global control was faster at reducing coffee rust in comparison to local control.



 §.§ Parameter estimation and uncertainty effects


In addition to the simulations, we analyzed the estimation of parameters in the model by estimating their probability density function with the Approximate Bayesian Calculation rejection algorithm in a similar way as the one developed in <cit.>. We analyzed the variables in R using the RNetLogo package to connect NetLogo with R and the EasyABC package to perform the ABC sequencing. We estimated the parameters of both global and local control strategies using the parameter values α=0.5,β=0.8,h=0.07,γ=0.3,ψ=1,μ = 0.6,ρ=0.6,δ=0.3. To sample the Approximate Bayesian Calculation, we implemented the Beaumont algorithm using the EasyABC package, based on <cit.>, with a prior uniform distribution on the specific range of each parameter as depicted in Table <ref>.

To analyze the uncertainty effect of each parameter, we report the mean, standard deviation, and 95% confidence interval of each of the accepted samples using the EasyABC package. We used the following statistic as an accepting measure for the sampling process:


    Λ(Ĥ) = √(∑_t=1^T∑_i=1^n∑_j=1^m(H_{i,j}(t)-H_{i,j}(t))^2),


where Ĥ is the matrix of simulated values of coffee rust with entries H_{i,j}(t) and H_{i,j}(t) is the original simulation keeping all the parameters fixed. The expected value of this summary statistic corresponds to the mean of 1000 simulations with the same parameter values as H_{i,j}(t) and the tolerance values correspond to σ̂ and 2σ̂, where σ̂ corresponds to the standard deviation of 1000 simulations with the same parameter values as H_{i,j}(t). The results of these simulations are found in Table <ref>.



These results let us understand how sensitive the parameters are when they are estimated and also the degree of identifiability of the parameters by analyzing how close the true values are with respect to the center of each corresponding confidence interval. Notice that in most cases the true value lies within a relative distance to this center of less than 10%, which tells us that these parameters can be successfully recovered from simulations. This can be appreciated in particular with the β and γ parameters, which have a low variance relative to their mean, which means that those parameters are not very sensitive to sample variation. This also let us notice that both control strategies are relatively robust, and the estimation process of their parameters with a real data set is feasible.

The only parameter which failed to be recovered is the h parameter, which could be explained by the possibility that for a relatively small value of h, the influence of other parameters to the overall behaviour of the system could be more important than the effect of h, which could be explained by the high variance relative to the mean values and the relatively large confidence intervals.



§ DISCUSSION

We have analyzed a stochastic model that describes the dispersal of coffee rust in a coffee plantation and its interaction with antagonistic bacteria. When introducing bacteria into the system, two strategies were explored, local and global control of coffee rust. With global control, we found at equilibrium that the proportion of infected tree leaves 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. Hence, it is viable to reduce coffee rust in coffee plantations using bacteria, which results in the prevention of problems related to the use of fungicides, such as damage to the soil health <cit.> and the reduced resistance of coffee trees to other plagues such as the green scale <cit.>.

With local control, since the irrigation rate of bacteria is not homogeneous throughout the coffee plantation, when analyzing equilibria of the model, there may be cases where we have different basic reproductive numbers for each tree. This could lead to the possibility of having ℛ_0>1 in some areas of the plantation, while also ℛ_0<1 in others, which can lead to an unstable behavior near its respective equilibrium point. Albeit an interesting possibility, we presume that additional hypotheses need to be made to explore this scenario.

Furthermore, qualitatively the local and global control strategies do not show a significant difference. In fact, the effect of the irrigation rates of bacteria (μ,ρ,δ) for the values presented in the framework of Table <ref> did not differ greatly, and the results mainly depended in a balance between the coffee rust growth rate (h) and the conversion rate of bacteria to coffee rust (γ). In other words, the coffee rust growth inside a tree depends mostly in its natural growth and the effect the bacteria population has on this population.

For small values of h (coffee rust growth rate), the coffee rust spores decrease rapidly, and increases when the growth rate decreases and the bacteria conversion rate increases. However, when the coffee rust growth rate is high, the coffee rust spores grow steadily even if the bacteria conversion rate is close to 1, therefore the bacteria can reduce the coffee rust of a completely infected leaf by itself. Moreover, the immigration and emigration rates only determined how fast the coffee rust dispersed through the plantation and had little to no effect on the infected leaves in the plantation.

There are several challenges when applying this model in a real-life context. One comes from the environmental dependence of the unknown parameters. Several environmental variables have to be put into consideration, such as amount of rainfall, temperature, and altitude to get an explicit value of the parameters h,α,β <cit.>. Although there is no reported value of the conversion rate of bacteria to coffee rust changing by environmental factors, this possibility cannot be discarded. Another situation comes from calculating the R_0 value, which faces the problem of determining the amount of coffee rust in equilibrium if the only factor considered is dispersal (Π), which depends on plantation size and the probability P_{i,j}. Therefore, these parameters should also be calculated in a case by case basis. However, this model defines a theoretical framework that can be used for developing biological control solutions to coffee producers and whose parameters have been shown to be possible to determine. The equilibrium could be even further analyzed if a method for empirically determining the value of the parameter Π is developed.

In conclusion, although we found no significant difference between the effect of global control and local control strategies, there are important differences in the application of each strategy. The global control strategy works like a safety net for the coffee plantation, where most of the plantation is being protected in a simple manner, however its application cost is higher than the application of the local control strategy, especially in bigger plantations, since it would require a larger labor force applying the solution throughout the coffee plantation. 

On the other hand, although the local control strategy is more cost-effective, since it works as a reactive strategy, the plantation has to be in constant monitoring to attack coffee rust as soon as it appears. This requires workers to diligently check the coffee plantation to detect coffee rust in time, which incurs in higher production costs. Although, we can argue these costs are lower compared to a frequent application of the bacteria solution in a global control strategy. Therefore, under this scenario, a global strategy is ideal for smaller plantations, while trying to optimize costs in larger coffee plantations using local strategy.



§ ACKNOWLEDGEMENTS
 
The authors acknowledge the feedback and orientation given by the attendees of the XIX National Meeting of Mathematical Biology at Universidad de Colima, Mexico, and Noel Molina and the agronomists met at the Field Trip of Coffee Producers held in Alajuela, in October 2017. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

99
alvarado2007
	M. Alvarado and S. Alvarado,
	Cultivo y beneficiado del café,
	Editorial Universidad Estatal a Distancia, San José, 2007.

avelino2015
	J. Avelino, M. Cristancho, S. Georgiou, P. Imbach, L. Aguilar, G. Bornemann, P. Laderach, F. Anzueto, A. Hruska, and C. Morales,
	The coffee rust crises in Colombia and Central America (2008–2013): impacts, plausible causes and proposed solutions,
	Food Sec, 7 (2015), 303–321.

beaumont2009
    M. A. Beaumont, J. M. Cornuet, J. M. Marin, and C. P. Robert,
    Adaptative approximate Bayesian computation,
    Biometrika, 96 (2009), 983–990.

brown2018
	G. D. Brown, A. T. Porter, J. J. Olsen, and J. A. Hinman,
	Approximate Bayesian computation for spatial SEIR(S) epidemic models,
	Spatial Spatio-temporal Epidemiol, 24 (2018), 27–37.

carvalho2011
	C. R. Carvalho, R. C. Fernandez, G. Mendes, R. W. Barreto, and H. C. Evans,
	Cryptosexuality and the Genetic Diversity Paradox in Coffee Rust, Hemileia vastatrix,
	PLOS ONE, 6 (2011).

ccc2002 Castillo-Chavez, C., Feng, Z., Huang, W. (2002). On the computation of R_0 and its role on global stability. Mathematical Approaches for Emerging and Reemerging Infectious Diseases, Springer. 

chatto2002
	J. Chattopadhayay, R. R. Sarkar, and S. Mandal,
	Toxin-producing plankton may act as a biological control for planktonic blooms: field study and mathematical modelling,
	J Theor Bol, 215 (2002), 333–344.	

clifford1985
	M. Clifford and K. Willson,
	Coffee: Botany, Biochemistry and Production of Beans and Beverage,
	The AVI Publishing Company, New York, 1985.

collinge2016
	D. Collinge,
	Plant Pathogen Resistance Biotechnology,
	John Wiley & Sons, New Jersey, 2016.

cressey2013
	D. Cressey,
	Coffee rust regains foothold,
	Nature, 493 (2013), 587.
	
dorighello2015
	D. V. Dorighello, W. Bettiol, N. Borlina, and R. M. Boas de Campos Leite,
	Controlling Asian soybean rust (Phakopsora pachyrhizi) with Bacillus spp. and coffee oil,
	Crop Prot, 67 (2015), 59–65.

haddad2009
	F. Haddad, L. Maffia, E. Mizubuti, and H. Texeira,
	Biological control of coffee rust by antagonistic bacteria under field conditions in Brazil,
	Biol Control, 49 (2009), 114–119.
	
haddadf2013
	F. Haddad, R. Saraiva, E. Mizubuti, R. Romeiro, and L. Maffia,
	Antifungal
compounds as a mechanism to control hemileia vastatrix by antagonistic bacteria,
	Trop Plant Pathol, 38 (2013), 398–405.

icafe2011
	Guía Técnica para el Cultivo del Café, Guide of Costa Rican Research Center of Coffee,
	2011. Available from <http://www.icafe.cr/wp-content/uploads/cicafe/documentos/GUIA-TECNICA-V10.pdf>.

imn2017
	Boletín Metereológico Mensual Mayo 2017, Report of Costa Rican National Institute of Metereology,
	2017. Available from <https://www.imn.ac.cr/boletin-meteorologico>.

jacksond2012
	D. Jackson, J. Skillman, and J. Vandermeer,
	Indirect biological control of the coffee leaf rust, Hemileia vastatrix, by the entomogenous fungus Lecanicillium lecanii in a complex coffee agroecosystem,
	Biol Control, 61 (2012), 89–97.
	
jie2013
	H. Jie, H. Feng, L. Huang, and Z. Kang,
	Biological control of wheat stripe rust by an endophytic Bacillus subtilis strain E1R-j in greenhouse and field trials,
	Crop Prot, 43 (2013), 201–206.
	
kranz1990
	J. Krans,
	Epidemics of Plant Diseases: Mathematical Analysis and Modelling,
	Springer, Berlin, 1990.

filho2010
	R. Lanna Filho, H. Ferro, and R. Pinho,
	Controle biológico mediado por Bacillus subtilis,
	Revista Trópica, 4 (2010), 12–20.

li2012
	X. Li, X. Ding, L. Xia, Y. Sun, C. Yuan, and J. Yin,
	Proteomic Analysis of Bacillus thuringiensis Strain 4.0718 at Different Growth Phases,
	Sci World J, 2012 (2012).

little1993
	I. Little, R. Cooper, W. Corden, and S. Rajapatriana,
	Boom, Crisis, and Adjustment: The Macroeconomic Experience of Developing Countries,
	Oxford University Press, Oxford, 1993.

loland2004
	J. Loland and B. Singh,
	Copper contamination of soil and vegetation in coffee orchards after long-term use of Cu fungicides,
	Nutr Cycl Agroecosys, 63 (2004), 203–211.
	
mills1996
	N. J. Mills and W. M. Getz,
	Modelling the biological control of insect pests: a review of host-parasitoid models,
	Ecol Model, 92 (1996), 121–143.

netlogo
	U. Wilensky,
	Netlogo, Center for Connected Learning and Computer-Based Modeling, Northwestern University.
	Available from <http://ccl.northwestern.edu/netlogo/>.
	
nutmanf1963
	F. Nutman and F. Roberts,
	Studies on the Biology of Hemileia vastatrix Berk. & Br.,
	Trans Brit Mycol Soc, 46 (1963), 27–44.
	
obadia2012
    T. Obadia, R. Haneef, and P.Y. Boelle,
    The R0 package: a toolbox to estimate reproduction numbers for epidemic outbreaks,
    BMC Med Inform Decis Mak, 12 (2012).

pendergrast2010
	M. Pendergrast,
	Uncommon Grounds: The History of Coffee and How it Transformed our World,
	Basic Books, New York, 2010.

pielaata1998
	A. Pielaat and F. Bosch,
	A model for dispersal of plant pathogens by rainsplash,
	Math Med Biol, 15 (1998), 117–134.

rafikov2008
	M. Rafikov, J. M. Balthazar, H. F. von Bremen,
	Mathematical modeling and control of population systems: Applications in biological pest control,
	Appl Math Comput, 200 (2008), 557–573.

rayner1961
	R. Rayner,
	Germination and Penetration studies on coffee rust (Hemileia vastatrix B. & Br.),
	Ann Appl Biol, 49 (1961), 497–505.

robinson1964
	J. B. Robinson,
	A handbook on Arabica coffee in Tanganyika,
	Tanganyika Coffee Board, Dar es Salaam, 1964.

scholl2012
	M. Schollenberger, K. Felczak, and K. Borowska,
	The effect of epiphytic bacteria on selected rust fungi,
	6 th International Plant Protection Symposium at University of Debrecen, 50 (2012), 98–100.

smith1998
    M.C. Smith, J. Holt, L. Kenyon, and C. Foot,
    Quantitative epidemiology of Banana Bunchy Top Virus Disease and its control,
    Plant Pathol, 47 (1998), 177–187.

vandermeer2010b
    J. Vandermeer and A. King,
    Consequential classes of resources: Subtle global bifurcation with dramatic ecological consequences in a simple population model,
    J Theor Biol, 263 (2010), 237–241.

vandermeerj2010b
	J. Vandermeer, I. Perfecto, and S. Philpott,
	Ecological Complexity and Pest Control in Organic Coffee Production: Uncovering an Autonomous Ecosystem Service,
	BioScience, 60 (2010), 527–537.
	
vandermeerj2014
	J. Vandermeer, D. Jackson, and I. Perfecto,
	Qualitative Dynamics of the Coffee Rust Epidemic: Educating Intuition with Theoretical Ecology,
	BioScience, 63 (2014), 210–218.
	
waller1982
	J. M. Waller,
	Coffee Rust-epidemiology and control,
	Crop Prot, 1 (1982), 385–404.




