Mathematical modelling and analysis of coffee berry disease dynamics on a coffee farm

1.   Introduction

Coffee is a plantation crop that is well adapted to different eco-physiological conditions of tropical and subtropical highlands. It is the favourite beverage and second most traded commodity next to crude oil in the world ^[1]. The coffee industry is estimated to be over fanxiexian_myfh100 billion worldwide with an average consumption of 500 billion cups per year. It is mostly consumed in the developed nations and produced by tropical countries, often less developed, being a big source of their economy ^[1,2]. However, coffee production has been challenged by diseases, weeds and pests. For instance, coffee berry disease (CBD) caused by Colletotrichum kahawae is a major challenge to Arabica coffee production in African countries ^[3]. It is a fungal plant pathogen ^[4,5]. The symptoms of CBD include black depressed wounds on coffee berries (see Figure 1) ^[6].

The first report of CBD (Colletotrichum kahawae) dates back to 1922 in western Kenya ^[7]. Soon after, the fungus quickly spread throughout most of the African countries (see Figure 2) ^[8]. The disease attacks flowers and fruits at all stages of growth, but it is more destructive to young berries especially during the expanding period 4–l6 weeks after flowering ^[9]. This disease decreases both the quality and yield of coffee berries. For instance, Ethiopia, Kenya and Cameroon annually lose up to 29.9, 75 and 60%, respectively. Losses may reach up to 100% ^[2,10] in high rainfall, humidity and altitude areas.

Because CBD can become very severe and there is a lack of effective control measures, there is much concern that the fungus could spread to coffee-growing areas on other continents. Currently, however, the disease is only prevalent in Africa at high altitudes with relatively high humidity ^[11]. There are various applications for CBD management such as chemical treatments, cultural practices, the use of resistant varieties and biological control ^{[12,13,14,15,16,17,18]}. Several chemical products (fungicides) have been evaluated and recommended for CBD control measures such as different copper fungicides, organic fungicides and mixtures of the two. For instance, fungicide spray against CBD starting six weeks after the main flowering for six rounds at a four-week interval during a crop season was recommended in Ethiopia. Cultural practices include a variety of management tactics, such as mixed cropping with shade plants, pruning infected branches, the destruction of infected material and the removal of mummified berries to create environmental conditions that limit CBD development ^[14]. Chemical control seems the most effective method but if incorrectly used, it causes ecological risks. On the contrary, the use of resistant varieties and cultural practices is cost-effective, biologically safe and environmentally friendly.

Mathematics is becoming an important tool for studying the evolution of plant pests and diseases. Some mathematical models have been formulated and analyzed to explain the dynamics of plant disease transmission and assess their controls. For example, Fotsa et al. ^[19] introduced the mathematical modelling and optimal control of Anthracnose disease. Anthracnose attacks a wide range of commercial crops, including coffee, mango, banana, blueberry, cherry, and strawberry. Cunniffe and Gilligan ^[20] investigated epidemiological models for plant pathogens. These models ensured maximum transferability across the widest range of host-pathogen systems by using the common currency of the field; they illustrated the results in a class of model that has been experimentally tested for plant disease. Fotso et al. ^[21] introduced and analyzed a mathematical model of coffee berry borer with optimal control strategies. Among coffee diseases, CBD, coffee wilt disease, and coffee leaf rust are the most destructive diseases threatening coffee production in developing countries like Ethiopia. However, according to the survey results, mostly CBD was more widespread than the other coffee diseases ^[22]. In all the previous studies, the mathematical model for CBD epidemics concerning fungal pathogen and vector population was not considered. In the present work, we developed a nonlinear deterministic mathematical model for the dynamics of CBD infestation on a coffee farm; and also we also applied their qualitative analysis.

2.   Materials and methods

2.1.   Biological background and pathosystem

The reproduction of the coffee tree involves flowers containing both females and males and their pollination occurs mainly by wind and to a lesser extent by insects ^[23]. CBD is locally spread between coffee trees and branches by wind and rain ^[11]. But, the common vectors of long- and medium-distance CBD dispersal are insects, birds, and coffee harvesters ^[6]. Hence, we assumed those common vectors that, after contacting a fungal pathogen from the environment or already CBD infected coffee berries, then contact healthy coffee berries as responsible vectors for the spread of CBD. The fungal pathogen in the environment may be transported to the coffee plant via infected vectors; once the coffee berry becomes infected, the fungus increases within the infected coffee and destroys it. Twig bark, flower cushions, and mummified berries are considered to be sources of primary inoculum ^[24]. Most diseased berries drop prematurely, but those that stick to the branches are the main sources of the secondary inoculum (new Colletotrichum kahawae) conidia that are dispersed to contaminate other healthy berries.

