# Modelling and optimal strategy to control coffee berry borer

Yves Fotso Fotso, Suzanne Touzeau, Berge Tsanou, Samuel Bowong, Frédéric Grognard

► To cite this version:

Yves Fotso Fotso, Suzanne Touzeau, Berge Tsanou, Samuel Bowong, Frédéric Grognard. Modelling and optimal strategy to control coffee berry borer. Mathematical Methods in the Applied Sciences, 2021, 44 (18), pp.14569-14592. 10.1002/mma.7726. hal-03334090

HAL Id: hal-03334090

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

Submitted on 23 Jan 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
## RESEARCH ARTICLE

# Modelling and optimal strategy to control coffee berry borer

Yves Fotso Fotso*<fnref n="1" />,<fnref n="2" />,<fnref n="6" /> | Suzanne Touzeau<fnref n="2" />,<fnref n="3" /> | Berge Tsanou<fnref n="1" />,<fnref n="4" />,<fnref n="6" /> | Samuel Bowong<fnref n="5" />,<fnref n="6" /> | Frédéric Grognard<fnref n="2" />

### Correspondence

*Yves Fotso Fotso, Department of Mathematics and Computer science, University of Dschang, Cameroon. Email: fotsofyves@yahoo.fr

### Summary

Coffee berry borer *Hypothenemus hampei* (Coleoptera: Scolytidae), denoted CBB, is the most important insect pest of coffee worldwide, with a high impact on the economy of coffee producing countries. The insect spends a great part of its life cycle inside the coffee berry and causes severe crop losses. Biological control based on the use of an entomopathogenic fungus is a major alternative to chemical pesticides in order to control CBB. The fungus is sprayed on the coffee berries to kill CBB when the insects drill holes to penetrate inside the berries. Our aim in this work is to optimise the fungus application, using a modelling approach. We first formulate a mathematical model describing the infestation dynamics of coffee berries by CBB. We analyse the model and show that the stability of the pest-free and coexistence equilibria depends on the basic reproduction number. To introduce a control variable corresponding to the application of an entomopathogenic fungus, we then extend the model and include the fungus dynamics. We formulate an optimal problem which consists in maximising the coffee yield, while minimising the control cost, as well as the CBB population for the next cropping season. The existence of the optimal control and the necessary optimality condition are established using Pontryagin's Maximum Principle. The optimal control problem is solved numerically using the BOCOP software and simulations are provided, showing that the use of entomopathogenic fungus effectively controls CBB.

### KEYWORDS:

population dynamics, plant-pest interaction, optimal control, ODE model, numerical simulations

MSC classification: 34D20; 49J15; 49K15; 92D30

## 1 | INTRODUCTION

For many years, the coffee economy has been a major issue in the world trade. Originating from Ethiopia, coffee is mainly produced by developing countries amongst which Brazil, Vietnam, Columbia are the top producers in the world in the dereasing order. As far as Africa is concerned, the top producers are Ethiopia, Ivory Coast, Uganda and Cameroon [1]. Coffee is mainly consumed by developed countries in America, Asia and Europe [1]. There are more than a hundred species of coffee, but only two of these are commercially traded: *coffea canephora* (predominantly known as robusta) and *coffea arabica*. Coffee production has increased in recent years. According to the United States Department of Agriculture, the world coffee production increased from 152.9 million bags (a bag weighing 60 kg) in 2015–2016 to 174.6 million bags in 2018–2019, and it is forecast at 169.3

<footnote n="1">¹Department of Mathematics and Computer science, University of Dschang, Cameroon</footnote>

<footnote n="2">²Université Côte d'Azur, INRIA, INRAE, CNRS, Sorbonne Université, BIOCORE, France</footnote>

<footnote n="3">³Université Côte d'Azur, INRAE, CNRS, ISA, France</footnote>

<footnote n="4">⁴Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria, South Africa</footnote>

<footnote n="5">⁵Department of Mathematics and Computer Science, University of Douala, Cameroon</footnote>

<footnote n="6">⁶IRD, Sorbonne Université, UMMISCO, Bondy, France</footnote>
million bags in 2019–2020 [1]. Coffee production is affected by several pests that cause considerable damages. The main pests encountered in coffee plantations are *Hypothenemus hampei*, *Leucoptera coffeella*, *Leucoptera caffeina*, *mealybugs* and nematodes. The most devastating is *Hypothenemus hampei* [2, 3], a coffee berry borer hereafter denoted by CBB. It is an insect that develops in coffee berries at different stages of ripening. Originally from Central Africa, CBB can now be found in most coffee growing countries throughout the world [2]. It can cause yield losses of 30–35% and it deteriorates the berry quality with up to 100% of berries perforated at harvest time [4]. Worldwide damages are estimated at more than USD 500 million each year [5] and threaten more than 20 million farming families [2, 6, 7].

Coffee berry borer control worldwide is mainly based on cultural practices, which consist in removing dropped berries and strip-picking trees after harvest or stump pruning [3, 8]. Chemical insecticides such as endosulfan can be used, but they have harmful consequences on the environment and induce the development of resistance [9, 10]. As a consequence, alternative control methods have been developed, such as trapping, using attractive traps to capture colonising CBB as they emerge from the berries [11, 12, 13], and biological control. The most commonly used biocontrol agents are parasitoids (*Prorops Nasuta*, *Phyrostichus coffea*), predators (*Cephalonomia stephanoderis*, *Cathartus quadricollis*) and entomopathogenic fungi (*Beauveria bassiana*, *Vuillemin*) [3, 7, 8, 14, 15]. *Beauveria bassiana* naturally infects CBB and was reported in less than 1% of CBB in Brazil, but up to 44% in Nicaragua, 60% in India and 71% in Cameroon [8].

To determine how to best apply a biocontrol agent such as this entomopathogenic fungus, field experiments are particularly time consuming and expensive [8], as several plots would need to be monitored during several months. Moreover, it is difficult to control all factors in the field, such as for instance the initial CBB population. Mathematical modelling and optimal control are hence relevant cost-effective alternatives to assess control strategies and long term dynamics of plant-pest interactions [16]. In the literature, various models have been developed to get insights into the evolution of CBB. A statistical model analyses the spatio-temporal CBB infestation in Colombia [17]. A tritrophic coffee–CBB–parasitoids simulation model was developed [18] and further extended [19], based on field data and including detailed biological processes. Due to its complexity, this model is analytically untractable, so neither can its qualitative behaviour be investigated, nor its data-dependent outcomes generalised. Recently, we have proposed a preliminary model to describe the infestation dynamics of coffee berries by CBB, which takes into account berry availability [20].

Our aims in this work are twofold: firstly, to extend this preliminary model and formulate a more comprehensive model, with a control that can persist over time; secondly, to propose and implement an optimal control strategy in order to maximise the yield while minimising the cost, as well as the CBB population at the end of the cropping season to prevent too high infestations for the next season. This control problem was first introduced in [20], but the control used was non persistent and a quadratic objective function was considered. In this work, we propose a more realistic optimal control problem, based on the release of the entomopathogenic fungus *B. bassiana* and a linear objective function, which represents the cost of fungus application and hence the profit better than a quadratic function.

This paper is organised as follows. The mathematical model is formulated in Section 2, based on the CBB life cycle; the positivity and asymptotic behaviour of the model are determined in Section 3. Section 4 presents the optimal control problem, which aims at maximising the profit while minimising the CBB population for the next cropping season. The problem is solved analytically using Pontryagin's maximum principle (PMP) [21] in Section 4 and numerically using a direct method implemented in the BOCOP software in Section 5. Finally, conclusions are drawn in Section 6.

# 2 | MODELLING COFFEE BERRIES-CBB INTERACTION

## 2.1 | Biological background

Coffee phenology depends on climate and geographical region with two and more flowering periods in some regions that allow more than one harvest per year. For example, there are two flowering periods in the central coffee region of Colombia [8]; in the main production areas of Brazil, two to four flowering periods per year are common [19]. This favours the CBB, because the insect finds food and shelter for reproduction all year long, making its control more difficult [3, 8].

The life cycle of CBB is hemimetabolous (incomplete) with four stages: egg, larva, nymph and adult [3, 8, 22]. The first three life stages are spent inside the coffee berry. Adult males remain in the berry, while fertilised female, which are able to fly, emerge from the berry to find new berries for oviposition. These females, usually called colonising females, are responsible for CBB dispersal in the plantation. They search for an uninfested berry, which they enter by drilling a hole. They lay 2 to 3 eggs per day during approximately 20 days and remain inside the berry until they die [8]. The adult female lifespan ranges from 87 to
282 days, but males are shorter lived with a longevity reported around 40 days [8]. Coffee berry borer has a female biased sex ratio of 10:1 females to males [23], with mating taking place inside the berry between mates of same generation [24]. Due to the extended oviposition period of the females, all CBB life stages can be found inside a berry at the same time. The life cycle from eclosion to adults takes between 28 and 30 days [3].

## 2.2 | Model formulation

The model describes the infestation dynamics of coffee berries by CBB during a cropping season in a plantation, taking into account the CBB life cycle and the berry availability. At each time unit $t$, the coffee berries in the plantation are subdivided into two compartments: the healthy berries represented by their number $s(t)$ and the infested coffee berries represented by $i(t)$. As only fertilised females are responsible for CBB dispersal, we only consider this development stage in the model. We subdivide these females into two compartments: the colonising females represented by their number $y(t)$, which correspond to the flying fertilised females looking for their host, and the infesting females represented by their number $z(t)$, which correspond to the females that are laying eggs inside the berries. Figure 1 illustrates the infestation dynamics of coffee berries by CBB.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Model diagram. Healthy berries are produced at a constant rate ($\Lambda$). CBB colonising females ($y$) enter healthy berries ($s$) with ratio-dependent rate (parameters $\alpha$, $\beta$ and berry-CBB conversion rate $\varepsilon$). Colonising CBB become infesting females ($z$) and healthy berries become infested ($i$). Infesting females lay eggs inside the berries, from which colonising females emerge (at rate $\phi$). All CBB and berry stages undergo mortality (at rates $\mu$, $\nu$, $\mu_y$, $\mu_z$).</figcaption></figure>

In accordance with Figure 1, we derive the following system of nonlinear ordinary differential equations:

$$ 
\left\{
\begin{aligned}
\dot{s} &= \Lambda - \beta \frac{sy}{y + \alpha s} - \mu s, \\
\dot{i} &= \beta \frac{sy}{y + \alpha s} - \nu i, \\
\dot{y} &= \phi z - \varepsilon \beta \frac{sy}{y + \alpha s} - \mu_y y, \\
\dot{z} &= \varepsilon \beta \frac{sy}{y + \alpha s} - \mu_z z.
\end{aligned}
\right.
$$

The model hypotheses are described below.

* We assume that flowering occurs throughout the season, which translates into a constant healthy berry production rate $\Lambda$.

* Although the sex ratio is largely female biased, we assume that there are enough males in each generation to fertilise young females (no fertility issues were raised in the field).

* We represent the infestation through a ratio-dependent function $f(s, y) = \beta \frac{y}{sy+\alpha s}$. Ratio-dependence [25, 26, Chapters 1 and 6] is particularly adapted to this pathosystem as there is a one-to-one interaction between CBB and berries. Ratio-dependent interactions are invariant to a proportional change of all abundances in the system. Moreover, classical Holling
type II functional response implements a saturation term on the resource (due to handling time or satiety), which is ill-adapted to this case, especially when CBB are abundant. With a ratio-dependent function, when healthy berries are largely outnumbered by colonising females ($s \ll y$), berries determine the infestation response ($f(s, y) \approx \beta s$), and vice versa. The definition of $f$ is extended to $f(0, 0) = 0$ [27]. Parameter $\beta$ represents the infestation rate per unit of time and $\alpha$ the infestation half saturation constant.

* Super-parasitism is rarely observed in a plantation: colonising females can usually detect infested coffee berries, thanks to the hole drilled in the berry by the previous CBB, and they rather choose another healthy berry to lay their eggs [18, 19]. Hence, a healthy coffee berry is usually infested by a single colonising female. The parameter $\varepsilon$, which is the number of colonising females that can infest a healthy coffee berry and become infesting females, is hence equal to 1 CBB per berry.

* The number $\phi$ of colonising females emerging from the berries per infesting female and per time unit is supposed constant, as infesting females lay eggs at a constant rate over their fertility period.

* The mortality rates of healthy and infested coffee berries are $\mu$ and $\nu$ respectively, with $\mu < \nu$ as CBB infestation damages the berries.

* $\mu_y$ is the natural mortality rate of colonising females; the average life expectancy $1/\mu_y$ of colonising females is much shorter than the average adult lifespan, as colonising females do not have access to coffee berries, which are their sole food source. Infesting females survive their egg-laying period for quite some time but, as they then do not contribute to the infestation dynamics anymore, they are not accounted for in the z-compartment. Hence, $1/\mu_z$ corresponds to the average egg-laying period of infesting females.

Since variable $i$ has no impact on the dynamics of the three remaining variables, we can reduce our analysis to the following model:

$$ 
\left\{ 
\begin{aligned}
\dot{s} &= \Lambda - \beta \frac{sy}{y + \alpha s} - \mu s \\
\dot{y} &= \phi z - \varepsilon \beta \frac{sy}{y + \alpha s} - \mu_y y \\
\dot{z} &= \varepsilon \beta \frac{sy}{y + \alpha s} - \mu_z z
\end{aligned}
\right.
\qquad (1) $$

# 3 | MODEL ANALYSIS

## 3.1 | Positivity

For system (1) to be mathematically and biologically meaningful, it is important to verify that all state variables remain nonnegative and are bounded for nonnegative initial conditions and for all $t \ge 0$. We have the following result.

**Proposition 1.** The compact domain

$$ \Gamma := \left\{ (s, y, z) \in \mathbb{R}_+^3 : s \le \frac{\Lambda}{\mu}, \varepsilon s + z \le \frac{\varepsilon \Lambda}{\zeta}, y \le \frac{\varepsilon \phi \Lambda}{\zeta \mu_y} \right\}, $$

where $\zeta = \min\{\mu, \mu_z\}$, is positively invariant for the flow of system (1).

*Proof.* Since the right-hand-side of system (1) is locally Lipschitz continuous on the non-negative orthant, this system has a unique solution for any non-negative initial condition. Hence, in order to prove positivity of the solution of this system, we can consider each state variable at 0 (with the other taking non-negative values) and check that it cannot become non positive. Since we have $\dot{s} = \Lambda > 0$ when $s = 0$, $\dot{y} = \phi z \ge 0$ when $y = 0$ and $\dot{z} = \varepsilon \beta \frac{sy}{y+\alpha s} \ge 0$ when $z = 0$, the first orthant is invariant, so that all solutions of system (1) remain in $\mathbb{R}_+^3$ for any initial condition $x_0 \in \mathbb{R}_+^3$.

Using the nonnegativity of the state variables of system (1), one first gets:

$$ \dot{s} \le \Lambda - \mu s \qquad (2) $$

so that, when $s = \frac{\Lambda}{\mu}$, $\dot{s} \le 0$. We also have:

$$ \varepsilon \dot{s} + \dot{z} \le \varepsilon \Lambda - \varepsilon \mu s - \mu_z z \le \varepsilon \Lambda - \zeta (\varepsilon s + z) \qquad (3) $$
so that, when $\epsilon s + z = \frac{\epsilon\Lambda}{\zeta}$, then $\epsilon\dot{s} + \dot{z} \le 0$. Finally, we have:

$$ \dot{y} \le \phi z - \mu_y y \le \phi \frac{\epsilon\Lambda}{\zeta} - \mu_y y, \quad (4) $$

when $z \le \frac{\epsilon\Lambda}{\zeta}$, so that, when $y = \frac{\epsilon\phi\Lambda}{\zeta\mu_y}$, one has $\dot{y} \le 0$. Therefore, the compact set $\Gamma$ is positively invariant. Thus, by the classical results of dynamical systems, for any initial condition in $\mathbb{R}_+^2 \setminus \{(0, 0)\} \times \mathbb{R}_+$ the Cauchy problem associated with system (1) has a unique solution in $\mathbb{R}_+$.

## 3.2 | Asymptotic behaviour

In this section, we study the existence and stability of the model equilibria. Although the model is only valid over a cropping season $[0, t_f]$, such an asymptotic analysis gives an idea of the system behaviour during this finite period, especially if the transient dynamics are fast, which can be expected with CBB that very efficiently colonise plantations. System (1) has a pest-free steady state $\mathcal{E}^0 = (s^0, 0, 0)$ with $s^0 = \frac{\Lambda}{\mu}$. We compute the basic reproduction number, $\mathcal{R}_0$, using the next generation operator approach [28]. Let $F = [\phi z, 0]$ be the emergence of CBB and $V = [\frac{\epsilon\beta_{sy}}{y+\alpha s} + \mu_y y, -\frac{\epsilon\beta_{sy}}{y+\alpha s} + \mu_z z]$ the outflows minus the inflows of the CBB compartments, emergence excepted. The Jacobian matrices of $F$ and $V$ at the pest-free steady state $\mathcal{E}^0$ are respectively given by:

$$ F = \begin{bmatrix} 0 & \phi \\ 0 & 0 \end{bmatrix} \quad \text{and} \quad V = \begin{bmatrix} \epsilon\beta\frac{1}{\alpha} + \mu_y & 0 \\ -\epsilon\beta\frac{1}{\alpha} & \mu_z \end{bmatrix}. $$

The basic reproduction number is defined as the spectral radius of the next generation matrix $FV^{-1}$:

$$ \mathcal{R}_0 = \rho(FV^{-1}) = \frac{\phi\epsilon\beta\frac{1}{\alpha}}{\mu_z(\epsilon\beta\frac{1}{\alpha} + \mu_y)}. \quad (5) $$

**Remark 1.** The biological interpretation of the basic reproduction number $\mathcal{R}_0$ is the following. The factor $\frac{\epsilon\beta\frac{1}{\alpha}}{\epsilon\beta\frac{1}{\alpha} + \mu_y}$ is the average number of infesting females produced by a colonising female during its lifespan, whereas the factor $\frac{\phi}{\mu_z}$ is the average number of colonising females produced by an infesting female. Then the basic reproduction number $\mathcal{R}_0$, corresponding to the product of these two factors, represents the average number of new infesting (or colonising) females originated from a single infesting (or colonising) female during its lifespan.

Based on the basic reproduction number, we obtain the following stability results.

**Theorem 1.** 1. The pest-free steady state $\mathcal{E}^0$ is globally asymptotically stable (GAS) if $\mathcal{R}_0 \le 1$.

2. When $\mathcal{R}_0 > 1$, there exists a unique coexistence steady state $\bar{\mathcal{E}} = (\bar{s}, \bar{y}, \bar{z})$ for system (1) where $\bar{s}, \bar{y}$ and $\bar{z}$ are defined as follows:

$$ \bar{s} = \frac{\Lambda\left(\frac{1}{\mu_y}(\epsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \alpha\right)}{\frac{1}{\mu_y}(\mu + \beta)(\epsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \mu\alpha}; \quad \bar{y} = \frac{1}{\mu_y}(\epsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1)\bar{s}; \quad \bar{z} = \frac{\epsilon\beta}{\mu_z}\frac{\bar{y}\bar{s}}{\bar{y} + \alpha\bar{s}}. \quad (6) $$

Moreover, this unique steady state $\bar{\mathcal{E}}$ is locally asymptotically stable.

*Proof.* See Appendix A.

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

Theorem 1 has a biological meaning: the CBB population disappears from the coffee plantation if the basic reproduction number $\mathcal{R}_0 \le 1$, so no control is needed, except to hasten the CBB decline. In contrast, when $\mathcal{R}_0 > 1$, the CBB population persists in the coffee plantation. Hence, control measures should aim at modifying model parameters such that $\mathcal{R}_0$ is brought below 1, in order to reduce the damages and the CBB population at the end of the cropping season.

## 3.3 | Time-dependent effective reproduction number

The basic reproduction number gives information on the onset of the epidemic in an entirely susceptible population. The status of an epidemic in a partially infected population is given by the effective reproduction number $\mathcal{R}_e(t)$, whose value compared
to 1 indicates if the epidemic is expanding, endemic or contracting. At time t, it is defined as the average number of secondary cases per infectious case present at time t. A number of formulas have been proposed to approximate $R_e(t)$, but they mostly require that the total population of susceptible and infected is constant, and are defined in the case where infection is directly transmitted from an infectious individual to a susceptible [29]. Neither is satisfied in our case, as the host population strongly increases during the season and as the infection is vector-transmitted.

Therefore, we propose an exact formula for $R_e(t)$. Within our model, we can compute the number of secondary colonising females produced in average by a colonising female present at time t (the same can be done for infesting females). We identify the number of colonising females present at time t and follow their decrease, through mortality and transfer to the infesting female compartment. These infesting females are then computed and used to calculate the secondary colonising females. Even though the actual season is limited to a length $t_f$, we evaluate $R_e(t)$ by computing all the secondary colonising females generated over an infinite horizon. Doing so, we include those that would have been produced after the end of the season, so that $R_e(t)$ is not influenced by the proximity to the end of the actual season. The detailed calculation is given in Appendix B and leads to the following expression:

$$R_e(t) = \int_{t}^{+\infty} \int_{t}^{\tau} \phi \varepsilon \beta \frac{s(\hat{t})}{y(\hat{t}) + \alpha s(\hat{t})} e^{-\varepsilon \beta \int_{t}^{\hat{t}} \frac{s(\tau)}{y(\tau) + \alpha s(\tau)} d\tau - \mu_y(\hat{t} - t) - \mu_z(\tau - \hat{t})} d\hat{t} d\tau \quad (7)$$

We note that, when $y(0) \approx 0$ and $s(0)$ is at the pest-free equilibrium, $\frac{s(t)}{y(t)+\alpha s(t)} \approx \frac{1}{\alpha}$ at all times, so that $R_e(0) = R_0$. Moreover, at the endemic equilibrium (6), $R_e = 1$.

An approximation of the time-dependent effective reproduction number may be obtained from the basic reproduction number (5), by replacing $\frac{1}{\alpha}$ by $\frac{s(t)}{y(t)+\alpha s(t)}$:

$$R(t) = \frac{\phi \varepsilon \beta \frac{s(t)}{y(t)+\alpha s(t)}}{\mu_z \left( \varepsilon \beta \frac{s(t)}{y(t)+\alpha s(t)} + \mu_y \right)} = \frac{\phi \varepsilon \beta}{\mu_z} \frac{s(t)/y(t)}{(\varepsilon \beta + \mu_y \alpha) s(t)/y(t) + \mu_y}, \quad (8)$$

where $s(t)$ and $y(t)$ are the solutions of system (1). Assuming that the population does not vary too much, which determines the validity of this approximation, the time-dependent approximated reproduction number can be interpreted as the basic reproduction number at time t, that is the product of the average number of infesting females produced by a colonising female during its lifespan $\frac{\varepsilon \beta \frac{s(t)}{y(t)+\alpha s(t)}}{\varepsilon \beta \frac{s(t)}{y(t)+\alpha s(t)} + \mu_y}$ and the average number of colonising females produced by an infesting female $\frac{\phi}{\mu_z}$. As for the exact effective reproduction number, at the pest-free equilibrium $R = R_0$ and at the endemic equilibrium $R = 1$.

# 4 | OPTIMAL CONTROL STRATEGY

## 4.1 | Problem statement

In this section, we formulate a continuous optimal control problem which consists in maximising the profit, by means of an entomopathogenic fungus.

Entomopathogenic fungi grow at the expense of insects and usually kill them. *Beauveria bassiana* is the most studied fungus to control coffee pests, especially CBB [30, 31, 32, 33, 34, 35, 36]. This fungus has the advantage of not being harmful to humans, of having no negative impact on coffee berries, and of being infective through a simple contact [8, 31, 32, 33, 34]. It is sprayed on the surface of the coffee berries and it infects the colonising CBB females when they begin to drill an entry hole into a coffee berry. The fungus then rapidly grows, feeding on the nutrients in the CBB body and producing toxins in the process. When the CBB dies after 3 to 5 days, *B. bassiana* covers the carcass in a layer of white mold [30]. Abiotic factors limit the efficacy of the fungus on the plant canopy, leading to a progressive decay of the inoculum [8].

Let us denote by $v(t)$ the fungus load at time t. The dynamics of the entomopathogenic fungus on the berries is:

$$\dot{v}(t) = -\gamma v(t) + h(t),$$

with initial condition $v(0) \in \mathbb{R}_+$. Parameter $\gamma$ denotes the natural decay rate of the fungus. The function $h(t)$ represents the fungus application rate at time t, with an upper bound $h_M > 0$. The effect of the fungus is limited and saturates at $\xi \in (0, 1)$:

$$\sigma(v) = \frac{\xi v}{v + k},$$
So we obtain the following controlled system:

$$ 
\left\{
\begin{aligned}
\dot{s} &= \Lambda - (1 - \sigma(v))\beta \frac{sy}{y + \alpha s} - \mu s \\
\dot{y} &= \phi z - \varepsilon \beta \frac{sy}{y + \alpha s} - \mu_y y \\
\dot{z} &= (1 - \sigma(v))\varepsilon \beta \frac{sy}{y + \alpha s} - \mu_z z \\
\dot{v} &= -\gamma v + h(t)
\end{aligned}
\right.
\qquad (9) $$

with initial conditions $(s(0), y(0), z(0), v(0)) \in \mathbb{R}_+^4$.

The set $\mathcal{U}$ of admissible control functions is defined as:

$$ \mathcal{U} = \{h \in L^{\infty}(0, t_f) : 0 \le h(t) \le h_M, \forall t \in [0, t_f]\}. \qquad (10) $$

The well-posedness of the controlled system, which is a biological requirement, is handled by the following proposition.

**Proposition 2.** For all nonnegative initial conditions $(s(0), y(0), z(0), v(0)) \in \mathbb{R}_+^4$, the solutions $(s(t), y(t), z(t), v(t))$ of system (9) remain nonnegative and are bounded for all time $t > 0$. Moreover, the bounded set

$$ \Omega = \Gamma \times \left\{ v \in \mathbb{R}_+ : v(t) \le \frac{h_M}{\gamma} \right\} $$

is positively invariant.

*Proof.* Since $h$ is a nonnegative function, $v$ remains nonnegative at all times for any nonnegative initial condition $v(0)$. Then, as $\sigma(v) \le 1$ and proceeding as in the proof of Proposition 1, the other variables also remain nonnegative and bounded. This implies that all trajectories stay in the first quadrant for all initial conditions in $\mathbb{R}_+^4$. Moreover, $\Gamma$ also remains positively invariant for the three first state variables. According to the standard comparison principle [37] and using the fact that $h(t) \le h_M$, $v(t)$ is positively invariant in $\{v \in \mathbb{R}_+ / v(t) \le \frac{h_M}{\gamma}\}$. Therefore, system (9) is positively invariant in $\Omega$. $\square$

Our primary goal is to maximise the profit generated by the berry yield, by applying an entomopathogenic fungus to control the CBB population. An option would be to bring the CBB population down as fast as possible in order to increase the yield, but we want to keep the control cost in check. Furthermore, as CBB survive between cropping seasons, we also want to reduce the CBB population at the end of the season considered, so as not to unduly burden the next season yield. So we aim at maximising the berry yield, while minimising the control cost as well as the CBB population at the end of the cropping season. We assume that the harvest takes place at the end of the cropping season, denoted by $t_f$. We also assume that almost all berries, and hence almost all infesting females $z$ that are inside berries, are picked at $t_f$, so that only colonising females $y$ remain in the plantation after harvest. The application cost $C$ of a fungus unit is considered constant and only healthy coffee berries are sold, at a constant berry price per berry $\xi_s$. The optimal control problem is then formulated as follows.

**Problem 1.** Find an admissible fungus application strategy $h^*(.) \in \mathcal{U}$ and the corresponding state variable $x^*(.) = (s^*(.), y^*(.), z^*(.), v^*(.))^T$ maximising the following objective functional:

$$ \mathcal{J}(x, h) = \xi_s s(t_f) - \int_0^{t_f} Ch(t)dt - \xi_y y(t_f), \qquad (11) $$

subject to the control system (9) and with initial condition in $\mathbb{R}_+^4$.

This type of problem is called a Bolza optimal control problem [21]. The first term of the objective functional represents the crop yield, that is the harvest of healthy coffee berries at the end of cropping season multiplied by the coffee berry price $\xi_s$. The second term represents the control cost, with $C$ the price of a fungus unit. The first two terms hence represent the financial gain $G$. Finally, the last term represents the CBB population that remains after harvest in the plantation, consisting of colonising females, weighted by parameter $\xi_y$. So $\mathcal{J}$ represents the profit, penalised by the final CBB population.
Let us consider the vectors $x(t) = (s(t), y(t), z(t), v(t))^T$ and $e_4 = (0, 0, 0, 1)^T$. We define $f: \mathbb{R}^4 \to \mathbb{R}^4$ as follows:

$$f(x) = \begin{pmatrix} \Lambda - (1 - \sigma(\nu))\beta \frac{sy}{y+\alpha s} - \mu s \\ \phi z - \varepsilon \beta \frac{sy}{y+\alpha s} - \mu_y y \\ (1 - \sigma(\nu))\varepsilon \beta \frac{sy}{y+\alpha s} - \mu_z z \\ -\gamma v \end{pmatrix}$$

Then the model (9) can be rewritten in the following compact form:

$$\dot{x}(t) = f(x(t)) + h(t)e_4.$$

## 4.2 | Necessary optimal conditions

Based on Corollary III.4.1 of Theorem III.4.1 in Fleming and Rishel [38] (see also [39]), we establish the existence of an optimal strategy that maximises the objective functional (11) in the following theorem.

**Theorem 2.** The optimal control (1) admits a solution, i.e. there exists an optimal strategy $h^*(.) \in \mathcal{U}$ with associated optimal solution $x^* = (s^*, y^*, z^*, v^*)$ which maximises the objective functional $\mathcal{J}(x, h)$.

*Proof.* Corollary III.4.1 in Fleming and Rishel [38] imposes three conditions:

1. $\mathcal{G}$, denoting the class of all $(x(0), h)$ with $x(0) \in \mathbb{R}_+^4$ and $h \in \mathcal{U}$, such that the corresponding solution of system (9) is defined over $[0, t_f]$, is nonempty.

2. The admissible control set $\mathcal{U}$ is closed and convex.

3. The integrand of objective functional $\mathcal{J}(x, h)$ is concave on $\mathcal{U}$ and bounded above by $-C_1|h|^\kappa + C_2$ for some $C_1 > 0$ and $\kappa > 1$.

Firstly, let us choose some $x(0) \in \mathbb{R}_+^4$ and $h \in \mathcal{U}$. Since system (9) has bounded coefficients and any solution is bounded on the finite interval time $[0, t_f]$, then by the Caratheodory theorem (see for instance in [40, p. 82]), there exists a unique solution of system (9) with initial condition $x(0)$ and control $h$. This implies that the domain $\mathcal{G}$ is nonempty. By definition, the admissible set control $\mathcal{U}$ is closed and convex, hence the second condition is satisfied. For the third condition, it is easy to see that the integrand function of the objective functional is concave since it is linear in control $h(t)$. Moreover, taking $C_1 = \frac{C}{h_M}$, $\kappa = 2$ and $C_2 = 0$, the last condition is satisfied. Therefore, there exists an optimal strategy. □

We now apply the Pontryagin's Maximum Principle [21] to establish the first-order necessary optimality conditions of problem 1. Let $p \in \mathbb{R}^4$ and $p^0 \in \mathbb{R}$. The Hamiltonian associated to problem 1 is defined as:

$$\begin{aligned} \mathcal{H}(x, h, p, p^0) &= \langle p, f(x) + h e_4 \rangle - p^0 Ch \\ &= p_1 \left[ \Lambda - (1 - \sigma(\nu))\beta \frac{sy}{y + \alpha s} - \mu s \right] + p_2 \left[ \phi z - \varepsilon \beta \frac{sy}{y + \alpha s} - \mu_y y \right] \\ &\quad + p_3 \left[ (1 - \sigma(\nu))\varepsilon \beta \frac{sy}{y + \alpha s} - \mu_z z \right] + p_4 [-\gamma v + h] - p^0 Ch, \end{aligned} \quad (12)$$

where $\langle \cdot, \cdot \rangle$ represents the scalar product in $\mathbb{R}^4$.

According to the Pontryagin's Maximum Principle, if $h^*(.) \in \mathcal{U}$ is an optimal strategy and $x^*(.)$ is its associated trajectory, there exists an absolutely continuous vector mapping $p: [0, t_f] \to \mathbb{R}^4$ and a scalar $p^0 \ge 0$ that satisfy the following conditions:

(i). The vector $(p(\cdot), p^0)$ is non trivial.
(ii). The adjoint vector $p$ satisfies $\dot{p}(t) = -\nabla_x \mathcal{H}(x(t), h(t), p(t), p^0)$ for a.e. $t \in [0, t_f]$, that is:

$$ 
\left\{
\begin{aligned}
\dot{p}_1 &= [(1 - \sigma(v))(p_1 - \epsilon p_3) + \epsilon p_2] \frac{\beta y^2}{(y + \alpha s)^2} + \mu p_1 \\
\dot{p}_2 &= [(1 - \sigma(v))(p_1 - \epsilon p_3) + \epsilon p_2] \frac{\beta \alpha s^2}{(y + \alpha s)^2} + \mu_y p_2 \\
\dot{p}_3 &= -\phi p_2 + \mu_z p_3 \\
\dot{p}_4 &= (-p_1 + \epsilon p_3) \frac{d\sigma}{dv} \frac{\beta sy}{y + \alpha s} + \gamma p_4.
\end{aligned}
\right.
\qquad (13) $$

(iii). For a.e. $t \in [0, t_f]$:

$$ h^*(t) \in \underset{\omega \in [0, h_M]}{\operatorname{argmax}} \mathcal{H}(x(t), \omega, p(t), p^0). \qquad (14) $$

(iv). The following transversality conditions, obtained by deriving the terminal terms of the objective functional (11) relative to each system variable and applied at final point $t_f$, are verified:

$$ p_1(t_f) = \zeta_s p^0, \quad p_2(t_f) = -\zeta_y p^0, \quad p_3(t_f) = p_4(t_f) = 0. \qquad (15) $$

By definition, an extremal trajectory of the optimal control problem is a quadruplet $(x(\cdot), h(\cdot), p(\cdot), p^0)$ satisfying equations (9) and (12)–(15). An optimal solution is necessarily an extremal. If $p^0 = 0$, then we have $p(t_f) = 0$, which contradicts condition ((i)) above. Hence $p^0 \neq 0$, which can be normalised to $p^0 = 1$. The $p^0$ argument is then dropped from the Hamiltonian expression.

We now study the structure of the optimal control. As $\mathcal{H}$ is linear with respect to control $h(t)$, maximising $\mathcal{H}$ depends on the sign of

$$ \partial_h \mathcal{H}(x(t), h(t), p(t)) = -C + p_4(t). $$

So we introduce the switching function:

$$ \Phi(t) = -C + p_4(t), \quad \forall t \in [0, t_f], $$

and characterise the optimal control $h^*$ by:

$$ h^*(t) = \begin{cases} 0 & \text{if } \Phi(t) < 0 \\ \text{undefined} & \text{if } \Phi(t) = 0 \\ h_M & \text{if } \Phi(t) > 0. \end{cases} $$

If $h^*(t)$ is equal to 0 or $h_M$ for a.e. $t \in [0, t_f]$, the optimal solution is called a *bang-bang* solution.

**Remark 2.** At the end of cropping season, according to equation (15), $\Phi(t_f) = -C < 0$, so the control vanishes, that is $h(t_f) = 0$.

## 4.3 | Singular control

A control is called singular if there exists a nonempty interval $[t_1, t_2] \subset [0, t_f]$ such that $\Phi(t) = 0$ for every $t \in [t_1, t_2]$.

**Theorem 3.** If there exists a singular control $h_{\text{sing}}^*(t)$ over $[t_1, t_2]$, it is of order 1 and locally maximising. Moreover, it is given as a function of the state and adjoint variables by:

$$ h_{\text{sing}}^*(t) = -\frac{S(x, p)}{Q(x)}, $$

where

$$ 
\begin{aligned}
S(x, p) &= \left[ -(1 - \sigma(v))\gamma C - \epsilon p_2 \frac{d\sigma}{dv} \frac{\beta sy}{y + \alpha s} \right] \frac{\beta y^2}{(y + \alpha s)^2} - \frac{2\gamma^2 C v}{v + k} - \mu \gamma C \\
&\quad + \frac{d\sigma}{dv} \frac{\beta sy}{y + \alpha s} \left[ -\epsilon\phi p_2 + \epsilon(\mu_z - \mu)p_3 \right] - \frac{\gamma C y}{s(y + \alpha s)} \dot{s} - \frac{\gamma C \alpha s}{y(y + \alpha s)} \dot{y}, \\
Q(x) &= \frac{2\gamma C}{v + k}.
\end{aligned}
$$
Proof. Assuming that a singular control exists for all $t \in [t_1, t_2]$, we have $\Phi(t) = 0$. Then differentiating the switching function with respect to time $t$ along the singular control and adjoint system (13), we obtain:

$$ \dot{\Phi}(t) = \dot{p}_4(t) = (-p_1 + \varepsilon p_3) \frac{d\sigma}{dv} \beta \frac{sy}{y + \alpha s} + \gamma p_4(t), $$

which does not depend explicitly on control $h$. So we differentiate $\Phi$ a second time:

$$ \begin{aligned} \ddot{\Phi}(t) = & \left[ (1 - \sigma(v))(-p_1 + \varepsilon p_3) \frac{d\sigma}{dv} \beta \frac{sy}{y + \alpha s} - \varepsilon p_2 \frac{d\sigma}{dv} \frac{sy}{y + \alpha s} \right] \beta \frac{y^2}{(y + \alpha s)^2} \\ & + [-\mu p_1 - \varepsilon \phi p_2 + \mu_z p_3] \frac{d\sigma}{dv} \beta \frac{sy}{y + \alpha s} - (-\gamma v + h) \frac{2}{v + k} (-p_1 + \varepsilon p_3) \frac{d\sigma}{dv} \beta \frac{sy}{y + \alpha s} \\ & + (-p_1 + \varepsilon p_3) \frac{d\sigma}{dv} \beta \frac{sy}{y + \alpha s} \left( \dot{s} \frac{y}{s(y + \alpha s)} + \dot{y} \frac{\alpha s}{y(y + \alpha s)} \right) + \gamma \dot{p}_4(t). \end{aligned} $$

Then, as $\Phi(t) = 0$ and $\dot{\Phi}(t) = 0$ for all $t \in [t_1, t_2]$, we obtain:

$$ p_4(t) = C, \quad (-p_1 + \varepsilon p_3) \frac{d\sigma}{dv} \beta \frac{sy}{y + \alpha s} = -\gamma C \quad \forall t \in [t_1, t_2]. $$

Substituting the equations above into the expression of $\ddot{\Phi}$, we obtain:

$$ \ddot{\Phi}(t) = S(x, p) + hQ(x). $$

As $\ddot{\Phi}(t) = 0$ along the singular trajectory, we obtain the singular control expression:

$$ h_{sing}^{*}(t) = - \frac{S(x, p)}{Q(x)}. $$

The number of derivatives of the switching function required to determine the singular control is 2, which implies that the singular control is of order 1. For the singular control to be maximising, the classical Legendre Clebsch condition [39, 41] requires that:

$$ (-1) \partial_h \left[ \frac{d^2}{dt^2} \partial_h H(x(\cdot), h(\cdot), p(\cdot)) \right] \le 0. $$

Since $\frac{d^2}{dt^2} \partial_h H(x(\cdot), h(\cdot), p(\cdot)) = \ddot{\Phi}(t) = S(x, p) + hQ(x)$, we have $(-1) \partial_h [\ddot{\Phi}(t)] = -Q(x) \le 0$, because the state variable $v$ is nonnegative for all time $t \ge 0$. Then the singular control $h_{sing}^*(\cdot)$ is a locally maximising strategy. $\square$

We can summarise the results above in the following theorem.

**Theorem 4.** Let $h^* \in \mathcal{V}$ be an optimal control that maximises the objective functional (11). Then, there exists an adjoint vector function $p(t) \in \mathbb{R}^4$ satisfying (13) with transversality conditions (15). Moreover, the structure of the optimal control is given by:

$$ h^*(t) = \begin{cases} 0 & \text{if } \Phi(t) < 0, \\ h_{sing}^*(t) & \text{if } \Phi(t) = 0 \text{ for } t \in [t_1, t_2], \\ h_M & \text{if } \Phi(t) > 0. \end{cases} \qquad (16) $$

# 5 | NUMERICAL SIMULATIONS

In this section, we investigate numerically the effect of the optimal strategy on the infestation process of coffee berries by CBB.

## 5.1 | Parameter values

Some parameters, particularly those associated with CBB, are estimated from the literature and expert opinions. For instance, the mortality of females is chosen such that $1/\mu_y + 1/\mu_z$ plus the larval development time of CBB lies within the interval 87–282 days, which represents the female lifespan [8]. The emergence rate of females is estimated from the average number of eggs laid by females, which is one to three egg per day [8]. The entomopathogenic fungus is supposed to be viable up to a few months. Fungus cost and application rate are derived from biocontrol products on the market. For the small family plantation considered in this study (0.12 ha), the recommended application rate ranges between 60 and 180 g. To translate this amount into a daily maximal rate, we make the hypothesis that there will be no more than one application per week, so the maximal daily
rate is a little less than the value we conservatively chose. The coffee cost is derived from market values. The duration of a coffee cropping season is taken from [8]. The remaining parameter values are assumed. All parameter values are given in Table 1. With these parameter values, the basic reproduction number $R_0 = 12$ for the model without control described in system (1), which is consistent with the fast CBB dynamics observed in the field. However, we illustrate the optimal control strategy when $R_0 < 1$ and the pest-free equilibrium is stable in Appendix C.1.

<table><caption>TABLE 1 Model and control parameters.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>Λ</td><td>Production rate of new coffee berries</td><td>1200 berries.day<sup>-1</sup></td></tr><tr><td>μ</td><td>Natural mortality rate of healthy coffee berries</td><td>0.002 day<sup>-1</sup></td></tr><tr><td>φ</td><td>Emergence rate of new colonising females</td><td>2 day<sup>-1</sup></td></tr><tr><td>β</td><td>Infestation rate</td><td>0.01 day<sup>-1</sup></td></tr><tr><td>α</td><td>Infestation half saturation constant</td><td>0.7 female.berry<sup>-1</sup></td></tr><tr><td>ε</td><td>Conversion rate from coffee berries to CBB</td><td>1 female.berry<sup>-1</sup></td></tr><tr><td>μ<sub>y</sub></td><td>Natural mortality rate of colonising females</td><td>1/20 day<sup>-1</sup></td></tr><tr><td>μ<sub>z</sub></td><td>Natural mortality rate of infesting females</td><td>1/27 day<sup>-1</sup></td></tr><tr><td>ξ</td><td>Fungus load maximal effectiveness</td><td>0.8</td></tr><tr><td>k</td><td>Fungus load half saturation constant</td><td>200 g.day<sup>-1</sup></td></tr><tr><td>γ</td><td>Fungus decay rate</td><td>1/50 day<sup>-1</sup></td></tr><tr><td>ζ<sub>s</sub></td><td>Price of coffee berries</td><td>2.10<sup>-3</sup> US$.berry<sup>-1</sup></td></tr><tr><td>ζ<sub>y</sub></td><td>Cost of remaining colonising females</td><td>10<sup>-4</sup> US$.female<sup>-1</sup></td></tr><tr><td>C</td><td>Cost of fungus</td><td>0.022 US$.g<sup>-1</sup></td></tr><tr><td>h<sub>M</sub></td><td>Maximal fungus application rate in the plantation</td><td>30 g.day<sup>-1</sup></td></tr><tr><td>t<sub>f</sub></td><td>Duration of a cropping season</td><td>250 days</td></tr></tbody></table>

## 5.2 | Optimal control

To solve our optimal control problem numerically, we use the BOCOP software [42, 43], an open source toolbox in which the optimal control problem is approximated by a finite dimensional optimisation problem. In BOCOP we choose the Lobatto discretisation method (implicit, 4 stages, order 6), 2000 time steps, and a $10^{-5}$ tolerance.

We use the following initial conditions. At the beginning of a cropping season, there are no berries, so $s(0) = 0$. Moreover, we assume that there are only colonising females, as infesting females were eliminated from the plantation when berries were manually removed during the preceding harvest, so $z(0) = 0$. We also assume that the entomopathogenic fungus is absent in the plantation, so the initial load $v(0) = 0$. We consider several initial conditions for the colonising females: $y(0) = 10^2$ in Figure 2, a fairly low number of females corresponding to a new infestation of the plantation; $y(0) = 10^6$ in Figure 4, which is close to the coexistence equilibrium value without control and hence corresponds to a worst case scenario, where all females have survived since the end of the preceding cropping season; and $y(0) = 10^4$ in Figure 3, which is an intermediate value. Moreover, initial conditions $y(0) = 10^3$ in Figure C2 and $y(0) = 10^5$ in Figure C3 are presented in Appendix C.2.

There are six panels on each figure. Panel (a) represents the evolution of the healthy coffee berries $s(t)$ during the cropping season, without CBB (black), with initial CBB in the plantation but no control (blue), and with initial CBB and optimal fungus application (red). Similarly, panels (b) and (c) represent the dynamics of the colonising $y(t)$ and infesting $z(t)$ CBB females, respectively, without and with optimal control. Panel (e) represents the optimal fungus application $h(t)$, the control that maximises the penalised profit $J$ defined in equation (11). Panel (d) represents the dynamics of the fungus load $v(t)$ in the plantation and panel (e) the effect of this load, that saturates through function $\sigma(v) = \frac{\xi v}{v+k}$. Finally in panels (a-c), the horizontal dotted blue lines represent the coexistence equilibrium without control.
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal control (16) (plain red curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^2$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve). The coexistence steady state without control is indicated in (a–c) (dotted blue lines).</figcaption></figure>

In Figures 2 (b) and 3 (b), the trajectories of the colonising females increase, both in the uncontrolled (dashed blue curve) and controlled (red curve) cases. In contrast, for higher initial values in Figure 4 (b) and to a lesser extent Figure C3 (b) in Appendix C.2, the colonising female trajectories first decrease, due to the lack of healthy coffee berries at the beginning of the season, and later increase.

For all initial conditions, at the end of the cropping season, the uncontrolled solution (dashed blue curves in subfigures (a–c)) is very close to the coexistence steady state: $\bar{s} = 1.06 \times 10^5$, $\bar{y} = 10.5 \times 10^5$, $\bar{z} = 0.27 \times 10^5$ (dotted blue horizontal lines). Since the basic reproduction number $R_0 = 12 > 1$, according to Theorem 1, the CBB population persists in the plantation and the coexistence steady state is locally stable. Without control, the berry production at the end of the cropping season is more than halved, compared to the pest-free case (black curves).

The optimal control depends on the initial condition. In Figures 2 (e) and 3 (e) the control structure is *bang(30)*–*singular-bang(0)*, whereas in Figure 4 (e) it is *bang(0)*–*singular-bang(0)*; both are consistent with the theoretical results in Remark 2 and Theorem 3. In the latter case, the colonising female population is initially higher and first decreases. As a consequence, the fungus is first applied later than for lower initial conditions. Interestingly, the switch between the *bang(30)* and the *singular* control occurs earlier as the initial number of colonising females increases, as shown in Figures 2 (e), C2 (e), and 3 (e). In parallel, the switch between the *bang(0)* and the *singular* control occurs later as the initial number of colonising females increases, as shown in Figures C3 (e) and 4 (e). In all cases, more effort is devoted earlier to control the borers when the initial population is lower, which allows to maintain the population at lower levels. The control pay-off is higher for lower pest populations, which justifies more important control efforts. When the CBB initial population is high, as in Figures 4 and C3, the CBB population first decreases, so there is initially no control effort. In all cases, once the fungus is applied, its effect (subfigures (f)) rapidly increases but remains below its maximum value ($\xi = 0.8$). The colonising (subfigures (b), red curves) and infesting (subfigures
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal control (16) (plain red curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^4$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve). The coexistence steady state without control is indicated in (a–c) (dotted blue lines).</figcaption></figure>

<table><caption>**TABLE 2** Penalised profit, financial gain, and colonising CBB females remaining in the plantation after harvest at the end of the cropping season ($t_f = 250$ days), without control and with optimal fungus application, for several initial abundances of colonising females.</caption><thead><tr><th rowspan="2">y(0)</th><th colspan="3">Penalised profit J (US$)</th><th colspan="3">Financial gain G (US$)</th><th colspan="3">Remaining CBB y(t<sub>f</sub>) (×10<sup>5</sup>)</th></tr><tr><th>No control</th><th>Optimal</th><th>Increase<sup>†</sup></th><th>No control</th><th>Optimal</th><th>Increase<sup>†</sup></th><th>No control</th><th>Optimal</th><th>Decrease<sup>†</sup></th></tr></thead><tbody><tr><td>10<sup>2</sup></td><td>108</td><td>296</td><td>174%</td><td>206</td><td>320</td><td>55%</td><td>9.79</td><td>2.38</td><td>76%</td></tr><tr><td>10<sup>3</sup></td><td>108</td><td>233</td><td>116%</td><td>204</td><td>275</td><td>35%</td><td>9.58</td><td>4.22</td><td>56%</td></tr><tr><td>10<sup>4</sup></td><td>108</td><td>207</td><td>92%</td><td>203</td><td>256</td><td>26%</td><td>9.48</td><td>4.85</td><td>49%</td></tr><tr><td>10<sup>5</sup></td><td>107</td><td>199</td><td>86%</td><td>201</td><td>248</td><td>23%</td><td>9.42</td><td>4.86</td><td>48%</td></tr><tr><td>10<sup>6</sup></td><td>107</td><td>195</td><td>82%</td><td>201</td><td>243</td><td>21%</td><td>9.37</td><td>4.81</td><td>49%</td></tr></tbody></table>

Penalised profit/financial gain in the pest-free case: $J = G = 472$ US$.

†Relative increase and decrease compared to the case without control.

(c), red curves) female trajectories (subfigures (c), red curves) are hence much lower than in the case without control (dashed blue curves).

Table 2 summarises the plantation state at the end of the cropping season (for $t_f = 250$ days), with and without optimal control, depending on the initial conditions. As we assume that all berries have been harvested, only colonising females (y) remain in the plantation. The financial gain $G$ corresponds to the coffee berry yield, minus the control cost (when applicable).
<figure><img src="image_6.png" /><figcaption>**FIGURE 4** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal control (16) (plain red curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^6$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve). The coexistence steady state without control is indicated in (a–c) (dotted blue lines).</figcaption></figure>

The penalised profit $J$, defined in equation (11), corresponds to the coffee berry yield, minus the control cost (when applicable) and the remaining CBB “cost”, as we also aim at reducing the CBB population at the end of the cropping season, so as not to unduly burden the next season yield. This partly explains why the penalised profit in the pest-free case is much higher than when CBB are present in the plantation: the remaining CBB “cost” is around 95 US$ without control, and 49 US$ with optimal control. The control cost is around 60–80 US$ and also impacts the profit. At the end of the cropping season, regardless the initial condition, the remaining CBB population, the financial gain, and the penalised profit are very similar in the case without control. It is expected, as the CBB dynamics are very fast and trajectories get very close to the coexistence equilibrium before the end of the cropping season. In the case with control, there are notable differences. The financial gain and the penalised profit decrease with the initial CBB population size, especially for the lower values ($y(0) = 10^2$ to $10^4$ females). For higher values ($y(0) = 10^4$ to $10^6$ females), the lack of susceptible berries at the beginning of the cropping season regulates the number of CBB, so the impact of the initial number of CBB is limited. Remaining CBB population, financial gain, and penalised profit vary when the entomopathogenic fungus is applied, compared to the case without control: the remaining CBB population is more or less half as abundant, the monetary gain is notably improved and the penalised profit is at worst almost twice as high (82% increase) and at best almost three times higher (174% increase).

## 5.3 | Effective reproduction numbers

We also compute the exact (7) and approximated (8) time-dependent effective reproduction numbers defined for system (1) without control. For system (9) with optimal control (16), the approximated time-dependent reproduction number is similarly
<figure><img src="image_7.png" /><figcaption>FIGURE 5 Exact (dashed curves) and approximated (plain curves) time-dependent reproduction numbers for system (1) without control (blue curves) and system (9) with optimal control (16) (red curves), for initial conditions $s(0) = z(0) = v(0) = 0$ and various values of $y(0)$. In the case without control, the exact reproduction number $R_c(t)$ (dashed blue line) is computed according to equation (7), the approximated number $R(t)$ (plain blue line) corresponds to equation (8). In the case with optimal control the approximated number $R^*(t)$ (plain red line), corresponds to equation (17).</figcaption></figure>

defined as follows:

$$R^*(t) = \frac{\phi \varepsilon \beta [1 - \sigma(v(t))] \frac{s(t)}{y(t) + \alpha s(t)}}{\mu_z \left( \varepsilon \beta \frac{s(t)}{y(t) + \alpha s(t)} + \mu_y \right)} = \frac{\phi \varepsilon \beta [1 - \sigma(v(t))] \frac{s(t)}{y(t)}}{(\varepsilon \beta + \mu_y \alpha) \frac{s(t)}{y(t)} + \mu_y}, \quad (17)$$

where $s(t)$, $y(t)$ and $v(t)$ are the solutions of system (9) with optimal control (16). Note that $R^*(t)$, just as $R(t)$ in equation (8), increases with the healthy berries per capita ratio $s(t)/y(t)$. The exact time-dependent reproduction number cannot be defined in the case with control. Indeed, to compute the average number of colonising females produced by a colonising female present at time $t$, we would need to know the optimal control strategy after the end of the cropping season (to estimate the offspring of the CBB that appear close to the end of the season), but the control is not defined. Results are shown in Figure 5.

In the case without control, the exact and approximated time-dependent reproduction numbers (blue curves) are reasonably close. As there are no berries at initial time $t = 0$, the approximated numbers (plain blue curves) are equal to zero, but not the exact numbers as colonising females are present. With few initial colonising females in Figure 5 (a), which corresponds to Figure 2, shortly after $t = 0$, both exact and approximated numbers are close to the basic reproduction number $R_0 = 12$. Indeed, as healthy berries (Figure 2 (a)) grow faster than colonising females (Figure 2 (b)), the plantation is close to a pest-free state, with few CBB compared to healthy berries. The time-dependent reproduction numbers then decrease as colonising females become abundant compared to healthy berries. This decrease is faster for the exact calculation, as the female offspring of a colonising female present at time $t$ will be borne at a later and less favourable time in terms of healthy berries per capita ($s(t)/y(t)$). Both time-dependent reproduction numbers tend towards 1, as the system gets close to the endemic equilibrium, in which the plantation can only sustain the replacement of a colonising female by a single colonising female.

With more initial colonising females in Figure 5 (b–c), which corresponds to Figures 3 and 4, both exact and approximated time-dependent reproduction numbers first increase but remain much lower than the basic reproduction number $R_0 = 12$. Indeed, the system is quite far from the pest-free state, especially for the higher value of initial colonising females (Figure 5 (c)). In this case, the exact time-dependent reproduction number starts below the threshold value 1, as there are not enough berries to sustain the very high CBB population, which first declines (Figure 4 (b)). For both initial conditions, the time-dependent reproduction numbers increase, as healthy berries become available and then decreases towards 1, as the system gets close to the endemic equilibrium.

The optimal control lowers the peaks of the approximated effective reproduction number (red curves). At the beginning of the simulation, the effective reproduction numbers with control are lower than without control, thanks to the control term $[1 - \sigma(v(t))]$ in equation (17) which represents the effect of fungus application. However, the population dynamics are fairly different with and without control: the healthy berries per capita ratio is much larger with control, so the effective reproduction numbers with control decrease less. With few initial colonising females, the effective reproduction number with control remains notably higher than 1 (Figure 5 (a)), as healthy berries per capita are still high thanks to the control efficiency; the system with control is fairly
close to the system without pest (Figure 2). With many initial colonising females, the healthy berries per capita are particularly low (Figure 4), hence the very low effective reproduction number with control (Figure 5 (c)); the increase at the end is due to the interruption of fungus application.

## 6 | CONCLUSIONS

In this paper, we have developed a basic mathematical model to study the infestation dynamics of coffee berries by CBB in a plantation during a cropping season. This model is based on the CBB life cycle and takes into account the availability of healthy coffee berries. Our theoretical analysis shows that the CBB dynamics varies according to a threshold parameter called the basic reproduction $R_0$: when $R_0 < 1$, the CBB population decreases until extinction; when $R_0 > 1$ the pest-free steady state becomes unstable and a locally stable coexistence steady state appears. In the latter case, the CBB population needs to be controlled.

To reduce the infestation in the plantation, we have designed an optimal control problem based on the use of an entomopathogenic fungus, which attacks colonising CBB as they drill into the berries. The objective of this control problem was to maximise the healthy coffee berry yield at the end of cropping season, while minimising the cost of control as well as the remaining CBB population for the next cropping season. We showed that an optimal control exists and that it can be characterised using the Pontryagin's Maximum Principle. In our numerical simulations, we obtained bang-singular-bang optimal controls that are consistent with our theoretical results. Moreover, we showed that applying the entomopathogenic fungus efficiently controls the CBB population, which is halved, and doubles the penalised profit, compared to the case without control. However, because of the control cost, the profit is still less than half of the pest-free profit (without control, the profit is only slightly more than a fifth of the pest-free profit).

Using synthetic pesticides instead of entomopathogenic fungi would reduce the control cost and might also entail a better short-term efficiency, but their extensive use leads to the development of resistance and has deleterious consequences of the environment and public health [9, 10]. Biopesticides such as *Beauveria bassiana* may be more costly to produce, but they do not cause such problems [8, 31, 32]. Moreover, coffee produced with biopesticides can be sold on the organic, where it would fetch premium prices that could compensate the higher entomopathogenic fungus cost. Finally, combining several control methods based on biological control and cultural practices [3, 8] could improve the efficiency and cost of CBB control.

In a future work, it would be interesting and more realistic to apply the entomopathogenic fungus at discrete times. However, we believe that for a start, the powerful tools of optimal control on the continuous system give us a first idea about the optimal strategy, depicted in Figures 2-4. To pursue this idea, we integrate on fixed time intervals the optimal fungus application obtained in the continuous case and apply the resulting cumulated fungus loads at the beginning of each corresponding interval. The total quantity of fungus applied and the control costs are hence similar in the continuous and discrete cases. With weekly fungus applications, the dynamics and hence the profit are very close for the discrete and optimal continuous controls (see Figure C4 in Appendix C.3). However, when the fungus is applied every two weeks, as shown in Figure 6, the discrete control (plain magenta curve) is more efficient than the optimal continuous control (plain red curve): the healthy berry production is higher (subfigure (a)) and the colonising female population is lower (subfigure (b)) at the end of the season. Therefore, the penalised profit with the discrete fungus application ($J = 242 US$) is greater than with the optimal continuous fungus application ($J = 207 US$). This is probably due to the higher loads present in the field when the fungus is applied (subfigure (d)). As no bounds were set on the discrete control, the fungus applications may be fairly large in this case (subfigure (e')). This suggests that discrete fungus applications may yield better results than optimal continuous control, provided that fairly large quantities can be applied at a time.

In this paper, we consider a single cropping season. Nevertheless, as CBB survive between cropping seasons, we incorporated in our optimisation criterion a term that penalises the CBB population at the end of the season, so as not to unduly burden the next season yield. The control we obtain is hence still relevant in a multi-seasonal perspective, as it does not completely disregard what happens after the cropping season. The next step is to extend our model beyond a single cropping season, using a semi-discrete formalism to account for discrete events such as harvest or control methods such as dropped berry removal between cropping seasons.
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal continuous control (16) (plain red curves) and with discrete control applied every 2 weeks (plain magenta curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^4$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve). The coexistence steady state without control is indicated in (a–b) (dotted blue lines).</figcaption></figure>

## ACKNOWLEDGMENTS

This work was supported by EPITAG, an Inria associated team part of the LIRIMA (https://team.inria.fr/epitag/), as well as the Collège doctoral régional de l'Afrique Centrale et des Grands Lacs "Mathématiques, Informatique, Biosciences et Géo-sciences de l'environnement" (AUF, French Embassy in Cameroon, IRD). We would like to thank Dr. Pierre Martinon for helpful discussions about the BOCOP software.

## Conflict of interest

The authors declare no potential conflict of interests.

**How to cite this article:** Y. Fotso Fotso, S. Touzeau, B. Tsanou, S. Bowong, and F. Grognard (2020), Modelling and optimal strategy to control coffee berry borer infestation, *Math Meth Appl Sci*.
# APPENDIX

## A PROOF OF THEOREM 1

### A.1 Global stability of pest-free steady state

In order to prove the global asymptotical stability of the pest-free equilibrium point, we first show attractivity in finite time of a set that is slightly larger than $\Gamma$ defined in Proposition 1. We will show that any solution initiated in $\mathbb{R}_+^3$ outside of

$$ \Gamma_{\delta} := \left\{ (s, y, z) \in \mathbb{R}_+^3 : s \le \frac{\Lambda}{\mu} + \delta, \varepsilon s + z \le \frac{\varepsilon \Lambda}{\zeta} + \delta, y \le \frac{\varepsilon \phi \Lambda}{\zeta \mu_y} + \frac{\phi \delta}{\mu_y} + \delta \right\} $$

for some $\delta > 0$, enters $\Gamma_{\delta}$ in finite time. In fact, if $(s(0), y(0), z(0))$ is outside $\Gamma_{\delta}$, it violates at least one of the inequalities defining $\Gamma_{\delta}$. If we first suppose that $s > \frac{\Lambda}{\mu} + \delta$, equation (2) yields:

$$ \dot{s} \le \Lambda - \mu s < \Lambda - \mu \left( \frac{\Lambda}{\mu} + \delta \right) = -\mu\delta < 0. \quad (A1) $$

$\dot{s} < -\mu\delta$ implies that:

$$ s(t) < s(0) - \mu\delta t. \quad (A2) $$

If we suppose that $s(t)$ never reaches $\Lambda/\mu + \delta$, we can write inequality (A2) at time $t_{\delta} = (s(0) - \Lambda/\mu - \delta)/(\mu\delta) \ge 0$. It then reads $s(t_{\delta}) < \Lambda/\mu + \delta$, which is in contradiction with the hypothesis that $s(t)$ never reaches $\Lambda/\mu + \delta$. Hence, $s(t)$ reaches $\Lambda/\mu + \delta$ in finite time and will stay below this value afterwards.

Similarly, if we now suppose that $\varepsilon s + z > \frac{\varepsilon\Lambda}{\zeta} + \delta$, equation (3) yields:

$$ \varepsilon\dot{s} + \dot{z} \le \varepsilon\Lambda - \zeta(\varepsilon s + z) < \varepsilon\Lambda - \zeta\left(\frac{\varepsilon\Lambda}{\zeta} + \delta\right) = -\zeta\delta $$

and $\varepsilon s + z$ will reach $\frac{\varepsilon\Lambda}{\zeta} + \delta$ in finite time and stay below this value afterwards.

Finally, supposing now that $y > \frac{\varepsilon\phi\Lambda}{\zeta\mu_y} + \frac{\phi\delta}{\mu_y} + \delta$ with $\varepsilon s + z \le \frac{\varepsilon\Lambda}{\zeta} + \delta$ (which may be valid at time 0 or after the previously described convergence in finite time), equation (4) yields:

$$ \dot{y} \le \phi z - \mu_y y < \phi \left( \frac{\varepsilon\Lambda}{\zeta} + \delta \right) - \mu_y \left( \frac{\varepsilon\phi\Lambda}{\zeta\mu_y} + \frac{\phi\delta}{\mu_y} + \delta \right) = -\delta\mu_y $$

and $y$ will reach $\frac{\varepsilon\phi\Lambda}{\zeta\mu_y} + \frac{\phi\delta}{\mu_y} + \delta$ in finite time and stay below this value afterwards. Hence, the compact set $\Gamma_{\delta}$ is reached in finite time and is positively invariant.

Local asymptotic stability of the pest-free steady state is readily obtained when $R_0 < 1$ as a consequence of the construction of $R_0$ through the next generation matrix approach of [28]. Also this equilibrium is unstable if $R_0 > 1$.

For global attractivity of the equilibrium, we will now distinguish two cases depending on the value taken by the threshold $R = \frac{\phi}{\mu_z}$.

If $R < 1$ (which implies $R_0 < 1$), let us consider the total population of female CBB $w = y + z$. Its dynamics is governed by:

$$ \dot{w} = -\mu_z(1 - R)z - \mu_y y \le -\kappa w, $$

where $\kappa = \min\{\mu_z(1 - R), \mu_y\}$. So using the comparison principle, we have $w(t) \le w(0)\exp(-\kappa t)$ which converges to 0, when $t \to +\infty$. As a consequence, the variable $s(t)$ will converge to $\frac{\Lambda}{\mu}$ when $t \to +\infty$, and the pest-free equilibrium is globally asymptotically stable.

If $R \ge 1$ and $R_0 \le 1$, let us consider the following Lyapunov candidate function $V(s, y, z) = \mu_y y + \phi z$, whose time derivative verifies:

$$ \begin{aligned} \dot{V}(s, y, z) &= -\mu_z \mu_y y + (-\mu_z + \phi) \frac{\varepsilon \beta s y}{y + \alpha s} \\ &= -\mu_z \mu_y y + (-\mu_z + \phi) \varepsilon \beta \frac{1}{\alpha} y - (-\mu_z + \phi) \varepsilon \beta \frac{1}{\alpha} y + (-\mu_z + \phi) \frac{\varepsilon \beta s y}{y + \alpha s} \\ &= (-\mu_z (\varepsilon \beta \frac{1}{\alpha} + \mu_y) + \varepsilon \beta \frac{1}{\alpha} \phi) y - (-\mu_z + \phi) \varepsilon \beta \frac{y}{\alpha} \left(1 - \frac{\alpha s}{y + \alpha s}\right) \\ &= -\left[\mu_z (\mu_y + \varepsilon \beta \frac{1}{\alpha})(1 - R_0) + (R - 1)\mu_z \varepsilon \beta \frac{y}{\alpha y + \alpha s}\right] y. \end{aligned} $$
As $\mathcal{R} \ge 1$ and $\mathcal{R}_0 \le 1$, $\dot{\mathcal{V}} \le 0$, so $\mathcal{V}$ is indeed a Lyapunov function for the pest-free steady state defined in $\Gamma_\delta$. This proves that $\mathcal{E}^0$ is stable. Its remains to prove that $\mathcal{E}^0$ is attractive. We shall use the LaSalle invariant principle. First of all notice that when $\mathcal{R}_0 = 1$, then $\mathcal{R} = 1 + \frac{\alpha\mu_y}{\varepsilon\beta} > 1$, such that it is impossible to have both $\mathcal{R}_0$ and $\mathcal{R}$ equal to 1 simultaneously. This implies that the expression inside the square brackets is always positive unless $y = 0$. Finally, the largest invariant set contained in $M = \{(s, y, z) \in \Gamma_\delta : \dot{\mathcal{V}} = 0\} = \{(s, 0, z) \in \Gamma_\delta\}$ is the pest-free steady state $\mathcal{E}^0$. Indeed, if $z \ne 0$, we have $\dot{y} > 0$, so that the corresponding solution leaves $M$. Also, if $s \ne \Lambda/\mu$, the $s$ dynamics in backward time ($y$ and $z$ remaining zero) bring the solution outside $\Gamma_\delta$. Thus, using the LaSalle invariance principle [44], we conclude that the pest-free steady state $\mathcal{E}^0$ is attractive inside $\Gamma_\delta$ which, combined with the global attractivity of $\Gamma_\delta$ and the local stability, ensures that $\mathcal{E}^0$ is globally asymptotically stable.

## A.2 Existence and stability of the coexistence steady state

Let $\bar{\mathcal{E}} = (\bar{s}, \bar{y}, \bar{z})$ be the positive coexistence steady state of system (1). Then, the components of $\bar{\mathcal{E}}$ solve for positive values the following system:

$$ 
\left\{ 
\begin{array}{l}
\Lambda - \beta \frac{\bar{s}\bar{y}}{\bar{y} + \alpha\bar{s}} - \mu\bar{s} = 0, \\
\phi\bar{z} - \varepsilon\beta \frac{\bar{s}\bar{y}}{\bar{y} + \alpha\bar{s}} - \mu_y\bar{y} = 0, \\
\varepsilon\beta \frac{\bar{s}\bar{y}}{\bar{y} + \alpha\bar{s}} - \mu_z\bar{z} = 0.
\end{array}
\right. 
$$

Multiplying the first equation by $\bar{y} + \alpha\bar{s}$ and eliminating $\bar{z}$ between the last two equations, we obtain, after some computations:

$$ 
\left\{ 
\begin{array}{l}
(\Lambda - \mu\bar{s})(\bar{y} + \alpha\bar{s}) - \beta\bar{s}\bar{y} = 0, \\
(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1)\bar{s} - \mu_y\bar{y} = 0.
\end{array}
\right. 
$$

Eliminating $\bar{y}$ between these two equations, we obtain:

$$ \bar{s} = \frac{\Lambda \left( \frac{1}{\mu_y} (\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \alpha \right)}{\frac{1}{\mu_y} (\mu + \beta)(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \mu\alpha}, $$

from which we build the other two components of $\bar{\mathcal{E}}$ as follows:

$$ \bar{y} = \frac{1}{\mu_y} (\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1)\bar{s}; \quad \bar{z} = \frac{\varepsilon\beta}{\mu_z} \frac{\bar{y}\bar{s}}{\bar{y} + \alpha\bar{s}}. $$

To investigate the stability of the coexistence steady state $\bar{\mathcal{E}}$, we compute the Jacobian matrix of system (1) evaluated at the steady state $\bar{\mathcal{E}}$:

$$ J = \begin{pmatrix} -\frac{\beta\bar{y}^2}{(\bar{y}+\alpha\bar{s})^2} - \mu & -\frac{\beta\alpha\bar{s}^2}{(\bar{y}+\alpha\bar{s})^2} & 0 \\ -\frac{\varepsilon\beta\bar{y}^2}{(\bar{y}+\alpha\bar{s})^2} & -\frac{\varepsilon\beta\alpha\bar{s}^2}{(\bar{y}+\alpha\bar{s})^2} - \mu_y & \phi \\ \frac{\varepsilon\beta\bar{y}^2}{(\bar{y}+\alpha\bar{s})^2} & \frac{\varepsilon\beta\alpha\bar{s}^2}{(\bar{y}+\alpha\bar{s})^2} & -\mu_z \end{pmatrix}. $$

Its characteristic polynomial is then given by:

$$ \lambda^3 + a_2 \lambda^2 + a_1 \lambda + a_0 = 0, $$
where:

$$a_2 = \mu_y + \mu + \mu_z + \frac{\varepsilon\beta\alpha + \beta\frac{1}{\mu_y^2}(\varepsilon\beta + \alpha\mu_y)^2(\mathcal{R}_0 - 1)^2}{\left[\frac{1}{\mu_y}(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \alpha\right]^2},$$

$$a_1 = \frac{\mu(\mu_z + \mu_z) \left( \frac{\varepsilon\beta(\phi - \mu_z)}{\mu_y \mu_z} \right)^2 + (\mu_y + \mu_z) \frac{\beta}{\mu_y^2} (\varepsilon\beta + \alpha\mu_y)^2 (\mathcal{R}_0 - 1)^2}{\left[ \frac{1}{\mu_y} (\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \alpha \right]^2}$$

$$+ \frac{\frac{\varepsilon\beta}{\mu_y}(\phi - \mu_z)(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \mu\varepsilon\beta\alpha}{\left[ \frac{1}{\mu_y}(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \alpha \right]^2},$$

$$a_0 = \frac{\mu\varepsilon\beta\frac{(\phi-\mu_z)}{\mu_y}(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \frac{\mu_z\beta}{\mu_y}(\varepsilon\beta + \alpha\mu_y)^2(\mathcal{R}_0 - 1)^2}{\left[ \frac{1}{\mu_y}(\varepsilon\beta + \alpha\mu_y)(\mathcal{R}_0 - 1) + \alpha \right]^2}.$$

The Routh-Hurwitz stability criterion imposes that $a_2 > 0$, $a_1 > 0$, $a_0 > 0$ and $a_2a_1 - a_0 > 0$. The first three are trivial when $\mathcal{N} > 1$, and the last one is verified as easy computations show term compensations between the negative terms of $-a_0$ and some terms within $a_1a_2$.

## B COMPUTATION OF THE TIME-DEPENDENT EFFECTIVE REPRODUCTION NUMBER

We compute the number of secondary colonising females produced by the colonising females present at time $t_r$. For that, we first follow the evolution of these colonising females that are present at time $t_r$; we define these by $\hat{y}$, whose dynamics are only decreasing, through mortality $\mu_y$ and transformation into infesting females (at rate $\varepsilon\beta \frac{s}{y+\alpha s}$):

$$\dot{\hat{y}} = - \left( \varepsilon\beta \frac{s(t)}{y(t) + \alpha s(t)} + \mu_y \right) \hat{y} \quad (B3)$$

and then, for $t \ge t_r$

$$\hat{y}(t) = y(t_r)e^{-\varepsilon\beta \int_{t_r}^{t} \frac{s(r)}{y(r)+\alpha s(r)} dr - \mu_y(t-t_r)}.$$

The dynamics of the infesting females produced by the latter colonising females $\hat{y}(t)$ are then:

$$\begin{aligned}
\dot{\hat{z}} &= \varepsilon\beta \frac{s}{y+\alpha s} \hat{y} - \mu_z \hat{z} \\
&= \varepsilon\beta \frac{s(t)}{y(t) + \alpha s(t)} e^{-\varepsilon\beta \int_{t_r}^{t} \frac{s(r)}{y(r)+\alpha s(r)} dr - \mu_y(t-t_r)} y(t_r) - \mu_z \hat{z} \quad (B4) \\
&= \eta(t)y(t_r) - \mu_z \hat{z}
\end{aligned}$$

whose solution is:

$$\hat{z}(t) = \hat{z}(t_r)e^{-\mu_z(t-t_r)} + \int_{t_r}^{t} \eta(\hat{t})y(t_r)e^{-\mu_z(t-\hat{t})}d\hat{t}.$$

Since we only consider infesting females that are initiated by the colonising females present at time $t_r$, we have $\hat{z}(t_r) = 0$, so that the secondary colonising females are given by:

$$\dot{\hat{y}} = \phi \hat{z} \quad (B5)$$
and

$$
\begin{align*}
\bar{y}(T) &= \int_{t_r}^{T} \phi \hat{z}(\tau) d\tau \\
&= \int_{t_r}^{T} \phi \int_{t_r}^{\tau} \eta(\hat{t}) e^{-\mu_z(\tau-\hat{t})} d\hat{t} d\tau y(t_r) \\
&= \int_{t_r}^{T} \phi \int_{t_r}^{\tau} \varepsilon \beta \frac{s(\hat{t})}{y(\hat{t}) + \alpha s(\hat{t})} e^{-\varepsilon \beta \int_{t_r}^{\hat{t}} \frac{s(r)}{y(r)+\alpha s(r)} dr - \mu_y(\hat{t}-t_r)} e^{-\mu_z(\tau-\hat{t})} d\hat{t} d\tau y(t_r).
\end{align*}
$$

If we pick $T = t_f$ in order to take into account the actual season length, $\bar{y}(T) \to 0$ as $t_r$ nears the end of the season, whatever the status of the infection. It hides a situation that might be troublesome for the onset of the next season. Hence, we opt to have $T \to +\infty$ so that all potential offspring are taken into account. We then obtain $\mathcal{R}_e(t)$ by renaming $t_r$ as $t$ and taking the limit for $T$ going to $+\infty$ of $\frac{\bar{y}(T)}{y(t)}$:

$$
\mathcal{R}_e(t) = \int_{t}^{+\infty} \int_{t}^{\tau} \phi \varepsilon \beta \frac{s(\hat{t})}{y(\hat{t}) + \alpha s(\hat{t})} e^{-\varepsilon \beta \int_{t}^{\hat{t}} \frac{s(r)}{y(r)+\alpha s(r)} dr - \mu_y(\hat{t}-t) - \mu_z(\tau-\hat{t})} d\hat{t} d\tau
$$

This closed-form formula has the advantage of being self-contained. In practice, in order to obtain $\mathcal{R}_e(t)$ for some $t$, we simulate system (1) until time $t$ and then simulate (1)-(B3)-(B4)-(B5) for a long duration from time $t$ (until some large time $T$). In the end, we obtain $\mathcal{R}_e(t)$ from the ratio $\frac{\bar{y}(T)}{y(t)}$.

C EXTRA SIMULATIONS

C.1 Simulation for a stable pest-free equilibrium

Coffee berry borers infest berries very efficiently, so the simulations we performed in Section 5.2, based on the best estimates
we could find for the model parameters, correspond to a situation in which the pest-free equilibrium is unstable and the basic
reproduction number $\mathcal{R}_0 \gg 1$.

To observe what happens with the optimal control when $\mathcal{R}_0 < 1$, that is when the pest-free equilibrium is stable, we drastically
reduced the infestation rate value to $\beta = 0.0005 \text{ day}^{-1}$ which corresponds to $\mathcal{R}_0 = 0.78 < 1$. We observe in Figure C1
that, with or without optimal control, the CBB population declines and the healthy coffee berry dynamics are very close to the
pest-free case. Moreover, the penalised profit $J = 472 \text{ US}$ is similar to the financial gain in the pest-free case. The optimal
control is singular and very little entomopathogenic fungus is applied on the coffee berries to accelerate the decline of the CBB
population. In this case, the optimal control is hence basically useless.

C.2 Simulations for intermediate initial conditions

Figure C2 illustrates the model dynamics for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^3$ colonising females,
Figure C3 for $y(0) = 10^5$ colonising females.

C.3 Simulation with discrete control

Figure C4 shows that with weekly fungus applications, the healthy berry production and the colonising female population at the end of the season are very close to those observed with the optimal continuous control. The penalised profit is hence only slightly lower in the discrete case ($J = 200 \text{ US}$) than in the optimal continuous case ($J = 207 \text{ US}$).
<figure><img src="image_9.png" /><figcaption>FIGURE C1 Simulation of system (1) without control (dashed blue curves) and system (9) with optimal control (16) (plain red curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^4$ colonising females, when the infestation rate $\beta = 0.0005 \text{ day}^{-1}$ which corresponds to $R_0 = 0.78 < 1$. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**FIGURE C2** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal control (16) (plain red curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^3$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted green curve). The coexistence steady state without control is indicated in (a–c) (dotted blue lines).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**FIGURE C3** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal control (16) (plain red curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^5$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted green curve). The coexistence steady state without control is indicated in (a–c) (dotted blue lines).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**FIGURE C4** Simulation of system (1) without control (dashed blue curves) and system (9) with optimal continuous control (16) (plain red curves) and weekly discrete control (plain magenta curves), for initial condition $s(0) = z(0) = v(0) = 0$ and $y(0) = 10^4$ colonising females. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve). The coexistence steady state without control is indicated in (a–b) (dotted blue lines).</figcaption></figure>
## References

[1] USDA. Coffee: World markets and trade, December 2019. https://www.fas.usda.gov/data/coffee-world-markets-and-trade (last accessed 2020-03-02).

[2] R. H. Le Pelley. *Pests of Coffee*. Tropical science series. Longmans, London, UK, 1968.

[3] A. Damon. A review of the biology and control of the coffee berry borer, Hypothenemus hampei (Coleoptera: Scolytidae). *Bull. Entomol. Res.*, 90(6):453-465, 2000. DOI: 10.1017/S0007485300000584.

[4] H. Ochoa Milián, O. Campos Almengor, B. Vidal S, and B. Decazy. Determinación de pérdidas en la cosecha por broca del fruto del cafeto Hypothenemus hampei Ferr., en función de diferentes porcentajes de infestación. *Revista Cafetalera – ANACAFE*, 303:23-27, 1989.

[5] Fernando E. Vega, Eric Rosenquist, and Wanda Collins. Global project needed to tackle coffee crisis. *Nature*, 425(6956):343, 2003. DOI: 10.1038/425343a.

[6] Elsie Burbano, Mark Wright, Donald E. Bright, and Fernando E. Vega. New record for the coffee berry borer, Hypothenemus hampei, in Hawaii. *J. Insect Sci.*, 11(1):117, 2011. DOI: 10.1673/031.011.11701.

[7] Fernando Vega, Francisco Infante, Alfredo Castillo, and Juliana Jaramillo. The coffee berry borer Hypothenemus hampei (Ferrari) (Coleoptera: Curculionidae): A short review, with recent findings and future research directions. *Terr. Arthropod Rev.*, 2:129-147, 2009. DOI: 10.1163/187498209X12525675906031.

[8] Luis F. Aristizábal, Alex E. Bustillo, and Steven P. Arthurs. Integrated pest management of coffee berry borer: Strategies from Latin America that could be useful for coffee farmers in Hawaii. *Insects*, 7(1):6, 2016. DOI: 10.3390/insects7010006.

[9] L. O. Brun, C. Marcillaud, V. Gaudichon, and D. M. Suckling. Cross resistance between insecticides in coffee berry borer, Hypothenemus hampei (Coleoptera: Scolytidae) from New Caledonia. *Bull. Entomol. Res.*, 84(2):175-178, 1994. DOI: 10.1017/S0007485300039651.

[10] L. O. Brun, C. Marcillaud, V. Gaudichon, and D. M. Suckling. Endosulfan resistance in Hypothenemus hampei (Coleoptera: Scolytidae) in New Caledonia. *J. Econ. Entomol.*, 82(5):1311-1316, 1989. DOI: 10.1093/jee/82.5.1311.

[11] B. P. Dufour and B. Frérot. Optimization of coffee berry borer, Hypothenemus hampei Ferrari (Col., Scolytidae), mass trapping with an attractant mixture. *J. Appl. Entomol.*, 132(7):591-600, 2008. DOI: 10.1111/j.1439-0418.2008.01291.x.

[12] A. Dufour, M. O González, J. J. Mauricio, B. A. Chavez, and R. Ramirez. Validation of coffee berry borer (CBB) trapping with the Brocap trap. In 20th International Conference on Coffee Science, pages 1243-1247, Bangalore, India, 1998.

[13] Russell Messing. The coffee berry borer (Hypothenemus hampei) invades Hawaii: Preliminary investigations on trap response and alternate hosts. *Insects*, 3(3):640-652, 2012. DOI: 10.3390/insects3030640.

[14] J. Jaramillo, C. Borgemeister, and P. Baker. Coffee berry borer Hypothenemus hampei (Coleoptera: Curculionidae): searching for sustainable control strategies. *Bull. Entomol. Res.*, 96(3):223-233, 2006. DOI: 10.1079/BER2006434.

[15] S. T. Murphy and D. Moore. Biological control of the coffee berry borer, *Hypothenemus hampei* (Ferrari) (Coleoptera, Scolytidae): previous programmes and possibilities for the future. *Biocontrol News and Information*, 11(2):107-117, 1990.

[16] Suzanne Lenhart and John T. Workman. *Optimal Control Applied to Biological Models*. Mathematical and Computational Biology Series. Chapman & Hall/CRC, 2007.

[17] Ramiro Ruiz-Cárdenas, Renato Martins Assunção, and Clarice Garcia Borges Demétrio. Spatio-temporal modelling of coffee berry borer infestation patterns accounting for inflation of zeroes and missing values. *Sci. Agric.*, 66(1):100-109, 2009. DOI: 10.1590/S0103-90162009000100014.

[18] Andrew P. Gutierrez, Amador Villacorta, Jose R. Cure, and C. Ken Ellis. Tritrophic analysis of the coffee (Coffea arabica) - coffee berry borer [Hypothenemus hampei (Ferrari)] - parasitoid system. *An. Soc. Entomol. Bras.*, 27(3):357-385, 1998. DOI: 10.1590/S0301-80591998000300005.
[19] Daniel Rodríguez, José Ricardo Cure, Andrew Paul Gutierrez, José Miguel Cotes, and Fernando Cantor. A coffee agroecosystem model: II Dynamics of coffee berry borer. Ecol. Model., 248:203–214, 2013. DOI: 10.1016/j.ecolmodel.2012.09.015.

[20] Yves Fotso, Frédéric Grognard, Berge Tsanou, and Suzanne Touzeau. Modelling and control of coffee berry borer infestation. In 14th African Conference on Research in Computer Science and Applied Mathematics, Stellenbosch, South Africa, October 2018.

[21] L. S. Pontryagin, V. G. Boltyanskii, E. F. Gamkrelize, and E. F. Mishchenko. The mathematical theory of optimal processes. Wiley, New York, 1962.

[22] J. F. Barrera. Dynamique des populations du scolyte des fruits du caféier, Hypothenemus hampei (Coleoptera: Scolytidae) et lutte biologique avec le parasitoïde Cephalonomia stephanoderis (Hymenoptera: Bethylidae), au Chiapas, Mexique. PhD thesis, Université Paul Sabatier, Toulouse, France, 1994.

[23] P. S. Baker, J. F. Barrera, and A. Rivas. life-history studies of the coffee berry borer (Hypothenemus hampei, Scolytidae) on coffee trees in Southern Mexico. J. Appl. Ecol., 29(3):656–662, 1992. DOI: 10.2307/2404473.

[24] Weliton D. Silva, Gabriel Moura Mascarin, Emiliana Manesco Romagnoli, and José Maurício S. Bento. Mating behavior of the coffee berry borer, Hypothenemus hampei (Ferrari) (Coleoptera: Curculionidae: Scolytinae). J. Insect. Behav., 25:408–417, 2012. DOI: 10.1007/s10905-011-9314-4.

[25] Roger Arditi and Lev R. Ginzburg. Coupling in predator-prey dynamics: ratio-dependence. J. Theor. Biol., 139(3):311–326, 1989. DOI: 10.1016/S0022-5193(89)80211-5.

[26] Roger Arditi and Lev R. Ginzburg. How species interact: altering the standard view on trophic ecology. Oxford University Press, 2012.

[27] Yang Kuang and Edoardo Beretta. Global qualitative analysis of a ratio-dependent predator-prey system. J. Math. Biol., 36(4):389–406, 1998. DOI: 10.1007/s002850050105.

[28] P. Van den Driessche and James Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci., 180(1–2):29–48, 2002. DOI: 10.1016/S0025-5564(02)00108-6.

[29] Katelyn M. Gostic, Lauren McGough, Edward B. Baskerville, Sam Abbott, Keya Joshi, Christine Tedijanto, Rebecca Kahn, Rene Niehus, James A. Hay, Pablo M. De Salazar, Joel Hellewell, Sophie Meakin, James D. Munday, Nikos I. Bosse, Katharine Sherrat, Robin N. Thompson, Laura F. White, Jana S. Huisman, Jérémie Scire, Sebastian Bonhoeffer, Tanja Stadler, Jacco Wallinga, Sebastian Funk, Marc Lipsitch, and Sarah Cobey. Practical considerations for measuring the effective reproductive number, $r_t$. PLoS Comput. Biol., 16(12):e1008409, 2020. DOI: 10.1371/journal.pcbi.1008409.

[30] J. Amang À. Mbang, P. Mounjouenpou, R. J. Mahob, M. Mbarga Amougo, J. Mouen Bedimo, S. Nyassé, L. Dibog, L. Bidzanga Nomo, I. R. Tchouamo, and R. Babin. Evaluation naturelle de l'impact de beauveria bassiana : champignon enthomopathogène dans la dynamique de population de *Hypothenemus hampei*, scolyte de baies des cerises de Coffea canefora. Afr. Crop Sci. J., 20(S2):443–451, 2012.

[31] Elsie Greco and Mark Wright. Efficacy of Beauveria bassiana applications on coffee berry borer across an elevation gradient in Hawaii. Biocontrol. Sci. Techn., 2018. DOI: 10.1080/09583157.2018.1493088.

[32] Richard Samuels, Rui Pereira, and Carlos Gava. Infection of the coffee berry borer Hypothenemus hampei (Coleoptera: Scolytidae) by Brazilian isolates of the entomopathogenic fungi Beauveria bassiana and Metarhizium anisopliae (Deuteromycotina: Hyphomycetes). Biocontrol. Sci. Techn., 12:631–635, 2002. DOI: 10.1080/0958315021000016298.

[33] N. Haraprasad, S. R. Niranjana, H. S. Prakash, H. S. Shetty, and Seema Wahab. Beauveria bassiana - A potential mycopathogenic agent for the efficient control of coffee berry borer, Hypothenemus hampei (Ferrari) in India. Biocontrol Sci. Techn., 11(2):251–260, 2001. DOI: 10.1080/09583150120035675.
[34] Pablo Benavides, Carmenza Góngora, and Alex Bustillo. IPM program to control coffee berry borer Hypothenemus hampei, with emphasis on highly pathogenic mixed strains of Beauveria bassiana, to overcome insecticide resistance in Colombia. In Farzana Khan Perveen, editor, *Insecticides: Advances in Integrated Pest Management*, chapter 23, pages 511–540. IntechOpen, Rijeka, Croatia, 2012. DOI: 10.5772/28740.

[35] W. De La Rosa, R. Alatorre, J. Trujillo, and J. F. Barrera. Virulence of Beauveria bassiana (Deuteromycetes) Strains Against the Coffee Berry Borer (Coleoptera: Scolytidae). *J. Econ. Entomol.*, 90(6):1534–1538, 1997. DOI: 10.1093/jee/90.6.1534.

[36] Alex E. Bustillo, Martha G. Bernal, Pablo Benavides, and Bernardo Chaves. Dynamics of Beauveria bassiana and Metarhizium anisopliae infecting Hypothenemus hampei (Coleoptera: Scolytidae) populations emerging from fallen coffee berries. *Fla. Entomol.*, 82(4):491–498, 1999.

[37] P. Waltman. *A second course in elementary differential equations*. Academic press, Orlando, 1968.

[38] W. H. Fleming and R. W. Rishel. *Deterministic and Stochastic Optimal Control*. Springer-Verlag, New York, 1975.

[39] L. G. De Pillis, W. Gu, K. R. Fister, K. Maples, A. Murugan, T. Neal, and K. Yoshida. Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls. *Math. Biosci.*, 209(1):292–315, 2007. DOI: 10.1016/j.mbs.2006.05.003.

[40] D. L. Luckes. *Differential Equations: Classical to Controlled*. Academic Press, New York, 1982.

[41] Arthur J. Krener. The high order maximal principle and its application to singular extremals. *SIAM J. Control Optim.*, 2(15):256–293, 1977. DOI: 10.1137/0315019.

[42] Team Commands, Inria Saclay. BOCOP: The optimal control solver, August 2019. http://bocop.org Version 2.2.1.

[43] Frédéric Bonnans, Pierre Martinon, D. Giorgi, V. Grélard, B. Heymann, J. Liu, S. Maindrault, and O. Tissot. BOCOP: A collection of examples, June 2019. http://bocop.org.

[44] J. P. LaSalle. Stability theory for ordinary differential equations. *J. Differ. Equ.*, 4(1):57–65, 1968. DOI: 10.1016/0022-0396(68)90048-X.