# Mathematical modelling and optimal control of production losses caused by Miridae

Myriam Djoukwe Tapi, André Nana Yakam, Roger Tagne Wafo, Samuel Bowong

► To cite this version:

Myriam Djoukwe Tapi, André Nana Yakam, Roger Tagne Wafo, Samuel Bowong. Mathematical modelling and optimal control of production losses caused by Miridae. Mathematical Modelling of Natural Phenomena, 2023, 18, pp.28. 10.1051/mmnp/2023030. hal-04323646

HAL Id: hal-04323646

https://hal.science/hal-04323646v1

Submitted on 5 Dec 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 CC BY 4.0 - Attribution - International License
# MATHEMATICAL MODELLING AND OPTIMAL CONTROL OF PRODUCTION LOSSES CAUSED BY MIRIDAE

MYRIAM DJOUKWE TAPI<fnref n="1" />,<fnref n="2" />✉, ANDRÉ NANA YAKAM<fnref n="3" />, ROGER TAGNE WAFO<fnref n="1" />✉ AND SAMUEL BOWONG<fnref n="1" />,<fnref n="2" />*✉

**Abstract.** Cocoa mirid, *Sahlbergella singularis*, is the major pest of cocoa (*Theobroma cacao*) responsible for several damage in plots in West Africa and particularly in Cameroon. Occasional damage accounts for 30–40% of pod losses. However, when miridae affect the foliage, gradual wilting occurs and eventually, tree death. A few studies have focused on describing the time evolution of Miridae in the plot in Cameroon, yet numerous questions remain. The aim of this paper is to estimate and control the losses of production caused by the bites of miridea. To do this, we will formulate and study a mathematical model for the dynamics of pods that takes into account the feeding and egg-laying of adults miridae on pods. We present the theoretical analysis of the model. More precisely, we compute equilibria and derive a threshold parameter that determines the presence or not of miridae in the plot. Throughout numerical simulations, we found that miridae can cause approximately 39.21% of production losses (which represents approximatively USD 1276.8 revenue losses) when initially, one has 1200 plants in the plot. After, we aim to increase cocoa production through optimal control. Optimal control consists in reducing density the number of nymphs and adults miridae in the plot. We studied the controlled model and we found that losses with control shrink to 20.58% which corresponds to USD 670.32 income revenue.

**Mathematics Subject Classification.** 35K57, 93D05, 65M06, 92D30.

Received January 11, 2023. Accepted August 19, 2023.

## 1. INTRODUCTION

Originally from Latin America, cocoa (*Theobroma cacao*) is one of the most important perennial crops worldwide [ICCO, 2016]. It is essential for the livelihood of millions of small producers in Africa especially in Cameroon and has a worldwide consumption. African countries supply over 75% of world cocoa [17]. Cocoa production was introduced in the coastal zones of Cameroon in the late 19th century from Latin America ([17], [22]). Cocoa plays an important role in Cameroon's economic development [3]. With its production, Cameroon is the world's

<footnote>Keywords and phrases: *Sahlbergella singularis*, Cocoa, mathematical models, stability, optimal control.</footnote>

<footnote n="1">¹ Department of Mathematics and Computer Science, Faculty of Science, The University of Douala, PO Box 24157 Douala-Cameroon.</footnote>

<footnote n="2">² IRD, Sorbonne University, UMMISCO, F-93143 Bondy, France.</footnote>

<footnote n="3">³ Faculty of Economics Science and Applied Management, The University of Douala, PO Box 17273, Douala, Cameroon.</footnote>

<footnote>* Corresponding author: sbowong@fs-univ-douala.cm/sbowong@gmail.cm</footnote>
fourth largest cocoa producer after Côte d'Ivoire, Ghana, and Indonesia [ICCO, 2014]. The cocoa sector is a source of employment for about four million individuals and it is the Cameroon's major agricultural export crop. The revenue generated from cocoa exports accounts for about 14% of non-oil exports in 2012, particularly to Europe [ICCO, 2014]. Compared to other agricultural activities, cocoa has been a leading sub-sector in Cameroon's economic growth and development and it remains the main cash crop to more than 75% of the population [28]. Cacao sector is very important for the economics of rural communities and consequently for that of their countries respectively.

Despite its real importance, cacao production faced many problems and damages. Wastes in cocoa production accounts for over 30% in annual losses and insect pests like capsids bugs (*Miridae* spp.) cause ripening of immature cocoa pods [26]. Other damages arise from parasites such as mistletoe (*Loranthus parasiticus*) which cause death of cocoa trees and tillage from cassava cultivators who cut the roots of cocoa trees. Damages are also known to result from overheating and stress when trees are exposed to direct sunlight. Cacao miridae feed on cacao by inserting their mouth parts into the plant tissues, injecting saliva into the lesion, and sucking the digestion products. This saliva has a marked histolytic effect probably due to the activity of the esterase [31]. Damages caused by bites and saliva action lead to the death of young shoots. In contrast, damage is usually less important on lignified tissues of hardened twigs and stems. Miridae and fungi attacks result in cankering or bark roughening, destruction of the flower cushions, severe dieback of twigs and branches, and sometimes partial degradation of cacao farms. High numbers of feeding punctures may cause a distortion of young pods during growth, or even death of severely damaged fruits [32]. Feeding lesions on shoots usually induce the drying up of buds and leaves of the terminal part of branches. Dry leaves stay on trees and they are easy to detect for several weeks. After a few months, feeding lesions on shoots develop into typical cankers, accumulating on branches and roughening the bark [15]. Yield losses attributed to miridae alone have not been precisely assessed in Cameroon. However, it is thought that cacao miridae may be responsible for yield losses that are higher than 25% in West Africa [2]. To better control miridae damage, development of pest management strategies is essential to prevent devastating impact on economy, food security, and biodiversity. Actually, in Cameroon, cacao miridae are generally controlled by chemical spraying. A wide range of spraying practices, sometimes very different from those recommended by agricultural extension has been reported ([4, 25]).

The interaction between miridae and cacao is essentially due to the action of miridae on cacao. Miridae use cacao for their feeding and their development. Miridae feed and lay on cacao: this action lead to damage observed on the tree: anecdotal damage on a pod and cumulated damage on the other parts of the tree. It is difficult to establish the link between Miridae and production losses because control is based on reducing miridae population. Moreover, the presence of miridae is often detected by observing damage in the plot. Bites on pods are check off on peduncle and according to the cortex size, damage on pods remain anecdotal. When the pod is already of considerable size, the harvest is not affected by miridae bites. The most harmful damages (damages which appear on the other parts of cacao) are cumulated over time and can lead to the destruction of the tree. Feeding and egg-laying of *S. singularis* lead to enormous damages which in turn lead to losses in production. The death on pods due to *S. singularis* occurs only in "cherelle" stage: this damage leads to the deformation of growing fruits or their drying in the case the peduncle is affected. The bites of *S. singularis* on mature pods don't affect harvesting but damages on the tree can lead to its destruction. In fact, the lesions on the tree due to miridae feeding are infected by a parasitical toadstool named *Calonectria rigidiuscula*. These lesions evolve into cankers which accumulate and weaken the branches and the trunk of the tree; this lead to the death of the entire branches, all the coronet and even the tree. This phenomenon is currently called "die-back". Miridae are very harmful pests for cocoa. Losses due to miridae are difficult to be estimated and can reach up to 30-40% of the potential production and their presence in the plot can lead to the destruction of cacao over time [1].

Mathematical approach appears as the most relevant tool to better estimate production losses due to miridae and also to find the way to increase production. Mathematical modelling has proven to be valuable in understanding the dynamics of many plant diseases and interaction between plants and bio-aggressors. In the special case of cacao, few models exist in the literature. In [12], the first authors and collaborators formulated several mathematical models to describe the time dynamics of the cocoa miridae. They also consider in a mathematical
model a biological control method based on mating-disrupting using artificial sex pheromone and trapping to limit the impact of miridae in the plot [13]. However, these two works do not take into account the interaction between the pods in several stages and miridae.

The aim of this paper is to estimate and control the production losses caused by miridae. We first formulate a coupled model between miridae and pods in order to evaluate the percentage of production losses. We consider several stages of development of the pods and we suppose that only adults miridae bit a pod. After the formulation of the model, we present the sensitivity analysis of the model in order to know the relative importance of some parameters that increase the production. We determine equilibria and their stability, the basic offspring number $R_0$ and its sensitivity analysis and estimate the production losses due to miridae. Next, optimal control is used to show the impact of miridae on pods production: optimal control consists in reducing the density number of nymphs and adults in the plot.

The paper is organized as follows. In Section 1, we formulate a mathematical model for the interaction between miridae population and pods development. Section 2 is devoted to the mathematical analysis of the mode. In the last section, we study the controlled model.

## 2. MATHEMATICAL MODELLING OF MIRIDAE AND PODS INTERACTION

In this section, we formulate a mathematical model for the interaction between miridae and pods population in the plot.

### 2.1. Biological backgrounds

Herein, we present biological backgrounds on the growth of pods and miridae which help us to formulate our mathematical model. The growth of pods and miridae need several stages. The life cycle of pods is made of three main stages: cherelle, young and mature pods and ripped pods. Cherelle is the first developmental stage of cocoa pods and pods less than 10 cm size enter this stage. After this stage, we have young and mature pods and the last stage is made of ripped pods which can be harvested and provide beans for commercialization. Pods are supported by a peduncle which proceeds from the development of thickness of the peduncle of flowers. Pods need from three to six months for their total development. This duration of time depends on the variety of cacao and it is necessary for the full development of fruits. After this, pod will ripe during one or two months, to undergo interior transformations and to change colours [18]. Sigmoid growth curve of fruits and their development stages is represented in Figure 1.

As cacao, the life cycle of *S. singularis* is dated of three stages: egg stage, nymph stage, and adult stage. These stage develop mainly either on pods either on shoots. The biological life cycle of *S. singularis* is given in Figure 2. Indeed, each individual needs to stay a certain time in each compartment to complete its stage development. The eggs are individually inserted into the host plant tissues [23] principally in the cortex of pods and sometimes under the bark of young shoots [14]. The incubation period of eggs is on average 15 days with a minimum of 9 days and a maximum of 21 days [19] before reaching nymph stage. Miridae *S. singularis* has a very long life cycle (eggs to adults). It is on average 40 days with a minimum of 36 days [7] and a maximum of 50 days [5].

### 2.2. Model formulation

Herein, we formulate a deterministic model for the interaction between the population dynamics of miridae and pods in the plot.

The pods population under consideration is grouped into different compartments according to their stage of development. We divide this population into two main classes or compartments representing their epidemiological status: Susceptible which are classified into three subclasses according to their development. At any time $t$, one has the following state:

(i) *Cherelle pods* $S_1$: Cherelles are the very first growth stage of pods. They are the pods resulting from flowering and having no more than 10 centimetres size.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Sigmoid growth curve of fruits and their development stages [29].</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Life cycle of S. singularis [12].</figcaption></figure>

(ii) *Young and mature pods S₂*: Young pods are the growing pods. After the cherelle stage, the pods will continue their development until the harvest.

(iii) *Ripped pods S₃*: This class contains pods that have completed their development, are ready to be harvested and whose seeds after fermentation and drying are commercialized.

(vi) *Infected pods I*: Infected pods are pods that have been impacted by miridae either through feeding or egg laying by females. We assume that an infected pod cannot be harvested. They are the pods for which the action of miridae is so strong that it prevents their development until harvest.

Thus, the total pod population at time t is

$$N(t) := S_1(t) + S_2(t) + S_3(t) + I(t). \qquad (2.1)$$

Based on its life cycle, miridae population is subdivided into three classes or compartments

(a) *Eggs E*: The eggs are individually inserted into the host plant tissues principally in the cortex of pods and sometimes under the bark of young shoots.
(b) Nymphs L: Nymphs come from the eggs. After on average 15 days, the eggs become nymphs.

(c) Adults A: Adults arise from nymphs that have completed their nymph development after on average of 25 days.