The CBD is highly dependent upon climatic factors: humidity, rainfall, and temperature ^[25,26,27]. These climatic conditions conducive to Colletotrichum kahawae typically occur at high elevations greater than 1600 m (i.e., the altitude at which coffee Arabica is grown); disease incidence is minimal below 1000 m ^[24]. The optimal temperatures for conidium germination and mycelium growth are in the range of 20 to 22℃ for CBD ^[28]. An increase of 0.61℃ in the global mean temperature has been recorded since the beginning of the twentieth century and the predicted warming of 2–6℃ by 2100 have direly increased the need to understand the impacts of climate change ^[29]. Most insects in temperate climates have optimum development at temperatures between 20 and 35℃. The total development time from egg to adult was 89.6 days at 20℃, 63.10 days at 25℃ and 55.81 days at 30℃, and no development at 35℃ ^[30]. Whereas at temperatures below 15℃ mating is limited, and movement such as flying becomes stalled ^[31]. Most conidia of Colletotrichum kahawae can be effectively dispersed by an optimum rainfall of 10 mm but heavy rainfall leads more to their leaching from the coffee tree canopy to the soil ^[32,33].

2.2.   Mathematical model

This study considers coffee berry and vector populations with the interaction of a fungal pathogen $(B)$. The coffee berry population is considered into two forms: the susceptible coffee berry $(S_{c})$ and infected coffee berry $(I_{c})$. Due to the limited size of a coffee plantation, we adopted a logistic growth function for the biomass density of coffee berries. The coffee berry population grows logistically ^[34] at a growth rate $r$ and an environmental carrying capacity $K$. Susceptible coffee berries move to the infected class following contacts with an infected vector at a per capita rate $\beta_{2}$; its mortality rate is $\theta$. The infected coffee berry's death rate is $\gamma$ due to the disease. The coffee berry, once becoming infected, never recovers, and gives no or very low yield of coffee. The vector population is divided into two form: the susceptible vector $(S_{v})$ and infected vector $(I_{v})$. The infected vectors (insect, birds) are assumed to transport the fungal pathogen to the coffee berry. The vectors can be infected from the fungal pathogen in the environment at rate $\beta_{1}$ or from the infected coffee berry at rate $\beta_{3}$. The susceptible vector is recruited at rate $\lambda$. The vectors' death rate is represented by $\delta$. The fungal pathogens in the environment may be transported to the coffee plant via the infected vectors, once the coffee berry becomes infected, the fungal increases within the infected coffee and destroy it. This invariably adds to the fungal pathogens in the environment at rate $\eta$ and the infected coffee berry fungal pathogen contribution to the environment is ($\eta I_c$). The fungal pathogen's decay rate is $m$. The flow diagram of CBD transmission is illustrated in Figure 3. For further information, the parameters and their biological meanings are given in Table 1.

The governing mathematical model is given by

Together with initial condition:

Note that in the Table 1, the unit of $\eta$ is to be ${\rm{Cells}}{\rm{(ml}}^{-1}){\rm{(Individual}}^{-1}){\rm{(day}}^{-1})$ = $\mho$.

3.   Results

3.1.   Model analysis

3.1.1.   Positivity and boundedness of solutions

In this section, we need to prove that the solutions of system (2.1) are nonnegative for all $t\geq 0.$ This will be stated as follows.

Theorem 1. Every solution of system (2.1) with initial condition (2.2) will remain positive in $\mathbb{R}_+^{5}$ for all $t > 0$.

Proof. From the system (2.1), we obtain

This implies that all the solutions with positive initial data remains nonnegative in $\mathbb{R}_+^{5}$ for all $t\geq 0$. This means, the region attracts all solutions of the governing system (2.1).

Moreover, it is easy to verify that there exist an invariant region $\Omega$ where a solution set for the system (2.1) is bounded.

Theorem 2. Every solution of system (2.1) initiating in $\mathbb{R}_+^{5}$ is uniformly bounded in the region $\Omega = \Omega_{c}\times\Omega_{v}\times\Omega_{p}$, where

with $M = \max\left\{S_{c}(0), \; K \right\}, \; h = \min\left\{1, \gamma+\eta \right\}$.

Proof. Let $N_{c}(t) = S_{c}(t)+ I_{c}(t)$ be the total population of coffee berry. Then differentiating $N_{c}(t)$ with respect to time $t$ and adding the first two equations of system (2.1), we get

The last inequality in Eq (3.1) can be rewritten as

Integrating Eq (3.2) by separation of variables and then as $t$ $\rightarrow$ $\infty$, the feasible region of the system (2.1) for coffee population is given by

We also consider the total population of vector as $N_{v}(t) = S_{v}(t)+ I_{v}(t)$. Then the derivative of $N_{v}(t)$ with respect to $t$ is given by

By solving Eq (3.3) and then as $t \rightarrow \infty$, we get the feasible region of the system (2.1) for vector population:

From the fifth equation of system (2.1) and by the comparison, we have

Solving the last inequality of Eq (3.5) and $t \rightarrow \infty$ results $0\leq B(t)\leq M(1+ r)\eta/mh$. Consequently, the feasible region of the system (2.1) is given by $\Omega = \Omega_{c}\times\Omega_{v}\times\Omega_{p}$ is positively invariant.

Hence, the solutions of the model are bounded. Therefore, the model is suitable to conduct the study, and the analysis of the system dynamics can be considered in the region $\Omega$.

3.1.2.   Existence and uniqueness of solutions

In this section, we provide the following results which guarantee that the CBD model governed by system (2.1) is epidemiologically and mathematically well posed.

Theorem 3. (Existence-uniqueness of solution). Let $\Theta$ be the domain:

and suppose that $f(t, x)$ satisfies the Lipschitz condition:

where $(t, x_1)$, $(t, x_2)$ $\in\Theta$, and $\kappa > 0$. Then, there exist a constant $\nu > 0$ such that there exists a unique continuous vector solution $x(t)$ of the system (2.1) in the interval $\left|t-t_0\right|\leq \nu$. From Eq (3.7), $f$ is satisfied by requirement that $\partial f_{i}/\partial x_{j}, \; i, j = 1, 2, ..., 5$ be continuous and bounded in $\Theta$.

If $f (t, x)$ has continuous partial derivative $\partial f_{i}/\partial x_{j}$ on a bounded closed convex domain $\mathbb{R}$, then it satisfies a Lipschitz condition in $\mathbb{R}\in(0, \infty)$. Our concern is in the domain:

Let $\Theta$ denote the region defined in Eq (3.6) such that Eqs (3.7) and (3.8) hold. Then, there exist a solution of system (2.1) which is bounded in $\Theta$.