Cherelle pods are recruited at constant rate $\Lambda$, which corresponds to the regular production along the year. Denote by $\sigma_S$ and $\sigma_A$ respectively, the maximum number of miridae's bites a pod can have per unit of time and the maximum number of times one miridae want to bite pods per unit time. Pods (cherelle, young and mature pods) become infected due to the feeding and egg-laying. We assumed that a fraction $r$ of the population of adults are females. We also assumed that the number of males adults is enough and can ensure the fertilization of all available females. Following [11], we assume that the infection rate from miridae to cherelle and miridae to young and mature pods are respectively modelled as follows:

$$ \lambda_1(A, N) = f(A, N) \beta_1 \frac{S_1}{N} \quad \text{and} \quad \lambda_2(A, N) = f(A, N) \beta_2 \frac{S_2}{N}, \qquad (2.2) $$

where

$$ f(A, N) = \frac{\sigma_S \sigma_A A}{\sigma_S + \sigma_A (A/N)}, \qquad (2.3) $$

represents the total number of miridae's bites on pods. Indeed, the force of infection from miridae to cherelle is defined, as the product of the maximum number of miridae's bites one pods has per unit of time $f(A, N)$, the probability that a bite leads to infectious cherelle $\beta_1$ and the probability that the pod is a cherelle $S_1/N$. The force of infection from miridae to young and mature pods is defined, as the product of the maximum number of miridae's bites on pods per unit of time $f(A, N)$, the probability that a bite leads to infectious young and mature pods $\beta_2$ and the probability that the pod is a young or mature $S_2/N$.

Since pods become infected through contacts with miridae (adult), the appearing of new infected pods per unit time is given by

$$ \lambda(N, A) (\beta_1 S_1 + \beta_2 S_2) A = \frac{\beta_1 \sigma_S \sigma_A A S_1}{\sigma_S N + \sigma_A A} + \frac{\beta_2 \sigma_S \sigma_A A S_2}{\sigma_S N + \sigma_A A}, \qquad (2.4) $$

where

$$ \lambda(N, A) = \frac{\sigma_S \sigma_A}{\sigma_S N + \sigma_A A}. \qquad (2.5) $$

The biological parameters $\mu_1$, $\mu_2$, $\mu_3$ and $\mu_I$ represent respectively the natural daily death rate of cherelle, young and mature, ripped and infected pods. $\alpha S_1$ is the additional death of cherelle due to wilt. $\gamma_1$ and $\gamma_2$ represent respectively the maturation rate of cherelle and young or mature pods.

Following [12], we denote by $\mu_E$, $\mu_L$ and $\mu_A$ the daily death rate of eggs, nymphs and adults, respectively. $1/\nu_E$ and $1/\nu_L$ represent, respectively the needed time for an egg to become nymph and for a nymph to become an adult. We suppose that, in the presence of pods, adult mortality is regulated by the consumption so that, adult mortality is modeled by $\frac{\mu_A}{1 + \theta_1 S_1 + \theta_2 S_2}$ where $\mu_A$, $\theta_1$ and $\theta_2$ are respectively the natural mortality of adults, the regulation death rate of adults by consumption on cherelle and young pods.

We now formulate a model based on the interaction between pods and miridae that is summarized in Figure 3. From the flowchart diagram in Figure 3, the dynamics of the interaction between miridae and pods is given by
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Structure of the model.</figcaption></figure>

the following system of ordinary differential equations:

$$
\left\{
\begin{aligned}
\dot{S}_1 &= \Lambda - \frac{\beta_1 \sigma_S \sigma_A A S_1}{\sigma_S N + \sigma_A A} - (\gamma_1 + \mu_1 + \alpha S_1)S_1, \\
\dot{S}_2 &= \gamma_1 S_1 - \frac{\beta_2 \sigma_S \sigma_A A S_2}{\sigma_S N + \sigma_A A} - (\gamma_2 + \mu_2) S_2, \\
\dot{S}_3 &= \gamma_2 S_2 - \mu_3 S_3, \\
\dot{I} &= \frac{\sigma_S \sigma_A A}{\sigma_S N + \sigma_A A} (\beta_1 S_1 + \beta_2 S_2) - \mu_I I, \\
\dot{E} &= r b A \left(1 - \frac{E}{K}\right) - (\nu_E + \mu_E) E, \\
\dot{L} &= \nu_E E - (\nu_L + \mu_L) L, \\
\dot{A} &= \nu_L L - \frac{\mu_A}{1 + \theta_1 S_1 + \theta_2 S_2} A,
\end{aligned}
\right.
\qquad (2.6)
$$
<table><caption>TABLE 1. Variables of system (2.6).</caption><thead><tr><th>Symbols</th><th>Biological meaning</th><th>Units</th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>Cherelle pods</td><td>Number</td></tr><tr><td>S<sub>2</sub></td><td>Young pods</td><td>Number</td></tr><tr><td>S<sub>3</sub></td><td>Ripped pods</td><td>Number</td></tr><tr><td>I</td><td>Infected pods</td><td>Number</td></tr><tr><td>E</td><td>Eggs</td><td>Number</td></tr><tr><td>L</td><td>Nymphs</td><td>Number</td></tr><tr><td>A</td><td>Adults miridae</td><td>Number</td></tr><tr><td>N</td><td>Total number of pods population</td><td>Number</td></tr></tbody></table>

<table><caption>TABLE 2. Parameter description of system (2.6).</caption><thead><tr><th>Parameters</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>Λ</td><td>Cherelle “recruitment” rate</td><td>Number × Days<sup>-1</sup></td></tr><tr><td>γ<sub>1</sub></td><td>Progression rate from cherelle to young pod</td><td>Days<sup>-1</sup></td></tr><tr><td>γ<sub>2</sub></td><td>Progression rate from young pod to ripe pod</td><td>Days<sup>-1</sup></td></tr><tr><td>α</td><td>Additional death rate due to wilt</td><td>Days<sup>-1</sup></td></tr><tr><td>σ<sub>S</sub></td><td>Number of bites supported by a pod per unit of time</td><td></td></tr><tr><td>σ<sub>A</sub></td><td>Number of bites done by Miridae per unit of time</td><td></td></tr><tr><td>β<sub>1</sub></td><td>Infection rate of cherelles stage</td><td>Days<sup>-1</sup></td></tr><tr><td>β<sub>2</sub></td><td>Infection rate of young pods stage</td><td>Days<sup>-1</sup></td></tr><tr><td>r</td><td>Proportion of female adults</td><td></td></tr><tr><td>b</td><td>Mean number of eggs laid by a mature female</td><td>Days<sup>-1</sup></td></tr><tr><td>1/ν<sub>E</sub></td><td>Time necessary for which an egg to become a nymph</td><td>Days</td></tr><tr><td>1/ν<sub>L</sub></td><td>Duration of the development of nymphs</td><td>Days</td></tr><tr><td>μ<sub>1</sub></td><td>Mortality rate of cherelles</td><td>Days<sup>-1</sup></td></tr><tr><td>μ<sub>2</sub></td><td>Mortality rate of young pods</td><td>Days<sup>-1</sup></td></tr><tr><td>μ<sub>3</sub></td><td>Mortality rate of ripped pods</td><td>Days<sup>-1</sup></td></tr><tr><td>μ<sub>I</sub></td><td>Mortality rate of infected pods</td><td>Days<sup>-1</sup></td></tr><tr><td>μ<sub>E</sub></td><td>Mortality of eggs</td><td>Days<sup>-1</sup></td></tr><tr><td>μ<sub>L</sub></td><td>Mortality of nymphs</td><td>Days<sup>-1</sup></td></tr><tr><td>μ<sub>A</sub></td><td>Mortality of adults</td><td>Days<sup>-1</sup></td></tr><tr><td>θ<sub>1</sub></td><td>Regulation death rate of adults by consumption cherelles</td><td>Days<sup>-1</sup></td></tr><tr><td>θ<sub>2</sub></td><td>Regulation death rate of adults by consumption young pods</td><td>Days<sup>-1</sup></td></tr></tbody></table>

with the following nonnegative initial conditions:

$$ \begin{cases} S_1(0) = S_1^0, \quad S_2(0) = S_2^0, \quad S_3(0) = S_3^0, \quad I(0) = I^0, \\ E(0) = E^0, \quad L(0) = L^0 \quad \text{and} \quad A(0) = A^0. \end{cases} \tag{2.7} $$

A complete list and description of all variables and parameters of system (2.6) is summarized in Tables 1 and 2, respectively. For biological reasons, all the parameters are nonnegative.

## 2.3. Sensitivity analysis

Sensitivity analysis aims to determine the optimal values of the model parameters that can increase the cocoa production. We consider the range of parameters values given in Table 3. For some parameters, we have a relatively good idea about the range of values, for others, this is more vague, but we choose values that, according to field experts, are the most reasonable. In any case, and whatever the methods used to make the
<table><caption>TABLE 3. Parameters values used for the sensitivity analysis of system (2.6).</caption><thead><tr><th>Parameters</th><th>Values for sensitivity analysis</th><th>Source</th></tr></thead><tbody><tr><td>Λ</td><td>[1000, 3000]</td><td>Estimated</td></tr><tr><td>γ<sub>1</sub></td><td>[0.0001, 0.9]</td><td>[29] and [27]</td></tr><tr><td>γ<sub>2</sub></td><td>[0.0001, 0.9]</td><td>Estimated</td></tr><tr><td>α</td><td>[0.0005, 0.5]</td><td>Estimated</td></tr><tr><td>σ<sub>S</sub></td><td>[1, 30]</td><td>Estimated</td></tr><tr><td>σ<sub>A</sub></td><td>[1, 30]</td><td>Estimated</td></tr><tr><td>β<sub>1</sub></td><td>[0.0001, 0.9]</td><td>Estimated</td></tr><tr><td>β<sub>2</sub></td><td>[0.0001, 0.9]</td><td>Estimated</td></tr><tr><td>r</td><td>[0.1, 0.8]</td><td>[5]</td></tr><tr><td>b</td><td>[1, 4]</td><td>Estimated</td></tr><tr><td>1/ν<sub>E</sub></td><td>[1, 100]</td><td>Estimated</td></tr><tr><td>1/ν<sub>L</sub></td><td>[1, 100]</td><td>Estimated</td></tr><tr><td>μ<sub>1</sub></td><td>[0.001, 0.8]</td><td>[27]</td></tr><tr><td>μ<sub>2</sub></td><td>[0.0001, 0.8]</td><td>[27]</td></tr><tr><td>μ<sub>3</sub></td><td>[0.0001, 0.8]</td><td>Estimated</td></tr><tr><td>μ<sub>I</sub></td><td>[0.0001, 0.8]</td><td>Estimated</td></tr><tr><td>μ<sub>E</sub></td><td>[0.005, 0.8]</td><td>Estimated</td></tr><tr><td>μ<sub>L</sub></td><td>[0.005, 0.8]</td><td>Estimated</td></tr><tr><td>μ<sub>A</sub></td><td>[0.005, 0.9]</td><td>Estimated</td></tr><tr><td>θ<sub>1</sub></td><td>[0.00001, 0.5]</td><td>Estimated</td></tr><tr><td>θ<sub>2</sub></td><td>[0.00001, 0.5]</td><td>Estimated</td></tr></tbody></table>

full sensitivity analysis, the results are discussed and valid only under the chosen intervals. We will use e-fast method which indicates which parameter uncertainty has the greatest impact on the output variability (see, for instance, Marino et al. [21] for further explanations). E-fast sensitivity method highlights the effects of the first order called main effects and the total effects that combine the main effects and all the interaction effects of the parameters on the outputs of system (2.6). It is a global sensitivity technique based on the decomposition or partitioning of the variance. The variance of the model output is decomposed into components resulting from the individual effects of parameters as well as their interactions.

Figure 4 presents the sensitivity analysis (e-fast method) of all variables of system (2.6). The white part in Figure 4 shows the first-order sensitivity index (main index), while the sum of two parts (white and gray) shows the total sensitivity index. It illustrates the sensitivity of twenty (20) parameters on the output variables $S_1$, $S_2$, $S_3$, $I$, $E$, $L$ and $A$. This method provides that parameters $\mu_3$ and $\gamma_1$ are the most sensitive parameters for ripped pods ($S_3$) and $\mu_L$ are the most sensitive parameters for adults population ($A$). For infected pods, all parameter have approximately the same sensitivity but $\mu_I$ is the most sensitive parameter. For miridae population, carrying capacity $K$ has a strong impact on all variables $E$, $L$ and $A$. $\nu_L$ and $\mu_L$ are the most sensitive parameters for nymph's population and $\nu_L$ is the most sensitive parameter for adults.

## 3. MODEL ANALYSIS

In this section, we study the long term dynamics of system (2.6).

### 3.1. Basic properties

It is important to show that our problem is well-posed so that solutions are always positive and bounded. Herein, we study the basics properties of system (2.6) which are essential in the proof of the stability results. We have the following result about the existence and positivity of solutions of system (2.6).

**Theorem 3.1.** For any nonnegative initial conditions, there exists a unique solutions of system (2.6) which is non-negative.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Efast sensitivity analysis of system (2.6).</figcaption></figure>
*Proof.* Since the right-hand side of system (2.6) is Lipschitz continuous, for any initial condition, there exists a unique maximal solution.

System (2.6) can be rewritten in the following compact form:

$$ egin{cases} \frac{dX}{dt} = A(X)X + F, \\ X(0) = X_0 \ge 0, \end{cases} \qquad (3.1) $$

where $X(t) = (S_1, S_2, S_3, I, E, L, A)^T$,

$$ A(X) = \begin{pmatrix} -B_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \gamma_1 & -B_2 & 0 & 0 & 0 & 0 & 0 \\ 0 & \gamma_2 & -\mu_3 & 0 & 0 & 0 & 0 \\ C_1 & C_2 & 0 & -\mu_I & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -B_3 & 0 & rb \\ 0 & 0 & 0 & 0 & \nu_E & -(\nu_L + \mu_L) & 0 \\ 0 & 0 & 0 & 0 & 0 & \nu_L & -C_3 \end{pmatrix} \quad \text{and} \quad F = \begin{pmatrix} \Lambda \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, $$

with

$$ B_1 = \gamma_1 + \mu_1 + \alpha S_1 + \frac{\sigma_S \sigma_A \beta_1 A}{\sigma_S N + \sigma_A A}, \quad B_2 = \gamma_2 + \mu_2 + \frac{\sigma_S \sigma_A \beta_2 A}{\sigma_S N + \sigma_A A}, \quad B_3 = (\nu_E + \mu_E) + \frac{rbA}{K}, $$

$$ C_1 = \frac{\sigma_S \sigma_A \beta_1 A}{\sigma_S N + \sigma_A A}, \quad C_2 = \frac{\sigma_S \sigma_A \beta_2 A}{\sigma_S N + \sigma_A A} \quad \text{and} \quad C_3 = \frac{\mu_A}{1 + \theta_1 S_1 + \theta_2 S_2}. $$

Since $A(X)$ is a Metzler matrix and $F \ge 0$, system (2.6) is positively invariant in $\mathbb{R}_+^7$ [8]. This means that any trajectory of the system starting from an initial state in the positive orthant $\mathbb{R}_+^7$ remains forever in $\mathbb{R}_+^7$. This achieves the proof. □

Now, let us prove that the state variables of system (2.6) are bounded. We have the following result.

**Theorem 3.2.** *The set*

$$ \Omega = \{(S_1, S_2, S_3, I, E, L, A) \in \mathbb{R}_+^7; S_1 \le S_1^0, S_2 \le S_2^0, S_3 \le S_3^0, I \le I^0, E \le K, L \le \frac{\nu_E K}{\mu_L + \nu_L} \text{ and } A \le \frac{\nu_E \nu_L K(1 + \theta_1 S_1^0 + \theta_2 S_2^0)}{\mu_A(\nu_L + \mu_L)}\}, \qquad (3.2) $$

is positively invariant by system (2.6) where

$$ S_1^0 = \frac{\sqrt{(\gamma_1 + \mu_1)^2 + 4\alpha\Lambda} - (\gamma_1 + \mu_1)}{2\alpha}, \quad S_2^0 = \frac{\gamma_1 S_1^0}{\gamma_2 + \mu_2}, \quad S_3^0 = \frac{\gamma_2 S_2^0}{\mu_3} \quad \text{and} \quad I^0 = \frac{\beta_1 A^0 S_1^0 + \beta_2 A^0 S_2^0}{\mu_I}. $$
Proof. From the first equation of system (2.6), one has that

$$ \dot{S}_1(t) \leq \Lambda - (\gamma_1 + \mu_1 + \alpha S_1(t)) S_1(t). $$

Using Gronwall's lemma yields

$$ S_1(t) \leq S_1^0 + \frac{1}{\exp((\gamma_1 + \mu_1 + 2\alpha S_1^0)t) \left(\frac{\gamma_1 + \mu_1 + 2\alpha S_1^0}{S_1(0) - S_1^0} + \alpha\right) - \frac{\alpha}{\gamma_1 + \mu_1 + 2\alpha S_1^0}}. $$

Thus, if $S_1(0) \leq S_1^0$, one has

$$ 0 \leq S_1(t) \leq S_1^0 \quad \text{for all } t \geq 0. $$

From the second equation of system (2.6), one has

$$ \dot{S}_2(t) \leq \gamma_1 S_1(t) - (\gamma_2 + \mu_2) S_2(t) $$

which leads to

$$ S_2(t) \leq S_2^0 + (S_2(0) - S_2^0)e^{-(\gamma_2+\mu_2)t}. $$

Then, if $S_2(0) \leq S_2^0$, one can deduce that

$$ 0 \leq S_2(t) \leq S_2^0 \quad \text{for all } t \geq 0. $$

Similarly, one proves that $S_3(t) \leq S_3^0$ and $I(t) \leq I^0$ for all $t \geq 0$.

From the fifth equation of system (2.6), one has

$$ \dot{E}(t) = r b A(t) \left(1 - \frac{E(t)}{K}\right) - (\nu_E + \mu_E) E(t). $$

Assume that there exists $\varepsilon > 0$ such that

$$ t_1 \leq t_1 + \varepsilon < T \quad \text{and} \quad E(t_1 + \varepsilon) > K. $$

Now define

$$ t_1^* = \inf\{t \geq t_1, L(t) \geq K\} $$

It then comes that

$$ E(t_1^*) = K. $$

Thus,

$$ E(t) = E(t_1^*) + \dot{E}(t_1^*)(t - t_1^*) + o(t - t_1^*) \quad \text{and} \quad \dot{E}(t_1^*) = -(\nu_E + \mu_E)E(t_1^*) < 0. $$
So, there exists $\varepsilon_1 > 0$ such that

$$t_1^* \le t < t_1^* + \varepsilon_1 \quad \text{and} \quad E(t_1^*) < K.$$

This is absurd because $t_1^* = \inf\{t \ge t_1, E(t) \ge K\}$. Consequently, $E(t) \le K_C$ for all $t \ge t_0$. By integrating equations in $L$ and $A$ separately in system (2.6) and according to Gronwall's lemma, one proves that

$$L(t) \le \frac{\nu_E K}{\nu_L + \mu_L} \quad \text{and} \quad A(t) \le \frac{\nu_E \nu_L K(1 + \theta_1 \tilde{S}_1 + \theta_2 \tilde{S}_2)}{\mu_A} \quad \text{for all } t \ge 0.$$

This concludes the proof.

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

## 3.2. Asymptotic behavior

Herein, we determine equilibria of system (2.6) and study their stabilities.

### 3.2.1. The pest free equilibrium (PFE) and its stability

The pest-free equilibrium point is obtained by solving all equations of system (2.6) equal to zero with $A = 0$. Note that $A = 0$ implies $E = L = I = 0$. With this in mind, we obtain the following system of equations:

$$\left\{
\begin{array}{l}
(\gamma_1 + \mu_1 + \alpha S_1) S_1 = \Lambda, \\
(\gamma_2 + \mu_2) S_2 = \gamma_1 S_1, \\
\gamma_2 S_2 = \mu_3 S_3.
\end{array}
\right.
\qquad (3.3)$$

Solving Eq. (3.3) gives

$$S_1^0 = \frac{(\gamma_1 + \mu_1) + \sqrt{(\gamma_1 + \mu_1)^2 + 4\alpha\Lambda}}{2\alpha}, \quad S_2^0 = \frac{\gamma_1 S_1^0}{(\gamma_2 + \mu_2)} \quad \text{and} \quad S_3^0 = \frac{\gamma_2 S_2^0}{\mu_3}. \qquad (3.4)$$

Then, the PFE of system (2.6) is

$$Q^0 = (S_1^0, S_2^0, S_3^0, 0, 0, 0, 0)^T. \qquad (3.5)$$

where $S_1^0, S_2^0$ and $S_3^0$ are defined as in equation (3.4).

Now, let us study the stability of the PFE. To do this, we need to compute the basic offspring number.

The basic offspring number $\mathcal{R}_0$ is one of the most important threshold for disease control. In this special case, it is defined as mean number of adults female produced by one adult female over its life span. We use the method of Van den Driessche and Watmough [30] to compute the basic offspring number of system (2.6).

Infected pods and the miridae population of system (2.6) can be rewritten in the following compact form:

$$\frac{dX}{dt} = \mathcal{F}(X) - \mathcal{V}(X),$$
where $X = (I, E, L, A)^T$, $\mathcal{F}(X)$ is the incidence rate of new infections, and $\mathcal{V}(X)$ is the transfer rate of individuals into, and out of, each sub-population defined as follows:

$$ \mathcal{F} = \begin{pmatrix} \lambda(N, A) A(\beta_1 S_1 + \beta_2 S_2) \\ r b A \left(1 - \frac{E}{K}\right) \\ 0 \\ 0 \end{pmatrix} \quad \text{and} \quad \mathcal{V} = \begin{pmatrix} \mu_I I \\ (\nu_E + \mu_E) E \\ -\nu_E E + (\nu_L + \mu_L) L \\ -\nu_L L + \frac{\mu_A A}{1 + \theta_1 S_1 + \theta_2 S_2} \end{pmatrix}. $$

The Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ at the pest-free equilibrium $Q^0$ are, respectively,

$$ F = \begin{pmatrix} 0 & 0 & 0 & \frac{\sigma_A \sigma_S}{\sigma_S (S_1^0 + S_2^0 + S_3^0)} (\beta_1 S_1^0 + \beta_2 S_2^0) \\ 0 & 0 & 0 & r b \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, $$

and

$$ V = \begin{pmatrix} \mu_I & 0 & 0 & 0 \\ 0 & \nu_E + \mu_E & 0 & 0 \\ 0 & -\nu_E & \nu_L + \mu_L & 0 \\ 0 & 0 & -\nu_L & \frac{\mu_A}{1 + \theta_1 S_1^0 + \theta_2 S_2^0} \end{pmatrix} $$

Then, the basic offspring number is the spectral radius of the next generation matrix, $FV^{-1}$ given by

$$ \mathcal{R}_0 = \rho(FV^{-1}) = \frac{rb\nu_L\nu_E(1 + \theta_1 S_1^0 + \theta_2 S_2^0)}{\mu_A(\nu_E + \mu_E)(\nu_L + \mu_L)}. \quad (3.6) $$

Following Van den Driessche and Watmough [30], we have the following result about the local stability of the pest free equilibrium $Q^0$.

**Proposition 3.3.** The pest free equilibrium $Q^0$ is locally asymptotically stable when $\mathcal{R}_0 \le 1$ and unstable when $\mathcal{R}_0 > 1$.

We point out that, it is important to always have $\mathcal{R}_0 \le 1$ so that the production of cocoa in the plot is not affected by the pest. To do so, it is important to determine the parameters that impact negatively and positively the basic offspring number $\mathcal{R}_0$ by performing the sensitivity analysis of $\mathcal{R}_0$. We stress that, the sensitivity analysis performed here is not to be confused with the sensitivity analysis presented in Section 2.3 where the analysis is performed to determine which parameters can be used to decrease the values of the state variables. Indeed, the sensitivity Analysis presented in Section 2.3 used the eFast method in order to detect the
<figure><img src="image_7.png" /><figcaption>FIGURE 5. PRCC sensitivity analysis of the basic offspring number $\mathcal{R}_0$.</figcaption></figure>

most sensitive parameters, that is the parameters that most influence the output variable of the model. This can help to predict the effect of each parameter on the model results and classify them according to their degree of sensitivity. Here, we are interested to determine the model parameters that significantly affect the basic reproduction number $\mathcal{R}_0$ since they are parameters that should be taken into consideration when considering an intervention strategy. To do so, we use the LHS-PRCC analysis. These two sensitivity analysis can help to identify the model parameters that are most influential in determining disease dynamics and have the advantage that the entire parameter is explored.

Figure 5 shows the LHS-PRCC analysis of the basic offspring number $\mathcal{R}_0$. One can observe that $\mathcal{R}_0$ is more impacted specifically by parameters $r$, $b$, $\Lambda$, $\theta_1$, $\nu_E$, $\nu_L$, $\alpha$ which have a positive effect on $\mathcal{R}_0$ and $\gamma_2$ and $\mu_A$ which have a negative effect.

The local asymptotic stability of the pest free equilibrium $Q^0$ does not guarantee the complete elimination of the pest in the plot. Only the global asymptotic stability of the pest-free equilibrium $Q^0$ ensures that the pest either dies out or persists within a cocoa plot. We are now going to study the global asymptotic stability of the pest-free equilibrium $Q^0$.

System (2.6) can be written in the following compact form:

$$ \begin{cases} \dot{x} = h(x, y), \\ \dot{y} = (\tilde{F} - \tilde{V})y, \end{cases} \qquad (3.7) $$

where $x = (S_1, S_2, S_3)^T$ represents the class of the non-infected pods, $y = (I, E, L, A)^T$ represents the class of infected pods and Miridae population,

$$ \tilde{F} = \begin{pmatrix} 0 & 0 & 0 & \lambda(N, A)(\beta_1 S_1 + \beta_2 S_2) \\ 0 & 0 & 0 & rb \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad \tilde{V} = \begin{pmatrix} \mu_I & 0 & 0 & 0 \\ 0 & \nu_E + \mu_E & 0 & 0 \\ 0 & -\nu_E & \nu_L + \mu_L & 0 \\ 0 & 0 & -\nu_L & \frac{\mu_A}{1 + \theta_1 S_1 + \theta_2 S_2} \end{pmatrix}, $$
and

$$h(x, y) = \begin{pmatrix} \Lambda - \lambda(N, A) \beta_1 A S_1 - (\gamma_1 + \mu_1 + \alpha S_1) S_1 \\ \gamma_1 S_1 - \lambda(N, A) \beta_2 A S_2 - (\gamma_2 + \mu_2) S_2 \\ \gamma_2 S_2 - \mu_3 S_3 \end{pmatrix} .$$

Since

$$\lambda(N, A) (\beta_1 S_1 + \beta_2 S_2) = \frac{\beta_1 \sigma_S \sigma_A S_1}{\sigma_S N + \sigma_A A} + \frac{\beta_2 \sigma_S \sigma_A S_2}{\sigma_S N + \sigma_A A} \leq \sigma_A \left( \frac{\beta_1 S_1}{N} + \frac{\beta_2 S_2}{N} \right) \leq \sigma_A (\beta_1 + \beta_2),$$

system (3.7) can be bounded by the following system of equations:

$$\begin{cases} \dot{\bar{x}} = h(\bar{x}, \bar{y}), \\ \dot{\bar{y}} = (\bar{F} - V) \bar{x}. \end{cases} \qquad (3.8)$$

where

$$\bar{F} = \begin{pmatrix} 0 & 0 & 0 & \sigma_A (\beta_1 + \beta_2) \\ 0 & 0 & 0 & r b \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} .$$

Then, the dynamics of system (3.7) is deduced by the dynamics of system (3.8). Moreover, system (3.8) is globally asymptotically stable if the spectral radius of matrix $\bar{F}V^{-1}$ is less than one, i.e $\rho(\bar{F}V^{-1}) \leq 1$. A single computation shows that $\rho(\bar{F}V^{-1}) \leq 1$ implies $\rho(FV^{-1}) \leq 1$ that is $R_0 \leq 1$. This means that the pest free equilibrium $Q^0$ is globally asymptotically stable when $R_0 \leq 1$. We have proved the following result.

**Theorem 3.4.** *The pest free equilibrium $Q^0$ is globally asymptotically stable in $\Omega$ when $R_0 \leq 1$.*

For numerical simulations, we use parameters values consigned in Table 3. Initial conditions are $S_1(0) = 1200$, $S_2(0) = S_3(0) = I(0) = E(0) = 400$, $L(0) = 360$ and $A(0) = 600$. Indeed, according to [6], the mean number of miridae in the plot is evaluated as 960. We consider for our numerical simulation 360 nymphs and 600 adults for initial conditions. The initial number of eggs and cherelle is chosen arbitrarily. Figure 6, next page presents the trajectories of system (2.6) when $\Lambda = 14, 400$; $\gamma_1 = 0.05$; $\gamma_2 = 0.027$; $\alpha = 0.001$; $\sigma_S = 20$; $\sigma_A = 26$; $\beta_1 = 0.09$; $\beta_2 = 0.02$; $r = 0.58$; $b = 3.28$; $\nu_E = 1/15$; $\nu_L = 1/25$; $\mu_1 = \mu_I = 0.05$; $\mu_2 = \mu_3 = 0.00469$; $\mu_E = 0.09$; $\mu_L = 0.19$; $\mu_A = 0.27$; $\theta_1 = \theta_2 = 0.00001$ (so that $R_0 = 0.6896 \leq 1$). From this figure, it clearly appears that the trajectories of system (2.6) converge to the pest free equilibrium $Q^0$ as shown in Theorem 3.4. This means that the pest disappears within a cocoa plot and we have an optimal production in the plot.

### 3.2.2. The pest equilibrium point (PEP) and its stability

Herein, we prove numerically the existence of a pest equilibrium point of system (2.6) and we study its stability.
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Numerical simulations for system (2.6) when $\Lambda = 14,400$; $\gamma_1 = 0.05$; $\gamma_2 = 0.027$; $\alpha = 0.001$; $\sigma_S = 20$; $\sigma_A = 26$; $\beta_1 = 0.09$; $\beta_2 = 0.02$; $r = 0.58$; $b = 3.28$; $\nu_E = 1/15$; $\nu_L = 1/25$; $\mu_1 = \mu_I = 0.05$; $\mu_2 = \mu_3 = 0.00469$; $\mu_E = 0.09$; $\mu_L = 0.19$; $\mu_A = 0.27$; $\theta_1 = \theta_2 = 0.00001$ (so that $R_0 = 0.6896 \le 1$). (a) Cherelle $S_1$; (b) Young pods $S_2$; (c) Ripped pods $S_3$; (d) Inefected pods $I$; (e) Eggs Miridae $E$; (f) Nymphs Miridae $L$ and (g) Adults Miridae $A$.</figcaption></figure>
Let $\overline{Q} = (\bar{S}_1, \bar{S}_2, \bar{S}_3, \bar{I}, \bar{E}, \bar{L}, \bar{A})$ be a pest equilibrium point (equilibrium with infested pods), with $\bar{S}_1, \bar{S}_2, \bar{S}_3, \bar{I}, \bar{E}, \bar{L}$, and $\bar{A}$ satisfy the following system of equations:

$$ 
\left\{
\begin{aligned}
& \bar{\lambda}(\bar{A}, \bar{N}) \beta_1 \bar{A} \bar{S}_1 + (\mu_1 + \gamma_1 + \alpha \bar{S}_1) \bar{S}_1 = \Lambda, \\
& \bar{\lambda}(\bar{A}, \bar{N}) \beta_2 \bar{A} \bar{S}_2 + \gamma_2 \bar{S}_2 + \mu_2 \bar{S}_2 = \gamma_1 \bar{S}_1, \\
& \gamma_2 \bar{S}_2 = \mu_3 \bar{S}_3, \\
& \bar{\lambda}(\bar{A}, \bar{N}) \bar{A} (\beta_1 \bar{S}_1 + \beta_2 \bar{S}_2) = \mu_I \bar{I}, \\
& r b \bar{A} \left(1 - \frac{\bar{E}}{K}\right) = (\nu_E + \mu_E) \bar{E}, \\
& \nu_E \bar{E} = (\nu_L + \mu_L) \bar{L}, \\
& \nu_L \bar{L} = \frac{\mu_A \bar{A}}{1 + \theta_1 \bar{S}_1 + \theta_2 \bar{S}_2}.
\end{aligned}
\right.
\qquad (3.9) $$

Solving system (3.9) in terms of $\bar{S}_1$ and $\bar{S}_2$ gives

$$ \bar{E} = K \left(1 - \frac{1}{\bar{R}_0}\right), \quad \bar{L} = \frac{\nu_E K (\bar{R}_0 - 1)}{\bar{R}_0 (\nu_L + \mu_L)}, \quad \bar{A} = \frac{K (\nu_E + \mu_E) (\bar{R}_0 - 1)}{rb} \qquad (3.10) $$

$$ \text{and} \quad \bar{I} = \frac{\Lambda + (\mu_1 + \alpha \bar{S}_1) \bar{S}_1 + (\gamma_2 + \mu_2) \bar{S}_2}{\mu_I}, $$

where

$$ \bar{R}_0 = \frac{r b \nu_L \nu_E (1 + \theta_1 \bar{S}_1 + \theta_2 \bar{S}_2)}{\mu_A (\nu_E + \mu_E) (\nu_L + \mu_L)}. \qquad (3.11) $$

Now, plugging equation (3.10) into the first and second equations of (3.9) gives

$$ \bar{S}_1 = \bar{S}_1 [f(\bar{S}_1, \bar{S}_2)] \quad \text{and} \quad \bar{S}_2 = \bar{S}_2 [g(\bar{S}_1, \bar{S}_2)] \qquad (3.12) $$

where

$$ \begin{aligned} f(\bar{S}_1, \bar{S}_2) &= \frac{\sigma_S \sigma_A \mu_3 \mu_I K (\nu_E + \mu_E) (\bar{R}_0 - 1) \beta_1 \bar{S}_1 + \Delta (\mu_1 + \gamma_1 + \alpha \bar{S}_1) \bar{S}_1}{\Lambda \Delta}, \\ g(\bar{S}_1, \bar{S}_2) &= \frac{\sigma_S \sigma_A \mu_3 \mu_I K (\nu_E + \mu_E) (\bar{R}_0 - 1) \beta_2 \bar{S}_2 + \Delta (\mu_2 + \gamma_2) \bar{S}_2}{\Delta \gamma_1 \bar{S}}. \end{aligned} \qquad (3.13) $$

with

$$ \Delta = \sigma_S r b [\bar{S}_1(\mu_I + \mu_1 + \alpha \bar{S}_1)\mu_3 + \bar{S}_2(\mu_I(\gamma_2 + \mu_3) + \mu_3(\gamma_2 + \mu_2)) + \Lambda \mu_3] + \sigma_A \mu_3 \mu_I K (\nu_E + \mu_E)(\bar{R}_0 - 1), $$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Existence of a pest equilibrium point for system (2.6). In (a) and (b), one can see there is a curve in the $(\bar{S}_1, \bar{S}_2)$ which plane along $f(\bar{S}_1, \bar{S}_2) = 1$ and $g(\bar{S}_1, \bar{S}_2) = 1$ respectively. (c) There is a unique point at which $f(\bar{S}_1, \bar{S}_2) = g(\bar{S}_1, \bar{S}_2) = 1$. This point determines a pest equilibrium point $\bar{Q}$. (d) The contour curves $f(\bar{S}_1, \bar{S}_2) = 1$ and $g(\bar{S}_1, \bar{S}_2) = 1$, and there is a unique intersection point. We choose for numerical simulations $\Lambda = 14,400$; $\gamma_1 = 0.05$; $\gamma_2 = 0.027$; $\alpha = 0.001$; $\sigma_S = 20$; $\sigma_A = 26$; $\beta_1 = 0.09$; $\beta_2 = 0.02$; $r = 0.58$; $b = 3.28$; $\nu_E = 1/15$; $\nu_L = 1/25$; $\mu_1 = \mu_I = 0.05$; $\mu_2 = \mu_3 = 0.00469$; $\mu_E = 0.001$; $\mu_L = 0.03$; $\mu_A = 0.07$; $\theta_1 = \theta_2 = 0.00001$ (so that $\mathcal{R}_0 = 20.2345 > 1$).</figcaption></figure>

Since we are looking for a pest equilibrium point $\bar{Q}$ such that $\bar{S}_1 \neq 0$ and $\bar{S}_2 \neq 0$, equation (3.12) can be simplified as

$$f(\bar{S}_1, \bar{S}_2) = 1 \quad \text{and} \quad g(\bar{S}_1, \bar{S}_2) = 1. \qquad (3.14)$$

From equations (3.13) and (3.14), the interior pest equilibrium point corresponds to the intersection point $(\bar{S}_1, \bar{S}_2)$ of the two curves $f(\bar{S}_1, \bar{S}_2) = 1$ and $g(\bar{S}_1, \bar{S}_2) = 1$ with $\bar{S}_1 > 0$ and $\bar{S}_2 > 0$. Since equation (3.14) is very difficult to solve analytically due to the high non linearity of $f$ and $g$, we can numerically plot these two curves and examine how the intersection points change with system.

**Remark 3.5.** It is straightforward to show that $\bar{\mathcal{R}}_0$ is less than $\mathcal{R}_0$, so that $\bar{\mathcal{R}}_0 > 1$ implies $\mathcal{R}_0 > 1$.
<table><caption>TABLE 4. Estimation of production losses.</caption><thead><tr><th></th><th>Prod. without Mirids</th><th>Prod. with Mirids</th><th>Losses</th></tr></thead><tbody><tr><td>Ripe pods</td><td>51 000</td><td>31 000</td><td>39.21%</td></tr><tr><td>Adults miridae</td><td>0</td><td>3150</td><td></td></tr></tbody></table>

Figure 7 illustrates the existence of an interior pest equilibrium point $\bar{Q}$ when $\mathcal{R}_0 > 1$. From this figure, the surfaces $f(\bar{S}_1, \bar{S}_2)$ and $g(\bar{S}_1, \bar{S}_2)$ are plotted and the curves $f(\bar{S}_1, \bar{S}_2) = 1$ and $g(\bar{S}_1, \bar{S}_2) = 1$ are shown as intersections of the surface with the plane (see Fig. 7c and d). Figure 7b illustrates that there is a unique point $(\bar{S}_1, \bar{S}_2)$ at which $f(\bar{S}_1, \bar{S}_2) = g(\bar{S}_1, \bar{S}_2) = 1$. So, we have numerically shown the following result.

**Lemma 3.6.** System (2.6) has a unique pest equilibrium point when $\mathcal{R}_0 > 1$.

Now, we study the local asymptotic stability of the pest equilibrium point $\bar{Q}$. To do this, we use the Central Manifold Theory [9] as described in Theorem of Castillo- Chavez and Song [10].

**Theorem 3.7.** System (2.6) exhibits a bifurcation at $\mathcal{R}_0 = 1$ if the quantity $a_1$ defined as in equation (A.10) is negative, otherwise there exists a unique pest equilibrium point $\bar{Q}$ which is locally asymptotically stable in $\Omega \setminus \{Q_0\}$ for $\mathcal{R}_0 > 1$, but close to 1 and miridae population reached the highest possible level.

The proof is given in Appendix A.

For the rest of numerical simulations, we choose $\Lambda = 14,400$; $\gamma_1 = 0.05$; $\gamma_2 = 0.027$; $\alpha = 0.001$; $\sigma_S = 20$; $\sigma_A = 26$; $\beta_1 = 0.09$; $\beta_2 = 0.02$; $r = 0.58$; $b = 3.28$; $\nu_E = 1/15$; $\nu_L = 1/25$; $\mu_1 = \mu_I = 0.05$; $\mu_2 = \mu_3 = 0.00469$; $\mu_E = 0.001$; $\mu_L = 0.03$; $\mu_A = 0.07$; $\theta_1 = \theta_2 = 0.00001$ (so that $\mathcal{R}_0 = 20.2345 > 1$). Initial conditions are chosen to be $S_1(0) = 1200$, $S_2(0) = S_3(0) = I(0)$, $E(0) = 300$, $L(0) = 360$ and $A(0) = 600$.

Figure 8 presents the trajectories of system (2.6) when $\mathcal{R}_0 > 1$. From this figure, it clearly appears that the trajectories of system (2.6) converge to the pest equilibrium point $\bar{Q}$ as shown in Theorem 3.7. This means that the pest persists within the cocoa plot which leads to the production losses.

### 3.3. Estimation of production losses caused by Miridae

Herein, we provide and discuss some numerical simulations to estimate the production losses caused by Miridae within a cocoa plot.

Figure 9 shows the time evolution of ripped pods and adults miridae population with and without miridae in the plot during a growth season. As expected, the number of harvested pods decreases considerably. We will therefore record these results in Table 4 in order to estimate the production losses of cocoa within a plot.

Table 4 recapitulates the number of ripped pods with and without the infestation by miridae using parameters values in Figs 6 and 8 and the number of adults miridae presents in the plot. From this table, it is evident that the production losses is approximately 39.21%. One can observe that about 3, 150 adults miridae are present in the cacao plot.

According to the estimation of production losses and the presence level of miridae within a cocoa plantation, it is then important to reduce the production losses caused by miridae and the presence level of miridae. This is the objective of the next section.

## 4. CONTROLLING PRODUCTION LOSSES CAUSED BY MIRIDAE USING
OPTIMAL CONTROL THEORY

In this section, we formulate a continuous optimal control problem which consists in minimizing the damages caused by miridae within a cocoa plot. The optimal control approach is pretty common to investigate the best control strategy to reduce the presence of miridae and thus increase the production within a cocoa plantation.
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Numerical simulations for system (2.6) when $\Lambda = 14,400$; $\gamma_1 = 0.05$; $\gamma_2 = 0.027$; $\alpha = 0.001$; $\sigma_S = 20$; $\sigma_A = 26$; $\beta_1 = 0.09$; $\beta_2 = 0.02$; $r = 0.58$; $b = 3.28$; $\nu_E = 1/15$; $\nu_L = 1/25$; $\mu_1 = \mu_I = 0.05$; $\mu_2 = \mu_3 = 0.00469$; $\mu_E = 0.001$; $\mu_L = 0.03$; $\mu_A = 0.07$; $\theta_1 = \theta_2 = 0.00001$ (so that $\mathcal{R}_0 = 20.2345 > 1$). (a) Time evolution of cherelle $S_1$; (b) Young pods $S_2$; (c) Ripped pods $S_3$; (d) Infected pods $I$; (e) Eggs Miridae $E$; (f) Nymphs Miridae $L$ and (g) Adults Miridae $A$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Ripped pods and adults Miridae with and without Miridae in the plot. Without miridae (blue line) and with miridae (magenta line). (a) Ripped pods $S_3$; (b) adults miridae A.</figcaption></figure>

## 4.1. Optimal Control model

Herein, the interaction of pods and miridae is formulated in an optimal control model where the control objective is to determine the rate of control such that the pest population keep under the injury level. There are many strategies of control like cultural management, varietal management, chemical management, semiochemical management. In this paper, the control u represents chemical management, based on the use of chemical insecticide. This control is the most widespread to reduce the miridae population in most countries producers of cocoa like Cameroon. Note that the synthetic insecticides like λ-cyhalothrine and imidacloprid have a long residual effect, but it depends on several environmental factors, like rainfall and is not usually used by cocoa farmers in most countries producers of cocoa. We stress that the control action on the populations of adults and nymphs should be different. In order to show the effects of the control action on the populations of adults and nymphs, we have introduced two positive parameters $\varepsilon_1$ and $\varepsilon_2$ ($\varepsilon_1 > \varepsilon_2$) which are both less than one in order to differentiate the effects of control on those two populations.

Using the same parameters and classes as in system (2.6), the system of differential equations describing the controlled model is

$$ 
\left\{
\begin{aligned}
\dot{S}_1 &= \Lambda - \beta_1 \lambda(N, A) A S_1 - (\gamma_1 + \mu_1 + \alpha S_1)S_1, \\
\dot{S}_2 &= \gamma_1 S_1 - \beta_2 \lambda(N, A) A S_2 - (\gamma_2 + \mu_2) S_2, \\
\dot{S}_3 &= \gamma_2 S_2 - \mu_3 S_3 \\
\dot{I} &= \lambda(N, A) A (\beta_1 S_1 + \beta_2 S_2) - \mu_I I, \\
\dot{E} &= r b A \left(1 - \frac{E}{K}\right) - (\nu_E + \mu_E) E, \\
\dot{L} &= \nu_E E - (\nu_L + \mu_L) L - \varepsilon_1 u L, \\
\dot{A} &= \nu_L L - \frac{\mu_A}{1 + \theta_1 S_1 + \theta_2 S_2} A - \varepsilon_2 u A,
\end{aligned}
\right.
\quad (4.1) $$
where $\lambda(N, A)$ is defined as in equation (2.5).

The function $0 \le u \le 1$ represents the vector control to reduce the number of nymphs and adults population in the plot. Products use for control are lethal for adults and nymphs. However, adults can still feed and lay on pods.

The given objective function is

$$C(u) = \int_{0}^{T} (B_{1} L(t) + B_{2} A(t) + B_{3} u^{2}(t)) dt = \int_{0}^{T} g(t, X, u)dt, \quad (4.2)$$

where $T$ is the final time and the coefficients $B_1, B_2$ and $B_3$ are positive weights to balance the factors. The terms $B_1 L$ and $B_2 A$ are the cost of Miridae attacks while $B_3 u^2$ is the cost of treatment efforts. Our aim is to minimize the number of nymphs $L(t)$ and adults $A(t)$, while minimizing the cost of control $u(t)$. Thus, we seek an optimal control $u^*$ such that

$$C(u^*) = \min_{u} \{C(u), u \in \mathcal{U}\}, \quad (4.3)$$

where $\mathcal{U}$ is the control set defined as:

$$\mathcal{U} = \{u \in L^1(0, T), 0 \le u \le 1\}. \quad (4.4)$$

The control is a quadratic form.

The existence of the optimal control can be obtained by using a result by Fleming and Rishel [16]. We have the following result.

**Theorem 4.1.** Consider the control problem with system (4.1), there exists $u^* \in \mathcal{U}$ such that

$$\min_{u \in \mathcal{U}} J(u) = J(u^*).$$

*Proof.* To use an existence result, we must check the following properties:

1. The set of controls and corresponding state variables is nonempty. By Theorem 3.2, we have the existence of solution of system (2.6) with bounded coefficients, which gives condition 1.

2. The set of control $\mathcal{U}$ is convex and closed. We note that the solutions are bounded. Our control set satisfies condition 2, so $\mathcal{U}$ is closed and convex set.

3. The right hand side of the state system is bounded by a linear function in the state and control variables. By definition, each right hand side of system (2.6) is continuous and can be written as a linear function of $u$ with coefficients depending on time and state.

4. The integrand of the objective functional is convex on $\mathcal{U}$.

5. There exist $c_1, c_2 > 0$ and $\beta > 1$ such that the integrand $g(t, X, u)$ of the objective functional satisfies

$$g(t, X, u) \ge -c_1 + c_2(|u|^2)^{\beta/2}.$$

Indeed,

$$\begin{align*}
g(t, X, u) &= B_1 L + B_2 A + B_3 u^2 \\
&\ge B_3 u^2 \\
&\ge -c_1 + B_3(|u^2|)^1
\end{align*}$$
Take $c_1 \in \mathbb{R}_+$, $c_2 = B_3$ and $\beta = 2$.

Thus, there exists and optimal control. This achieves the proof.

The necessary conditions that an optimal control must satisfy come from Pontryagin's Minimum principle [24] converts into a problem of minimizing point-wise a Hamiltonian $H$, with respect to $u$:

$$ 
\begin{aligned}
H = B_1 L + B_2 A + B_3 u^2 &+ \lambda_1 \left( \Lambda - \frac{\beta_1 \sigma_S \sigma_A A S_1}{\sigma_S N + \sigma_A A} - (\gamma_1 + \mu_1 + \alpha S_1) S_1 \right) \\
&+ \lambda_2 \left( \gamma_1 S_1 - \frac{\beta_2 \sigma_S \sigma_A A S_2}{\sigma_S N + \sigma_A A} - (\gamma_2 + \mu_2) S_2 \right) + \lambda_3 (\gamma_2 S_2 - \mu_3 S_3) \\
&+ \lambda_4 \left( \frac{\sigma_S \sigma_A A (\beta_1 S_1 + \beta_2 S_2)}{\sigma_S N + \sigma_A A} - \mu_I I \right) + \lambda_5 \left( r b A \left( 1 - \frac{E}{K} \right) - (\nu_E + \mu_E) E \right) \\
&+ \lambda_6 (\nu_E E - (\nu_L + \mu_L) L - \varepsilon_1 u L) + \lambda_7 \left( \nu_L L - \frac{\mu_A A}{1 + \theta_1 S_1 + \theta_2 S_2} - \varepsilon_2 u A \right),
\end{aligned} 
\quad (4.5) $$

where $\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6$ and $\lambda_7$ are the adjoint functions associated with their respective states. Note that in $H$, each adjoint function multiplies the right-hand side of the differential equation of its corresponding state function. The first terms in $H$ comes from the integrand of the objective functional. Thus, the adjoint variable $\lambda_j; j \in \{S_1, S_2, S_3, I, E, L, A\}$ together with our state system determine our optimality system. We have the following result.

**Theorem 4.2.** Given an optimal control $u^*$ and the corresponding states $S_1, S_2, S_3, I, E, L$ and $A$ that minimizes $C(u)$ over $U$, there exists adjoint variables $\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6$ and $\lambda_7$ satisfying:

$$ 
\left\{
\begin{aligned}
\dot{\lambda}_1 &= \lambda_1(\mu_1 + 2\alpha S_1) + \gamma_1(\lambda_1 - \lambda_2) - (\lambda_1 - \lambda_4) \left( \frac{\beta_1 \sigma_S \sigma_A A}{\sigma_S N + \sigma_A A} - \frac{\beta_1 \sigma_S^2 \sigma_A A S_1}{(\sigma_S^2 N + \sigma_A A)^2} \right) \\
&\quad - (\lambda_2 - \lambda_4) \frac{\beta_2 \sigma_S^2 \sigma_A A S_2}{(\sigma_S N + \sigma_A A)^2} - \lambda_7 \frac{\mu_A \theta_1 A}{(1 + \theta_1 S_1 + \theta_2 S_2)^2}, \\
\dot{\lambda}_2 &= \lambda_2 \mu_2 + \gamma_2(\lambda_2 - \lambda_3) - (\lambda_2 - \lambda_4) \left( \frac{\beta_2 \sigma_S^2 \sigma_A A}{\sigma_S N + \sigma_A A S_2} - \frac{\beta_2 \sigma_S \sigma_A A S_2}{(\sigma_S^2 N + \sigma_A A)^2} \right) \\
&\quad - (\lambda_1 - \lambda_4) \frac{\beta_1 \sigma_S^2 \sigma_A A S_1}{(\sigma_S N + \sigma_A A)^2} - \lambda_7 \frac{\mu_A \theta_2 A}{(1 + \theta_1 S_1 + \theta_2 S_2)^2}, \\
\dot{\lambda}_3 &= \lambda_3 \mu_3 - \frac{\sigma_S^2 \sigma_A \beta_1 A S_1}{(\sigma_S N + \sigma_A A)^2} (\lambda_1 - \lambda_4) - \frac{\sigma_S^2 \sigma_A \beta_2 A S_2}{(\sigma_S N + \sigma_A A)^2} (\lambda_2 - \lambda_4), \\
\dot{\lambda}_4 &= \lambda_4 \mu_I - \frac{\sigma_S^2 \sigma_A \beta_1 A S_1}{(\sigma_S N + \sigma_A A)^2} (\lambda_1 - \lambda_4) - \frac{\sigma_S^2 \sigma_A \beta_2 A S_2}{(\sigma_S N + \sigma_A A)^2} (\lambda_2 - \lambda_4), \\
\dot{\lambda}_5 &= \lambda_5 \left( \frac{r b A}{K} + \mu_E \right) + (\lambda_5 - \lambda_6) \nu_E, \\
\dot{\lambda}_6 &= -B_1 + \lambda_6 (\varepsilon_1 u + \mu_L) + (\lambda_6 - \lambda_7) \nu_L, \\
\dot{\lambda}_7 &= -B_2 + \lambda_7 \left( \varepsilon_2 u + \frac{\mu_A}{1 + \theta_1 S_1 + \theta_2 S_2} \right) + \left( \frac{\beta_1 \sigma_S \sigma_A A}{\sigma_S N + \sigma_A A} - \frac{\beta_1 \sigma_S^2 \sigma_A A S_1}{(\sigma_S^2 N + \sigma_A A)^2} \right) (\lambda_1 - \lambda_4) \\
&\quad + \left( \frac{\beta_2 \sigma_S \sigma_A A}{\sigma_S N + \sigma_A A} - \frac{\beta_2 \sigma_S^2 \sigma_A A S_2}{(\sigma_S^2 N + \sigma_A A)^2} \right) (\lambda_2 - \lambda_4) - \lambda_5 r b \left( 1 - \frac{E}{K} \right),
\end{aligned}
\right. 
\quad (4.6) $$
with the following transversality conditions:

$$ \lambda_1(T) = \lambda_2(T) = \lambda_3(T) = \lambda_4(T) = \lambda_5(T) = \lambda_6(T) = \lambda_7(T) = 0, \quad (4.7) $$

where $T = 50$ days is the free terminal time. Furthermore, the optimal control is characterized by

$$ u^* = \max \left\{ 0, \min \left\{ 1, \frac{\lambda_6 L^* + \lambda_7 A^*}{2 B_3} \right\} \right\}. \quad (4.8) $$

*Proof.* The differential equations governing the adjoint variables are obtained by differentiation of the Hamiltonian function, evaluated at the optimal control:

$$ \dot{\lambda}^*(t) = -\frac{\partial H}{\partial X}(t, X^*, \lambda^*(t), u^*) $$

$$ \frac{\partial H}{\partial u}(t, X^*, \lambda^*(t), u^*) = 0 $$

The control characterization for $u^*$ comes from $\frac{\partial H}{\partial u} = 0$ whenever $0 < u^* < 1$ and taking bounds into account.

$$ \frac{\partial H}{\partial u}(t, X^*, \lambda(t)^*, u^*) = 0 \implies u^* = \begin{cases} \frac{\lambda_6 L^* + \lambda_7 A^*}{2 B_3} & \text{if } 0 < \frac{\lambda_6 L^* + \lambda_7 A^*}{2 B_3} < 1, \\ 0 & \text{if } \frac{\lambda_6 L^* + \lambda_7 A^*}{2 B_3} \le 0, \\ 1 & \text{if } \frac{\lambda_6 L^* + \lambda_7 A^*}{2 B_3} \ge 1. \end{cases} $$

Thus

$$ u^* = \max \left\{ 0, \min \left\{ 1, \frac{\lambda_6 L^* + \lambda_7 A^*}{2 B_3} \right\} \right\}. $$

This concludes the proof.

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

Next, we discuss the numerical solutions of the optimality system and the corresponding optimal control, the parameter choices and the interpretations from various cases.

## 4.2. Optimal control numerical simulations

Herein, we study numerically an optimal vector parameter control for system (4.1). Optimal control is obtained by solving optimality system, consisting of seven (07) differential equations from the state and the adjoint. An iterative scheme is used for solving the optimality system. We start to solve the state equations with a guess for the controls over the simulated time using Euler scheme. Because of the transversality conditions (4.7), the adjoint equations are solved by a backward Euler scheme using the current iterations solutions of the state equation. Then, the control is updated by using a convex combination of the previous control and the value from the characterization (4.8). This process is repeated and iterations are stopped if the values of the unknowns at the previous iterations are very close to the ones at the present iterations [20].

For the numerical simulations, we used the following weight factors $B_1 = 0.15$, $B_2 = 0.1$ and $B_3 = 0.015$. We have the following initial conditions. At the beginning of the season, there are no pods in the plot so $S_1(0) = S_2(0) = S_3(0) = I(0) = 0$, $E(0) = L(0) = 0$ and $A(0) = 960$ corresponding to the presence of miridae
<figure><img src="image_13.png" /><figcaption>FIGURE 10. Simulations of system (4.1) without control (red line) and with control (blue line).</figcaption></figure>

(a) Cherelles $S_1$; (b) Young pods $S_2$; (c) infected pods $I$; (d) ripped pods $S_3$.

in the plot. The parameter values used are $\beta_1 = 0.09$, $\beta_2 = 0.02$, $\gamma_1 = 0.05$, $\gamma_2 = 0.027$, $\sigma_S = 20$, $\sigma_A = 26$, $\mu_1 = 0.05$, $\mu_2 = \mu_3 = 0.00469$, $\mu_I = 0.05$, $\alpha = 0.001$, $r = 0.58$, $b = 3.28$, $\nu_L = 1/25$, $\nu_E = 1/15$, $\mu_L = 0.03$, $\mu_E = 0.001$, $\mu_A = 0.07$, $\theta_1 = 0.0001$, $\theta_2 = 0.0001$, $\Lambda = 14400$, $K = 5000$, $\varepsilon_1 = 0.75$ and $\varepsilon_2 = 0.25$ (so that $R_0 = 20.2365$) to illustrate the optimal control strategy.

Figure 10 shows the time evolution of cherelles, young pods, infected pods and ripped pods with and without control. During the control period, it is evident that the number of infected pods remains at zero and begins to increase when the control is stopped. However, throughout the growth period, it remains lower than the number of infected pods present in the plot in the absence of control. Throughout the production period, the number of harvested pods and young pods is greater when considering the control in the plantation even if the control is only carried out over a short period.

The time evolution of miridae population with and without control is depicted in Figure 11. From this figure, one can observe that eggs, nymphs and adults, miridae population decreases to zero during the control period and begins to grow when the control is stopped. However, throughout the growth period, each population remains lower or equal than the level in the absence of control. One can also observe that the effect of the control on the population of nymphs is more significant than on the population of adults. However, a decrease of the population of nymphs will indirectly induce a decrease on the population of adults. Thus, the optimal control strategy is very efficient to increase pods production and reduce the miridae population in the plot.
<figure><img src="image_14.png" /><figcaption>FIGURE 11. Simulations of system (4.1) without control (red line) and with control (blue line).</figcaption></figure>

(a) Eggs Miridae *E*; (b) nymphs *L*; (c) adults *A*.

Figure 12 presents the time evolution of vector control function *u*(*t*). Vector control is applied during 50 days and after it is supposed to be equal to zero.

Numerical simulations were made over 180 days corresponding to the growth of the pods: from the appearance of cherelle to the harvest. The control was done over 50 days, the first days after the appearance of cherelle. After this time, the control is assumed to be equal to zero.

4.2.1. *Estimation of production losses*

Herein, we estimate the production losses in the presence under the control action.

Figure 13, page 27 presents the number of ripped pods and adults miridae in the plot with and without miridae under the vector control in the plot. As expected, the number of ripped pods in the plot without miridae is greater than the number of ripped pods with miridae, under or not the control. We recapitulate these results in Table 5 in order to estimate the production. From Table 5, it is evident that the production losses with control in the plot is approximately 20.58% whereas the production losses without control is 39.21%.

Figure 14, page 28 presents the effect of control on nymphs and adults miridae. As expected, under the vector control, nymphs and adults population decreases rapidly toward zero. When control is stopped, nymphs and adults miridae population gradually grow and converge towards the uncontrolled dynamics (see Fig. 11).
<figure><img src="image_15.png" /><figcaption>FIGURE 12. Optimal control $u(t)$ for system (4.1).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 13. Ripped pods and adults miridae control numerical simulations for system (4.1) without Miridae (green line), with Miridae and control (blue line) and without Miridae and without control (red line). (a) Ripped pods $S_3$, (b) adults miridae $A$.</figcaption></figure>

<table><caption>TABLE 5. Estimation of production losses under the control action.</caption><thead><tr><th></th><th>Without<br>Miridae</th><th>Production<br>without control</th><th>Production<br>with control</th><th>Losses<br>without control</th><th>Losses<br>with control</th></tr></thead><tbody><tr><td>Ripe pods</td><td>51 000</td><td>31 000</td><td>40 500</td><td>39.21%</td><td>20.58%</td></tr></tbody></table>

Figures 15 presents the production $S_3(t)$ during the season without miridae (green, colour), for three different initial values of miridae in the plot but no control (red colour), and with initial miridae and optimal control (blue) and the curve of optimal control for several initial miridae in the plot. It seems that the initial value of miridae in the plot do not impact the efficiency of vector control. Indeed, for these different values of initial miridae in the plot, the production with control seems to be the same at the end of the season.
<figure><img src="image_17.png" /><figcaption>FIGURE 14. Effect of vector control on nymphs (a) and adults miridae population (b).</figcaption></figure>

<table><caption>TABLE 6. Estimation of losses of gain.</caption><thead><tr><th></th><th>Without<br>Miridae</th><th>Without<br>control</th><th>With<br>control</th><th>Losses<br>without control</th><th>Losses<br>with control</th></tr></thead><tbody><tr><td>Ripe pods</td><td>51 000</td><td>31 000</td><td>40 500</td><td>39.21%</td><td>20.58%</td></tr><tr><td>Commercial<br>cocoa</td><td>2142 kg</td><td>1302 kg</td><td>1701 kg</td><td>840 kg</td><td>441 kg</td></tr><tr><td>Farmer's<br>gain</td><td>3255.84 USD</td><td>1979.04 USD</td><td>2585.52 USD</td><td>1276.8 USD</td><td>670.32 USD</td></tr></tbody></table>

### 4.2.2. Estimation of losses of gain

Annual production is estimated at 51 000 pods harvested in the absence of mirids and 31 000 pods in the presence of mirids. A pod contains approximately 30–40 beans. The weight of a bean after elimination of the pulp and husk varies between 1.3 g and 2.3 g; after drying, the weight of the bean varies between 0.9 g and 1.5 g. Fermentation and drying of the beans of the pods lead to obtaining of commercial cacao [18]. If we assume that a pod contains on average 35 beans and that a commercial cocoa bean weighs 1.2 g on average, a pod makes it possible to obtain 42g of commercial cocoa and therefore the production in a plot of 1200 plants is evaluated at 2,142 kg of market cocoa in the absence of mirids and 1302 kg in the presence of mirids. A the end of the season, the financial gain for the farmer is evaluated by the number of kilogramme of beans obtained. The price of a kilogram of commercial cocoa is estimated to 1025 FCFA or 1.52 USD in Cameroon on September 28, 2022 (ONCC).

Table 6 recapitulates the losses of gain obtained by a farmer. The losses of gain with control in the plot is approximately 670.32 USD whereas the losses of gain obtained without control is 1276.8 USD for 1,200 plants in the plot.

## 5. CONCLUSION AND DISCUSSIONS

In this paper, we have formulated a mathematical model for the interaction between pods and miridae in which we have incorporated an optimal control. The objective was to evaluate and control the losses of production due to miridae. To do this, we formulated a controlled model based on vector control to reduce pests population and increase production. A qualitative analysis of the models has been presented and our main findings can be summarized as follows:
<figure><img src="image_18.png" /><figcaption>FIGURE 15. Simulation of system (2.6) ripped pods and adults Miridae with optimal control for initial conditions of adults Miridae populations. Without control (green colour) and system (4.1) (a) $A(0) = 100$; (b) $u(t)$ when $A(0) = 100$; (c) $A(0) = 2000$; (d) $u(t)$ when $A(0) = 2000$; (e) $A(0) = 4000$; (f) $u(t)$ when $A(0) = 4000$.</figcaption></figure>
1. The system is well posed and we proved the basic properties: positivity and boundedness of solutions, the existence and uniqueness of solutions.

2. For the model without control, we have two equilibria $Q^0$ and $\bar{Q}$. we show that there exists a threshold parameter, $\mathcal{R}_0$, also called the basic offspring number, that summarizes the dynamics of the system. Equilibrium $Q^0$ corresponds to the case when there is no pest in the plot and equilibrium $\bar{Q}$ corresponds to the case where there are pest and pods in the plantation. We proved that equilibrium $Q^0$ is GAS whenever $\mathcal{R}_0$ is less than one and unstable when $\mathcal{R}_0$ is greater than one. Numerical simulations show that equilibrium $\bar{Q}$ is unique and we proved that it is LAS when $\mathcal{R}_0$ is greater than one.

3. As an application, another objective of our work was to study optimal control model. Optimal control $u$ consists in using vector control to reduce nymphs and adults and increase production in the plot. Optimal control model was studied and we determined an optimal control strategy $u^*$.

4. Numerical simulations have been presented to illustrate the obtained theoretical results. Through numerical simulations, we found that

(a) The sensitivity analysis of the model and sensitivity analysis of the basic offspring number $\mathcal{R}_0$ have been investigated

(b) The presence of miridae in the plot can cause about 39.21% of losses of ripped pods within a plantation where 1200 plants.

(c) Optimal control help us to reduce the losses of production to 20.58% and is independent of initial miridae values in the plot.

(d) Numerical simulations tend to show that the effect of control action on the population of nymphs is more significant than on the population of adults.

(e) The initial value of miridae in the plot do not impact the efficiency of vector control. Indeed, for these different values of initial miridae in the plot, the production with control seems to be the same at the end of the season.

(f) For a plantation of 1200 cocoa trees, the farmer loses an average of 1276.8 USD when he leaves his plot under miridae effects. When he applies optimal control (based on the elimination of nymphs and adults miridae), the loss is 670.32 USD. Optimal control therefore contributes to the increase of production and consequently to the increase of gain.

Different improvements and extensions of this work can include introducing time-dependent parameters in order to integrate the fluctuation of environmental factors due to periodic variations of climate; to formulate and study a multi seasonal model with optimal control in order to better estimate the production. This multi seasonal model with optimal control can take into account all the recommended periods of treatment of the plot along the year.

## APPENDIX A. PROOF OF THEOREM 3.7

In this appendix, we give the proof of Theorem 3.7 on the local asymptotic stability of the pest equilibrium point $\bar{Q}$ of system (2.6). We first recall the theorem of Castillo-Chavez and Song [10].

**Theorem A.1.** (Castillo-Chavez and Song [10]) Let

$$ \frac{dx}{dt} = f(x, \phi), \quad f: \mathbb{R}^n \times \mathbb{R} \leftrightarrow \mathbb{R}^n \text{ and } f \in C^2(\mathbb{R}^n \times \mathbb{R}). \qquad (A.1) $$

Without loss of generality, it is assumed that $0$ is an trivial equilibrium for system (A.1) for all values of the parameter $\phi$, that is

$$ f(0, \phi) \equiv 0 \text{ for all } \phi. \qquad (A.2) $$

Assume
$A_1$ A = $D_x f(\mathbf{0}) = \left(\frac{\partial f_i}{\partial x_j}(\mathbf{0})\right)$ is the linearisation matrix of system (A.1) around the equilibrium $\mathbf{0}$ with $\phi$ evaluated at $\mathbf{0}$. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts;

$A_2$ Matrix A has a non-negative right eigenvector $w$ and a left eigenvector $v$ corresponding to the zero eigenvalue.

Let $f_k$ be the kth component of $f$ and

$$a_1 = \sum_{k,i,j=1}^{n} v_k w_i w_j \frac{\partial^2 f_k}{\partial x_i \partial x_j}(\mathbf{0}) \qquad (A.3)$$

$$b_1 = \sum_{k,i=1}^{n} v_k w_i \frac{\partial^2 f_k}{\partial x_i \partial \phi}(\mathbf{0}) \qquad (A.4)$$

The local dynamics of system (A.1) around $\mathbf{0}$ are totally determined by a and b.

(i) $a_1 > 0, b_1 > 0$. When $\phi < 0$ with $|\phi| \ll 1$, $\mathbf{0}$ is locally asymptotically stable, and there exists a positive unstable equilibrium; when $0 < \phi \ll 1$, $\mathbf{0}$ is unstable and there exists a negative and locally asymptotically stable equilibrium;

(ii) $a_1 < 0, b_1 < 0$. When $\phi < 0$ with $|\phi| \ll 1$, $\mathbf{0}$ is unstable; when $0 < \phi \ll 1$, $\mathbf{0}$ is locally asymptotically stable and there exists a positive unstable equilibrium;

(iii) $a_1 > 0, b_1 < 0$. When $\phi < 0$ with $|\phi| \ll 1$, $\mathbf{0}$ is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0 < \phi \ll 1$, $\mathbf{0}$ is stable, and a positive unstable equilibrium appears;

(iv) $a_1 < 0, b_1 > 0$. When $\phi$ changes from negative to positive, $\mathbf{0}$ changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

In order to apply Castillo-Chavez and Song theorem [10] on system (2.6), the following simplification and change of variables are first of all made.

Let $x_1 = S_1, x_2 = S_2, x_3 = S_3, x_4 = I, x_5 = 4, x_6 = L$ and $x_7 = A$. Further, by using the vector notation $x = (x_1, x_2, x_3, x_4, x_5)$, system (2.6) can be written in the form $\dot{x} = f(x)$, with $f = (f_1, f_2, f_3, f_4, f_5, f_6, f_7)$, as follows:

$$\left\{
\begin{aligned}
\dot{x}_1 &= f_1 = \Lambda - \lambda(A, N) \beta_1 x_1 x_7 - (\gamma_1 + \mu_1 + \alpha x_1)x_1, \\
\dot{x}_2 &= f_2 = \gamma_1 x_1 - \lambda(A, N) \beta_2 x_7 x_2 - (\gamma_2 + \mu_2)x_2, \\
\dot{x}_3 &= f_3 = \gamma_2 x_1 - \mu_3 x_3, \\
\dot{x}_4 &= f_4 = \lambda(A, N) x_7(\beta_1 x_1 + \beta_2 x_2) - \mu_I x_4, \\
\dot{x}_5 &= f_5 = r b x_7 \left(1 - \frac{x_5}{K}\right) - (\nu_E + \mu_E)x_5, \\
\dot{x}_6 &= f_6 = \nu_E x_5 - (\nu_L + \mu_L)x_6, \\
\dot{x}_7 &= f_7 = \nu_L x_6 - \frac{\mu_A}{1 + \theta_1 x_1 + \theta_2 x_2} x_7.
\end{aligned}
\right.
\qquad (A.5)$$
Consider $\mathcal{R}_0 = 1$. Suppose, further, that $\nu_E = \nu_E^*$ is chosen as a bifurcation parameter. Solving $\mathcal{R}_0 = 1$ give

$$ \nu_E = \nu_E^* = \frac{\mu_A \mu_E (\nu_L + \mu_L)}{rb(1 + \theta_1 S_1^0 + \theta_2 S_2^0) - \mu_A (\nu_L + \mu_L)}. $$

The Jacobian matrix around the pest free equilibrium $Q^0$ when $\nu_E = \nu_E^*$ is given by

$$ J(Q^0) = \begin{pmatrix} -(\gamma_1 + \mu_1 + 2\alpha S_1^0) & 0 & 0 & 0 & 0 & 0 & \lambda_1 \beta_1 S_1^0 \\ \gamma_1 & -(\gamma_2 + \mu_2) & 0 & 0 & 0 & 0 & \lambda_1 \beta_2 S_2^0 \\ 0 & \gamma_2 & -\mu_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\mu_I & 0 & 0 & \lambda_1 (\beta_1 S_1^0 + \beta_2 S_2^0) \\ 0 & 0 & 0 & 0 & -(\nu_E + \mu_E) & 0 & rb \\ 0 & 0 & 0 & 0 & \nu_E & -(\nu_L + \mu_L) & 0 \\ 0 & 0 & 0 & 0 & 0 & \nu_L & -\frac{\mu_A}{1 + \theta_1 S_1^0 + \theta_2 S_2^0} \end{pmatrix}, $$

where

$$ \lambda_1 = \frac{\sigma_S \sigma_A}{\sigma_S(S_1^0 + S_2^0 + S_3^0)}. $$

When $\nu_E = \nu_E^*$, it is straightforward to show that 0 is a simple eigenvalue of $J(Q_0)$ (all other eigenvalues have negative real parts). In order to apply Castillo-Chavez and Song theorem [10], we need to compute:

* **Eigenvectors of J under the condition $\nu_E = \nu_E^*$**. When $\mathcal{R}_0 = 1$, it can be shown that the Jacobian matrix $J(Q^0)$ has a right eigenvector (corresponding to the zero eigenvalue), given by $u = (u_1, u_2, u_3, u_4, u_5, u_6, u_7)$, where,

$$ \begin{align} u_1 &= \frac{\lambda_1 \beta_1 S_1^0}{\gamma_1 + \mu_1 + 2\alpha S_1^0}, & u_2 &= \frac{\lambda_1 \gamma_1 \beta_1 S_1^0}{(\gamma_1 + \mu_1 + 2\alpha S_1^0)(\gamma_2 + \mu_2)} + \frac{\lambda_1 \beta_2 S_2^0}{\gamma_2 + \mu_2}, & u_3 &= \frac{\gamma_2}{\mu_3} u_2, \tag{A.6} \\ u_4 &= \frac{\lambda_1(\beta_1 S_1^0 + \beta_2 S_2^0)}{\mu_I}, & u_5 &= \frac{rb}{\nu_E^* + \mu_E}, & u_6 &= \frac{\mu_A}{\nu_L(1 + \theta_1 S_1^0 + \theta_2 S_2^0)} & u_7 &= 1. \notag \end{align} $$

Similarly, the components of the left eigenvector (corresponding to the zero eigenvalue for $J^T(Q_0)$) denoted $v = (v_1, v_2, v_3, v_4, v_5, v_6, v_7)$ are given by

$$ v_1 = v_2 = v_3 = v_4 = 0, \quad v_5 = 1, \quad v_6 = \frac{\nu_E^* + \mu_E}{\nu_E^*} \quad \text{and} \quad v_7 = \frac{(\nu_E^* + \mu_E)(\nu_L + \mu_L)}{\nu_E^* \nu_L}. \tag{A.7} $$

* **Computation of $b_1$:** For system (A.5), the associated non-zero partial derivatives $\partial f$ at $Q^0$ are given by

$$ \frac{\partial^2 f_5}{\partial x_5 \partial \nu_E} = -1 \quad \text{and} \quad \frac{\partial^2 f_6}{\partial x_5 \partial \nu_E} = 1. \tag{A.8} $$
Then, it follows that

$$
\begin{align*}
b_1 &= \sum_{i,k=1}^{7} v_k u_i \frac{\partial^2 f_k}{\partial x_i \partial \nu_E}(Q^0), \\
&= v_5 u_5 \frac{\partial^2 f_5}{\partial x_5 \partial \nu_E}(Q^0) + v_6 u_5 \frac{\partial^2 f_6}{\partial x_5 \partial \nu_E}(Q^0), \\
&= (v_6 - v_5) u_5, \\
&= \frac{r b \mu_E}{\nu_E^* + \mu_E} > 0.
\end{align*}
$$

• Computation of $a_1$: Since $v_1 = v_2 = v_3 = v_4 = 0$,

$$
a_1 = \sum_{i,j=1}^{7} v_5 u_i u_j \frac{\partial f_5}{\partial x_i \partial x_j} + v_6 u_i u_j \frac{\partial f_6}{\partial x_i \partial x_j} + v_7 u_i u_j \frac{\partial f_7}{\partial x_i \partial x_j}.
$$

For system (A.5), the associated non-zero partial derivatives of $f$ at $Q^0$ are given by:

$$
\begin{align}
\frac{\partial^2 f_5}{\partial x_5 \partial x_7} &= \frac{\partial^2 f_5}{\partial x_7 \partial x_5} = -\frac{rb}{K}, \quad \frac{\partial^2 f_7}{\partial x_1 \partial x_7} = \frac{\partial^2 f_5}{\partial x_7 \partial x_1} = \frac{\mu_A \theta_1}{(1 + \theta_1 x_1 + \theta_2 x_2)^2}, \tag{A.9} \\
\text{and } \quad \frac{\partial^2 f_7}{\partial x_2 \partial x_7} &= \frac{\partial^2 f_5}{\partial x_7 \partial x_2} = \frac{\mu_A \theta_2}{(1 + \theta_1 x_1 + \theta_2 x_2)^2}. \nonumber
\end{align}
$$

Then, it follows that

$$
\begin{align}
a_1 &= 2 v_5 u_5 u_7 \frac{\partial^2 f_5}{\partial x_5 \partial x_7} + 2 v_7 u_1 u_7 \frac{\partial^2 f_7}{\partial x_7 \partial x_1} + \dots + 2 v_7 u_2 u_7 \frac{\partial^2 f_7}{\partial x_7 \partial x_2}, \tag{A.10} \\
&= -2 \frac{r^2 b^2}{K(\nu_E^* + \mu_E)} + \frac{2\lambda_1 \beta_2 S_1^0 \mu_A \theta_1 (\nu_E^* + \mu_E)(\nu_L + \mu_L)}{\nu_E^* \nu_L (\gamma_1 + \mu_1 + 2\alpha S_1^0)(1 + \theta_1 x_1 + \theta_2 x_2)^2}, \nonumber \\
&\quad + 2 \left( \frac{\lambda_1 \gamma_1 \beta_1 S_1^0}{(\gamma_1 + \mu_1 + 2\alpha S_1^0)(\gamma_2 + \mu_2)} + \frac{\lambda_1 \beta_2 S_2^0}{\gamma_2 + \mu_2} \right) \left( \frac{(\nu_E^* + \mu_E)(\nu_L + \mu_L)}{\nu_E^* \nu_L} \right) \frac{\mu_A \theta_2}{(1 + \theta_1 S_1 + \theta_2 S_2)^2}. \nonumber
\end{align}
$$

Thus, depending on the values of the parameters of system (A.5), the value of $a_1$ can be positive or negative.
Since $b_1 > 0$, the conclusion follows from Theorem A.1 items (i) and (iv). This achieves the proof of Theorem 3.7.

Funding Information. The first author is grateful to the APSA foundation for financial support during the preparation and the finalization of this manuscript.

REFERENCES

[1] R. Adu-Acheampong, J. Jiggins, A. van Huis, A.R. Cudjoe, V. Johnson, O. Sakyi-Dawson, K. Ofori-Frimpong, P. Osei-Fosu, E. Tei-Quartey, W. Jonfia-Essien et al., *Int. J. Trop. Insect Sci.* **34** (2014) 58–71.

[2] J.C. Anikwe, A.A. Omoloye, P.O. Aikpokpodion, F.A. Okelana and A.B. Eskes, *Crop Protect.* **28** (2009) 350–355.

[3] J. Armathé Amougou, M. Tchindjang, H. Unusa and R.A. Batha, Soleil. *J. Cameroon Acad. Sci.* **11** (2013).

[4] R. Babin, *Contribution à l'amélioration de la lutte contre le miride du cacaoyer Sahlbergella singularis Hagl. (Hemiptera: Miridae).* *Influence des facteurs agro-écologiques sur la dynamique des populations du ravageur*. Ph.D. thesis, Université Paul Valéry-Montpellier III (2009).

[5] R. Babin, D.H. Bisseleua, L. Dibog and J.P. Lumaret, *J. Appl. Entomol.* **132** (2008) 366–374.

[6] R. Babin, M. Gerben ten Hoopen, C. Cilas, F. Enjalric, Yede, P. Gendre and J.P. Lumaret, *Agric. Forest Entomol.* **12** (2010) 69–79.
[7] R. Babin, J.C. Anikwé, L. Dibog and J.P. Lumaret, *Entomol. Exp. Appl.* **141** (2011) 25–34.

[8] A. Berman and R.J. Plemmons, Nonnegative matrices in the mathematical sciences. *SIAM* (1994).

[9] J. Carr, Heidelberg, Berlin, 1981.

[10] C. Castillo-Chavez and B. Song, *Math. Biosci. Eng.* **1** (2004) 361.

[11] N. Chitnis, J.M. Cushing and J. Hyman, *SIAM J. Appl. Math.* **67** (2006) 24–45.

[12] M. Djoukwe Tapi, L. Bagny Beilhe, S. Bowong and Y. Dumont, *Math. Methods Appl. Sci.* **41** (2018) 8673–8696.

[13] M. Djoukwe Tapi, L. Bagny-Beilhe and Y. Dumont, *Nonlinear Anal. Real World Appl.* **54** (2020) 103082.

[14] P.F. Entwiste, Insects and cocoa, in *Cocoa, edited by G.A.R. Wood and R.A. Lass. Longman, London* (1985) 366–443.

[15] P.F. Entwistle *et al.*, *Pests of cocoa*, 1972.

[16] H. Fleming Wendell and W. Rishel Raymond, *Deterministic and Stochastic Optimal Control, Vol. 1*. Springer Science & Business Media (2012).

[17] B. Fule Chi, 2014.

[18] P. Jagoret, *Analyse et évaluation de systèmes agroforestiers complexes sur le long terme : application aux systèmes de culture à base de cacaoyer au Centre Cameroun*. Thèse soutenue pour l'obtention du titre de Docteur de Montpellier Supagro (2011).

[19] R. Kumar and A.K. Ansari, *Zool. J. Linnean Soc.* **54** (1974) 1–29.

[20] S. Lenhart and J.T. Workman, *Optimal Control Applied to Biological Models*. Chapman and Hall/CRC (2007).

[21] S. Marino, B. Hogue Ian, J. Ray Christian and D.E. Kirschner, *J. Theor. Biol.* **254** (2008) 178–196.

[22] N. Nkelle, The sustainability of Cameron cocoa economy, 2016.

[23] P. Petithuguenin, *Plantations Recherche Dev.* **5** (1998) 393–411.

[24] S. Pontryagin Lev, *Mathematical Theory of Optimal Processes*. CRC Press (1987).

[25] D.J. Sonwa, O. Coulibaly, S.F. Weise, A. Adesina Akinwumi and J.J. Janssens Marc, *Crop Protect.* **27** (2008) 1159–1164.

[26] E. Tamanjong and R. Neena, Insect threat to Cameroon cocoa crop. Sap-sucking bug infests African nation as demand for chocolate ingredient is set to outstrip supply, 2016.

[27] P. Takam Soh, E.P. Nguéma Ndong, H. Gwet and M. Ndoumbè-Nkeng, *J. Roy. Stat. Soc. C (Appl. Stat.)* **62** (2013) 741–760.

[28] N. Tcharbuahbokengo, Cocoa production in Cameroon, in *AFTA 2005 conference proceedings*, 2005.

[29] G.M. ten Hoopen, P. Deberdt, M. Mbenoun and C. Cilas, *Ann. Appl. Biol.* **160** (2012) 260–272.

[30] P. Van den Driessche and J. Watmough, *Math. Biosci.* **180** (2002) 29–48.

[31] G. Williams, *Bul. Entomol. Res.* **44** (1953) 101–119.

[32] R. Yede, Babin, C. Djieto-Lordon, C. Cilas, L. Dibog, R. Mahob and C.F. Bilong Bilong, *J. Econ. Entomol.* **105** (2012) 1285–1292.

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the *Subscribe to Open model* (S2O). We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting [subscribers@edpsciences.org](mailto:subscribers@edpsciences.org).

More information, including a list of supporters and financial transparency reports, is available at [https://edpsciences.org/en/subscribe-to-open-s2o](https://edpsciences.org/en/subscribe-to-open-s2o).