Proof. {Suppose that the right side of system (2.1) is re-written as:

Since the functions in Eqs (3.9)–(3.13) are polynomial, then they are infinitely differentiable. Thus, for the detailed proof, we can follow the proof of classical Cauchy-Lypschitz theorem ^[43]. Therefore, all the partial derivatives exist and are finite, then there exists a solution for the model and hence, we can say that there exist a unique solution of system (2.1) in the domain $\Theta$.

3.1.3.   Disease-free equilibrium point (DFE)

At DFE, there is no infections (i.e., $I_{c} = I_{v} = B = 0$) in the populations. To find DFE, we equate the right hand side of system (2.1) to zero and solving for the noninfected state variables we get $S_{c} = 0$ or $S_{c} = K-K\theta/r, \; S_{v} = \lambda/\delta$. Thus, the equilibrium points of the system (2.1) are $(0, \; 0, \; \lambda/\delta, \; 0, \; 0)$ and $(K(r-\theta)/r, \; 0, \; \lambda/\delta, \; 0, \; 0)$. Hence, the DFE is given by

Remark 1. The disease-free equilibrium poin $E_{0}$ in Eq (3.14) is biologically feasible when $r > \theta$.

3.1.4.   Basic reproduction number ($R _{0}$)

The basic reproduction number is the number of secondary infections that one infectious individual would create over the duration of the infectious period ^[37]. The expression of $R _{0}$ for the system (2.1) can be derived using the next generation matrix method ^[38]. The first step to find $R_{0}$ is rewriting the model equations starting with newly infective classes:

The right hand side of system (3.15) is written as $f-g$, where

Then by linearization approach, the associated matrices of $f$ and $g$ at $E_{0}$ are given by

The basic reproduction number, $R _{0} = \rho(\mathbb{FG}^{-1})$ is the spectral radius of the product $\mathbb{FG}^{-1}$. Thus, it is given by

3.1.5.   Local stability of disease-free equilibrium

In this section, we shall investigate the local stability of disease-free equilibrium point based on the basic reproduction number, $R _{0}$.

Theorem 4. If $R_{0} < 1$, then $E_{0}$ of the system (2.1) is locally asymptotically stable in $\Omega$.

Proof. The Jacobian matrix of the system (2.1) at $E_{0}$ is given by

From the Jacobian matrix (3.18), we obtain characteristic polynomial:

From Eq (3.19), we obtain that $\Lambda_{1} = -\delta < 0, \;$ $\Lambda_{2} = -(r-\theta) < 0$ (see Remark 1, $r > \theta$) and

where

Using Routh-Hurwitz criteria ^[38], Eq (3.20) has eigenvalues which are negative real part whenever $A_{0} > 0$, $A_{1} > 0$, $A_{2} > 0$, $\; A_{1}A_{2}-A_{0} > 0$ and $A_{0}(A_{1}A_{2}-A_{0}) > 0$. Therefore, $E_{0}$ is locally asymptotically stable for $R_{0} < 1$.

3.1.6.   Global stability of disease-free equilibrium

Theorem 5. If $R_{0} < 1$, then $E_{0}$ of the system (2.1) is globally asymptotically stable in $\Omega$.

Proof. To proof this theorem, we first define the linear Lyapunov function:

where $\epsilon_{1}, \epsilon_{2}$ and $\epsilon_{3}$ are positive constants to be computed.

Then differentiating $U$ with respect to $t$, we obtain

with $S_{c}^{0} = K(r-\theta)/r$ and $S_{v}^{0} = \lambda/\delta$. The $\epsilon_{1}, \epsilon_{2}, \epsilon_{3}$ can be chosen such that

If we choose $\epsilon_{1} = 1$, then $\epsilon_{2} = \beta_{2}K(r-\theta)/r\delta$ and $\epsilon_{3} = \beta_{1}\beta_{2}K(r-\theta)\lambda/r\delta^{2} m$. Substituting $\epsilon_{1}, \epsilon_{2}$ and $\epsilon_{3}$ into Eq (3.21) leads to

Since $R_0 < 1$ and $R_0\geq0$ imply that $R_{0}^{2} < 1$, we have $dU/dt\leq0$.

Moreover, the largest compact invariant set in $\left\{(S_{c}, \; I_{c}, \; S_{v}, \; I_{v}, \; B)\in \Omega: dU/dt = 0\right\}$ is the singleton $\left\lbrace E_{0}\right\rbrace$. By LaSalle ^[39], it then implies that $E_{0}$ is globally asymptotically stable in $\Omega$.

3.1.7.   Endemic equilibrium point (EEP)

Endemic equilibrium point ($E^{*}$) exists when CBD persist in the populations. To obtain EEP, we set each right hand equation in the system (2.1) equal to zero. That is,

From the first equation of Eq (3.24) we find $S_{c}^{*} = 0$ or $r(1-(S_{c}^{*}+I_{c}^{*})/K)-\beta_{2}I_{v}^{*}-\theta = 0$.

Case i: $S_{c}^{*} = 0$, we obtain $I_{c}^{*} = \delta m/(\beta_1\eta+m\beta_3), \; S_{v}^{*} = \lambda/2\delta, \; I_{v}^{*} = \lambda/2\delta$, $B^{*} = \eta \delta/(\beta_1\eta+m\beta_3)$. This endemic point is a point where the disease kills all coffee berries.

Case ii: $r(1-(S_{c}^{*}+I_{c}^{*})/K)-\beta_{2}I_{v}^{*}-\theta = 0$, then we get $S_{c}^{*}, \; S_{v}^{*}, I_{v}^{*}$ and $B^{*}$:

where $I_{c}^{*}$ is the positive root of

with

The feasibility of the endemic equilibrium is determined by the next proposition.

Proposition 1. (Feasibility of the endemic equilibrium).

Let $E^{*} = (S_{c}^{*}, I_{c}^{*}, S_{v}^{*}, I_{v}^{*}, B^{*})$ be as Eq (3.25).

(i) If $B > 0$ and $C < 0$ or $B < 0$ and $C < 0$, then there is a unique endemic equilibrium $E^{*}$.

(ii) If $B < 0$ or $C > 0$, then there exist two endemic equilibrium $E^{*}$ of the system.

Proof. The result follows Descartes' rule of signs to $f(I_{c}^{*})$ given in Eq (3.26).

3.1.8.   Local stability of endemic equilibrium point

Theorem 6. If $R_{0} > 1$, $B_{i} > 0, i = 1, 2, ..., 5$, $B_{1}B_{2}-B_{3} > 0, B_{1}(B_{2}B_{3}+B_{5})-B_{3}^{2}-B_{1}^{2}B_{4} > 0, (B_{2}B_{3}$ $-B_{1}(B_{2}^{2}-2B_{4}))B_{5}-B_{4}^{2}(-B_{1}B_{2}B_{3}+B_{3}^{2}+B_{1}^{2}B_{4})-$ $B_{5}^{2} > 0, B_{5}(-B_{4}(-B_{1}B_{2}B_{3}+$ $B_{3}^{2}+B_{1}^{2}B_{4})+(B_{2}B_{3}-B_{1}(B_{2}^{2}-2B_{4}))B_{5}-B_{5}^{2}) > 0$, then $E^{*}$ of system (2.1) is locally asymptotically stable in $\Omega$.

Proof. The Jacobian matrix of system (2.1) at $E^{*}$ is given by

The eigenvalues of matrix (3.27) are computed from

Let us consider the non-zero entities of matrix (3.28) as:

$b_{11} = -\theta-r(I_{c}^{*}+2S_{c}^{*})/K-\beta_{2}I_{v}^{*}+r, \; b_{12} = -(r-\theta), \;$$b_{21} = \beta_{2}I_{v}^{*}, \;$$b_{22} = -\gamma-\eta, \\ $$b_{32} = -\beta_{3}S_{v}^{*}, \;$$b_{42} = \beta_{3}S_{v}^{*}, \;$$b_{52} = \eta, \;$$b_{33} = -\beta_{1}B^{*}-\beta_{3}I_{c}^{*}-\delta, \;$ $b_{43} = \beta_{1}B^{*}+\beta_{3}I_{c}^{*}, \\ $$b_{14} = -\beta_{2}S_{c}^{*}, \;$$b_{24} = \beta_{2}S_{c}^{*}, \;$$b_{44} = -\delta, \;$$b_{35} = -\beta_{1}S_{v}^{*}, \;$$b_{45} = \beta_{1}S_{v}^{*}, \;$$b_{55} = -m$.

Then, the characteristic polynomial of Eq (3.28) is given by

where

Using the Routh Hurwitz criterion ^[39], the endemic equilibrium $E^{*}$ is locally asymptotically stable for $R_{0} > 1$ if $B_{i} > 0, i = 1, 2, ..., 5$ and

Hence, all eigenvalues of Eq (3.29) evaluated at $E^{*}$ have negative real parts whenever $R_{0} > 1$ and the equilibrium $E^{*}$ is locally asymptotically stable.

3.1.9.   Global stability of endemic equilibrium

Theorem 7. If $R_{0} > 1$, then $E^{*}$ of the model (2.1) is globally asymptotically stable in $\Omega$.

Proof. Let us consider a Volterra-type Lyapunov function:

where $\zeta_1, \zeta_2, \zeta_3, \zeta_4$ and $\zeta_5$ are non-negative constants. Then by taking the derivative of $V$ with respect to $t$, we obtain

Since $E^{*}$ is an equilibrium point, Eq (3.24) can be rewritten as:

Substituting the relations in Eq (3.32) into Eq (3.31), we obtain

Expanding and grouping Eq (3.33), we have

Choose $\zeta_1, \zeta_2, \zeta_3, \zeta_4$ and $\zeta_5$ such that the expressions in the brackets in Eq (3.34) vanish. That is

Fixing $\zeta_1 = \delta$, then Eq (3.35) leads to

Replacing Eq (3.36) in Eq (3.34), we obtain

or equivalently

From Eq (3.38), since $\dfrac{\eta I_c^{*}S_c^{*}}{(\beta_1+\beta_3m)S_v^{*}B^{*}}\left(\dfrac{\beta_2I_v^{*}}{I_c^{*}}-\dfrac{r}{K}\right) = \zeta_3\geq0$, we have

With this in mind, and if we let

then $(dV/dt)(S_{c}, \; I_{c}, \; S_{v}, \; I_{v}, \; B)\leq 0$ if and only if $S_{c} = S_{c}^{*}$, $\; I_{c} = I_{c}^{*}$, $\; S_{v} = S_{v}^{*}$, $\; I_{v} = I_{v}^{*}, \; B = B^{*}$ and $Z+(-W)\leq 0$.

Therefore, the largest compact invariant set in $\left\{(S_{c}, \; I_{c}, \; S_{v}, \; I_{v}, \; B)\in \Omega: dV/dt = 0\right\}$ is the singleton $\left\lbrace E^{*}\right\rbrace$. By LaSalle ^[39], it then implies that $E^{*}$ is globally asymptotically stable in $\Omega$.

3.1.10.   Bifurcation analysis

To investigate the possibility of existence of the bifurcation analysis of the system (2.1) at $R_{0} = 1$, we use the center manifold theory ^[40]. So using this approach, the following theorem can be established.

Theorem 8. If $R_{0} < 1$, then the system (2.1) exhibits forward bifurcation at $R_{0} = 1$.

Proof. Using center manifold theory ^[40], we carry out bifurcation analysis of system (2.1) at $R_{0} = 1$. Let us consider the change of variables $S_{c} = x_{1}, \;$ $I_{c} = x_{2}, \;$ $S_{v} = x_{3}, \;$ $I_{v} = x_{4}, \;$ $B = x_{5}$. Furthermore, using the vector notation: $X = (x_{1}, x_{2}, x_{3}, x_{4}, x_{5})^{T}$. Then the system (2.1) can be written in the form $dX/dt = (g_{1}, \; g_{2}, \; g_{3}, \; g_{4}, \; g_{5})^{T}$ as:

We consider $\beta_{2}$ as bifurcation parameter so that $R_{0} = 1$ if

The Jacobian matrix of Eq (3.39) at disease-free equilibrium $E_{0}$ is given by

The right eigenvector, $u = (u_{1}, \; u_{2}, \; u_{3}, \; u_{4}, \; u_{5})^{T}$ is computed from $J(E_{0})u = 0$ so that

where $u_{4} = u_{4} > 0.$ Similarly, the left eigenvector, $z = (z_{1}, \; z_{2}, \; z_{3}, \; z_{4}, \; z_{5})$ is computed from $zJ(E_{0}) = 0$ so that $z_{1} = z_{3} = 0, \; z_{2} = \dfrac{\delta r}{K(r-\theta)\beta_{2}^{*}}z_{4}, \; z_{5} = \dfrac{\beta_{1}\lambda}{\delta m}z_{4};\; \; z_{4} = z_{4} > 0.$ The nonzero second partial derivatives of $g_{1}, \; g_{2}, \; g_{3}$ and $g_{4}$ are given as follows:

All the others second partial derivatives of $g_{i}, \; i = 1, ..., 5$ are zero. Based on center manifold theorem ^[40], the coefficients $a$ and $b$ are given by

Since $a$ is negative and $b$ is positive, then the system (2.1) exhibits forward bifurcation at $R_{0} = 1$ and there exists at least one stable endemic equilibrium when $R_{0} > 1$ (See Figure 4).

3.1.11.   Sensitivity of the basic reproduction number

In this section, we compute the sensitivity index of basic reproduction number, $R_{0}$, with respect to the main parameters ^[41]. This help us to check and identify parameters which highly affect $R_{0}$ on the eradication or spread of CBD. For example, sensitivity index of $R_{0}$ with respect to $K$ is given by

In the same way, the sensitivity index of $R_{0}$ are also computed and the corresponding sensitivity indeces are given in Table 2.

3.1.12.   Interpretation of the sensitivity indices

The parameters: $K, \; \beta_{1}, \; \beta_{2}, \; \beta_{3}, \; r$ and $\lambda$ have positive sensitivity indices show that the more impact on expanding the disease in the population if their values are increasing since $R_{0}$ increases. Besides, the parameters: $\theta, \; \eta, \; \delta, \; \gamma$ and $m$ have negative sensitivity indices show that less impact on expanding the disease as their values increase. The bar diagram of the sensitivity indices of $R_{0}$ in Table 2 is depicted in Figure 5. Hence, the parameters: $K, \; \beta_{2}$ and $\lambda$ are the most in influential in expanding the disease.

Furthermore, in Figure 6, we have seen that the contour plot of basic reproduction number $R_{0}$ with respect to the parameters $K$, $\beta_{2}$ and $K$, $\lambda$, respectively. It can be observed from Figure 6(a) that for decreasing value of $K$ and $\beta_{2}$, $R_{0}$ decreases. For increasing $K$ and $\lambda$ the basic reproduction number $R_{0}$ increases, which has been shown in Figure 6(b).

3.2.   Numerical simulations

In order to illustrate the analytical results, the numerical simulations of the system (2.1) were performed. Some parameter values were assumed and some of them were obtained from literature. In addition to parameter values in Table 1, we assumed the initial data: $S_{c}(0) = 100, \; I_{c}(0) = 50, \; S_{v}(0) = 50, \; I_{v}(0) = 10, \; B(0) = 2$.

Figures 7(a) and 8(a) depict the global asymptotic behaviour of the infected coffee and infected vector at the disease-free equilibrium $E_{0}$ when $R_{0} < 1$ where the disease dies out regardless of large initial values of the infectious population. On the other hand, Figures 7(b) and 8(b) show the global asymptotic behaviour of the infected coffee and infected vector at the endemic equilibrium $E^{*}$ when $R_{0} > 1$ where the disease persists regardless of the initial sizes of the infectious population.

From the Figure 9(a) susceptible coffee berries decreases while the infected coffee berries increases with $R_{0} > 1$ due to there is very favorable condition for breeding of pathogen and vector that transmit the disease. The susceptible vector population in the Figure 9(c) decreases while infected vector population increases. The CBD pathogen population in the Figure 9(b) increases to a maximum point and then proceed decrease with minimum at unfavorable climate.

In the Figure 10(a), (b), the effects of the modification factor for asymptomatic coefficient at different values: $K = 150$ and $\beta_{2} = 0.000795455$ are displayed. The projection shows that the cumulative number of individuals becoming infected coffee berries is high when $K$ and $\beta_{2}$ become high and becoming decreasing with low values of $K$ and $\beta_{2}$, which minimize also the values of basic reproduction number less one. As shown in Figure 10(c), (d), when the rates: $K$ and $\beta_{2}$ increase, the number of individuals becoming infected vector is high and decreasing the value of decreases the number of individuals becoming infected. As shown in Figure 11(a), (b), when the rates: $K$ and $\beta_{2}$ increase, the concentration of the pathogen becomes high, and decreasing the value of decreases the pathogen.

4.   Discussion and conclusions

In this paper, we formulated and analyzed a nonlinear deterministic mathematical model for CBD transmission dynamics in a coffee farm. A compartmental mathematical model by including coffee plants and vectors of disease with the interaction of fungal pathogen was investigated. We obtained the feasible region where the model is epidemiologically and mathematically meaningful. The positivity, existence and uniqueness of model solutions are also demonstrated. The equilibria (disease-free and endemic equilibrium) of the model were computed. We calculated the basic reproduction number at the disease-free equilibrium point by using the next-generation matrix method. The local stability of disease-free and endemic equilibria is shown using the Routh-Hurwitz criteria. Besides, the global stability of equilibria is proved by using the Lyapunov function. Moreover, bifurcation analysis is proved by the Center Manifold theory. The sensitivity indices of basic reproduction number with respect to the main parameters are computed. Further, sensitivity indices are graphically shown and the most influential parameters in expanding the disease are identified. Despite coffee playing a dominant role in the social, cultural, and national economy, the country's coffee industry is potentially at risk due to climate changes. We discussed the impact of Colletotrichum kahawae fungus and vectors of disease in relation to climatic variables. More vectors and CBD pathogens are produced at the favourable climate which indirectly increases the rate of infected coffee berries. Thus, the frequency and severity of climate extremes are increasing and making adaptation an absolute imperative necessity through using current information on climate variability to develop long-term plans for managing CBD. Since CBD is highly dependent upon climatic factors, we are doing for the future with a paper applying the climatic variability in this model. Numerical simulations of the model suggested we apply effective control interventions on some parameters to control the disease. This will be done using chemical, cultural, and biological control strategies. We conclude that according to our mathematical model we have an endemic equilibrium point and the coffee disease remains endemic provided that the basic reproduction number is greater than one. Thus, we are developing the model with optimal control for future work to eradicate the disease.

Declaration of Competing Interest

The authors declare that there is no conflict of interests.

Acknowledgements

The authors wish to express their special gratitude to the editor and the reviewers for the helpful comments given for this paper. This work was supported by Adama Science and Technology University (grant number: ASTU/SP-R/001/19). We would like to express our appreciation for their support.