# Self-Financing Model for Cabbage Crops with Pest Management

Aurelien Kambeu Youmbi, Suzanne Touzeau, Frédéric Grognard, Berge Tsanou

► To cite this version:

Aurelien Kambeu Youmbi, Suzanne Touzeau, Frédéric Grognard, Berge Tsanou. Self-Financing Model for Cabbage Crops with Pest Management. 2024. hal-04589904

HAL Id: hal-04589904

https://hal.inrae.fr/hal-04589904

Preprint submitted on 27 May 2024

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
# Self-Financing Model for Cabbage Crops with Pest Management

Aurelien Kambeu Youmbi, Suzanne Touzeau, Frédéric Grognard, Berge Tsanou

► To cite this version:

Aurelien Kambeu Youmbi, Suzanne Touzeau, Frédéric Grognard, Berge Tsanou. Self-Financing Model for Cabbage Crops with Pest Management. Mathematical Biosciences, 2024, 378, pp.109332. 10.2139/ssrn.4793810 . hal-04589904v2

HAL Id: hal-04589904

https://hal.inrae.fr/hal-04589904v2

Submitted on 8 Nov 2024

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
# Self-Financing Model for Cabbage Crops
with Pest Management

Aurelien KAMBEU YOUMBI<fnref n="1" />,<fnref n="2" />,*, Suzanne TOUZEAU<fnref n="2" />,<fnref n="3" />,
Frédéric GROGNARD<fnref n="3" />, Berge TSANOU<fnref n="1" />,<fnref n="4" />

March 2024

## Abstract

Smallholder farmers rely on their farm earnings to cover operating costs and generate income. That is not an easy task because of the pests, which reduce yields and generate plant protection costs. The farm yield and plant protection depend on the budget capacity of the farmer. In this work, we want to explore conditions for a sustainable and self-financing cabbage farm. We propose then a non-linear mathematical model for cabbage crops by considering the current account of the plantation as a dynamic variable. We assume that this variable increases due to the sale of cabbages, and provides for the seedling purchase, the plant protection costs, and the grower's income. In the first part, we analyze the model without pest management. We determine how the budget must be spent and we show the existence of a double transcritical bifurcation. We quantify the seasonal yield and income, and estimate the damage due to pest herbivory. In the second part, we analyze a slightly simplified version of our model and obtain the existence of a backward bifurcation. Furthermore, we show that botanical pesticides can be used to prevent and eradicate the pest spread with relatively low plant protection costs.

**Keywords:** Dynamic system, Reproductive number, Self-financing smallholder, Vegetable crop, Botanical pesticide, *Plutella xylostella*

## Highlights:

* Bio-economic modeling is more adapted for understanding agricultural issues.

* Cabbage is a robustly sustainable and self-financing crop.

* Using botanical pesticides is an affordable pest management technique.

# 1 Introduction

Agriculture occupies a crucial place in the development and research programs of many countries. The added value generated by agriculture, forestry, and fishing increased by 73% between 2000 and 2019 with a peak of around USD 3.5 trillion in 2018 [1]. In Africa, this increase is around USD 234 billion. In 2020, primary crops production was around 9.3 billion tonnes. Vegetables represent 12% of primary crops production, that is 1,128 million tonnes [2], cabbage being the fourth most widely cultivated vegetable, with a 6% share of the overall vegetable production [2]. Cabbage (*Brassica oleracea varieta capitata*) is used for human consumption [3, 4] and features in several diet programs [5, 6]. Furthermore, cabbage crops constitute a substantial source of income for African smallholder farmers [7, 8].

FAO estimates that annually up to 40% of the global crop production is lost due to pests, which results in costs over USD 220 billion, including at least USD 70 billion due to invasive insects. As for the diamondback moth (*Plutella xylostella*), it attacks cruciferous crops [9]; it is the most destructive pest of cabbage crops worldwide [10, 11, 12, 13]. The insect can be found in Senegal [10], Benin [14], Ghana [15], South Africa [16], Togo [17], Brazil [18], India [12], Canada [19]. Its development cycle in tropical regions such as India allows for up to 16 generations in a single year [12]. The most significant damage is observed during the third larval stage. The larvae feed on the soft parts of cabbage leaves [10, 9] and can cause the death of the host in the case of heavy infestation [20]. Controlling this pest is hence a major issue.

<footnote n="1">¹ Department of Mathematics and Computer Science, URMA, University of Dschang, P.O. Box 67 Dschang, Cameroon</footnote>

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

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

<footnote n="4">⁴ IRD, UMMISCO, Sorbonne University, Bondy, France</footnote>

<footnote>* Corresponding author: :aurelienkambeu@gmail.com</footnote>
# Highlights

## Self-Financing Model for Cabbage Crops with Pest Management

Aurelien Kambeu Youmbi, Suzanne Touzeau, Frédéric Grognard, Berge Tsanou

* Bio-economic modeling is more adapted for understanding agricultural issues.

* Cabbage is a robustly sustainable and self-financing crop.

* Using botanical pesticides is an affordable pest management technique.
Pest management based on synthetic pesticides has its limits, as the pest developed resistance to most of the active substances in commonly used pesticides [21, 9, 22]. In addition, an excessive use of synthetic pesticides is harmful to the environment and causes significant health problems for human [17, 23]. In this context, botanical pesticides based on crude extracts from locally-growing plants are a viable alternative [24, 25, 12, 15]. The use of these pesticides is safe for humans and their environment [17, 26, 27]. In addition to their biocide properties, botanical pesticides have repellent and antifeedant effects [28, 15, 27], and attract natural enemies of plutella xylostella [29]. Furthermore, botanical pesticides can be produced locally using simple methods. This has the direct effect of reducing treatment costs and therefore, substantially increasing profits for smallholder farmers [30, 31]. In general, smallholder farms must be self-financing; that is, the income generated by the farm must cover all expenses incurred by the farm, as well as the owner's living expenses.

Bio-economic models have been extensively studied (see [32, 33] and references therein). In these models, the net economic revenue is given by the Gordon's algebraic equation [32] and resources growth following natural laws [34]. This classic approach cannot adequately address our self-financing concern. In fact, the farm yield and plant protection depend on the budget capacity of the farmer. A model that incorporates this budgetary constraint as a dynamic variable is therefore more appropriate [35, 36].

Thus, we developed a bio-economic model with a threefold purpose: to quantify the damage caused by the attack of plutella larvae in a cabbage field; to assess the effectiveness of using botanical pesticides in controlling pest; and to evaluate the sustainability and the profitability of such systems. While taking into account the specific features of the interaction between cabbage and its pest, our model presents how a smallholder farmer can distribute the income from his farm to not only sustain his farming activity but also make a living from it, with or without pest management based on botanical pesticides.

This paper is organized as follows. Section 2 describes the self-financing model formulation and establishes the positivity and boundedness properties of its solutions. In Section 3, we study the uncontrolled model. Different equilibria and their stability conditions are obtained together with some numerical simulations. In Section 4, we study a slightly simplified controlled model, which has the same dynamical properties as the complete model presented in Section 2. Pest management is based on the use of botanical pesticides, both as treatment and prophylaxis. In this setting, the model exhibits a backward phenomenon. Theoretical results are then illustrated by numerical experiments. Finally, the conclusion is drawn in Section 5.

# 2 Mathematical model and basic properties

In this section, we build the mathematical model and establish some of its basic properties.

## 2.1 Model formulation

As young cabbages are more susceptible than adult ones [10], the biomass in our model is split into a young susceptible biomass and a resistant adult biomass, denoted by $B_y$ and $B_a$ respectively. We assume a logistic growth for each biomass with growth rates $r_{0y}$ and $r_{0a}$, and an aging rate $\gamma$, where $1/\gamma$ represents the average susceptibility duration. We assume that the cultivated soil land contains enough nutrients to allow normal growth of plants. This situation is observable with volcanic soils which are very fertile, as its is the case of in the West region in Cameroon [37, 38]. Adult biomass is harvested at rate $h$. We assume that, without investment in replanting, the plantation is inevitably doomed to disappear. In mathematical language, this natural assumption is described by the following conditions:

$$ \alpha_y := \gamma - r_{0y} > 0 \text{ and } \alpha_a := h - r_{0a} > 0. \qquad (1) $$

In order to evaluate the plantation sustainability, we introduce the plantation current account as a dynamic variable $M$. This account is credited as a result of sales of mature biomass harvested and sold at unit price $q$. The money is withdrawn from this account at rate $k$ to invest to buy new plants, to protect the plantation from pests and finally to cover the farmer's living expenses. The young biomass plantation rate is represented by the bounded function $k_1 M/(1 + \delta_M M)$, where $k_1$ and $\delta_M$ represent the investment rate in young biomass and the fear cost of biomass competition, respectively. We chose a bounded function for the plantation rate as a large account balance does not necessarily directly translate into larger fields and manpower. The plant protection cost is $\alpha M$ where $\alpha$ represents the investment rate for pesticide acquisition such that:

$$ k > k_1 + \alpha. \qquad (2) $$

We denote the smallholder farmer's total income by a variable $T$, and its dynamics is modeled by the following equation:

$$ \dot{T} = kM - \frac{k_1 M}{1 + \delta_M M} - \alpha M. $$
# Self-Financing Model for Cabbage Crops with Pest Management

Aurelien Kambeu Youmbi<sup>a,b,*</sup>, Suzanne Touzeau<sup>b,c</sup>, Frédéric Grognard<sup>c</sup> and Berge Tsanou<sup>a,d,e</sup>

## ARTICLE INFO

**Keywords:**

Dynamic system

Reproductive number

Self-financing smallholder

Vegetable crop

Botanical pesticide

*Plutella xylostella*

## ABSTRACT

Smallholder farmers rely on their farm earnings to cover operating costs and generate income. That is not an easy task because of the pests, which reduce yields and generate plant protection costs. The farm yield and plant protection depend on the budget capacity of the farmer. In this work, we want to explore conditions for a sustainable and self-financing cabbage farm. We propose then a non-linear mathematical model for cabbage crops by considering the current account of the plantation as a dynamic variable. We assume that this variable increases due to the sale of cabbages, and provides for the seedling purchase, the plant protection costs, and the grower's income. In the first part, we analyze the model without pest management. We determine how the budget must be spent and we show the existence of a double transcritical bifurcation. We quantify the seasonal yield and income, and estimate the damage due to pest herbivory. In the second part, we analyze a slightly simplified version of our model and obtain the existence of a backward bifurcation. Furthermore, we show that botanical pesticides can be used to prevent and eliminate the pest spread with relatively low plant protection costs.

## 1. Introduction

Agriculture occupies a crucial place in the development and research programs of many countries. The added value generated by agriculture, forestry, and fishing increased by 73% between 2000 and 2019 with a peak of around USD 3.5 trillion in 2018 [12]. In Africa, this increase is around USD 234 billion. In 2020, primary crop production was around 9.3 billion tonnes. Vegetables represent 12% of primary crops production, that is 1,128 million tonnes [13], cabbage being the fourth most widely cultivated vegetable, with a 6% share of the overall vegetable production [13]. Cabbage (*Brassica oleracea varieta capitata*) is used for human consumption [5, 6] and features in several diet programs [24, 45]. Furthermore, cabbage crops constitute a substantial source of income for African smallholder farmers [33, 47].

Food and Agriculture Organization (FAO) estimates that annually up to 40% of the global crop production is lost due to pests, which results in costs of over USD 220 billion, including at least USD 70 billion due to invasive insects. As for the diamondback moth (*Plutella xylostella*), it attacks cruciferous crops [32]; it is the most destructive pest of cabbage crops worldwide [44, 18, 15, 35]. The insect can be found in Senegal [44], Benin [39], Ghana [2], South Africa [26], Togo [1], Brazil [10], India [15], and Canada [28]. Its development cycle in tropical regions such as India allows for up to 16 generations in a single year [15]. The most

significant damage is observed during the third larval stage. The larvae feed on the soft parts of cabbage leaves [44, 32] and can cause the death of the host in the case of heavy infestation [49]. Controlling this pest is hence a major issue.

Pest management based on synthetic pesticides has its limits, as pests develop resistance to most of the active substances in commonly used pesticides [46, 32, 14]. In addition, an excessive use of synthetic pesticides is harmful to the environment and causes significant health problems for humans [1, 21]. In this context, botanical pesticides based on crude extracts from locally-growing plants are a viable alternative [36, 50, 15, 2]. The use of these pesticides is safe for humans and their environment [1, 23, 34]. In addition to their biocide properties, botanical pesticides have repellent and antifeedant effects [19, 2, 34], and attract natural enemies of *plutella xylostella* [9]. Furthermore, botanical pesticides can be produced locally using simple methods. This has the direct effect of reducing treatment costs and therefore, substantially increasing profits for smallholder farmers [3, 31]. In general, smallholder farms must be self-financing; that is, the income generated by the farm must cover all expenses incurred by the farm, as well as the owner's living expenses.

Bio-economic models have been extensively studied (see [16, 20] and references therein). In these models, the net economic revenue is given by the Gordon's algebraic equation [16] and resources growth following natural laws [43]. This classic approach cannot adequately address our self-financing concern. In fact, the farm yield and plant protection depend on the budget capacity of the farmer. A model that incorporates this budgetary constraint as a dynamic variable is therefore more appropriate [30, 38].

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

<footnote><sup>a</sup>Department of Mathematics and Computer Science, URMA, University of Dschang, P.O. Box 67 Dschang, Cameroon</footnote>

<footnote><sup>b</sup>Université Côte d'Azur, INRAE, ISA, France</footnote>

<footnote><sup>c</sup>Université Côte d'Azur, Inria, INRAE, CNRS, MACBES, France</footnote>

<footnote><sup>d</sup>IRD, UMMISCO, Sorbonne University, Bondy, France</footnote>

<footnote><sup>e</sup>Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria 0002, South Africa</footnote>

<footnote>*Corresponding author</footnote>

<footnote>a aurelienkambeu@gmail.com (A. Kambeu Youmbi); suzanne.touzeau@inrae.fr (S. Touzeau); frederic.grognard@inria.fr (F. Grognard); bergetsanou@gmail.com (B. Tsanou)
ORCID(s): 0009-0009-5725-8102 (A. Kambeu Youmbi); 0000-0001-5255-8960 (S. Touzeau); 0000-0003-0805-775X (F. Grognard); 0000-0002-2298-7972 (B. Tsanou)</footnote>

<footnote>1</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic diagram for system (3). Parameter values are given in Table 1. Young ($B_y$) and Adult ($B_a$) biomass growth at the rate $r_{0y}$ and $r_{0a}$. Young biomass age at the rate $\gamma$. Adult biomass harvested at the rate $h$ and sold (at the unit price $q$). The money from the sale is transferred to the current account ($M$). This money is withdrawn at the rate $k$ and used for young biomass acquisition (parameters $k_1$ and $\delta_M$), pesticide acquisition (parameter $\alpha$) and farmer daily income ($T$). Larvae ($L$) feed young biomass (parameter $\psi$). Botanical pesticide ($P$) is used for pest management with an antifeedant effect (parameter $b$). Pesticide is absorbed by larvae (parameter $\phi$) which leads to a lethal effect (parameter $\lambda$). All biomass, larvae and pesticide undergo mortality ($\mu_y$, $\mu_a$, $\mu_L$, $\phi_L$, $\phi_P$).</figcaption></figure>

We represent the pest compartment by the variable $L$. The larvae functional response is modeled by the Holling type I function $\psi B_y$, where $\psi$ is the consumption rate. The growth of larvae is assumed to be logistic-like with the modified growth rate $c\psi B_y L$ and the competition term $\mu_L L^2$ where $c$ and $\mu_L$ represent the conversion efficiency and the per capita competition death rate of larvae respectively. In addition, a linear mortality is taken in account with the rate $\phi_L$.

The most remarkable features of botanical pesticides are repellent and antifeedant effects. We modified the larvae functional response of cabbage attack by adding a $bP+1$ denominator to represent the antifeedant effect of botanical pesticides. Botanical pesticides purchased and sprayed in the plantation are also absorbed by larvae following a functional response $\phi_P$. This absorption is responsible for an additional larvae mortality rate $\lambda\phi_P$ where $\lambda$ represents the efficacy of pesticides. Furthermore, botanical pesticides decay naturally at a rate $\phi_P$.

According to the above considerations, we obtain the following model which describes our self-financing system:

$$ 
\begin{aligned}
\dot{B}_y &= r_{0y}B_y - \mu_y B_y^2 + \frac{r_1k_1M}{1+\delta_M M} - \frac{\psi L B_y}{bP+1} - \gamma B_y, \\
\dot{B}_a &= \gamma B_y + r_{0a}B_a - \mu_a B_a^2 - hB_a, \\
\dot{L} &= c \frac{\psi L B_y}{bP+1} - \lambda \phi_P L - \mu_L L^2 - \phi_L L, \\
\dot{P} &= a_P \alpha M - \phi_P L - \phi_P P, \\
\dot{M} &= qhB_a - kM.
\end{aligned}
\quad (3) $$

For model (3), the initial conditions are given by:

$$ B_y(0) \ge 0, B_a(0) \ge 0, L(0) \ge 0, P(0) \ge 0, \text{ and } M(0) \ge 0. \quad (4) $$

A complete description of model parameters is given in Table 1 and model schematic diagram in Fig. 1.

## 2.2 Positivity and boundedness of solutions

Local existence, uniqueness and non-negativity of the solutions of the Cauchy problem (3–4) are direct consequences of Cauchy-Lipschitz’s theorem and a simple analysis of the vector-field on the boundary of the positive orthant $\mathbb{R}_+^5$ [39]. The following theorem establishes the boundedness and global existence of the unique solution of system (3).
Thus, we developed a bio-economic model with a three-fold purpose: to quantify the damage caused by the attack of plutella larvae in a cabbage field; to assess the effectiveness of using botanical pesticides in controlling pests; and to evaluate the sustainability and the profitability of such systems. While taking into account the specific features of the interaction between cabbage and its pests, our model presents how a smallholder farmer can distribute the income from his farm to not only sustain his farming activity but also make a living from it, with or without pest management based on botanical pesticides.

This paper is organized as follows. Section 2 describes the self-financing model formulation and establishes the positivity and boundedness properties of its solutions. In Section 3, we study the uncontrolled model. Different equilibria and their stability conditions are obtained together with some numerical simulations. In Section 4, we study a slightly simplified controlled model, which has the same dynamical properties as the complete model presented in Section 2. Pest management is based on the use of botanical pesticides, both as treatment and prophylaxis. In this setting, the model exhibits a backward bifurcation. Theoretical results are then illustrated by numerical experiments. Finally, the conclusion is drawn in Section 5.

## 2. Mathematical model and basic properties

In this section, we build the mathematical model and establish some of its basic properties.

### 2.1. Model formulation

As young cabbages are more susceptible than adult ones [44], the biomass in our model (i.e. the total mass of the cabbage plants) is split into a young susceptible biomass and a resistant adult biomass, denoted by $B_y$ and $B_a$ respectively. We assume a logistic growth for each biomass with intrinsic growth rates $r_{0y}$ and $r_{0a}$, and an aging rate $\gamma$, where $1/\gamma$ represents the average susceptibility duration. We assume that the cultivated soil land contains enough nutrients to allow normal growth of plants. This situation is observable with volcanic soils which are very fertile, as is the case of the West region in Cameroon [22, 48]. Adult biomass is harvested at rate $h$. We assume that, without investment in replanting, the plantation will inevitably disappear. In mathematical language, this natural assumption is described by the following conditions:

$$ \alpha_y := \gamma - r_{0y} > 0 \text{ and } \alpha_a := h - r_{0a} > 0. \quad (1) $$

We assume that the smallholder farmer has sufficient fertile land and that his main constraint to plant new cabbages or spray pesticide is the availability of operating funds. In order to evaluate the plantation sustainability, we introduce the plantation current account as a dynamic variable $M$. This account is credited as a result of sales of mature biomass harvested proportionally at the rate $h$ and sold at unit price $q$. The money is withdrawn from this account at rate $k$ as investments to buy new plants, to protect the plantation from

pests and finally to cover the farmer's living expenses. The young biomass plantation rate is represented by the bounded function $k_1 M/(1 + \delta_M M)$, where $k_1$ and $\delta_M$ represent the investment rate in young biomass and the saturation coefficient for young biomass acquisition, respectively. We chose a bounded function for the plantation rate as a large account balance does not necessarily directly translate into larger fields and manpower. The plant protection cost is $\alpha M$ where $\alpha$ represents the investment rate for pesticide acquisition such that:

$$ k > k_1 + \alpha. \quad (2) $$

The farmer's cumulative income, represented by the variable $T$, is what remains of the money withdrawn by the smallholder farmer from the current account ($kM$) after investing in young biomass and pesticide. Its dynamic is then modeled as follows

$$ \dot{T} = kM - \frac{k_1 M}{1 + \delta_M M} - \alpha M. $$

Note that the expected link between pesticides sprayed and cabbage biomass in the field is made indirectly through $M$. Indeed, the farmer invests in both seedlings and treatment depending on the balance of the current account.

We represent the pest compartment by the variable $L$. The larvae functional response is modeled by the Holling type I function $\psi B_y$, where $\psi$ is the consumption rate. The growth of larvae is assumed to be logistic-like with the modified growth rate $c\psi B_y L$ and the competition term $\mu_L L^2$ where $c$ and $\mu_L$ represent the conversion efficiency and the per capita competition death rate of larvae respectively. In addition, a linear mortality is taken in account with the rate $\phi_L$.

We denote the volume of botanical pesticides used for plant protection by the variable $P$. Botanical pesticides are purchased at unit price $1/a_P$ and sprayed in the plantation. The most remarkable features of botanical pesticides are the repellent and antifeedant effects which make the plant less attractive or toxic to the pest. They hence decrease herbivory and increase larvae mortality. The biomass consumption rate $\psi L$ is multiplied by "1/(1 + $bP$)", a simple decreasing function of $P$ accounting for the pesticide antifeedant effect. Botanical pesticides are absorbed by larvae following a functional response $\phi P$. This absorption is responsible for an additional larvae mortality rate $\lambda\phi P$, where $\lambda$ represents the toxicity of pesticides. Furthermore, botanical pesticides decay naturally at a rate $\phi_P$.

According to the above considerations, we obtain the following model which describes our self-financing system:
<table><caption>Table 1: Description of parameters in system 3. The estimation procedure is explained in Appendix A</caption><thead><tr><th>Parameters</th><th>Descriptions</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>r<sub>0y</sub></td><td>Growth rate of young biomass</td><td>1.6 × 10<sup>-2</sup></td><td>day<sup>-1</sup></td></tr><tr><td>μ<sub>y</sub></td><td>Competition rate of young biomass</td><td>1.07 × 10<sup>-6</sup></td><td>kg<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>r<sub>1</sub></td><td>Inverse of the unit purchase price of young biomass</td><td>2.17</td><td>USD.kg<sup>-1</sup></td></tr><tr><td>k<sub>1</sub></td><td>Investment rate for young biomass acquisition</td><td>0.8</td><td>day<sup>-1</sup></td></tr><tr><td>δ<sub>M</sub></td><td>Fear cost of biomass competition</td><td>10<sup>-5</sup> or 1.17 × 10<sup>-2‡</sup></td><td>USD<sup>-1</sup></td></tr><tr><td>γ</td><td>Aging rate of young biomass</td><td>0.02</td><td>day<sup>-1</sup></td></tr><tr><td>r<sub>0a</sub></td><td>Growth rate of adult biomass</td><td>1.6 × 10<sup>-2</sup></td><td>kg.day<sup>-1</sup></td></tr><tr><td>μ<sub>a</sub></td><td>Competition rate of adult biomass</td><td>1.68 × 10<sup>-6</sup></td><td>kg<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>h</td><td>Harvest rate</td><td>0.03</td><td>day<sup>-1</sup></td></tr><tr><td>c</td><td>Conversion efficiency</td><td>0.6</td><td>larvae.kg<sup>-1</sup></td></tr><tr><td>ψ</td><td>Biomass consumption rate by larvae</td><td>4.46 × 10<sup>-5</sup></td><td>larvae<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>b</td><td>Pesticide antifeedant effect</td><td>0.01</td><td>l<sup>-1</sup></td></tr><tr><td>λ</td><td>Efficacy of botanical pesticide</td><td>0.04</td><td>larvae.l<sup>-1</sup></td></tr><tr><td>φ</td><td>Uptake rate of botanical pesticide by larvae</td><td>0.05</td><td>larvae<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>μ<sub>L</sub></td><td>Competition rate of larvae</td><td>2 × 10<sup>-6</sup></td><td>larvae<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>φ<sub>L</sub></td><td>Natural mortality rate of larvae</td><td>0.07</td><td>day<sup>-1</sup></td></tr><tr><td>a<sub>P</sub></td><td>Inverse of the pesticide unit price</td><td>2.16</td><td>l.USD<sup>-1</sup></td></tr><tr><td>α</td><td>Investment rate for pesticide acquisition</td><td>1.35 × 10<sup>-2</sup></td><td>day<sup>-1</sup></td></tr><tr><td>φ<sub>P</sub></td><td>Natural decay rate of botanical pesticide</td><td>0.1</td><td>day<sup>-1</sup></td></tr><tr><td>q</td><td>Unit selling price of adult biomass</td><td>0.69</td><td>USD.kg<sup>-1</sup></td></tr><tr><td>k</td><td>Plantation self-financing rate</td><td>0.9</td><td>day<sup>-1</sup></td></tr></tbody></table>

**Theorem 1.** The domain $\Omega$ defined by

$$ \Omega = \left\{ (B_y; B_a; L; P; M) \in \mathbb{R}_+^5 | \begin{array}{l} B_y \le K_y, B_a \le \frac{\gamma K_y}{\alpha_a}, L \le \frac{c\psi K_y}{\mu_L}, \\ M \le \frac{qh\gamma K_y}{k\alpha_a}, P \le \frac{a_p\alpha qh\gamma K_y}{k\alpha_a\phi_P} \end{array} \right\}, \quad \text{where} \quad K_y = \frac{r_1 k_1}{\delta_M \alpha_y} \tag{5} $$

is positively invariant for the flow of the system (3) and attractive.

Furthermore, the solution $(B_y(t); B_a(t); L(t); P(t); M(t))$ of system (3) with initial conditions satisfying assumption (4) is uniquely defined for all $t \ge 0$.

*Proof.* See Appendix B. $\square$

# 3 First case study: self-financing uncontrolled model

In this section, we assume that the smallholder farmer does not apply any pest management strategy. Under these assumptions, we set $\alpha = 0$ and $P = 0$ in system (3) and obtain the following self-financed uncontrolled model:

$$ \begin{aligned} \dot{B}_y &= r_{0y}B_y - \mu_y B_y^2 + \frac{r_1 k_1 M}{1 + \delta_M M} - \psi L B_y - \gamma B_y, \\ \dot{B}_a &= \gamma B_y + r_{0a} B_a - \mu_a B_a^2 - h B_a, \\ \dot{L} &= c\psi L B_y - \mu_L L^2 - \phi_L L, \\ \dot{M} &= qh B_a - kM. \end{aligned} \tag{6} $$

## 3.1 Equilibrium and stability analysis

### 3.1.1 Equilibria

Sub-model (6) potentially admits three equilibria:

* the extinction equilibrium $E_0 = (0; 0; 0; 0)$;

* the pest-free equilibrium (PFE) $E_1 = (B_{y1}; B_{a1}; 0; M_1)$ when $R_B > 1$;

* the endemic equilibrium $E^* = (B_y^*; B_a^*; L^*; M^*)$ when $R_B > 1$ and $R_0 > 1$;
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Schematic diagram for system (3). Parameter values are given in Table 1. Young ($B_y$) and Adult ($B_a$) biomass grow at the rate $r_{0y}$ and $r_{0a}$. Young biomass age at the rate $\gamma$. Adult biomass are harvested at the rate $h$ and sold (at the unit price $q$). The money from the sale is transferred to the current account ($M$). This money is withdrawn at the rate $k$ and used for young biomass acquisition (parameters $k_1$ and $\delta_M$), pesticide acquisition (parameter $\alpha$) and farmer daily income ($T$). Larvae ($L$) feed on young biomass (parameter $\psi$). Botanical pesticide ($P$) is used for pest management with an antifeedant effect (parameter $b$). Pesticide is absorbed by larvae (parameter $\phi$) which leads to a lethal effect (parameter $\lambda$). All biomass, larvae and pesticide undergo mortality ($\mu_y, \mu_a, \mu_L, \phi_L, \phi_P$).</figcaption></figure>

$$
\begin{align*}
\dot{B}_y &= r_{0y}B_y + \frac{r_1k_1M}{1+\delta_M M} - \mu_y B_y^2 - \frac{\psi L B_y}{bP+1} - \gamma B_y, \\
\dot{B}_a &= \gamma B_y + r_{0a}B_a - \mu_a B_a^2 - hB_a, \\
\dot{L} &= c \frac{\psi L B_y}{bP+1} - \lambda \phi P L - \mu_L L^2 - \phi_L L, \\
\dot{P} &= a_P \alpha M - \phi P L - \phi_P P, \\
\dot{M} &= qhB_a - kM.
\end{align*}
$$

For model (3), the initial conditions are given by:

$$B_y(0) \ge 0, B_a(0) \ge 0, L(0) \ge 0, P(0) \ge 0, \text{ and } M(0) \ge 0. \quad (4)$$

Note that all the parameters and state variables in model (3) are nonnegative. A complete description of model parameters is given in Table 1 and model schematic diagram in Fig. 1.

## 2.2. Positivity and boundedness of solutions

Local existence, uniqueness and non-negativity of the solutions of the Cauchy problem (3-4) are direct consequences of Cauchy-Lipschitz's theorem and a simple analysis of the vector-field on the boundary of the nonnegative orthant $\mathbb{R}_+^5$ [25]. The following theorem establishes the boundedness and global existence of the unique solution of system (3) subject to the natural assumption (1).

**Theorem 1.** The domain $\Omega$ defined by

$$\Omega = \left\{ (B_y; B_a; L; P; M) \in \mathbb{R}_+^5 | B_y \le K_y, B_a \le \frac{\gamma K_y}{\alpha_a}, L \le \frac{c\psi K_y}{\mu_L}, M \le \frac{qh\gamma K_y}{k\alpha_a}, P \le \frac{a_p\alpha qh\gamma K_y}{k\alpha_a\phi_P} \right\}, \text{ where } K_y = \frac{r_1k_1}{\delta_M\alpha_y},$$

(5)

is positively invariant for the flow of the system (3) and attractive. Furthermore, the solution $(B_y(t); B_a(t); L(t); P(t); M(t))$ of system (3) with initial conditions satisfying assumption (4) is uniquely defined for all $t \ge 0$.

PROOF. See Appendix B.

# 3. First case study: self-financing uncontrolled model

In this section, we assume that the smallholder farmer does not apply any pest management strategy. Under the later assumption, which amounts in setting $\alpha = 0$ and $P = 0$ in system (3), we obtain the following self-financed uncontrolled model supplemented by the natural assumption (1):

$$
\begin{align}
\dot{B}_y &= \frac{r_1 k_1 M}{1 + \delta_M M} - \mu_y B_y^2 - \psi L B_y - \alpha_y B_y, \nonumber \\
\dot{B}_a &= \gamma B_y - \mu_a B_a^2 - \alpha_a B_a, \nonumber \\
\dot{L} &= c \psi L B_y - \mu_L L^2 - \phi_L L, \nonumber \\
\dot{M} &= qh B_a - k M. \nonumber
\end{align}
$$

## 3.1. Equilibrium and stability analysis

### 3.1.1. Equilibria

Sub-model (6) potentially admits three admissible non-negative equilibria:

* the extinction equilibrium $E_0 = (0; 0; 0; 0)$;
where

$$ \mathcal{R}_B = \frac{k_1 \gamma r_1 q h}{k \alpha_y \alpha_a} \quad \text{and} \quad \mathcal{R}_0 = \frac{c \psi (\alpha_a + \mu_a B_{a_1}) B_{a_1}}{\gamma \phi_L}. \qquad (7) $$

The biomass reproductive number $\mathcal{R}_B$ indicates if cabbage production is financially sustainable for the farm given a cabbage price and the investment policy. The basic reproductive number $\mathcal{R}_0$ gives the number of larvae produced by one larvae during its lifetime when introduced in a healthy plantation. $\mathcal{R}_0$ can be obtained directly by applying the standard next-generation matrix approach [40].

The existence of $E_0$ is trivial. Thereafter, we show the feasibility of equilibria $E_1$ and $E^*$.

**Pest-free equilibrium** We obtain the coordinates of $E_1$ by solving the following algebraic equations:

$$ r_1 \frac{k_1 M}{1 + \delta_M M} - \mu_y B_y^2 - \alpha_y B_y = 0, \qquad (8) $$

$$ \gamma B_y - \mu_a B_a^2 - \alpha_a B_a = 0, \qquad (9) $$

$$ qhB_a - kM = 0. \qquad (10) $$

Using (9) and (10) in (8) and the fact that $B_a \neq 0$, we obtain the following equation in $B_a$:

$$ G_1(B_a) = \frac{r_1 k_1 q h \gamma}{k + \delta_M q h B_a} - \left( \alpha_y (\alpha_a + \mu_a B_a) + \frac{\mu_y}{\gamma} (\alpha_a + \mu_a B_a)^2 B_a \right) = 0. \qquad (11) $$

From Equation (11), we have:

$$ (i) \quad G_1(B_a)' = - \frac{r_1 k_1 \delta_M q^2 h^2 \gamma}{(k + \delta_M q h B_a)^2} - \alpha_y \mu_a - \frac{\mu_y}{\gamma} \left( (\alpha_a + \mu_a B_a)^2 + 2\mu_a (\alpha_a + \mu_a B_a) Ba \right) < 0 \quad \text{for all } B_a \ge 0; $$

$$ (ii) \quad G_1(0) = \alpha_y \alpha_a (\mathcal{R}_B - 1) > 0, \text{ if and only if } \mathcal{R}_B > 1; $$

$$ (iii) \quad G_1 \left( \frac{\gamma K_y}{\alpha_a} \right) < \frac{r_1 k_1 \alpha_a}{\delta_M K_y} - \alpha_a \alpha_y = 0. $$

Thus, equation $G_1(B_a) = 0$ has a unique positive solution denoted by $B_{a_1}$ in the interval $(0; \frac{\gamma K_y}{\alpha_a})$ if and only if $\mathcal{R}_B > 1$. The remaining components of $E_1$ are obtained as follows: $B_{y_1} = \frac{(\alpha_a + \mu_a B_{a_1}) B_{a_1}}{\gamma}$ and $M_1 = \frac{qhB_{a_1}}{k}$.

**Endemic equilibrium** The components of $E^*$ are obtained by solving the following algebraic equations:

$$ r_1 \frac{k_1 M}{1 + \delta_M M} - \mu_y B_y^2 - \psi L B_y - \alpha_y B_y = 0, \qquad (12) $$

$$ \gamma B_y - \mu_a B_a^2 - \alpha_a B_a = 0, \qquad (13) $$

$$ c \psi B_y - \mu_L L - \phi_L = 0, \qquad (14) $$

$$ qhB_a - kM = 0. \qquad (15) $$

Plugging (13) in (12) we obtain for all $B_a \in (0; B_{a_1})$:

$$ L(B_a) = \frac{G_1(B_a)}{\psi(\alpha_a + \mu_a B_a)}, \qquad (16) $$

where $G_1$ is given by (11).

Similarly, plugging (13) and (15) in (14), we obtain for $B_a \ge \beta$, where $\beta$ denotes the positive root of the right hand-side of equation (17):

$$ L(B_a) = \frac{c \psi (\alpha_a + \mu_a B_a) B_a - \gamma \phi_L}{\mu_L \gamma}. \qquad (17) $$

Equating (16) and (17), we have $G_1(B_a) = G_2(B_a)$ where

$$ G_2(B_a) = \psi(\alpha_a + \mu_a B_a) \left[ \frac{c \psi (\alpha_a + \mu_a B_a) B_a - \gamma \phi_L}{\mu_L \gamma} \right]. $$

It is clear that $G_2(B_a)' > 0$ for all $B_a \ge \beta$, as it ensures that both factors are positive and increasing.

Since $G_1(B_a) > 0$ if and only if $B_a < B_{a_1}$ and $G_2(B_a) > 0$ if and only if $B_a > \beta$ it comes that the two curves $G_1$ and $G_2$ intersect with $L(B_a) > 0$ if and only if $\beta < B_{a_1}$. This intersection is unique as $G'_1(\cdot) < 0$ and $G'_2(\cdot) > 0$ in the $[\beta, B_{a_1}]$ interval. Using Equation (7), the condition $\beta < B_{a_1}$ is equivalent to $\mathcal{R}_0 > 1$. Let $B_a^*$ be the interception point of $G_1$ and $G_2$. Substituting $B_a$ by $B_a^*$ in (13), (15), and (16), we obtain positive values of $B_y^*, M^*$ and $L^*$, respectively. Thus, sub-model (6) has a unique endemic equilibrium $E^*$ if and only if $\mathcal{R}_B > 1$ and $\mathcal{R}_0 > 1$.
<table><caption>**Table 1**

Description of parameters in system (3). The estimation procedure is explained in Appendix A.</caption><thead><tr><th>Parameters</th><th>Descriptions</th><th>Values</th><th>Units</th><th>Sources</th></tr></thead><tbody><tr><td>r<sub>0y</sub></td><td>Growth rate of young biomass</td><td>1.6 &times; 10<sup>-2</sup></td><td>day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>&mu;<sub>y</sub></td><td>Competition rate of young biomass</td><td>1.07 &times; 10<sup>-6</sup></td><td>kg<sup>-1</sup>&middot;day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>r<sub>1</sub></td><td>Inverse of the unit purchase price of young biomass</td><td>0.46</td><td>USD.kg<sup>-1</sup></td><td>[41]</td></tr><tr><td>k<sub>1</sub></td><td>Investment rate for young biomass acquisition</td><td>0.85</td><td>day<sup>-1</sup></td><td>-</td></tr><tr><td>&delta;<sub>M</sub></td><td>Saturation coefficient for young biomass acquisition</td><td>10<sup>-5</sup> or 1.17 &times; 10<sup>-2</sup><sup>‡</sup></td><td>USD<sup>-1</sup></td><td>-</td></tr><tr><td>&gamma;</td><td>Aging rate of young biomass</td><td>0.02</td><td>day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>r<sub>0a</sub></td><td>Growth rate of adult biomass</td><td>1.6 &times; 10<sup>-2</sup></td><td>day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>&mu;<sub>a</sub></td><td>Competition rate of adult biomass</td><td>1.68 &times; 10<sup>-6</sup></td><td>kg<sup>-1</sup>&middot;day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>h</td><td>Harvest rate</td><td>0.03</td><td>day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>c</td><td>Conversion efficiency</td><td>0.6</td><td>larvae.kg<sup>-1</sup></td><td>[29]</td></tr><tr><td>&psi;</td><td>Biomass consumption rate by larvae</td><td>4.46 &times; 10<sup>-5</sup></td><td>larvae<sup>-1</sup>&middot;day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>b</td><td>Pesticide antifeedant effect</td><td>0.01</td><td>liter<sup>-1</sup></td><td>-</td></tr><tr><td>&lambda;</td><td>Toxicity of botanical pesticide</td><td>0.04</td><td>larvae.liter<sup>-1</sup></td><td>-</td></tr><tr><td>&phi;</td><td>Uptake rate of botanical pesticide by larvae</td><td>0.05</td><td>larvae<sup>-1</sup>&middot;day<sup>-1</sup></td><td>[29]</td></tr><tr><td>&mu;<sub>L</sub></td><td>Competition rate of larvae</td><td>2 &times; 10<sup>-6</sup></td><td>larvae<sup>-1</sup>&middot;day<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>&phi;<sub>L</sub></td><td>Natural mortality rate of larvae</td><td>0.07</td><td>day<sup>-1</sup></td><td>[15, 18]</td></tr><tr><td>a<sub>P</sub></td><td>Inverse of the pesticide unit price</td><td>2.16</td><td>liters.USD<sup>-1</sup></td><td>Appendix A</td></tr><tr><td>&alpha;</td><td>Investment rate for pesticide acquisition</td><td>1.35 &times; 10<sup>-2</sup></td><td>day<sup>-1</sup></td><td>-</td></tr><tr><td>&phi;<sub>P</sub></td><td>Natural decay rate of botanical pesticide</td><td>0.1</td><td>day<sup>-1</sup></td><td>[3]</td></tr><tr><td>q</td><td>Unit selling price of adult biomass</td><td>0.69</td><td>USD.kg<sup>-1</sup></td><td>[3]</td></tr><tr><td>k</td><td>Plantation self-financing rate</td><td>0.9</td><td>day<sup>-1</sup></td><td>-</td></tr></tbody></table>

‡We consider $\delta_M = 10^{-5}$ in Section 3 and $\delta_M = 1.17 \times 10^{-2}$ in Section 4.

* the pest-free equilibrium (PFE) $E_1 = (B_{y_1}; B_{a_1}; 0; M_1)$ when $\mathcal{R}_B > 1$;

* the endemic equilibrium $E^* = (B_y^*; B_a^*; L^*; M^*)$ when $\mathcal{R}_B > 1$ and $\mathcal{R}_0 > 1$;

where

$$ \mathcal{R}_B = \frac{k_1 \gamma r_1 q h}{k \alpha_y \alpha_a} \quad \text{and} \quad \mathcal{R}_0 = \frac{c \psi (\alpha_a + \mu_a B_{a_1}) B_{a_1}}{\gamma \phi_L}. \quad (7) $$

The biomass reproductive number $\mathcal{R}_B$ indicates if cabbage production is financially sustainable for the farm given a cabbage price and the investment policy. The basic reproductive number $\mathcal{R}_0$ gives the number of larvae produced by one larvae during its lifetime when introduced in a healthy plantation. $\mathcal{R}_0$ can be obtained directly by applying the standard next-generation matrix approach [11].

The existence of $E_0$ is trivial. Thereafter, we show the feasibility of equilibria $E_1$ and $E^*$.

**Pest-free equilibrium** We obtain the components of $E_1$ by solving the following algebraic equations:

$$ r_1 \frac{k_1 M}{1 + \delta_M M} - \alpha_y B_y - \mu_y B_y^2 = 0, \quad (8) $$

$$ \gamma B_y - \mu_a B_a^2 - \alpha_a B_a = 0, \quad (9) $$

$$ qhB_a - kM = 0. \quad (10) $$

Using (9) and (10) in (8) and the fact that $B_a \neq 0$ (and multiplying by $1/\gamma$), we obtain the following equation in $B_a$:

$$ G_1(B_a) = \frac{r_1 k_1 q h \gamma}{k + \delta_M q h B_a} - (\alpha_y (\alpha_a + \mu_a B_a) + \frac{\mu_y}{\gamma} (\alpha_a + \mu_a B_a)^2 B_a) = 0. \quad (11) $$

From Equation (11), we have:

$$ \begin{array}{l} (i) \quad G'_1(B_a) = -\frac{r_1 k_1 \delta_M q^2 h^2 \gamma}{(k + \delta_M q h B_a)^2} - \alpha_y \mu_a - \\ \qquad \frac{\mu_y}{\gamma} \left( (\alpha_a + \mu_a B_a)^2 + 2\mu_a (\alpha_a + \mu_a B_a) Ba \right) < 0 \text{ for all} \\ \qquad B_a \ge 0; \end{array} $$

$$ (ii) \quad G_1(0) = \alpha_y \alpha_a (\mathcal{R}_B - 1) > 0, \text{ if and only if } \mathcal{R}_B > 1; $$

$$ (iii) \quad G_1 \left( \frac{\gamma K_y}{\alpha_a} \right) < \frac{r_1 k_1 \alpha_a}{\delta_M K_y} - \alpha_a \alpha_y = 0. $$

Thus, equation $G_1(B_a) = 0$ has a unique positive solution denoted by $B_{a_1}$ in the interval $(0, \frac{\gamma K_y}{\alpha_a})$ if and only if $\mathcal{R}_B > 1$. The remaining components of $E_1$ are obtained as follows:

$$ B_{y_1} = \frac{(\alpha_a + \mu_a B_{a_1}) B_{a_1}}{\gamma} \quad \text{and} \quad M_1 = \frac{qh B_{a_1}}{k}. $$

**Endemic equilibrium** The components of $E^*$ are obtained by solving the following algebraic equations:

$$ r_1 \frac{k_1 M}{1 + \delta_M M} - \mu_y B_y^2 - \psi L B_y - \alpha_y B_y = 0, \quad (12) $$

$$ \gamma B_y - \mu_a B_a^2 - \alpha_a B_a = 0, \quad (13) $$

$$ c\psi B_y - \mu_L L - \phi_L = 0, \quad (14) $$

$$ qhB_a - kM = 0. \quad (15) $$
### 3.1.2 Stability

Regarding the local stability of the equilibria of sub-model (6), we have the following main result.

**Theorem 2.** • The extinction equilibrium $E_0$ is always feasible. It is locally asymptotically stable provided $\mathcal{R}_B < 1$ and unstable whenever $\mathcal{R}_B > 1$.

• The PFE $E_1$ is feasible if and only if $\mathcal{R}_B > 1$; it is then unique. It is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

• The endemic equilibrium $E^*$ is feasible if and only if $\mathcal{R}_B > 1$ and $\mathcal{R}_0 > 1$; it is then unique. Moreover, whenever this equilibrium exists, it is locally asymptotically stable.

*Proof.* See Appendix C. □

Moreover, we have a global stability result for the extinction equilibrium, stated in the following theorem.

**Theorem 3.** The extinction equilibrium $E_0$ is globally asymptotically stable in $\Omega$ if and only if $\mathcal{R}_B \le 1$.

*Proof.* We consider the following positive definite function:

$$V = k\alpha_a B_y + qhr_1k_1B_a + \frac{k\alpha_a}{c}L + r_1k_1\alpha_aM.$$

Differentiating the above function $V$ with respect to time, $t$, along the solutions of sub-model (6), after some algebraic simplifications, we obtain

$$\dot{V} = k\alpha_a\alpha_y(\mathcal{R}_B - 1)B_y - k\alpha_a\mu_yB_y^2 - qhr_1k_1\mu_aB_a^2 - \frac{k\alpha_a}{c}(\mu_L L^2 + \phi_L L) - \frac{r_1k_1\alpha_a\delta_M M^2}{1 + \delta_M M}.$$

Thus, $V$ is a strict Lyapunov function for equilibrium $E_0$ if $\mathcal{R}_B \le 1$. We conclude that $E_0$ is globally asymptotically stable in $\Omega$ if and only if $\mathcal{R}_B \le 1$. □

Let us now investigate the bifurcations of sub-model (6). The following result holds.

**Theorem 4.** Sub-model (6) presents a double forward transcritical bifurcation: the first one at $\mathcal{R}_B = 1$ and the second one at $\mathcal{R}_0 = 1$.

*Proof.* The forward transcritical bifurcation at $\mathcal{R}_B = 1$ follows from Theorem 2 and Theorem 3. Indeed, equilibrium $E_0$ changes from stable ($\mathcal{R}_B \le 1$) to unstable ($\mathcal{R}_B > 1$), while the pest-free equilibrium $E_1$ becomes positive from $E_0$ and asymptotically stable.

Similarly, when $\mathcal{R}_B > 1$, equilibrium $E_1$ changes from stable ($\mathcal{R}_0 < 1$) to unstable ($\mathcal{R}_0 > 1$), while the endemic equilibrium $E^*$ becomes positive from $E_1$ and asymptotically stable. □

The typical bifurcation diagrams near bifurcation points $\mathcal{R}_B = 1$ and $\mathcal{R}_0 = 1$ are illustrated in Fig. 2(a) and Fig. 2(b), respectively.

We establish the following main result.

**Theorem 5.** The domain

$$\Omega' = \{(B_y; B_a; L; M) \in \mathbb{R}_+^4 | B_y \le B_{y1}, B_a \le B_{a1}, M \le M_1\} \quad (18)$$

is positively invariant. Furthermore, the pest-free equilibrium $E_1 = (B_{y1}; B_{a1}; 0; M_1)$ is asymptotically stable in $\Omega'$ if and only if $\mathcal{R}_0 \le 1 < \mathcal{R}_B$.

*Proof.* See Appendix D □

## 3.2 Simulations

In this section, we perform some numerical simulations to substantiate the theoretical results. All parameter values are given in Table 1. With these values, the biomass reproductive number $\mathcal{R}_B = 3.21$ and the basic reproductive number $\mathcal{R}_0 = 1.64$.

The stability results obtained show that the asymptotic behavior of sub-model (6) is entirely defined by thresholds $\mathcal{R}_B$ and $\mathcal{R}_0$. The condition for a sustainable agroecosystem is $\mathcal{R}_B > 1$. Fig. 3(a) shows, for an array of purchasing prices of young plants from the nurseries, how the smallholder farmer should fix the unit selling price of its cabbages and what level of sustainability it ensures. Starting from the reference parameter values in Table 1, the selling price can drop or the purchase price can increase, but still the production generates a
Plugging (13) in (12) we obtain for all $B_a \in (0, B_{a_{c_1}})$:

$$L(B_a) = \frac{G_1(B_a)}{\psi(\alpha_a + \mu_a B_a)}, \quad (16)$$

where $G_1$ is given by (11).

Similarly, plugging (13) and (15) in (14), we obtain for $B_a \ge \beta$, where $\beta$ denotes the positive root of the right hand-side of equation (17):

$$L(B_a) = \frac{c\psi(\alpha_a + \mu_a B_a) B_a - \gamma\phi_L}{\mu_L \gamma}. \quad (17)$$

Equating (16) and (17), we have $G_1(B_a) = G_2(B_a)$ where

$$G_2(B_a) = \psi(\alpha_a + \mu_a B_a) \left[ \frac{c\psi(\alpha_a + \mu_a B_a) B_a - \gamma\phi_L}{\mu_L \gamma} \right].$$

It is not difficult to show that $G'_2(B_a) > 0$ for all $B_a \ge \beta$, as it ensures that both factors are positive and increasing.

Since $G_1(B_a) > 0$ if and only if $B_a < B_{a_1}$ and $G_2(B_a) > 0$ if and only if $B_a > \beta$ it comes that the two curves $G_1$ and $G_2$ intersect with $L(B_a) > 0$ if and only if $\beta < B_{a_1}$. This intersection is unique as $G'_1(\cdot) < 0$ and $G'_2(\cdot) > 0$ in the $[\beta, B_{a_1}]$ interval. Using Equation (7), the condition $\beta < B_{a_1}$ is equivalent to $\mathbb{R}_0 > 1$. Let $B_a^*$ be the intersection point of $G_1$ and $G_2$. Substituting $B_a$ by $B_a^*$ in (13), (15), and (16), we obtain positive values of $B_a^*$, $M^*$ and $L^*$, respectively. Thus, sub-model (6) has a unique endemic equilibrium $E^*$ if and only if $\mathbb{R}_B > 1$ and $\mathbb{R}_0 > 1$.

### 3.1.2. Stability

This section is devoted to the asymptotic stability properties of the admissible nonnegative equilibria of sub-model (6). We start with the local stability results, which recall the existence conditions of each equilibria and also establish the instability conditions. This leads to more general stability results of the equilibria $E_0$ and $E_1$. We end this section with a bifurcation result of our dynamical system (6).

Regarding the local stability of the equilibria of sub-model (6), we have the following main result.

**Theorem 2 (Local stability).**

* The extinction equilibrium $E_0$ is always feasible. It is locally asymptotically stable provided $\mathbb{R}_B < 1$ and unstable whenever $\mathbb{R}_B > 1$.

* The pest-free equilibrium $E_1$ is feasible if and only if $\mathbb{R}_B > 1$; it is then unique. It is locally asymptotically stable if $\mathbb{R}_0 < 1$ and unstable if $\mathbb{R}_0 > 1$.

* The endemic equilibrium $E^*$ is feasible if and only if $\mathbb{R}_B > 1$ and $\mathbb{R}_0 > 1$; it is then unique. Moreover, whenever this equilibrium exists, it is locally asymptotically stable.

**PROOF.** See Appendix C.

Moreover, we have a global stability result for the extinction equilibrium, stated in the following theorem.

**Theorem 3 (Global stability of $E_0$).** The extinction equilibrium $E_0$ is globally asymptotically stable in $\Omega$ if and only if $\mathbb{R}_B \le 1$.

**PROOF.** We consider the following positive definite function:

$$V = k\alpha_a B_y + qhr_1k_1B_a + \frac{k\alpha_a}{c}L + r_1k_1\alpha_a M.$$

Differentiating the above function $V$ with respect to time, $t$, along the solutions of sub-model (6), and after some algebraic simplifications, we obtain

$$\dot{V} = k\alpha_a\alpha_y (\mathbb{R}_B - 1) B_y - k\alpha_a\mu_y B_y^2 - qhr_1k_1\mu_a B_a^2 - \frac{k\alpha_a}{c} (\mu_L L^2 + \phi_L L) - \frac{r_1k_1\alpha_a\delta_M M^2}{1 + \delta_M M}.$$

Thus, $V$ is a strict Lyapunov function for equilibrium $E_0$ if $\mathbb{R}_B \le 1$. In addition, Theorem 3 ensures that $E_0$ is unstable when $\mathbb{R}_B > 1$. We conclude that $E_0$ is globally asymptotically stable in $\Omega$ if and only if $\mathbb{R}_B \le 1$.

We establish the following main result.

**Theorem 4 (Stability of $E_1$).** The domain

$$\Omega' = \{(B_y; B_a; L; M) \in \mathbb{R}_+^4 | B_y \le B_{y_1}, B_a \le B_{a_1}, M \le M_1\} \quad (18)$$

is positively invariant. Furthermore, the pest-free equilibrium $E_1 = (B_{y_1}; B_{a_1}; 0; M_1)$ is asymptotically stable in $\Omega'$ if and only if $\mathbb{R}_0 \le 1 < \mathbb{R}_B$.

**PROOF.** See Appendix D

Let us now investigate the bifurcations of sub-model (6). The following result holds.

**Theorem 5 (Bifurcation result).** Sub-model (6) presents a double forward transcritical bifurcation: the first one at $\mathbb{R}_B = 1$ and the second one at $\mathbb{R}_0 = 1$.

**PROOF.** The forward transcritical bifurcation at $\mathbb{R}_B = 1$ follows from Theorem 2 and Theorem 3. Indeed, equilibrium $E_0$ changes from stable ($\mathbb{R}_B \le 1$) to unstable ($\mathbb{R}_B > 1$), while the pest-free equilibrium $E_1$ becomes positive from $E_0$ and asymptotically stable.

Similarly, when $\mathbb{R}_B > 1$, equilibrium $E_1$ changes from stable ($\mathbb{R}_0 < 1$) to unstable ($\mathbb{R}_0 > 1$), while the endemic equilibrium $E^*$ becomes positive from $E_1$ and asymptotically stable.

The typical bifurcation diagrams near bifurcation points $\mathbb{R}_B = 1$ and $\mathbb{R}_0 = 1$ are illustrated in Fig. 2(a) and Fig. 2(b), respectively.

The stability results are summarized in Table 2.
<figure><img src="image_3.png" /><figcaption>Figure 2: Bifurcation diagrams for sub-model (6), according to the biomass reproductive number $\mathcal{R}_B$ and the basic reproductive number $\mathcal{R}_0$ defined in (7). (a) $\mathcal{R}_B$-bifurcation diagram, obtained by varying parameter $q$ and (b) $\mathcal{R}_0$-bifurcation diagram, obtained by varying parameter $\phi_L$, with $\mathcal{R}_B = 3.21 > 1$. All remaining parameter values are given in Table 1. The red, blue and orange lines represent the extinction $E_0$, pest-free $E_1$ and endemic $E^*$ equilibria, respectively. A solid line corresponds to a stable equilibrium, a dashed line to an unstable equilibrium.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Variation of the biomass reproductive number $\mathcal{R}_B$ as a function of: (a) the unit purchase price of young biomass $r_1^{-1}$ and the unit selling price of adult biomass $q$; (b) the plantation self-financing and living expense rate $k$ and the investment rate for young biomass $k_1$. The values of the remaining parameters are given in Table 1. The purple squares represent the parameter values Table 1. The red line represents $\mathcal{R}_B = 1$. The price ranges in sub-figure (a) are realistic values extracted from [30] for $q$ and [41] for $r_1^{-1}$. The white zone on sub-figure (b) represents the region where $k_1 > k$ which contradicts our modeling hypothesis (2).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2:** Bifurcation diagrams for sub-model (6), according to the biomass reproductive number $R_B$ and the basic reproductive number $R_0$ defined in (7). (a) $R_B$-bifurcation diagram, obtained by varying parameter $q$ and (b) $R_0$-bifurcation diagram, obtained by varying parameter $\phi_L$, with $R_B = 3.21 > 1$. All remaining parameter values are given in Table 1. The red, blue and orange lines represent the extinction $E_0$, pest-free $E_1$ and endemic $E^*$ equilibria, respectively. A solid line corresponds to a stable equilibrium, a dashed line corresponds to an unstable equilibrium.</figcaption></figure>

<caption>**Table 2**</caption>

Stability of the admissible equilibria of system (6). Note that if $R_B \ge 1$ then $R_0 \ge 0$.

<table><thead><tr><th>Equilibria</th><th>Existence</th><th>Asymptotic stability</th><th>Instability</th></tr></thead><tbody><tr><td>Extinction equilibrium E<sub>0</sub></td><td>Always exists</td><td>R<sub>B</sub> ≤ 1 (Theorem 3, global)</td><td>R<sub>B</sub> > 1 (Theorem 2)</td></tr><tr><td>Pest-free equilibrium (PFE) E<sub>1</sub></td><td>R<sub>B</sub> > 1 (Theorem 2)</td><td>R<sub>0</sub> ≤ 1 (Theorem 4, local)</td><td>R<sub>0</sub> > 1 (Theorem 2)</td></tr><tr><td>Endemic equilibrium E<sup>*</sup></td><td>R<sub>0</sub> > 1 (Theorem 2)</td><td>R<sub>0</sub> > 1 (Theorem 2, local)</td><td></td></tr></tbody></table>

Theorem 5 describes the double forward bifurcation which explains the stability changes of the equilibria.

## 3.2. Simulations

In this section, we perform some numerical simulations to substantiate the theoretical results. All parameter values are given in Table 1. With these values, the biomass reproductive number is $R_B = 3.21$ and the basic reproductive number is $R_0 = 1.64$.

The stability results obtained show that the asymptotic behavior of sub-model (6) is entirely defined by thresholds $R_B$ and $R_0$. The condition for a sustainable agroecosystem is $R_B > 1$. Fig. 3(a) shows, for an array of purchasing prices of young plants from the nurseries, how the smallholder farmer should fix the unit selling price of its cabbages and what level of sustainability it ensures. In fact, starting from

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

<figure><img src="image_7.png" /><figcaption>**Figure 3:** Variation of the biomass reproductive number $R_B$ as a function of: (a) the unit purchase price of young biomass $r_1^{-1}$ and the unit selling price of adult biomass $q$; (b) the plantation self-financing and living expense rate $k$ and the investment rate for young biomass $k_1$. The values of the remaining parameters are given in Table 1. The purple squares represent the points with the parameter values Table 1. The red line represents $R_B = 1$. The price ranges in sub-figure (a) are realistic values extracted from Amoabeng et al. [3] for $q$ and jardins de la Vallée de la Siagne [41] for $r_1^{-1}$. The white zone on sub-figure (b) represents the region where $k_1 > k$ which contradicts our modeling hypothesis (2).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Variation of the basic reproductive number $\mathcal{R}_0$ as a function of: (a) the mortality rate of larvae $\phi_L$ and the biomass consumption rate of larvae $\psi$; (b) the self-financing rate $k$ and the investment rate for young biomass $k_1$. The values of the remaining parameters are given in Table 1. The purple squares represent the parameter values in Table 1. The red line represents $\mathcal{R}_0 = 1$. The price ranges in sub-figure (b) are realistic values extracted from [30] for $q$ and [41] for $r_1^{-1}$. The white zone on sub-figure (b) represents the region where $k_1 > k$ which contradicts our modeling hypothesis (2); the gray area represents the region where $\mathcal{R}_B \le 1$ corresponding to the situation where the system converges globally towards the trivial equilibrium.</figcaption></figure>

positive result. To get below the sustainability threshold, both need to change drastically within the realistic price ranges. Cabbage is then a robustly sustainable crop.

In addition, Fig. 3(b) reveals that, with stable market prices (the parameters of cabbage and pest dynamics being set to their reference values in Table 1), the purchase of new seedlings must represent more than 21% (slope of the red line in Fig. 3(b)) of the total expense in order to satisfy the sustainability condition. Finally, in the absence of pests, cabbage is an easy crop to produce, with excellent yields per hectare for a very affordable investment.

Maintaining a healthy plantation without any pest management strategy is a difficult task. Fig. 4(a) shows that the basic reproduction number $\mathcal{R}_0$ is greater than 1 for most realistic values of the pest traits. This leads to the persistence of infestation. Nevertheless, it may be possible to control the pest by under-exploiting the farmland (low $k_1$ values to drop below the red line, but high enough to remain above the gray area) so that larvae disappear progressively, as shown in Fig. 4(b). Unfortunately, this solution is not feasible from an economic and food security point of view, as the yield and the revenues then drastically decrease.

Fig. 5 illustrates the stability results of sub-model (6), with the reference parameter values in Table 1. This situation corresponds to a farm without any pest management. As $\mathcal{R}_0 = 1.64 > 1$, the endemic equilibrium (orange star) is stable, while the pest-free equilibrium (red dot) is unstable. The convergence towards the endemic equilibrium occurs even with an initial condition close to the pest-free equilibrium (yellow curve).

Let us define the seasonal yield and income at the equilibrium by, respectively,

$$ \mathcal{Y} = h \times \overline{B_a} \times A_t \quad \text{and} \quad \mathcal{I} = q \times h \times \overline{B_a} \times A_t, \qquad (19) $$

where $A_t$ denotes the average duration of the cropping season and $\overline{B_a}$ is the level of adult biomass at the stable equilibrium. Considering that $A_t = 90$ days, the seasonal yield is around 7.6 tonnes for an income of 5,220 USD, which is consistent with the values per hectare reported by [30].

In Fig. 6, we present the dynamics of state variables $B_y$ (Fig. 6(a)) and $L$ (Fig. 6(b)). The convergence towards the endemic equilibrium is slow and requires at least 1500 days, with oscillations that are typical for a prey-predator system. In addition, the green curve describes a typical behavior of a prey-predator model: as the host biomass starts at a very low level, the pest population decreases, until the biomass reaches a satisfactory level and then the pest outbreak begins.

According to Fig. 4(a), an increase of the pest mortality rate $\phi_L$ implies a decrease of the basic reproductive number $\mathcal{R}_0$. By considering an increased value $\phi_L = 0.12$ compared to reference value in Table 1, $\mathcal{R}_0 = 0.96 < 1$, so the pest-free equilibrium becomes stable. This is illustrated in Fig. 7 with the same initial conditions as in Fig. 5. At the pest-free equilibrium, the seasonal yield is around 11.1 tonnes, corresponding to an income of 7,750 USD, which is consistent with [30] and [42].
<figure><img src="image_8.png" /><figcaption>Figure 4: Variation of the basic reproductive number $R_0$ as a function of: (a) the mortality rate of larvae $\phi_L$ and the biomass consumption rate of larvae $\psi$; (b) the self-financing rate $k$ and the investment rate for young biomass $k_1$. The values of the remaining parameters are given in Table 1. The purple squares the points with represent the parameter values in Table 1. The red line represents $R_0 = 1$. The price ranges in sub-figure (b) are realistic values extracted from Amoabeng et al. [3] for $q$ and jardins de la Vallée de la Siagne [41] for $r_1^{-1}$. The white zone on sub-figure (b) represents the region where $k_1 > k$ which contradicts our modeling hypothesis (2); the gray area represents the region where $R_B \le 1$ corresponding to the situation where the system converges globally towards the trivial equilibrium.</figcaption></figure>

the reference parameter values in Table 1, the selling price can drop or the purchase price can increase, but still the production generates a positive result. To get below the sustainability threshold, both prices need to change drastically within the realistic price ranges. Cabbage is then a robustly sustainable crop.

In addition, Fig. 3(b) reveals that, with stable market prices (the parameters of cabbage and pest dynamics being set to their reference values in Table 1), the purchase of new seedlings must represent more than 21% (slope of the red line in Fig. 3(b)) of the total expense in order to satisfy the sustainability condition. Finally, in the absence of pests, cabbage is an easy crop to produce, with excellent yields per hectare for a very affordable investment.

Maintaining a healthy plantation without any pest management strategy is a difficult task. Fig. 4(a) shows that the basic reproduction number $R_0$ is greater than 1 for most realistic values of the pest traits. This leads to the persistence of infestation. Nevertheless, it may be possible to control the pest by under-exploiting the farmland (low $k_1$ values to drop below the red line, but high enough to remain above the gray area) so that larvae disappear progressively, as shown in Fig. 4(b). Unfortunately, this solution is not feasible from an economic and food security point of view, as the yield and the revenues then drastically decrease.

Fig. 5 illustrates the stability results of sub-model (6), with the reference parameter values in Table 1. This situation corresponds to a farm without any pest management. As $R_0 = 1.64 > 1$, the endemic equilibrium (orange star) is stable, while the pest-free equilibrium (red dot) is unstable. The convergence towards the endemic equilibrium occurs even with an initial condition close to the pest-free equilibrium (yellow curve), suggesting the global asymptotic stability of

the endemic equilibrium in the domain $\Omega$ minus the pest-free equilibrium.

Let us define the seasonal yield and income at the equilibrium respectively by:

$$ \mathcal{Y} = h \overline{B_a} A_t \quad \text{and} \quad \mathcal{I} = q h \overline{B_a} A_t, \qquad (19) $$

where $A_t$ denotes the average duration of the cropping season and $\overline{B_a}$ is the level of adult biomass at the stable equilibrium. Considering that $A_t = 90$ days, the seasonal yield is around 7.6 tonnes for an income of 5,220 USD, which is consistent with the values per hectare reported by Amoabeng et al. [3].

In Fig. 6, we present the dynamics of state variables $B_y$ (Fig. 6(a)) and $L$ (Fig. 6(b)). The convergence towards the endemic equilibrium is slow and requires at least 1500 days, with damned oscillations that are typical for a prey-predator system. In addition, the green curve describes a typical behavior of a prey-predator model: as the host biomass starts at a very low level, the pest population decreases, until the biomass reaches a satisfactory level and then the pest outbreak begins.

According to Fig. 4(a), an increase of the pest mortality rate $\phi_L$ implies a decrease of the basic reproductive number $R_0$. By considering an increased value $\phi_L = 0.12$ compared to reference value in Table 1, $R_0 = 0.96 < 1$, so the pest-free equilibrium becomes stable. This is illustrated in Fig. 7 with the same initial conditions as in Fig. 5. At the pest-free equilibrium, the seasonal yield is around 11.1 tonnes, corresponding to an income of 7,750 USD, which is consistent with Amoabeng et al. [3] and Guvenc and Yildirim [17].
(a)

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

(b)

<figure><img src="image_7.png" /><figcaption>Figure 5: Phase portrait of sub-model (6) in sub-spaces (a) $B_y - B_a - M$ and (b) $B_y - B_a - L$. The orange star is the stable endemic equilibrium $E^*$ and the red dot is the unstable pest-free equilibrium $E_1$. Parameter values are given in Table 1. The initial conditions ($B_y(0)$; $B_a(0)$; $L(0)$; $M(0)$) are: (1000; 1000; 200; 300) in green, (4195.29; 4009.54; 0.88; 84.52) in yellow, (2000; 3000; 500; 200) in blue and (3000; 1000; 0.8; 0) in cyan.</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>Figure 6: Dynamics of state variables (a) $B_y$ and (b) $L$ of sub-model (6). Parameter values and initial conditions are the same as in Fig. 5.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 5: Phase portrait of sub-model (6) in sub-spaces (a) $B_y - B_a - M$ and (b) $B_y - B_a - L$. The orange star is the stable endemic equilibrium $E^*$ and the red dot is the unstable pest-free equilibrium $E_1$. Parameter values are given in Table 1. The initial conditions ($B_y(0)$; $B_a(0)$; $L(0)$; $M(0)$) are: (1000; 1000; 200; 300) in green, (4195.29; 4009.54; 0.88; 84.52) in yellow, (2000; 3000; 500; 200) in blue and (3000; 1000; 0.8; 0) in cyan.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6: Dynamics of state variables (a) $B_y$ and (b) $L$ of sub-model (6). Parameter values and initial conditions are the same as in Fig. 5.</figcaption></figure>

## 4. Second case study: self-financing controlled model

In this section, we perform a mathematical analysis of a slightly simplified version of model (3), in which we neglect the antifeedant action of the botanical pesticides ($b = 0$). Moreover, we assume that the farmer has enough experience to space the plants properly, such that one we can neglect the competition in the two biomass compartments ($\mu_y = \mu_a = 0$). The other values of parameters remain unchanged, as given in Table 1, except $\delta_M = 1.17 \times 10^{-2}$, which is adjusted to obtain reasonable estimations of the cabbage biomass.

Under the above considerations, we obtain from system (3) supplemented by the natural assumption (1) the following

slightly modified model:

$$ 
\begin{aligned}
\dot{B}_y &= \frac{r_1 k_1 M}{1 + \delta_M M} - \psi L B_y - \alpha_y B_y, \\
\dot{B}_a &= \gamma B_y - \alpha_a B_a, \\
\dot{L} &= c \psi L B_y - \lambda \phi P L - \mu_L L^2 - \phi_L L, \\
\dot{P} &= a_p \alpha M - \phi P L - \phi_p P, \\
\dot{M} &= q h B_a - k M.
\end{aligned}
\qquad (20) $$

### 4.1. Equilibrium and stability analysis

#### 4.1.1. Equilibria

We first determine the admissible nonnegative equilibria of sub-model (20). The extinction equilibrium $E_{c_0} =$
<figure><img src="image_10.png" /><figcaption>Figure 7: Phase portrait of sub-model (6) in sub-spaces (a) $B_y - B_a - M$ and (b) $B_y - B_a - L$. The red dot is the stable pest-free equilibrium $E_1$. Parameter values are given in Table 1 excepted $\phi_L = 0.12$ leading to $\mathcal{R}_0 = 0.96 < 1$. The initial conditions are the same as Fig. 5.</figcaption></figure>

# 4 Second case study: self-financing controlled model

In this section, we perform a mathematical analysis of a slightly simplified version of model (3), in which we neglect the antifeedant action of the botanical pesticides ($b = 0$). Moreover, we assume that the farmer has enough experience to space the plants properly, such that one we can neglect the competition in the two biomass compartments ($\mu_y = \mu_a = 0$). The other values of parameters remain unchanged, as given in Table 1, except $\delta_M = 1.17 \times 10^{-2}$, which is adjusted to obtain reasonable estimations of the cabbage biomass.

Under the above considerations, we obtain the following slightly modified model:

$$ 
\begin{aligned}
\dot{B}_y &= r_{0y}B_y + \frac{r_1k_1M}{1 + \delta_M M} - \psi L B_y - \gamma B_y, \\
\dot{B}_a &= \gamma B_y + r_{0a}B_a - hB_a, \\
\dot{L} &= c\psi L B_y - \lambda\phi P L - \mu_L L^2 - \phi_L L, \\
\dot{P} &= a_P \alpha M - \phi_P L - \phi_P P, \\
\dot{M} &= qhB_a - kM.
\end{aligned}
\quad (20) $$

## 4.1 Equilibrium and stability analysis

### 4.1.1 Equilibrium

We first determine the equilibria of sub-model (20). The extinction equilibrium $E_{c0} = (0; 0; 0; 0; 0)$ always exists. The unique pest-free equilibrium

$$ 
\begin{aligned}
E_{c1} &= \left( \frac{\alpha_a k (\mathcal{R}_B - 1)}{\gamma \delta_M q h}, \frac{k (\mathcal{R}_B - 1)}{\delta_M q h}, 0; \frac{a_P \alpha (\mathcal{R}_B - 1)}{\phi_P \delta_M}, \frac{\mathcal{R}_B - 1}{\delta_M} \right) \\
&= (B_{y_{c1}}; B_{a_{c1}}; 0; P_{c1}; M_{c1})
\end{aligned}
$$

exists whenever $\mathcal{R}_B > 1$, with $\mathcal{R}_B$ given by (7).

**Endemic equilibria** In order to study the existence of the endemic equilibrium, we recall that $\alpha_y = \gamma - r_{0y}$ and $\alpha_a = h - r_{0a}$, and introduce the new threshold

$$ \mathcal{R}_e = \frac{c\psi\alpha_a B_{a_{c1}}}{\gamma(\lambda\phi P_{c1} + \phi_L)} \quad (21) $$

representing the effective reproductive number.
(a)

<figure><img src="image_11.png" /><figcaption>**Figure 7:** Phase portrait of sub-model (6) in sub-spaces (a) $B_y - B_a - M$ and (b) $B_y - B_a - L$. The red dot is the stable pest-free equilibrium $E_1$. Parameter values are given in Table 1 excepted $\phi_L = 0.12$ leading to $R_0 = 0.96 < 1$. The initial conditions are the same as Fig. 5.</figcaption></figure>

(b)

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

(0; 0; 0; 0; 0) always exists. The unique pest-free equilibrium

$$E_{c_1} = \left( \frac{\alpha_a k (R_B - 1)}{\gamma \delta_M q h}; \frac{k (R_B - 1)}{\delta_M q h}; 0; \frac{a_p \alpha (R_B - 1)}{\phi_p \delta_M}; \frac{R_B - 1}{\delta_M} \right) = (B_{y_{c_1}}; B_{a_{c_1}}; 0; P_{c_1}; M_{c_1})$$

exists whenever $R_B > 1$, with $R_B$ given by (7).

**Endemic equilibria** In order to study the existence of the endemic equilibrium, we introduce the new threshold

$$R_e = \frac{c\psi\alpha_a B_{a_{c_1}}}{\gamma (\lambda\phi P_{c_1} + \phi_L)} \quad (21)$$

representing the effective reproductive number.

We want to solve the system:

$$r_1 \frac{k_1 M}{1 + \delta_M M} - \psi L B_y - \alpha_y B_y = 0, \quad (22)$$

$$\gamma B_y - \alpha_a B_a = 0, \quad (23)$$

$$c\psi B_y - \lambda\phi P - \mu_L L - \phi_L = 0, \quad (24)$$

$$a_p \alpha M - \phi_p P - \phi_p L = 0, \quad (25)$$

$$qh B_a - k M = 0. \quad (26)$$

Relations (23) and (26) give:

$$B_y = \frac{\alpha_a B_a}{\gamma} \quad \text{and} \quad M = \frac{qh B_a}{k}. \quad (27)$$

Introducing (27) in (22), we have:

$$L(B_a) = \frac{1}{\psi \alpha_a} \left( \frac{r_1 k_1 q h \gamma}{k + \delta_M q h B_a} - \alpha_a \alpha_y \right) = \frac{\alpha_y k \delta_M q h}{\psi (k + \delta_M q h B_a)} (B_{a_{c_1}} - B_a), \quad (28)$$

which is positive if and only if $B_a < B_{a_{c_1}}$.

Using both (28) and (26) in (25), we get:

$$P(B_a) = \frac{a_p \alpha q h B_a}{k (\phi_p + \phi \psi \alpha_a L(B_a))} = \frac{D_1 (k + \delta_M q h B_a) B_a}{D_2 + D_3 B_a}, \quad (29)$$

with

$$D_1 = \frac{a_p \alpha \psi \alpha_a q h}{k} > 0,$$

$$D_2 = k \phi \alpha_a \alpha_y (R_B - 1) + \phi_p \psi \alpha_a k > 0,$$

$$D_3 = \delta_M q h (\phi_p \psi \alpha_a - \alpha_a \alpha_y \phi).$$

The equivalence between the two definitions of $P$ implies that $D_2 + D_3 B_a > 0$ for all $B_a \in (0, B_{a_{c_1}})$.

Finally, using (27), (28) and (29) to replace $B_y$, $L$ and $P$ by functions of $B_a$ in (24), one has

$$G_3(B_a) = \frac{c\psi\alpha_a}{\gamma} B_a - \lambda\phi P(B_a) - \mu_L L(B_a) - \phi_L = 0 \quad (30)$$

that we will study on the biologically feasible interval $(0, B_{a_{c_1}})$ to determine the endemic equilibria.

The derivative of $L$ in (28) with respect to $B_a$ gives

$$L'(B_a) = - \frac{r_1 k_1 q^2 h^2 \gamma \delta_M}{\psi \alpha_a (k + \delta_M q h B_a)^2} < 0,$$

and its second derivative yields

$$L''(B_a) = \frac{r_1 k_1 q^3 h^3 \gamma \delta_M^2}{\psi \alpha_a (k + \delta_M q h B_a)^3} > 0. \quad (31)$$
We want to solve the system:

$$r_1 \frac{k_1 M}{1 + \delta_M M} - \psi L B_y - \alpha_y B_y = 0, \quad (22)$$

$$\gamma B_y - \alpha_a B_a = 0, \quad (23)$$

$$c \psi B_y - \lambda \phi P - \mu_L L - \phi_L = 0, \quad (24)$$

$$a_P \alpha M - \phi_P P - \phi_P L = 0, \quad (25)$$

$$q h B_a - k M = 0. \quad (26)$$

Relations (23) and (26) give:

$$B_y = \frac{\alpha_a B_a}{\gamma} \quad \text{and} \quad M = \frac{q h B_a}{k}. \quad (27)$$

Introducing (27) in (22), we have:

$$L(B_a) = \frac{1}{\psi \alpha_a} \left( \frac{r_1 k_1 q h \gamma}{k + \delta_M q h B_a} - \alpha_a \alpha_y \right), \quad (28)$$

which is positive if and only if $B_a < B_{a_1}$.

Using both (28) and (26) in (25), we get:

$$P(B_a) = \frac{a_p \alpha q h B_a}{k (\phi_p + \phi \psi \alpha_a L(B_a))} = \frac{D_1 (k + \delta_M q h B_a) B_a}{D_2 + D_3 B_a}, \quad (29)$$

with

$$D_1 = \frac{a_p \alpha \psi \alpha_a q h}{k} > 0,$$

$$D_2 = k \phi \alpha_a \alpha_y (\mathcal{R}_B - 1) + \phi_p \psi \alpha_a k > 0,$$

$$D_3 = \delta_M q h (\phi_P \psi \alpha_a - \alpha_a \alpha_y \phi).$$

The equivalence between the two definitions of $P$ implies that $D_2 + D_3 B_a > 0$ for all $B_a \in (0; B_{a_1})$.

Finally, using (27), (28) and (29) to replace $B_y$, $L$ and $P$ by functions of $B_a$ in (24), one has

$$G_3(B_a) = \frac{c \psi \alpha_a}{\gamma} B_a - \lambda \phi P(B_a) - \mu_L L(B_a) - \phi_L = 0 \quad (30)$$

that we will study on the biologically feasible interval $(0; B_{a_1})$ to determine the endemic equilibria.

The derivative of $L$ in (28) with respect to $B_a$ gives

$$L'(B_a) = - \frac{r_1 k_1 q^2 h^2 \gamma \delta_M}{\psi \alpha_a (k + \delta_M q h B_a)^2} < 0,$$

and its second derivative yields

$$L''(B_a) = \frac{r_1 k_1 q^3 h^3 \gamma \delta_M^2}{\psi \alpha_a (k + \delta_M q h B_a)^3} > 0. \quad (31)$$

Similarly, the first derivative of $P$ in (29) gives

$$P'(B_a) = D_1 \frac{\delta_M q h D_3 B_a^2 + 2 \delta_M q h D_2 B_a + D_2 k}{(D_2 + D_3 B_a)^2},$$

and

$$P''(B_a) = \frac{2 D_1 D_2 \delta_M r_1 k_1 \gamma \phi q^2 h^2}{(D_2 + D_3 B_a)^3}.$$

Because $D_2 + D_3 B_a > 0$ for all $B_a \in (0; B_{a_1})$, we then have $P''(B_a) > 0$.

Thus,

$$G_3''(B_a) = -\lambda \phi P''(B_a) - \mu_L L''(B_a) < 0. \quad (32)$$

It means that $G_3' : B_a \mapsto G_3'(B_a)$ is a decreasing function on $(0; B_{a_1})$. Remarking that $G_3(0) < 0$ since $P(0) = 0$ and $L(0) > 0$ when $\mathcal{R}_B > 1$, we can now distinguish two cases.

1. If $G_3'(0) \le 0$, the function $G_3$ defined by (30) does not have a positive solution on $(0; B_{a_1})$.

2. If $G_3'(0) > 0$, we distinguish two sub-cases.
Similarly, the first derivative of P in (29) gives

$$P'(B_a) = D_1 \frac{\delta_M q h D_3 B_a^2 + 2 \delta_M q h D_2 B_a + D_2 k}{(D_2 + D_3 B_a)^2},$$

and

$$P''(B_a) = \frac{2D_1 D_2 \delta_M r_1 k_1 \gamma \phi q^2 h^2}{(D_2 + D_3 B_a)^3}.$$

Because $D_2 + D_3 B_a > 0$ for all $B_a \in (0, B_{a_{c_1}})$, we then have $P''(B_a) > 0$.

Thus,

$$G_3''(B_a) = -\lambda \phi P''(B_a) - \mu_L L''(B_a) < 0. \quad (32)$$

It means that $G_3' : B_a \mapsto G_3'(B_a)$ is a decreasing function on $(0, B_{a_{c_1}})$. Remarking that $G_3(0) < 0$ since $P(0) = 0$ and $L(0) > 0$ when $R_B > 1$, we can now distinguish two cases.

1. If $G_3'(0) \le 0$, the function $G_3$ defined by (30) does not have a positive solution on $(0, B_{a_{c_1}})$.

2. If $G_3'(0) > 0$, we distinguish two sub-cases.

(a) If $G_3'(B_{a_{c_1}}) \ge 0$ then $G_3$ admits a unique positive root on $(0, B_{a_{c_1}})$ if and only if $G_3(B_{a_{c_1}}) > 0$.

(b) If $G_3'(B_{a_{c_1}}) < 0$ then, there exists a unique $\zeta \in (0, B_{a_{c_1}})$ such that $G_3'(\zeta) = 0$. By remarking that $G_3$ is a decreasing function on $(\zeta, B_{a_{c_1}})$, we always have $G_3(\zeta) > G_3(B_{a_{c_1}})$. And then:

* if $0 > G_3(\zeta) > G_3(B_{a_{c_1}})$, then $G_3$ has no positive root on $(0, B_{a_{c_1}})$;

* if $G_3(\zeta) > 0 > G_3(B_{a_{c_1}})$, then $G_3$ has two positive roots on $(0, B_{a_{c_1}})$;

* if $G_3(\zeta) > G_3(B_{a_{c_1}}) \ge 0$ or $G_3(\zeta) = 0$, then the function $G_3$ has a unique positive root on $(0, B_{a_{c_1}})$.

We can express the above conditions on $G_3$ in terms of the effective reproductive number $R_e$ given by (21) as follows.

Since

$$G_3(B_{a_{c_1}}) = \frac{c\psi\alpha_a}{\gamma} B_{a_{c_1}} - \lambda\phi P(B_{a_{c_1}}) - \mu_L L(B_{a_{c_1}}) - \phi_L \\ = \frac{c\psi\alpha_a}{\gamma} B_{a_{c_1}} - \lambda\phi P_{c_1} - \phi_L,$$

it is obvious that

$$G_3(B_{a_{c_1}}) > 0 \Leftrightarrow R_e > 1. \quad (33)$$

We will now analyze the threshold situations between the different sub-cases of case 2.(b). These situations are characterized by $G_3(\zeta) = 0$ and $G_3(B_{a_{c_1}}) = 0$. Due to the

concavity of the function $G_3$, these cannot be satisfied simultaneously for a unique set of parameters and we have $G_3(\zeta) > G_3(B_{a_{c_1}})$, so that

$$\frac{c\psi\alpha_a}{\gamma}\zeta - \mu_L L(\zeta) - \lambda\phi P(\zeta) > \frac{c\psi\alpha_a}{\gamma}B_{a_{c_1}} - \mu_L L(B_{a_{c_1}}) - \lambda\phi P(B_{a_{c_1}}) \quad (34)$$

In the following, we will analyze for which $\phi_L$ values $G_3(\zeta) = 0$ or $G_3(B_{a_{c_1}}) = 0$ are satisfied, all other parameters staying unchanged. Among other things, the choice of $\phi_L$ as variable parameter has the advantage of impacting neither the value of $\zeta$ nor that of $B_{a_{c_1}}$ so that inequality (34) still holds.

<figure><img src="image_13.png" /><figcaption>**Figure 8:** Plot of function $G_3(B_a)$ for different values of $\phi_L$: $\phi_L = 0.01 \text{ day}^{-1}$ (orange curve), $\phi_L^0 \approx 0.04$ (blue curve), $0.09 \text{ day}^{-1}$ (green curve), $\phi_L^s \approx 0.12 \text{ day}^{-1}$ (red curve) and $0.15 \text{ day}^{-1}$ (cyan curve). They correspond to hypothesis ($H_3$) of Theorem 6: $G_3'(0) = 1.86 \times 10^{-5} > 0$ and $G_3'(B_{a_{c_1}}) = -7.73 \times 10^{-4} < 0$. All remaining parameter values are given in Table 1. The stars are the zeros of $G_3$ in $(0, B_{a_{c_1}})$ and correspond to the endemic equilibria $E_c^*$ and $E_c^{**}$, whenever they exist. The vertical black line represents the line $B_a = B_{a_{c_1}}$ (PFE).</figcaption></figure>

If we consider that $G_3(\zeta) = 0$ holds for some $\phi_L^s$, we get

$$\phi_L^s = \frac{c\psi\alpha_a}{\gamma}\zeta - \mu_L L(\zeta) - \lambda\phi P(\zeta) \quad (35)$$

Similarly, $G_3(B_{a_{c_1}}) = 0$ for some $\phi_L^0$ gives

$$\phi_L^0 = \frac{c\psi\alpha_a}{\gamma} B_{a_{c_1}} - \lambda\phi P_{c_1}. \quad (36)$$

Injecting (35) and (36) into (34) then yields

$$\phi_L^s > \phi_L^0.$$

Noting in (21) that $\phi_L \mapsto R_e(\phi_L)$ is a decreasing function, we obtain $R_e(\phi_L^s) < R_e(\phi_L^0)$ with, as we have seen in (33)

$$R_e(\phi_L^0) = 1 \text{ and } R_e(\phi_L^s) = \frac{c\psi\alpha_a B_{a_{c_1}}}{\gamma(\lambda\phi_P P_{c_1} + \phi_L^s)}.$$
(a) If $G_3'(B_{a_{c_1}}) \ge 0$ then $G_3$ admits a unique positive root on $(0; B_{a_{c_1}})$ if and only if $G_3(B_{a_{c_1}}) > 0$.

(b) If $G_3'(B_{a_{c_1}}) < 0$ then, there exists a unique $\zeta \in (0; B_{a_{c_1}})$ such that $G_3''(\zeta) = 0$. By remarking that $G_3$ is a decreasing function on $(\zeta; B_{a_{c_1}})$, we always have $G_3(\zeta) > G_3(B_{a_{c_1}})$. And then:

* if $0 > G_3(\zeta) > G_3(B_{a_{c_1}})$, then $G_3$ has no positive root on $(0; B_{a_{c_1}})$;

* if $G_3(\zeta) > 0 > G_3(B_{a_{c_1}})$, then $G_3$ has two positive roots on $(0; B_{a_{c_1}})$;

* if $G_3(\zeta) > G_3(B_{a_{c_1}}) \ge 0$ or $G_3(\zeta) = 0$, then the function $G_3$ has a unique positive root on $(0; B_{a_{c_1}})$.

We can express the above conditions on $G_3$ in terms of the effective reproductive number $\mathcal{R}_e$ given by (21) as follows.

Since

$$ 
\begin{aligned}
G_3(B_{a_{c_1}}) &= \frac{c\psi\alpha_a}{\gamma}B_{a_{c_1}} - \lambda\phi P(B_{a_{c_1}}) - \mu_L L(B_{a_{c_1}}) - \phi_L \\
&= \frac{c\psi\alpha_a}{\gamma}B_{a_{c_1}} - \lambda\phi P_{c_1} - \phi_L,
\end{aligned}
 $$

it is obvious that

$$ G_3(B_{a_{c_1}}) > 0 \iff \mathcal{R}_e > 1. \qquad (33) $$

We will now analyze the threshold situations between the different sub-cases of case 2.(b). These situations are characterized by $G_3(\zeta) = 0$ and $G_3(B_{a_{c_1}}) = 0$. Due to concavity, these cannot be satisfied simultaneously for a unique set of parameters and we have $G_3(\zeta) > G_3(B_{a_{c_1}})$, so that

$$ \frac{c\psi\alpha_a}{\gamma}\zeta - \mu_L L(\zeta) - \lambda\phi P(\zeta) > \frac{c\psi\alpha_a}{\gamma}B_{a_{c_1}} - \mu_L L(B_{a_{c_1}}) - \lambda\phi P(B_{a_{c_1}}) \qquad (34) $$

In the following, we will analyze for which $\phi_L$ values $G_3(\zeta) = 0$ or $G_3(B_{a_{c_1}}) = 0$ are satisfied, all other parameters staying unchanged. Among other things, the choice of $\phi_L$ as varying parameter has the advantage of impacting neither the value of $\zeta$ nor that of $B_{a_{c_1}}$ so that inequality (34) still holds.

If we consider that $G_3(\zeta) = 0$ holds for some $\phi_L^s$, we get

$$ \phi_L^s = \frac{c\psi\alpha_a}{\gamma}\zeta - \mu_L L(\zeta) - \lambda\phi P(\zeta) \qquad (35) $$

Similarly, $G_3(B_{a_{c_1}}) = 0$ for some $\phi_L^0$ gives

$$ \phi_L^0 = \frac{c\psi\alpha_a}{\gamma}B_{a_{c_1}} - \lambda\phi P_{c_1}. \qquad (36) $$

Injecting (35) and (36) into (34) then yields

$$ \phi_L^s > \phi_L^0. $$

Noting in (21) that $\phi_L \mapsto \mathcal{R}_e(\phi_L)$ is a decreasing function, we obtain $\mathcal{R}_e(\phi_L^s) < \mathcal{R}_e(\phi_L^0)$ with, as we have seen in (33)

$$ \mathcal{R}_e(\phi_L^0) = 1 \text{ and } \mathcal{R}_e(\phi_L^s) = \frac{c\psi\alpha_a B_{a_{c_1}}}{\gamma(\lambda\phi_P P_{c_1} + \phi_L^s)}. $$

Thus, considering a given $\phi_L$ value, and since $G_3$ is a decreasing function of $\phi_L$, we make multiple use of

$$ 0 > G_3(\zeta) \Leftrightarrow G_3(\zeta)|_{\phi_L^s} > G_3(\zeta)|_{\phi_L} \Leftrightarrow \phi_L > \phi_L^s $$

and

$$ 0 > G_3(B_{a_{c_1}}) \Leftrightarrow G_3(B_{a_{c_1}})|_{\phi_L^0} > G_3(B_{a_{c_1}})|_{\phi_L} \Leftrightarrow \phi_L > \phi_L^0. $$

And then,

* $0 > G_3(\zeta) > G_3(B_{a_{c_1}}) \Leftrightarrow \phi_L > \phi_L^s > \phi_L^0 \Leftrightarrow \mathcal{R}_e(\phi_L) < \mathcal{R}_e(\phi_L^s) < 1$ (region above blue curve in Fig. 8).

* $G_3(\zeta) > 0 > G_3(B_{a_{c_1}}) \Leftrightarrow \phi_L^s > \phi_L > \phi_L^0 \Leftrightarrow \mathcal{R}_e(\phi_L^s) < \mathcal{R}_e(\phi_L) < 1$ (region between the blue and red curves in Fig. 8).

* $G_3(\zeta) > G_3(B_{a_{c_1}}) > 0 \Leftrightarrow \phi_L^s > \phi_L^0 > \phi_L \Leftrightarrow \mathcal{R}_e(\phi_L^s) < 1 < \mathcal{R}_e(\phi_L)$ (region below the red curve in Fig. 8).

In the next, $\mathcal{R}_e$ is exactly the threshold given by (21). The existence of endemic equilibria can be summarized as follows.
Thus, considering a given $\phi_L$ value, and since $G_3$ is a decreasing function of $\phi_L$, we make multiple use of

$$0 > G_3(\zeta) \Leftrightarrow G_3(\zeta), \phi_L^s > G_3(\zeta), \phi_L \Leftrightarrow \phi_L > \phi_L^s$$

and

$$0 > G_3(B_{a_{c_1}}) \Leftrightarrow G_3(B_{a_{c_1}}), \phi_L^0 > G_3(B_{a_{c_1}}), \phi_L \Leftrightarrow \phi_L > \phi_L^0.$$

Moreover,

* $0 > G_3(\zeta) > G_3(B_{a_{c_1}}) \Leftrightarrow \phi_L > \phi_L^s > \phi_L^0 \Leftrightarrow R_e(\phi_L) < R_e(\phi_L^s) < 1$ (region above blue curve in Fig. 8).

* $G_3(\zeta) > 0 > G_3(B_{a_{c_1}}) \Leftrightarrow \phi_L^s > \phi_L > \phi_L^0 \Leftrightarrow R_e(\phi_L^s) < R_e(\phi_L) < 1$ (region between the blue and red curves in Fig. 8).

* $G_3(\zeta) > G_3(B_{a_{c_1}}) > 0 \Leftrightarrow \phi_L^s > \phi_L^0 > \phi_L \Leftrightarrow R_e(\phi_L^s) < 1 < R_e(\phi_L)$ (region below the red curve in Fig. 8).

In what follows, $R_e$ is exactly the threshold given by (21). The existence of endemic equilibria can be summarized as follows.

**Theorem 6.** Assume that $R_B > 1$. Then the existence of endemic equilibria of sub-model (20) depends on $G'_3(0)$ and $G'(B_{a_{c_1}})$, with $G_3$ defined in (30):

(H1) $G'_3(0) \le 0$: sub-model (20) has no endemic equilibrium.

(H2) $G'_3(0) > 0$ and $G'_3(B_{a_{c_1}}) \ge 0$: sub-model (20) has a unique endemic equilibrium if and only if $R_e > 1$.

(H3) $G'_3(0) > 0$ and $G'_3(B_{a_{c_1}}) < 0$:

(i) If $R_e < R_e(\phi_L^s)$, then sub-model (20) has no endemic equilibrium.

(ii) If $R_e(\phi_L^s) < R_e < 1$, then sub-model (20) has two endemic equilibria, $E_c^*$ and $E_c^{**}$ such that $B_{a_c}^* < B_{a_c}^{**}$.

(iii) If $R_e = R_e(\phi_L^s)$ or $R_e \ge 1$, then sub-model (20) has a unique endemic equilibrium $E_c^*$.

Fig 8 illustrates the situation described above for different values of $\phi_L$. The red curve represents the function $G_3$ at the critical value $\phi_L^s \approx 0.12$ day$^{-1}$. Below this value, there is no root (cyan curve). Above this value, $G_3$ can have two (green curve) or one (orange curve) roots depending on the value of $\phi_L$.

### 4.1.2. Stability

Following a similar approach as in Section 3, we establish the following results for the trivial and pest-free equilibria.

**Theorem 7.** The extinction equilibrium, $E_{c_0}$, of sub-model (20), is globally asymptotically stable if and only if $R_B \le 1$. Otherwise, if $R_B > 1$, the pest-free equilibrium $E_{c_1}$ exists. It is locally asymptotically stable if $R_e < 1$ and unstable if $R_e > 1$.

Moreover, the system presents a forward transcritical bifurcation at $R_B = 1$.

**PROOF.** We suppose $R_B \le 1$ and we consider the following positive definite function,

$$V = qh\gamma B_y + qh\alpha_y B_a + \frac{qh\gamma}{c} L + \frac{k\alpha_a \alpha_y (1 - R_B)}{a_p \alpha} P + \alpha_a \alpha_y M.$$

Differentiating the above function $V$ with respect to time, $t$, along the solutions of sub-model (20), after some algebraic simplifications, we obtain

$$\dot{V} = -\frac{qh\gamma}{c}(\mu_L L^2 + \phi\lambda PL + \phi_L L) - \frac{k\alpha_a \alpha_y (1 - R_B)}{a_p \alpha} (\phi P L + \phi_p P) + \left[k\alpha_a \alpha_y (1 - R_B) + \frac{qh\gamma r_1 k_1}{1 + \delta_M M} - k\alpha_a \alpha_y\right] M.$$

And,

$$k\alpha_a \alpha_y (1 - R_B) + \frac{qh\gamma r_1 k_1}{1 + \delta_M M} - k\alpha_a \alpha_y < k\alpha_a \alpha_y (1 - R_B) + qh\gamma r_1 k_1 - k\alpha_y \alpha_a = 0.$$

Thus, $V$ is a strict Lyapunov function for equilibrium $E_{c_0}$ if $R_B \le 1$. We conclude that $E_{c_0}$ is globally asymptotically stable in $\Omega$ if $R_B \le 1$.

We assume that $R_B > 1$. The characteristic polynomial evaluated at the PFE $E_{c_1}$ is given by:

$$Q_{c_1}(\chi) = (\chi - (c\psi B y_{c_1} - \lambda\phi P_{c_1} - \phi_L))(\chi + \phi_P)(\chi^3 + B_2\chi^2 + B_1\chi + B_0),$$

where

$$B_2 = k + \alpha_y + \alpha_a, \\ B_1 = \alpha_y \alpha_a + k(\alpha_a + \alpha_y), \\ B_0 = k\alpha_y \alpha_a - \frac{k_1 qh\gamma r_1}{(1 + \delta_M M_{c_1})^2}.$$

We conclude as in Appendix (C).

The proof of forward transcritical bifurcation at $R_B = 1$ is trivial.

**Remark 1.** Under the hypothesis (H1) or (H2) of Theorem 6, we can prove, by adapting the proof of Theorem 4, that all the solutions starting in $\Omega''$ tend to the pest-free equilibrium $E_{c_1}$ if $R_e < 1 < R_B$, with $\Omega''$

$$\Omega'' = \{(B_y; B_a; L; P; M) \in \mathbb{R}_+^5 | B_y \le B_{y_{c_1}}, B_a \le B_{a_{c_1}}, M \le M_{c_1}\}. \quad (37)$$

If hypothesis (H2) holds, the system presents a transcritical bifurcation at the point $R_e = 1$.
<figure><img src="image_11.png" /><figcaption>Figure 8: Plots of function $G_3(B_a)$ for different values of $\phi_L$: $\phi_L = 0.01$ day$^{-1}$ (orange curve), $\phi_L^0 \approx 0.04$ (blue curve), 0.09 day$^{-1}$ (green curve), $\phi_L^s \approx 0.12$ day$^{-1}$ (red curve) and 0.15 day$^{-1}$ (cyan curve). They correspond to hypothesis ($H_3$) of Theorem 6: $G'_3(0) = 1.86 \times 10^{-5} > 0$ and $G'_3(B_{a_{c_1}}) = -7.73 \times 10^{-4} < 0$. All remaining parameter values are given in Table 1. The stars are the zeros of $G_3$ in $(0; B_{a_{c_1}})$ and the correspond to the endemic equilibria $E_c^*$ and $E_c^{**}$, whenever they exist. The vertical black line represents the line $B_a = B_{a_{c_1}}$ (PFE).</figcaption></figure>

**Theorem 6.** Assume that $\mathcal{R}_B > 1$. Then the existence of endemic equilibria of sub-model (20) depends on $G'_3(0)$ and $G'(B_{a_{c_1}})$, with $G_3$ defined in (30):

($H_1$) $G'_3(0) \le 0$: sub-model (20) has no endemic equilibrium.

($H_2$) $G'_3(0) > 0$ and $G'_3(B_{a_{c_1}}) \ge 0$: sub-model (20) has a unique endemic equilibrium if and only if $\mathcal{R}_e > 1$.

($H_3$) $G'_3(0) > 0$ and $G'_3(B_{a_{c_1}}) < 0$:

(i) If $\mathcal{R}_e < \mathcal{R}_e(\phi_L^s)$, then sub-model (20) has no endemic equilibrium.

(ii) If $\mathcal{R}_e(\phi_L^s) < \mathcal{R}_e < 1$, then sub-model (20) has two endemic equilibria, $E_c^*$ and $E_c^{**}$ such that $B_{a_c}^* < B_{a_c}^{**}$.

(iii) If $\mathcal{R}_e = \mathcal{R}_e(\phi_L^s)$ or $\mathcal{R}_e \ge 1$, then sub-model (20) has a unique endemic equilibrium $E_c^*$.

Fig 8 illustrates the situation described above for different values of $\phi_L$. The red curve represents the function $G_3$ at the critical value $\phi_L^s \approx 0.12$ day$^{-1}$. Below this value, there is no root (cyan curve). Above this value, $G_3$ can have two (green curve) or one (orange curve) roots depending on the value of $\phi_L$.

### 4.1.2 Stability

Following a similar approach as in Section 3, we establish the following results for the trivial and pest-free equilibria.

**Theorem 7.** The extinction equilibrium, $E_{c_0}$, of sub-model (20), is globally asymptotically stable if $\mathcal{R}_B < 1$. Otherwise, if $\mathcal{R}_B > 1$, the pest-free equilibrium $E_{c_1}$ exists. It is locally asymptotically stable if $\mathcal{R}_e < 1$ and unstable if $\mathcal{R}_e > 1$.

Moreover, the system presents a forward transcritical bifurcation at $\mathcal{R}_B = 1$.

*Proof.* We suppose $\mathcal{R}_B < 1$ and we consider the following positive definite function,

$$V = qh\gamma B_y + qh\alpha_y B_a + \frac{qh\gamma}{c} L + \frac{k\alpha_a\alpha_y(1 - \mathcal{R}_B)}{a_P\alpha} P + \alpha_a\alpha_y M.$$

Differentiating the above function $V$ with respect to time, $t$, along the solutions of sub-model (20), after some algebraic simplifications, we obtain

$$\begin{aligned} \dot{V} = & -\frac{qh\gamma}{c}(\mu_L L^2 + \phi\lambda P L + \phi_L L) \\ & -\frac{k\alpha_a\alpha_y(1 - \mathcal{R}_B)}{a_P\alpha}(\phi P L + \phi_P P) \\ & + \left[k\alpha_a\alpha_y(1 - \mathcal{R}_B) + \frac{qh\gamma r_1 k_1}{1 + \delta_M M} - k\alpha_a\alpha_y\right] M. \end{aligned}$$
<figure><img src="image_14.png" /><figcaption>Figure 9: Backward bifurcation diagram for sub-model (20), obtained by varying parameter: (a) $\phi_L$, with $R_e(\phi_L^*) \approx 0.68$ and (b) $\alpha$, with $R_e(\alpha^*) \approx 0.2$. All remaining parameter values are given in Table 1. The blue line represents pest-free $E_{c1}$. The solid and dashed red lines represent endemic equilibria $E_c^*$ and $E_c^{**}$. A solid line corresponds to a stable equilibrium, a dashed line to an unstable equilibrium.</figcaption></figure>

We explore the stability of the endemic equilibria, whenever they exist, numerically. Based on the parameter values in Table 1 (with $\delta_M = 1.17 \times 10^{-2}$ per USD), we vary parameter $\phi_L$. We find that $E_c^*$ is stable and $E_c^{**}$ is unstable. Thus, a bistability behavior occurs when $R_e(\phi_L^*) < R_e < 1$ (corresponding to $\phi_L \in (\phi_L^0, \phi_L^s)$).

The backward bifurcation diagram is shown in Fig. 9(a). In practice, the smallholder farmer cannot influence the natural mortality of larvae $\phi_L$.

Varying $\alpha$, the investment rate for pesticide acquisition, is more interesting from a practical point of view. Determining the threshold $\alpha^s$ corresponding to the lower bound for the existence of endemic equilibria, is not as straightforward as with parameter $\phi_L$. So we proceed numerically and obtained $\alpha^s \approx 0.086 \text{ day}^{-1}$ and $R_e(\alpha^s) \approx 0.2$. The threshold $\alpha^0 = \phi_P(c\psi B_{y_{c1}} - \phi_L)/(\lambda\phi a_p M_{c1})$, obtained by solving $R_e = 1$ in terms of $\alpha$. As previously, bi-stability behavior occurs when $R_e(\alpha^s) < R_e < 1$. The backward bifurcation diagram is shown in Fig. 9(b). We observe that $\alpha$ has more significant effect to the equilibria values than $\phi_L$.

## 4.2. Simulations

In this section, we carry out numerical simulations of sub-model (20). We recall that all parameter values are given in Table 1 (with $\delta_M = 1.17 \times 10^{-2}$ per USD). The biomass reproductive number $R_B = 3.21$ and the effective reproductive number $R_e = 0.85$.

Fig. 10(a) reveals that a high consumption or a low natural mortality rate makes pest management difficult. Fig. 10(b) reveals that the success of pest control depends on both the toxicity of the pesticide and the financial resources allocated. For a botanical pesticide with relatively low toxicity to eliminate larvae from the plantation, massive (and fairly unreasonable) investment is required. Conversely, a highly effective pesticide will require very little financial resources to fight the infestation. Thus, the smallholder

farmer needs to choose its products wisely for a cost-efficient plant protection.

Hypothesis (H1) of Theorem 6 corresponds to a situation where the botanical pesticides is extremely powerful ($\lambda \ge 7.5$ larvae/l). In practice, this situation cannot happen. Hypothesis (H2) of the same theorem corresponds to a situation where the smallholder farmer invests almost nothing in botanical pesticides (for example $\alpha = 5 \times 10^{-4} \text{ day}^{-1}$). Thus, the larvae act like there is no pesticides and sub-model (20) has the same dynamics as sub-model (6). Therefore, simulations were done under hypothesis (H3) of Theorem 6.

One can compute the cost of plant protection using the metric

$$ \alpha \overline{M} A_t, \qquad (38) $$

where $\overline{M}$ is the value of the component $M$ of the stable equilibrium and $A_t = 90$ days still denotes the average time of the agricultural season.

In our case (parameter values given in Table 1), when the pest management strategy works nicely (i.e. we are in the basin of attraction of the pest-free equilibrium $E_1$), using (19), we have a global yield around 22.2 tonnes for a total income of 15,320 USD. Using (38), the cost of plant protection is around 230 USD. These two metrics agree with the work in Amoabeng et al. [3]. Conversely, if the pest strategy does not work nicely (i.e. we are in the basin of attraction of the endemic equilibrium $E^*$), the expected yield will be around 10.4 tonnes for an income of 7,180 USD and a cost protection of only 107 USD. As shown in Fig. 11, in the bistability region, one cannot predict the behavior of the system. Indeed, even with a low initial level of larvae, the system can tend towards an endemic equilibrium (cyan curves), while it tends towards the pest-free equilibrium with a higher level (blue curves), all remaining initial conditions being the same.
And,

$$k\alpha_a\alpha_y(1 - \mathcal{R}_B) + \frac{qh\gamma r_1 k_1}{1 + \delta_M M} - k\alpha_a\alpha_y < k\alpha_a\alpha_y(1 - \mathcal{R}_B) + qh\gamma r_1 k_1 - k\alpha_y\alpha_a = 0.$$

Thus, $V$ is a strict Lyapunov function for equilibrium $E_{c_0}$ if $\mathcal{R}_B < 1$. We conclude that $E_{c_0}$ is globally asymptotically stable in $\Omega$ if $\mathcal{R}_B < 1$.

We assume that $\mathcal{R}_B > 1$. The characteristic polynomial evaluated at the PFE $E_{c_1}$ is given by:

$$Q_{c_1}(\chi) = (\chi - (c\psi By_{c_1} - \lambda\phi P_{c_1} - \phi_L))(\chi + \phi_P) (\chi^3 + B_2\chi^2 + B_1\chi + B_0),$$

where

$$B_2 = k + \alpha_y + \alpha_a, \\ B_1 = \alpha_y\alpha_a + k(\alpha_a + \alpha_y), \\ B_0 = k\alpha_y\alpha_a - \frac{k_1qh\gamma r_1}{(1 + \delta_M M_{c_1})^2}.$$

We conclude as in Appendix (C).

The proof of forward transcritical bifurcation at $\mathcal{R}_B = 1$ is trivial. $\square$

**Remark 1.** Under the hypothesis $(H_1)$ or $(H_2)$ of Theorem 6, we can prove, by adapting the proof of Theorem 5, that all the solutions starting in $\Omega''$ tend to the pest-free equilibrium $E_{c_1}$ if $\mathcal{R}_e \le 1 < \mathcal{R}_B$, with $\Omega''$

$$\Omega'' = \{(B_y; B_a; L; P; M) \in \mathbb{R}_+^5 | B_y \le B_{y_{c_1}}, B_a \le B_{a_{c_1}}, M \le M_{c_1}\}. \quad (37)$$

If hypothesis $(H_2)$ holds, the system presents a transcritical bifurcation at the point $\mathcal{R}_e = 1$.

We explored the stability of the endemic equilibria, whenever they existed, numerically. Based on the parameter values in Table 1 (with $\delta_M = 1.17 \times 10^{-2}$ per USD), we varied parameter $\phi_L$. We found that $E_c^*$ is stable and $E_c^{**}$ is unstable. Thus, a bistability behavior occurs when $\mathcal{R}_e(\phi_L^s) < \mathcal{R}_e < 1$ (corresponding to $\phi_L \in (\phi_L^0; \phi_L^s)$). The backward bifurcation diagram is shown in Fig. 9(a). In practice, the smallholder farmer cannot influence the natural mortality of larvae $\phi_L$.

Varying $\alpha$, the investment rate for pesticide acquisition, is more interesting from a practical point of view. Determining the threshold $\alpha^s$ corresponding to the lower bound for the existence of endemic equilibria, is not as straightforward as with parameter $\phi_L$. So we proceeded numerically and obtained $\alpha^s \approx 0.086$ day$^{-1}$ and $\mathcal{R}_e(\alpha^s) \approx 0.2$. The threshold $\alpha^0 = \phi_P(c\psi By_{c_1} - \phi_L)/(\lambda\phi_a M_{c_1})$, obtained by solving $\mathcal{R}_e = 1$ in terms of $\alpha$. As previously, bi-stability behavior occurs when $\mathcal{R}_e(\alpha^s) < \mathcal{R}_e < 1$. The backward bifurcation diagram is shown in Fig. 9(b). We observe that $\alpha$ has more significant effect to the equilibria values than $\phi_L$.

<figure><img src="image_12.png" /><figcaption>Figure 9: Backward bifurcation diagram for sub-model (20), obtained by varying parameter: (a) $\phi_L$, with $\mathcal{R}_e(\phi_L^s) \approx 0.68$ and (b) $\alpha$, with $\mathcal{R}_e(\alpha^s) \approx 0.2$. All remaining parameter values are given in Table 1. The blue line represents pest-free $E_{c_1}$. The solid and dashed red lines represent endemic equilibria $E_c^*$ and $E_c^{**}$. A solid line corresponds to a stable equilibrium, a dashed line to an unstable equilibrium.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 10: Variation of the effective reproductive number $R_e$ as a function of: (a) the natural mortality of larvae $\phi_L$ and the biomass consumption rate by larvae $\psi$; (b) the investment rate for plant protection $\alpha$ and the efficacy of pesticide $\lambda$. The values of the remaining parameters are given in Table 1. The purple squares represent the parameter values in Table 1. The red line represents the level 1</figcaption></figure>

## 4.2 Simulations

In this section, we carry out numerical simulations of sub-model (20). We recall that all parameters values are given in Table 1 (with $\delta_M = 1.17 \times 10^{-2}$ per USD). The biomass reproductive number $R_B = 3.21$ and the effective reproductive number $R_e = 0.85$.

Fig. 10(a) reveals that a high consumption or a low natural mortality rate makes pest management difficult. Fig. 10(b) reveals that the success of pest control depends on both the efficacy of the pesticide and the financial resources allocated. For a botanical pesticide with relatively low efficacy to eradicate larvae from the plantation, massive (and fairly unreasonable) investment is required. Conversely, a highly effective pesticide will require very little financial resources to fight the infestation. Thus, the smallholder farmer needs to choose its products wisely for a cost-efficient plant protection.

Hypothesis ($H_1$) of Theorem 6 corresponds to a situation where the botanical pesticides is extremely powerful ($\lambda \ge 7.5$ larvae/l). In practice, this situation cannot happen. Hypothesis ($H_2$) of the same theorem corresponds to a situation where the smallholder farmer invests almost nothing in botanical pesticides (for example $\alpha = 5 \times 10^{-4}$ day$^{-1}$). Thus, the larvae act like there is no pesticides and sub-model (20) has the same dynamics as sub-model (6). Therefore, simulations were done under hypothesis ($H_3$) of Theorem 6.

One can compute the cost of plant protection using the metric

$$ \alpha \times \overline{M} \times A_t, \qquad (38) $$

where $\overline{M}$ is the value of the component $M$ of the stable equilibrium and $A_t = 90$ days still denotes the average time of the agricultural season.

When the pest management strategy works nicely (i.e. the pest-free equilibrium is stable), using (19), we have a global yield around 22.2 tonnes for a total income of 15,320 USD. Using (38), the cost of plant protection is around 230 USD. These two metrics agree with the work in [30]. Conversely, if the pest strategy does not work nicely (i.e. an endemic equilibrium is stable), the expected yield will be around 10.4 tonnes for an income of 7,180 USD and a cost protection of only 107 USD. As shown in Fig. 11, in the bistabilty region, one cannot predict the behavior of the system. In fact, with a low level of larvae, the system tends towards endemic equilibrium, while with a high level of larvae, it tends towards pest-free equilibrium (blue and cyan curves).

We observe numerically that the original (3) and simplified (20) models present the same dynamics but with different equilibrium values, as shown in Appendix E.

## 5 Conclusion

With food security a growing problem, many governments have put agriculture and rural development at the heart of their public policies. This action is hampered by the presence of pests, which destroy crops and cause huge economic losses.

The aim of this study was to quantify the damage caused by the attack of *Plutella* larvae in a cabbage plantation, assess the effectiveness of using botanical pesticides as a pest management strategy, and evaluate the sustainability and profitability of such agroecosystem. To achieve our goals, we built an original mathematical
<figure><img src="image_15.png" /></figure>

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

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

<figure><img src="image_18.png" /><figcaption>**Figure 11:** Dynamics of state variables (a) $B_y$, (b) $L$, (c) $P$ and (d) $M$ of sub-model (20). Parameter values are given in Table 1. The initial conditions ($B_y(0)$; $B_a(0)$; $L(0)$; $M(0)$) are: (2791.24; 0; 51.89; 0; 0) in green, (4500; 0; 51.89; 0; 0) in blue, (2691.24; 0; 20; 0; 0) in red and (4500; 0; 0.65; 0; 0) in cyan. The effective reproductive number $R_e$ verifies $R_e(\alpha^s) \approx 0.2 < R_e = 0.85 < 1$. Bistability occurs: depending on the initial conditions, the system either converges towards the pest-free (red and blue) or towards the endemic (cyan and green) equilibrium.</figcaption></figure>

In the first case study, we considered a self-financing uncontrolled model, that is a model without any pest management corresponding to sub-model (6). We showed that when the condition $R_B > 1$ is satisfied, the threshold $R_0$ governs the dynamics of the resulting system. The system presents a transcritical bifurcation at $R_0 = 1$. When $R_0 < 1$, the larvae disappear naturally without any control action and the production is around 11.1 tonnes per hectare. If $R_0 > 1$, the presence of larvae leads to losses greater than 30% of the overall production.

In the second case study, we studied a slightly modified controlled model (20), obtained from the main model (3) by neglecting the antifeedant effect of botanical pesticide and the competition in the biomass compartments. We computed the effective reproductive number $R_e$ and noticed that it is smaller than $R_0$, suggesting that integrated pest management acts efficiently against the infestation. A botanical pesticide is a judicious choice because it is effective for both prophylaxis and treatment (as reported by Pavela [36]) with a reasonable treatment cost for the smallholder farmer. The presence of larvae leads to losses greater than 45% of the overall production. This percentage is higher than in the

uncontrolled case, which may seem contradictory, but can be explained by the absence of biomass competition in sub-model (20). This is consistent with the sensitivity analysis, which highlights the negative effect of having too much young biomass on the income. Unfortunately, the threshold $R_e$ does not solely govern the dynamics of the system.

Our model (both the original model and the slightly modified model) exhibits a backward bifurcation revealing that having $R_e$ less than one is not enough to guarantee the elimination of larvae in the plantation. Eradication occurs when the control is strong enough (i.e. $R_e < R_e(\alpha^s)$) which is not necessarily the best option for the smallholder farmer, because important investment in plant protection will reduce his net profit. This bistability linked to the backward bifurcation is dangerous because, depending on the initial condition, the system can converge either towards the pest-free equilibrium or towards an endemic equilibrium. However, it is possible to avoid this bi-stability by investing substantially in plant protection ($\alpha \ge 0.09 \text{ day}^{-1}$).

In this work, we assumed that the smallholder farmer has sufficient fertile land. Unfortunately, not all farmers have access to large amounts of arable and properly irrigated
<figure><img src="image_14.png" /><figcaption>Figure 11: Dynamics of state variables (a) $B_y$, (b) $L$, (c) $P$ and (d) $M$ of sub-model (20). Parameter values are given in Table 1. The initial conditions are: (2791.24; 0; 51.89; 0; 0) in green, (4500; 0; 51.89; 0; 0) in blue, (2691.24; 0; 20; 0; 0) in red and (4500; 0; 0.65; 0; 0) in cyan. The effective reproductive number $R_e$ verifies $R_e(\alpha^s) \approx 0.2 < R_e = 0.85 < 1$. Bistability occurs: according to the initial conditions, the system either converges towards the pest-free (red and blue) or the endemic (cyan and green) equilibrium.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 12:** Sensitivity analysis for economic and pesticide-related parameters on the seasonal income. All remaining parameter values are given in Table 1. The sample size for the latin hypercube sampling is N=3000. Partial rank correlation coefficients (PRCC) are computed according to [27]. Asterisks indicate significant PRCC (p-values below 0.01).</figcaption></figure>

Taking into account the limitation of fertile land could therefore be the focus of another case study, as well as including irrigation, soil fertilization and labor costs. Also, the functioning of the current account is very basic, and does not allow human actions to be properly described. Thus, another case study could be to consider discrete events to represent harvesting, cash withdrawals and investments. Furthermore, simplifying the pest's life cycle to its larval phase does not take into account the inherent delays in this cycle. Finally, another approach to modeling pesticide application in the plantation is to use a feedback control strategy (as proposed by Misra et al. [29]) or to combine prophylaxis and feedback control. All these extensions give rise to some modeling challenges which represent an avenue of our further works.

## A. Parameter estimation

Some parameters were estimated from the literature. The values of the conversion rate c and uptake rate of botanical pesticide by larvae were taken in Misra et al. [29]. Amoabeng et al. [3] reported that undamaged heads fetched 0.56 USD and 0.83 USD per kg for the major and minor seasons respectively. We chose an average cabbage price q = 0.69 USD per kg.

When replanted, the cabbage consists of 3 to 5 leaves [4, 2]. The weight at transplanting is therefore assumed to be 70g. We assume a growth period for the cabbage of 90 days (three months) for an average weight of 0.7 kg per head [2]. Sow et al. [44] report that there is a strong correlation between the age of the plant and the age of the diamondback moth larvae. In fact, plants less than 55 days old are favored by immature larvae (stages 2 and 3). Knowing also that larvae populations are more voracious at stage 3, we set the aging rate of young biomass $\gamma = 1/55 \approx 0.02 \text{ day}^{-1}$ and the harvest rate $h = 1/(90 - 55) \approx 0.03 \text{ day}^{-1}$.

One hectare can support around 35000 cabbage plants [3]. Then, the carrying capacity of a 1-hectare field is $35,000 \times 0.7 = 24,500$ kg, split into $K_a = 24,500 \times 35/90 = 9,528$ kg of adult biomass and $K_y = 24,500 \times 55/90 = 14,972$ kg of young biomass. Assuming that the weight of a young cabbage at 55 days after replanting is 0.4 kg, we compute the growth rate of adult biomass $r_{0a} = \ln(0.7/0.4)/35 \approx 1.6 \times 10^{-2} \text{ day}^{-1}$. Moreover, we assume that $r_{0y} \approx r_{0a} = 1.6 \times 10^{-2} \text{ day}^{-1}$. Thus, the competition rates are given by $\mu_a = 1.6 \times 10^{-2}/9528 \approx 1.68 \times 10^{-6} \text{ kg}^{-1}\text{day}^{-1}$, and, $\mu_y = 1.6 \times 10^{-2}/14972 \approx 1.07 \times 10^{-6} \text{ kg}^{-1}\text{day}^{-1}$.

According to Gautam et al. [15], the infestation is noticeable at more than 20 larvae per plant. In addition, Prasad et al. [37] reported that an increase of 5 larvae leads to an estimated increase in consumption of 19 g and 32 g for an exposure of 4 and 7 weeks respectively. Hence, we assume that the biomass consumption rate by larvae $\psi$ has values between $0.032/((20 + 5) \times 7 \times 7) \approx 2.27 \times 10^{-5}$ and $0.019/((20 + 5) \times 7 \times 4) \approx 4.57 \times 10^{-5} \text{ larvae}^{-1}\text{day}^{-1}$. The larvae development stage lasts about 9 to 32 days [15, 18]. Hence, we assume that the mortality rate of larvae $\phi_L$ is between 0.03 and $0.1 \text{ day}^{-1}$. Assuming a carrying capacity of 30 larvae per plant [15, 35], the carrying capacity of the plantation is $35,000 \times 30 = 1,050,000$ larvae. The competition of larvae $\mu_L$ is given by formula $c\psi K_y/1,050,000$. Thus, by considering the conversion efficiency $c = 0.6 \text{ larvae.kg}^{-1}$ [29], the range of $\mu_L$ is $[1.94 \times 10^{-6}, 3.91 \times 10^{-6}]$.

According to Amoabeng et al. [3], a hectare is treated between 7 and 8 times over a growing season. It is hence assumed that the natural depletion of botanical pesticide is $\phi_p = 8/90 \approx 0.1 \text{ day}^{-1}$. In addition, around 500 liters of pesticide per hectare are used [36] for a total cost around 306 USD [3]. Thus, the inverse of the pesticide price is $a_p \approx 500/231 \approx 2.16 \text{ liters}\text{USD}^{-1}$.

According to jardins de la Vallée de la Siagne [41], the unit price of a young cabbage plant is around 0.16 USD per plant, which weighs 70g. Thus, the inverse of the price of young biomass is estimated at $r_1 = 0.07/0.16 \approx 0.46 \text{ kg}\text{USD}^{-1}$.

## B. Proof of Theorem 1

Let $(B_y(t); B_a(t); L(t); P(t); M(t))$ be a local solution of model (3) defined on $(0; T_+)$ where $T_+ > 0$. Using the first equation of model (3), we have

$$ \begin{aligned} \dot{B}_y &= r_{0y}B_y - \mu_y B_y^2 + \frac{r_1 k_1 M}{1 + \delta_M M} - \frac{\psi L B_y}{bP + 1} - \gamma B_y \\ &\leq \frac{r_1 k_1}{\delta_M} - \alpha_y B_y. \end{aligned} $$

Applying Gronwall's lemma (Proposition 1.39 in [7]), we obtain

$$ B_y(t) \leq \frac{r_1 k_1}{\delta_M \alpha_y} + \left(B_y(0) - \frac{r_1 k_1}{\delta_M \alpha_y}\right) \exp(-\alpha_y t). \quad (B.1) $$
model that describes the agro-economic environment associated with cabbage culture. Indeed, we divided the
cabbage biomass into two groups to account for the larvae preference for young cabbage and natural resistance
of adult cabbage due to the firmness of leaves acquired with age. In order to take into account the self-financing
constraints of the farm, we introduced the plantation current account as a dynamical variable in the model.
This budget is allocated for purchasing new plants, botanical pesticides, and paying daily expenditures of the
smallholder farmer.

In this work, we showed that the farm sustainability depends on the biomass reproductive number $R_B$. If $R_B$ is smaller than one, the system will inevitably disappear. Conversely, $R_B$ values greater than one guarantee the sustainability of the entire ecosystem. Except in extreme market situations (i.e. very high purchase price of young plants $r_1^{-1} > 3$ USD/kg and very low selling price of mature cabbages $q < 0.4$ USD/kg) or when the farmer investment in the farm is particularly low (i.e. $k_1/k < 0.21$), analysis of this threshold shows that cabbage is an very sustainable crop. The mathematical analysis of two variants of the model proposed was carried out in details.

In the first case study, we considered a self-financing uncontrolled model; that is a model without any pest
management. We showed that when the condition $R_B > 1$ is satisfied, the threshold $R_0$ governs the dynamics
of the resulting system. The system presents a transcritical bifurcation at $R_0 = 1$. When $R_0 < 1$, the larvae
disappear naturally without any control action and the production is around 11.1 tonnes per hectare. If $R_0 > 1$,
the presence of larvae leads to losses greater than 30% of the overall production.

In the second case study, we studied a slightly modified controlled model (20), obtained from the main
model (3) by neglecting the antifeedant effect of botanical pesticide and the competition in the biomass com-
partments. We computed the effective reproductive number $R_e$ and noticed that it is smaller than $R_0$, suggesting
that integrated pest management acts efficiently against the infestation. A botanical pesticide is a judicious
choice because it is effective in both for prophylaxis and treatment(as reported by [24]) with a reasonable treat-
ment cost for the smallholder farmer. Unfortunately, the threshold $R_e$ does not solely govern the dynamics of
the system.

Our model (both the original model and the slightly modified model) exhibits a backward bifurcation re-
vealing that having $R_e$ less than one is not enough to guarantee the elimination of larvae in the plantation.
Eradication occurs when the control is strong enough (i.e. $R_e < R_e(\alpha^s)$) which is not necessarily the best
option for the smallholder farmer, because important investment in plant protection will reduce his net profit.
This bistability phenomenon is dangerous because, depending on the initial condition, the system can converge
either towards the pest-free equilibrium or towards an endemic equilibrium. However, it is possible to avoid
this bi-stability by investing substantially in plant protection ($\alpha \ge 0.09 \text{ day}^{-1}$).

Another approach of modeling pesticide application in the plantation is to use feedback control strategy (as proposed by [43]) or combine prophylaxis and feedback control. These extensions give rise to some optimization challenges which represent an avenue of our further works.

References

[1] Food and Agriculture Organization (FAO), World Food and Agriculture - Statistical yearbook 2021, rome.
https://doi.org/10.4060/cb4477en (2021).

[2] Food and Agriculture Organization (FAO), Agricultural production statistics.2000-2020, fAOSTAT Analytical Brief Series No. 41. Rome. https://www.fao.org/documents/card/en/c/cb9180en (2022).

[3] P. K. Baidoo, J. I. Adam, The Effects of Extracts of Lantana camara (L.) and Azadirachta indica (A. Juss) on the Population Dynamics of Plutella xylostella, Brevicoryne brassicae and Hellula undalis on Cabbage, Sustainable Agriculture Research 1 (2) (2012). doi:10.5539/sar.v1n2p229.

[4] P. K. Baidoo, M. B. Mochiah, Comparing the effectiveness of garlic (Allium sativum L.) and hot pepper (Capsicum frutescens L.) in the management of the major pests of cabbage Brassica oleracea (L.), Sustainable Agriculture Research 5 (2) (2016). doi:10.5539/sar.v5n2p83.

[5] B. Kusznierewicz, A. Bartoszek, L. Wolska, J. Drzewiecki, S. Gorinstein, J. Namieśnik, Partial characterization of white cabbages (brassica oleracea var. capitata f. alba) from different regions by glucosinolates, bioactive compounds, total antioxidant activities and proteins, LWT-Food Science and Technology 41 (1) (2008) 1-9. doi:10.1016/j.lwt.2007.02.007.

[6] A. Steinbrecher, J. Linseisen, Dietary intake of individual glucosinolates in participants of the EPIC-Heidelberg cohort study, Annals of Nutrition and Metabolism 54 (2) (2009) 87-96. doi:10.1159/000209266.
Thus, as the right-hand side is monotonic from $B_y(0)$ to $\frac{r_1 k_1}{\delta_M \alpha_y}$,

$$B_y(t) \leq \max(K_y; B_y(0)) \equiv B_{y\max}, \quad (B.2)$$

where $K_y = \frac{r_1 k_1}{\delta_M \alpha_y}$.

Using the remaining equations of model (3), we prove similarly that,

$$B_a(t) \leq \max \left( \frac{\gamma B_{y\max}}{\alpha_a}; B_a(0) \right) \equiv B_{a\max}, \\ M(t) \leq \max \left( \frac{qh B_{a\max}}{k}; M(0) \right) \equiv M_{\max}, \quad (B.3)$$

$$P(t) \leq \max \left( \frac{\alpha a_p M_{\max}}{\phi_p}; P(0) \right) \equiv P_{\max}.$$

Furthermore,

$$\dot{L} = c \frac{\psi L B_y}{bP + 1} - \lambda \phi P L - \mu_L L^2 - \phi_L L \\ \leq c \psi B_{y\max} L - \mu_L L^2 \quad (B.4)$$

$$L(t) \leq \frac{L(0) \frac{c \psi B_{y\max}}{\mu_L}}{L(0) + \left( \frac{c \psi B_{y\max}}{\mu_L} - L(0) \right) \exp(-c \psi B_{y\max} t)}.$$

Then, we obtain

$$L(t) \leq \max \left( L(0); \frac{c \psi B_{y\max}}{\mu_L} \right) \equiv L_{\max}. \quad (B.5)$$

We assume that the solution $(B_y(t); B_a(t); L(t); M(t))$ starts in the domain $\Omega$ which is described in Theorem 2. we deduce easily that $B_{y\max} = K_y$. Then, it follows from similar developments around relations (B.3) and (B.5) that $\Omega$ is positively invariant.

We focus our attention on the attractiveness of $\Omega$. Let us show that for any solution $(B_y(t); B_a(t); L(t); P(t); M(t))$ started outside of

$$\Omega_{\epsilon} := \left\{ (B_y; B_a; L; P; M) \in \mathbb{R}_+^5 \mid \begin{array}{l} B_y \le K_y + \epsilon, B_a \le K_a + \epsilon, \\ L \le K_L + \epsilon, M \le K_M + \epsilon, P \le K_P + \epsilon \end{array} \right\} \quad (B.6)$$

with $K_y = \frac{r_1 k_1}{\delta_M \alpha_y}$, $K_a = \frac{\gamma (K_y + \epsilon)}{\alpha_a}$, $K_M = \frac{qh(K_a + \epsilon)}{k}$, $K_P = \frac{a_p \alpha (K_M + \epsilon)}{\phi_P}$, $K_L = \frac{c \psi (K_y + \epsilon)}{\mu_L}$, and for any $\epsilon > 0$, enters $\Omega_\epsilon$ in finite time. If we first suppose that $B_y(0) > K_y + \epsilon$, we have

$$\dot{B}_y = r_{0y} B_y - \mu_y B_y^2 + \frac{r_1 k_1 M}{1 + \delta_M M} - \frac{\psi L B_y}{bP + 1} - \gamma B_y \\ \leq \frac{r_1 k_1}{\delta_M} - \alpha_y B_y < \frac{r_1 k_1}{\delta_M} - \alpha_y (K_y + \epsilon) < -\alpha_y \epsilon$$

which implies that

$$B_y(t) < B_y(0) - \alpha_y \epsilon t.$$

Writing the previous inequality at the time

$$t_{\epsilon} = (B_y(0) - K_y - \epsilon) / (\alpha_y \epsilon),$$

we get $B_y(t_\epsilon) < K_y + \epsilon$. Hence, $B_y(t)$ reaches $K_y + \epsilon$ in finite time $t_y < t_\epsilon$ and will stay below this value afterwards.

Similarly, if we assume that $B_a(0) > K_a + \epsilon$ with $B_y(t) \le K_y + \epsilon$ (which may be valid at time 0 or after the previously described convergence time $t_y$), the second equation of model (3) gives

$$\dot{B}_a \leq \gamma (K_y + \epsilon) - \alpha_a B_a < -\alpha_a \epsilon.$$

And $B_a$ will reach $K_a + \epsilon$ in finite $t_a$ and stay below this value afterwards.

Now, if we suppose that $M(0) > K_M + \epsilon$ with $B_a(t) \le K_a + \epsilon$ (which may be valid at time 0 or after the previously described convergence time $t_a$) the dynamic of $M$ in model (3) gives $\dot{M} \le qh(K_a + \epsilon) - kM < -ke$. Then $M$ will reach $K_M + \epsilon$ in finite $t_m$ and stay below this value afterwards. Repeating the same approach for $P$, we can conclude at the existence of convergence time $t_p$ at which $P$ reaches $K_P + \epsilon$ and stays below this value afterwards. Finally, we suppose $L(0) > c\psi(K_y + \epsilon)/\mu_L$ with $B_y(t) \le K_y + \epsilon$. Using the third equation of model (3), we have

$$\dot{L} \leq (c\psi(K_y + \epsilon) - \mu_L L) L < -\mu_L \epsilon L.$$

Hence, $L(t) < L(0) \exp(-\mu_L \epsilon t)$. At the time

$$t_{\epsilon} = \ln \left( \frac{L(0)}{K_L + \epsilon} \right) / (\mu_L \epsilon),$$

we have $L(t_\epsilon) < K_L + \epsilon$. We deduce that $L$ reaches $K_L + \epsilon$ in finite $t_l < t_\epsilon$ and stay below this value afterwards. We conclude that $\Omega_\epsilon$ is attractive for any $\epsilon > 0$. Hence, $\Omega$ is attractive.

We have just shown that the solution $(B_y(t); B_a(t); L(t); P(t); M(t))$, for any initial satisfying (4) is contained in a compact subset of $\mathbb{R}^5$. Using a classical result of dynamical system, we conclude that the solution $(B_y(t); B_a(t); L(t); P(t); M(t))$ is globally defined in $\mathbb{R}_+$.

### C. Proof of Theorem 2

The jacobian matrix of sub-model (6) is given by

$$\left( \begin{array}{ccc} -\alpha_y - 2\mu_y B_y - \psi L & 0 & -\psi B_y \\ \gamma & -\alpha_a - 2\mu_a B_a & 0 \\ c\psi L & 0 & c\psi B_y - 2\mu_L L - \phi_L \\ 0 & qh & 0 \end{array} \frac{r_1 k_1}{(1 + \delta_M M)^2} \right). \quad (C.1)$$

The characteristic polynomial evaluated at the extinction equilibrium $E_0$ is

$$Q_0(\chi) = (\chi + \phi_L) (\chi^3 + A_2 \chi^2 + A_1 \chi + A_0),$$

with

$$A_2 = k + \alpha_y + \alpha_a, \quad A_1 = \alpha_y \alpha_a + k (\alpha_a + \alpha_y),$$
[7] S. Mudzingwa, S. Muzemu, J. Chitamba, Pesticidal efficacy of crude aqueous extracts of Tephrosia vogelii L., Allium sativum L. and Solanum incanum L. in controlling aphids (Brevicoryne brassicae L.) in rape (Brassica napus L.) n controlling aphids (Brevicoryne brassicae L.) in rape (Brassica napus L.), Journal of Research in Agriculture 2 (1) (2013) 157-163.

[8] B. T. Tekle, N. B. Tesfu, Analysis of Head Cabbage Value Chain in Guji Zone, Southern Oromia, Ethiopia, International Journal of Agricultural Economics 8 (1) (2023) 14-26. doi:10.11648/j.ijae.20230801.13.

[9] N. Mpumi, R. S. Machunda, K. M. Mtei, P. A. Ndakidemi, Selected insect pests of economic importance to Brassica oleracea, their control strategies and the potential threat to environmental pollution in Africa, Sustainability 12 (9) (2020) 3824. doi:10.3390/su12093824.

[10] G. Sow, K. Diarra, L. Arvanitakis, D. Bordat, The relationship between the diamondback moth, climatic factors, cabbage crops and natural enemies in a tropical area, Folia Horticulturae 25 (1) (2013) 3-12. doi:10.2478/fhort-2013-0001.

[11] C. Huaripata, G. Sánchez, et al., Life Cycle of the Diamondback Moth Plutella xylostella L. (Lepidoptera: Plutelli-dae), in Broccoli and Cauliflower under Laboratory Conditions, Peruvian Journal of Agronomy 3 (1) (2019) 1-5. doi:10.21704/pja.v3i1.1132.

[12] M. Gautam, H. Singh, S. Kumar, V. Kumar, G. Singh, S. Singh, Diamondback moth, Plutella xylostella (Linnaeus) (Insecta: Lepidoptera: Plutellidae) a major insect of cabbage in India: A review, Journal of Entomology and Zoology Studies 6 (4) (2018) 1394-1399.

[13] R. Nofemela, R. Kfir, The role of parasitoids in suppressing diamondback moth, Plutella xylostella (L.)(Lepidoptera: Plutellidae), populations on unsprayed cabbage in the North West Province of South Africa, African Entomology 13 (1) (2005) 71-83. doi:10520/EJC32624.

[14] M. Sarfraz, A. Keddie, L. M. Dosdall, Biological control of the diamondback moth, Plutella xylostella: A review, Biocontrol Science and Technology 15 (8) (2005) 763-789. doi:10.1080/09583150500136956.

[15] B. W. Amoabeng, G. M. Gurr, C. W. Gitau, H. I. Nicol, L. Munyakazi, P. C. Stevenson, Tri-trophic insecticidal effects of african plants against cabbage pests, PloS one 8 (10) (2013) e78651. doi:10.1371/annotation/f0351003-b6f8-4249-ace5-bcd84dead916.

[16] H. Machekano, B. M. Mvumi, C. Nyamukondiwa, Diamondback moth, Plutella xylostella (L.) in Southern Africa: research trends, challenges and insights on sustainable management options, Sustainability 9 (2) (2017) 91. doi:10.3390/su9020091.

[17] L. K. Agboyi, K. M. Djade, K. M. Ahadji-Dabla, G. K. Ketoh, Y. Nuto, I. A. Glitho, Vegetable production in Togo and potential impact of pesticide use practices on the environment, International Journal of Biological and Chemical Sciences 9 (2) (2015) 723-736. doi:10.4314/ijbcs.v9i2.13.

[18] A. De Oliveira, J. Abreu de Siqueira, J. De Oliveira, J. Da Silva, M. Michereff Filho, Resistance of Brazilian diamondback moth populations to insecticides, Scientia Agricola 68 (2) (2011) 154-159. doi:10.1590/S0103-90162011000200004.

[19] P. G. Mason, T. Dancau, P. K. Abram, C. Noronha, P. L. Dixon, C. K. Parsons, M. Bahar, A. M. Bennett, J. Fernández-Triana, A. M. Brauner, et al., The parasitoid complex of diamondback moth, Plutella xylostella (Linnaeus)(Lepidoptera: Plutellidae), in Canada: impact and status, The Canadian Entomologist 154 (1) (2022) e12. doi:10.4039/tce.2021.51.

[20] J. Timbilla, K. Nyarko, A survey of cabbage production and constraints in Ghana, Ghana Journal of Agricultural Science 37 (1) (2004) 93-101. doi:10.4314/gjas.v37i1.2084.

[21] N. S. Talekar, A. M. Shelton, Biology, Ecology, and Management of the Diamondback moth, Annual Review of Entomology 38 (1) (1993) 275-301. doi:10.1146/annurev.en.38.010193.001423.

[22] M. J. Furlong, D. J. Wright, L. M. Dosdall, Diamondback Moth Ecology and Management: Problems, Progress, and Prospects, Annual Review of Entomology 58 (1) (2013) 517-541, pMID: 23020617. doi:10.1146/annurev-ento-120811-153605.

[23] K.-H. Kim, E. Kabir, S. A. Jahan, Exposure to pesticides and the associated human health effects, Science of the total environment 575 (2017) 525-535. doi:10.1016/j.scitotenv.2016.09.009.

[24] R. Pavela, Efficacy of three newly developed botanical insecticides based on pongam oil against Plutella xylostella L. larvae, Journal of Biopesticides 5 (1) (2012) 62-70.
$$R_B = \frac{k_1 \gamma r_1 q h}{k \alpha_y \alpha_a}, \quad A_0 = k \alpha_a \alpha_y (1 - R_B).$$

We remark that $-\phi_L$ is an eigenvalue of $Q_0$. The other eigenvalues are solutions of the polynomial function $\chi \mapsto \chi^3 + A_2\chi^2 + A_1\chi + A_0$.

If $R_B > 1$, Descartes rule reveals that $Q_0$ has a unique positive root. Using the Hartman-Grobman theorem (Theorem 4.7 in [7]), we conclude that $E_0$ is unstable. Assuming $R_B < 1$. We use the well-known Routh-Hurwitz stability's criterion for stability by verifying conditions: $A_2, A_1, A_0 > 0$ and $A_2A_1 - A_0 > 0$.

The positivity of $A_2, A_1$ and $A_0$ are obvious. A direct computation gives

$$A_2 A_1 - A_0 = (k + \alpha_y + \alpha_a) (\alpha_y \alpha_a + k (\alpha_a + \alpha_y)) - k \alpha_a \alpha_y (1 - R_B) = k^2 (\alpha_y + \alpha_a) + (\alpha_y + \alpha_a) (\alpha_y \alpha_a + k (\alpha_a + \alpha_y)) + r_1 k_1 q h \gamma > 0.$$

We conclude that $E_0$ is locally asymptotically stable.

We assume that $R_B > 1$. The characteristic polynomial evaluated at the PFE $E_1$ is given by

$$Q_1(\chi) = \left(\chi - \left(c\psi B_{y1} - \phi_L\right)\right) \left(\chi^3 + B_2\chi^2 + B_1\chi + B_0\right)$$

where

$$\begin{align*} B_2 &= k + \alpha_y + \alpha_a + 2\mu_a B_{a1} + 2\mu_y B_{y1} \\ B_1 &= (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) \\ &\quad + k (\alpha_a + \alpha_y + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \\ B_0 &= k (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) - \frac{k_1 q h \gamma r_1}{(1 + \delta_M M_1)^2}. \end{align*}$$

It is clear that $B_2, B_1 > 0$. Using Equation (11) and the fact that $1 + \delta_M M_1 > 1$, we have

$$\begin{align*} B_0 &> \frac{k}{1 + \delta_M M_1} \left( (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) \right. \\ &\qquad \left. - (\alpha_y (\alpha_a + \mu_a B_{a1}) + \mu (\alpha_a + \mu_a B_{a1})^2 B_{a1}) \right) \\ &> \frac{k}{1 + \delta_M M_1} \left( (\alpha_y + \mu (\alpha_a + \mu_a B_{a1}) B_a) (\alpha_a + \mu_a B_{a1}) \right. \\ &\qquad \left. - (\alpha_y (\alpha_a + \mu_a B_{a1}) + \mu (\alpha_a + \mu_a B_{a1})^2 B_{a1}) \right) = 0. \end{align*}$$

Moreover,

$$\begin{align*} B_2 B_1 - B_0 &= k^2 (\alpha_a + \alpha_y + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \\ &\quad + (\alpha_y + \alpha_a + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \\ &\quad \left( (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) + k (\alpha_a + \alpha_y + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \right) > 0. \end{align*}$$

According to Routh-Hurwitz stability's criterion, all the eigenvalues of the polynomial function $\chi \mapsto Q_1(\chi)$ have negative real parts. Thus, the local stability of $E_1$ depends only of the sign of $c\psi B_{y1} - \phi_L$. By noticing that $c\psi B_{y1} - \phi_L = \phi_L(R_0 - 1)$, it comes from Hartman-Grobman theorem that $E_1$ is locally asymptotically stable provided $R_0 < 1$ and unstable if $R_0 > 1$.

Now supposing $R_B > 1$ and $R_0 > 1$. The characteristic polynomial evaluated at $E^*$ is given by

$$Q_2(\chi) = \chi^4 + C_3 \chi^3 + C_2 \chi^2 + C_1 \chi + C_0 \quad (C.2)$$

where

$$\begin{align*} J_{11} &= \alpha_y + 2\mu_y B_y^* + \psi L^* \\ J_{22} &= \alpha_a + 2\mu_a B_a^* \\ C_3 &= \mu_L L^* + J_{11} + k + J_{22} \\ C_2 &= c\psi^2 L^* B_y^* + \mu_L L^* J_{11} + k J_{22} \\ &\quad + (k + J_{22}) (\mu_L L^* + J_{11}) \\ C_1 &= (k + J_{22}) (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) \\ &\quad + k J_{22} (\mu_L L^* + J_{11}) - \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\ C_0 &= k J_{22} (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) - \mu_L L^* \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2}. \end{align*}$$

We want to use again the Routh-Hurwitz stability's criterion. It consists in verifying the conditions: $C_3, C_2, C_1, C_0 > 0, C_3C_2 - C_0 > 0$ and $(C_3C_2 - C_1)C_1 - C_3^2C_0 > 0$.

It is clear that $C_3, C_2 > 0$. In addition, using Equation (12)

$$\begin{align*} C_1 &= (k + J_{22}) (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) \\ &\quad + k J_{22} (\mu_L L^* + J_{11}) - \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\ &> k J_{11} J_{22} - \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\ &> \frac{k}{1 + \delta_M M^*} \left( \left( \alpha_y + \frac{\mu_y}{\gamma} (\alpha_a + \mu_a B_a^*) B_a^* + \psi L^* \right) (\alpha_a + \mu_a B_a^*) - (\alpha_y + \psi L^*) (\alpha_a + \mu_a B_a^*) \right) \\ &\quad - \frac{\mu_y}{\gamma} (\alpha_a + \mu_a B_a^*)^2 B_a^* \right) = 0 \end{align*}$$

and

$$\begin{align*} C_0 &= k J_{22} (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) - \mu_L L^* \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\ &> \mu_L L^* \left( k J_{11} J_{22} - \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2} \right) > 0. \end{align*}$$

Furthermore,

$$\begin{align*} C_3 C_2 - C_1 &= (k + J_{22}) ((k + J_{22}) (\mu_y L^* + J_{11}) + k J_{22}) \\ &\quad + (\mu_y L^* + J_{11}) (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) \\ &\quad + (k_1 + J_{22}) (\mu_L L^* + J_{11}) + \frac{q h \gamma r_1 k_1}{(1 + \delta_M M^*)^2} > 0 \end{align*}$$

and,

$$(C_3 C_2 - C_1) C_1 - C_3^2 C_0 > (k^2 + J_{22}^2) (\mu_L^* + J_{11}) (k + J_{22}) \\ + (k J_{22} + k^2 + J_{22}^2) (\mu_L L^* + J_{11})^2 > 0.$$

We conclude by using Routh-Hurwitz stability criterion that the endemic equilibrium $E^*$, whenever it exists, is always locally asymptotically stable.
[25] M. M. Waiganjo, C. N. Waturu, J. M. Mureithi, J. Muriuki, J. Kamau, R. Munene, Use Of Entomopathogenic Fungi And Neem Bio-Pesticides For Brassica Pest Control And Conservation Of Their Natural Enemies, East African Agricultural and Forestry Journal 1& 2 (77) (2011).

[26] S. Kunbhar, L. B. Rajput, A. A. Gilal, G. A. Channa, J. G. M. Sahito, Impact of botanical pesticides against sucking insect pests and their insect predators in brinjal crop, J. Entomol. Zool. Stud 6 (2) (2018) 83-87.

[27] A. Nazeer, A. Mukhtar, S. Muhammad, U. Hidayat, I. Toheed, A. A.-M. Khalid, S. Kiran, U. Rafi, A. Saeed, A. A. H. A. Nibal, F. K. Hanem, S. Muhammad, *Botanical Insecticides Are a Non-Toxic Alternative to Conventional Pesticides in the Control of Insects and Pests*, in: H. A. F. El-Shafie (Ed.), Global Decline of Insects, IntechOpen, Rijeka, 2021, Ch. 6. doi:10.5772/intechopen.100416.
URL https://doi.org/10.5772/intechopen.100416

[28] M. B. Isman, Botanical insecticides, deterrents, and repellents in modern agriculture and an increasingly regulated world, Annual Review of Entomology 51 (2006) 45-66. doi:10.1146/annuev.ento.51.110104.151146.

[29] D. S. Charleston, R. Kfir, M. Dicke, L. E. M. Vet, Impact of botanical pesticides derived from Melia azedarach and Azadirachta indica on the biology of two parasitoid species of the diamondback moth, Biological Control 33 (2) (2005) 131-142. doi:10.1016/j.biocontrol.2005.02.007.

[30] B. W. Amoabeng, G. M. Gurr, C. W. Gitau, P. C. Stevenson, Cost:benefit analysis of botanical insecticide use in cabbage: implications for smallholder farmers in developing countries, Crop Protection 57 (2014) 71-76. doi:10.1016/j.cropro.2013.11.019.

[31] P. Mkenda, R. Mwanauta, P. C. Stevenson, P. Ndakidemi, K. Mtei, S. R. Belmain, Extracts from Field Margin Weeds Provide Economically Viable and Environmentally Benign Pest Control Compared to Synthetic Pesticides, PloS one 10 (11) (2015) e0143530. doi:10.1371/journal.pone.0143530.

[32] H. Gordon, Economic Theory of a Common-Property Resource: the Fishery, Journal of Political Economy 62 (2) (1954) 124-142. doi:10.1086/257497.

[33] J. Jiang, X. Li, X. Wu, The Dynamics of a Bioeconomic Model with Michaelis-Menten Type Prey Harvesting, Bull. Malays. Math. Sci. Soc. 46 (2) (2023) 57. doi:10.1007/s40840-022-01452-4.

[34] V. L. Smith, On Models of Commercial Fishing, Journal of Political Economy 77 (2) (1969) 181-198.
URL http://www.jstor.org/stable/1829764

[35] A. K. Misra, R. K. Rai, Y. Takeuchi, Modeling the effect of time delay in budget allocation to control an epidemic through awareness, International Journal of Biomathematics 11 (02) (2018) 1850027. doi:10.1142/S1793524518500274.

[36] R. K. Rai, A. K. Misra, Y. Takeuchi, Modeling the impact of sanitation and awareness on the spread of infectious diseases, Mathematical Biosciences and Engineering 16 (2) (2019) 667-700. doi:10.3934/mbe.2019032.

[37] G. Kome, R. Enang, B. Yerima, Knowledge and management of soil fertility by farmers in western Cameroon, Geoderma Regional 13 (2018) 43-51. doi:10.1016/j.geodrs.2018.02.001.

[38] S. Tetsopgang, F. Fonyuy, Enhancing growth quality and yield of cabbage (Brassica oleracea) while increasing soil pH, chemicals and organic carbon with the application of fines from volcanic pyroclastic materials on a tropical soil in wum, Northwest Cameroon, Africa, Scientific African 6 (2019) e00199. doi:10.1016/j.sciaf.2019.e00199.

[39] D. G. Luenberger, Introduction to dynamic systems: theory, models, and applications, Vol. 1, Wiley New York, 1979.

[40] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences 180 (1-2) (2002) 29-48. doi:10.1016/S0025-5564(02)00108-6.

[41] L. jardins de la Vallée de la Siagne, *Jeunes plants maraichers bio*, Website (2022).
URL https://www.associationjvs.com

[42] I. Guvenc, E. Yildirim, Increasing productivity with intercropping systems in cabbage production, Journal of Sustainable Agriculture 28 (4) (2006) 29-44. doi:10.1300/J064v28n04_04.
## D. Proof of Theorem 4

We assume that the solution $(B_y(t); B_a(t); L(t); M(t))$ starts in the domain $\Omega'$ which is described in Theorem 4.

If $B_a = B_{a_1}$, we have $\dot{B}_a = \gamma B_y - \mu_a B_{a_1}^2 - \alpha_a B_{a_1} = \gamma (B_y - B_{y_1}) \le 0$. then $B_a(t) \le B_{a_1}$ for all $t \ge t_0$. Similarly, if $M = M_1$, we obtain $\dot{M} = qhB_a - kM_1 = qh(B_a - B_{a_1})$. Hence, $M(t) \le M_1$ for all $t \ge t_0$. Finally, if $B_y = B_{y_1}$, we obtain

$$ 
\begin{aligned}
\frac{dB_y}{dt} &= r_{0y}B_{y_1} - \mu_y B_{y_1}^2 - \psi LB_{y_1} - \gamma B_{y_1} + r_1 \frac{k_1 M}{1 + \delta_M M} \\
&\le r_{0y}B_{y_1} - \mu_y B_{y_1}^2 - \gamma B_y + r_1 \frac{k_1 M}{1 + \delta_M M} \\
&\le r_1 k_1 \left( -\frac{M_1}{1 + \delta_M M_1} + \frac{M}{1 + \delta_M M} \right) \\
&\le r_1 k_1 \frac{M - M_1}{(1 + \delta_M M)(1 + \delta_M M_1)} \\
&\le 0.
\end{aligned}
$$

Then $B_y(t) \le B_{y_1}$. We conclude that $\Omega'$ is positively invariant.

Before dealing with the stability in $\Omega'$, we first prove the following half-way result.

**Theorem D.1.** The pest-free equilibrium $E_1$ is globally asymptotically stable in the sub-region $\mathbb{R}_+^3 \setminus \{0\}$ of hyper-plan $L = 0$ if and only if $1 < \mathcal{R}_B$.

**PROOF.** This proof is adapted from a result given in [42].

We consider

$$X_1 = g(X_1) \qquad (D.1)$$

with

$$X_1 = \begin{pmatrix} B_y \\ B_a \\ M \end{pmatrix}, \quad \text{and} \quad g(X_1) = \begin{pmatrix} -\mu_y B_y^2 - \alpha_y B_y + r_1 \frac{k_1 M}{1+\delta_M M} \\ \gamma B_y - \mu_a B_a^2 - \alpha_a B_a \\ qhB_a - kM \end{pmatrix}.$$

System (D.1) has two equilibria: 0 and $X_1^E = (B_{y_1}; B_{a_1}; M_1)$ on $\mathbb{R}_+^3$ if and only if $1 < \mathcal{R}_B$.

$\mathbb{R}_+^3$ is positively invariant. Furthermore, the Jacobian matrix of $g$ evaluated at $X_1$ is given by

$$ \mathcal{J}_g(X_1) = \begin{pmatrix} -\alpha_y - 2\mu_y B_y & 0 & \frac{r_1 k_1}{(1+\delta_M M)^2} \\ \gamma & -\alpha_a - 2\mu_a B_a & 0 \\ 0 & qh & -k \end{pmatrix}. \qquad (D.2) $$

$\mathcal{J}_g(X_1)$ is a cooperative irreducible matrix. Then, the flow $\phi_t$ of system (D.1) is strongly monotone in $\mathbb{R}_+^3$. In addition, for every $\epsilon > 0$, the set $Z = [0; z_\epsilon]$ with $z_\epsilon = (K_y + \epsilon; K_a + \epsilon; K_M + \epsilon)^T$, is positively invariant and attractive for the flow $\phi_t$ as $Z$ is the projection of $\Omega_\epsilon$ on this subspace (see (B.6)).

According to Perron-Frobenius theorem, the stability modulus $s$ of $\mathcal{J}_g(0)$, is a positive real eigenvalue of $\mathcal{J}_g(0)$ and there is a positive vector $v$ such that $\mathcal{J}_g(0)v = sv$.

Let $x$ be a positive vector of $\mathbb{R}_+^3$. Using attractiveness property of $Z$, we have $y = \phi_{t_1}(x) \in Z$ for large $t_1$.

Then, there is $\eta := \eta(y) > 0$ such that $\phi_t(y) \in [\eta v; z_\epsilon]$ which is a positive invariant order preserving interval and we have $g(z_\epsilon) < 0$ and $g(\eta v) \approx \mathcal{J}_g(0)\eta v = s\eta v > 0$. Using Proposition 3.2.1 in [42], $\phi_t(z_\epsilon)$ and $\phi_t(\eta v)$ converge monotonically both to $X_1^E$.

By monotonicity,

$$\eta v < \phi_t(\eta v) < \phi_t(y) < \phi_t(z_\epsilon) < z_\epsilon, \quad t > 0.$$

Letting $t \to +\infty$ leads to the $w$-limit set of $y$ is $\{X_1^E\}$. Then $\phi_t(x)$ converges to $X_1^E$ for all positive vector $x$ of $\mathbb{R}_+^3$.

Furthermore, $B_y$ can not stay equal to 0. Indeed, if $B_y = 0$ and $M > 0$ at a time $t_0$, we have $\dot{B}_y = r_1 k_1 M / (1 + \delta_M M) > 0$. In addition, if $B_y = M = 0$ and $B_a = 0$ at a time $t_0$, one has $\dot{M} = qhB_a > 0$ so that $M$ instantaneously becomes positive and it follows that $\dot{B}_y = r_1 k_1 M / (1 + \delta_M M) > 0$. One can reproduce the same approach for variables $B_a$ and $M$. Thus, except if $x(t_0) = 0$, $x$ instantaneously enters the positive orthant $\mathbb{R}_+^3$ and is amenable to the previous developments.

Thus, $X_1^E$ is asymptotically stable everywhere in $\mathbb{R}_+^3 \setminus \{0\}$.

For convenience, we set

$$\dot{X}_1 = F(X_1; L)$$

where

$$X_1 = \begin{pmatrix} B_y \\ B_a \\ M \end{pmatrix}, \quad \text{and}$$

$$F(X_1; L) = \begin{pmatrix} r_{0y} B_y - \mu_y B_y^2 - \psi L B_y - \gamma B_y + r_1 \frac{k_1 M}{1+\delta_M M} \\ \gamma B_y + r_{0a} B_a - \mu_a B_a^2 - h B_a \\ qh B_a - k M \end{pmatrix}.$$

We remark that

$$ 
\begin{aligned}
\dot{L} &= c\psi L B_y - \mu_L L^2 - \phi_L L \\
&= AL - \hat{G}(X_1; L)
\end{aligned}
$$

where, for all variables in $\Omega'$, one has

$$A = \phi_L (\mathcal{R}_0 - 1)$$

$$\text{and } \hat{G}(X_1; L) = c\psi(B_{y_1} - B_y)L + \mu_L L^2 \ge 0.$$

If $\mathcal{R}_0 \le 1$, one has $A \le 0$ and then $\dot{L} \le -\mu_L L^2$. Thus, $\lim_{t \to +\infty} L(t) = 0$.

According to Theorem D.1, $X_1^E$ is a globally asymptotically stable equilibrium of $\dot{X}_1 = F(X_1; 0)$, a limiting system of $\dot{X}_1 = F(X_1; L)$. Thus, $\lim_{t \to +\infty} X_1(t) = X_1^E$ (see Castillo-Chavez et al. [8]).
[43] A. Misra, N. Jha, R. Patel, Modeling the effects of insects and insecticides on agricultural crops with NSFD method, Journal of Applied Mathematics and Computing 63 (2020) 197–215. doi:10.1007/s12190-019-01314-6.

[44] J. Andaloro, K. Rose, A. Shelton, C. Hoy, R. Becker, Cabbage growth stages, New York's Food and Life Sciences Bulletin (101) (1983).

[45] S. Prasad, et al., Quantitative estimation of damage to crucifers caused by cabbageworm, cabbage looper, diatnonback moth and cabbage aphid., Indian Journal of Entomology 25 (3) (1963).

[46] L. Barreira, C. Valls, Ordinary differential equations: Qualitative theory, Vol. 137, American Mathematical Society, 2012.

[47] H. L. Smith, Irreducible Cooperative Systems, in: Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems, Vol. 41 of 0-8218-0393-X, American mathematical Society, Providence, RI, 1995, pp. 55–74.

[48] C. Castillo-Chavez, Z. Feng, W. Huang, On the computation of R0 and its role on global stability, Mathematical Approaches for Emerging and Re-emerging Infectious Diseases: An introduction 125 (2002) 31–65.

[49] Z. Shuai, P. Van den Driessche, Global stability of infectious disease models using lyapunov functions, SIAM Journal on Applied Mathematics 73 (4) (2013) 1513–1532.
URL https://www.jstor.org/stable/24510691

A Parameter estimation

Some parameters were estimated from the literature. The values of the conversion rate *c* and uptake rate of botanical pesticide by larvae were taken in [43]. [30] reported that undamaged heads fetched 0.56 USD and 0.83 USD per kg for the major and minor seasons respectively. We chose an average cabbage price *q* = 0.69 USD per kg [30].

When replanted, the cabbage consists of 3 to 5 leaves [44, 15]. The weight at transplanting is therefore assumed to be 70g. We assume a growth period for the cabbage of 90 days (three months) for an average weight of 0.7 kg per head [15]. [10] report that there is a strong correlation between the age of the plant and the age of the diamondback moth larvae. In fact, plants less than 55 days old are favored by immature larvae (stages 2 and 3). Knowing also that larvae populations are more voracious at stage 3, we set the aging rate of young biomass $\gamma = 1/55 \approx 0.02 \text{ day}^{-1}$ and the harvest rate $h = 1/(90 - 55) \approx 0.03 \text{ day}^{-1}$.

One hectare can support around 35000 cabbage plants [30]. Then, the carrying capacity of a 1-hectare field is 35,000 × 0.7 = 24,500 kg, split into $K_a = 24,500 \times 35/90 = 9,528$ kg of adult biomass and $K_y = 24,500 \times 55/90 = 14,972$ kg of young biomass. Assuming that the weight of a young cabbage at 55 days after replanting is 0.4 kg, we compute the growth rate of adult biomass $r_{0a} = \ln(0.7/0.4)/35 \approx 1.6 \times 10^{-2} \text{ day}^{-1}$. Moreover, we assume that $r_{0y} \approx r_{0a} = 1.6 \times 10^{-2} \text{ day}^{-1}$. Thus, the competition rates are given by $\mu_a = 1.6 \times 10^{-2}/9528 \approx 1.68 \times 10^{-6} \text{ kg}^{-1}\text{day}^{-1}$, and, $\mu_y = 1.6 \times 10^{-2}/14972 \approx 1.07 \times 10^{-6} \text{ kg}^{-1}\text{day}^{-1}$.

According to [12], the infestation is noticeable at more than 20 larvae per plant. In addition, [45] reported that an increase of 5 larvae leads to an estimated increase in consumption of 19 g and 32 g for an exposure of 4 and 7 weeks respectively. Hence, we assume that the biomass consumption rate by larvae $\psi$ has values between $0.032/((20+5) \times 7 \times 7) \approx 2.27 \times 10^{-5}$ and $0.019/((20+5) \times 7 \times 4) \approx 4.57 \times 10^{-5}$ larvae$^{-1}$.day$^{-1}$. The larvae development stage lasts about 9 to 32 days [12, 11]. Hence, we assume that the mortality rate of larvae $\phi_L$ is between 0.03 and $0.1 \text{ day}^{-1}$. Assuming a carrying capacity of 30 larvae per plant [12, 13], the carrying capacity of the plantation is 35,000 × 30 = 1,050,000 larvae. The competition of larvae $\mu_L$ is given by formula $c\psi K_y/1,050,000$. Thus, the range of $\mu_L$ is $[1.94 \times 10^{-6}; 3.91 \times 10^{-6}]$.

According to [30], a hectare is treated between 7 and 8 times over a growing season. It is hence assumed that the natural depletion of botanical pesticide is $\phi_p = 8/90 \approx 0.1 \text{ day}^{-1}$. In addition, around 500 liters of pesticide per hectare are used [24] for a total cost around 306 USD [30]. Thus, the inverse of the pesticide price is $a_p \approx 500/231 \approx 2.16 \text{ l.USD}^{-1}$.

According to [41], the unit price of a young cabbage plant is around 0.16 USD per plant, which weighs 70g. Thus, the inverse of price of young biomass is estimated at $r_1 = 0.07/0.16 \approx 0.46 \text{ kg.USD}^{-1}$.
<figure><img src="image_20.png" /><figcaption>**Figure E.1:** Dynamics of state variables (a) $B_y$, (b) $L$, (c) $P$ and (d) $M$ of sub-model (20). The initial conditions and parameter values are the same as in Fig. 12 except for $\alpha = 0.09 \text{ day}^{-1}$ leading to an effective reproductive number $R_e = 0.19 < R_e(\alpha^s) \approx 0.2 < 1$. For all initial conditions, the system converges to the pest-free equilibrium.</figcaption></figure>

In addition, Theorem 2 ensures that $E_1$ is unstable whenever $R_0 > 1$. Thus, the PFE $E_1$ is asymptotically stable in $\Omega'$ if and only if $R_0 \le 1 < R_B$.

We comment that the uniform persistence of sub-model (6) can be obtained in $\Omega'$ when $R_0 \le 1 < R_B$ [40].

## E. Other simulations

In this appendix, we simulate sub-model (20) and model (3) for different values of parameter $\alpha$, which represents the investment rate in plant protection. All remaining parameters are given in Table 1.

Figures Fig. E.1 and Fig. E.2 are similar to Fig. 12 (in which $\alpha = 1.35 \times 10^{-2} \text{ day}^{-1}$) and depict the dynamics of sub-model (20). Fig. E.1 presents a situation where the pest-free equilibrium is always stable and the system does not have an endemic equilibrium. It corresponds to a high investment in plant protection ($\alpha = 0.09 \text{ day}^{-1}$), associated with a cost of 1,530 USD computed according to (38). It generates a seasonal yield of 22.2 tonnes for an income of 15,320 USD, calculated from (19) at the pest-free equilibrium.

Fig. E.2 illustrates the situation where there is a unique endemic stable equilibrium and an unstable pest-free equilibrium. It corresponds to a low investment in plant protection ($\alpha = 0.01 \text{ day}^{-1}$), associated with a cost of 79 USD, and it generates a seasonal yield of 10.3 tonnes for an income of 7,111 USD at the endemic equilibrium.

The effective reproductive number of original model (3) is given by

$$R_e' = \frac{c\psi\alpha_a B_{a1}}{\gamma(\lambda\phi P_1 + \phi_L)(1 + bP_1)},$$

with $P_1 = a_p\alpha M_1/\phi_P$. The equivalence of the threshold $\alpha^0$ for the original model (3), obtained by solving $R_e'(\alpha) = 1$, is given by $\alpha^{0'} = \phi_P\Delta/(2a_P M_1\lambda b\phi)$, with $\Delta = - (\lambda\phi + b\phi_L) + \sqrt{(\lambda\phi + b\phi_L)^2 + 4\lambda b\phi (c\psi B_{y1} - \phi_L)}$.

Using parameter values in Table 1, we have $\alpha^{0'} \approx 7.34 \times 10^{-3} \text{ day}^{-1}$. The lower bound for the existence of endemic equilibria $\alpha^s \approx 0.022 \text{ day}^{-1}$. Fig. E.3 shows model (3), for different values of parameter $\alpha$. Comparing these dynamics with those of sub-model (20) depicted in panel (a) of Fig. 12, E.1 and E.2, we observe that both models have similar dynamics. In fact, we observe a bistability
B Proof of Theorem 1

Let $(B_y(t); B_a(t); L(t); P(t); M(t))$ be a local solution of model (3) defined on $(0; T_+)$ where $T_+ > 0$. Using the first equation of model (3), we have

$$
\begin{aligned}
\dot{B}_y &= r_{0y}B_y - \mu_y B_y^2 + \frac{r_1 k_1 M}{1 + \delta_M M} - \frac{\psi L B_y}{bP + 1} - \gamma B_y \\
&\leq \frac{r_1 k_1}{\delta_M} - \alpha_y B_y.
\end{aligned}
$$

Applying Gronwall’s lemma [46], we obtain

$$
B_y(t) \le \frac{r_1 k_1}{\delta_M \alpha_y} + \left( B_y(0) - \frac{r_1 k_1}{\delta_M \alpha_y} \right) \exp(-\alpha_y t). \quad (39)
$$

Thus, as the right-hand side is monotic from $B_y(0)$ to $\frac{r_1 k_1}{\delta_M \alpha_y}$,

$$
B_y(t) \le \max(K_y; B_y(0)) \equiv B_{y,\max}, \quad (40)
$$

where $K_y = \frac{r_1 k_1}{\delta_M \alpha_y}$.

Using the remaining equations of model (3), we prove similarly that,

$$
\begin{align}
B_a(t) &\le \max\left(\frac{\gamma B_{y_{\max}}}{\alpha_a}; B_a(0)\right) \equiv B_{a_{\max}}, \nonumber \\
M(t) &\le \max\left(\frac{qhB_{a_{\max}}}{k}; M(0)\right) \equiv M_{\max}, \tag{41} \\
P(t) &\le \max\left(\frac{\alpha a_p M_{\max}}{\phi_P}; P(0)\right) \equiv P_{\max}. \nonumber
\end{align}
$$

Furthermore,

$$
\begin{align}
\dot{L} &= c \frac{\psi L B_y}{bP + 1} - \lambda \phi P L - \mu_L L^2 - \phi_L L \nonumber \\
&\le c \psi B_{y_{\max}} L - \mu_L L^2 \tag{42} \\
L(t) &\le \frac{L(0) \frac{c \psi B_{y_{\max}}}{\mu_L}}{L(0) + \left( \frac{c \psi B_{y_{\max}}}{\mu_L} - L(0) \right) \exp(-c \psi B_{y_{\max}} t)}. \nonumber
\end{align}
$$

Then, we obtain

$$
L(t) \le \max \left( L(0); \frac{c \psi B_{y_{\max}}}{\mu_L} \right) \equiv L_{\max}. \quad (43)
$$

We assume that the solution $(B_y(t); B_a(t); L(t); M(t))$ starts in the domain $\Omega$ which is described in Theorem 2. we deduce easily that $B_{y_{\max}} = K_y$. Then, it follows from similar developments around relations (41) and (43) that $\Omega$ is positively invariant.

We focus our attention on the attractiveness of $\Omega$. Let us show that for any solution $(B_y(t); B_a(t); L(t); P(t); M(t))$ started outside of

$$
\Omega_{\epsilon} := \left\{ (B_y; B_a; L; P; M) \in \mathbb{R}_+^5 \mid \begin{array}{l} B_y \le K_y + \epsilon, B_a \le K_a + \epsilon, \\ L \le K_L + \epsilon, M \le K_M + \epsilon, P \le K_P + \epsilon \end{array} \right\} \quad (44)
$$

with $K_y = \frac{r_1 k_1}{\delta_M \alpha_y}$, $K_a = \frac{\gamma(K_y + \epsilon)}{\alpha_a}$, $K_M = \frac{qh(K_a + \epsilon)}{k}$, $K_P = \frac{a_p \alpha(K_M + \epsilon)}{\phi_P}$, $K_L = \frac{c \psi(K_y + \epsilon)}{\mu_L}$, and for any $\epsilon > 0$, enters $\Omega_\epsilon$ in finite time. If we first suppose that $B_y(0) > K_y + \epsilon$, we have

$$
\begin{aligned}
\dot{B}_y &= r_{0y}B_y - \mu_y B_y^2 + \frac{r_1 k_1 M}{1 + \delta_M M} - \frac{\psi L B_y}{bP + 1} - \gamma B_y \\
&\leq \frac{r_1 k_1}{\delta_M} - \alpha_y B_y < \frac{r_1 k_1}{\delta_M} - \alpha_y (K_y + \epsilon) < -\alpha_y \epsilon
\end{aligned}
$$

which implies that

$$
B_y(t) < B_y(0) - \alpha_y \epsilon t.
$$

Writing the previous inequality at the time

$$
t_{\epsilon} = (B_{y}(0) - K_{y} - \epsilon) / (\alpha_{y}\epsilon),
$$
<figure><img src="image_21.png" /></figure>

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

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

<figure><img src="image_24.png" /><figcaption>**Figure E.2:** Dynamics of state variables (a) $B_y$, (b) $L$, (c) $P$ and (d) $M$ of sub-model (20). The initial conditions and parameter values are the same as in Fig. 12 except for $\alpha = 0.01$ day$^{-1}$ leading to an effective reproductive number $R_e = 1.02 > 1$. For all initial conditions, the system converges towards the endemic equilibrium.</figcaption></figure>

for $\alpha \in (7.35 \times 10^{-3}; 0.022)$ day$^{-1}$ (see Fig. E.3(b) and Fig. E.3(c)). When $\alpha \ge 0.022$ day$^{-1}$, the solutions tend to the pest-free equilibrium of (3) (see Fig. E.3(d)) with a seasonal yield around 11.1 tonnes, and when $\alpha \le 7.34 \times 10^{-3}$ day$^{-1}$, they converge towards the endemic equilibrium of (3) (see Fig. E.3(a)) with a seasonal yield around 7.6 tonnes. Furthermore, $\alpha^0 \gg \alpha^{0'}$ means that a pesticide with a double action is more effective and costless than a classical one. Such results encourage vulgarizing the use of botanical pesticides.

## CRediT authorship contribution statement

**Aurelien Kambeu Youmbi:** Conceptualization, Modeling, Analysis, Simulation, Writing—Original Draft. **Suzanne Touzeau:** Modeling, Analysis, Simulation, Writing—Review & editing, Funding acquisition. **Frédéric Grognard:** Modeling, Analysis, Simulation, Writing—Review & editing. **Berge Tsanou:** Conceptualization, Modeling, Analysis, Simulation, Writing—Review & editing.

## Acknowledgments

This work was carried out within EPITAG (https://team.inria.fr/epitag/), an Inria Associate team part of the

LIRIMA, which is supported by INRAE via the TSARA program and by the French Embassy in Cameroon.

## Declaration of competing interest

The authors declare that they have no competing interests that are directly or indirectly related to this work.

## References

[1] Agboyi, L.K., Djade, K.M., Ahadji-Dabla, K.M., Ketoh, G.K., Nuto, Y., Glitho, I.A., 2015. Vegetable production in Togo and potential impact of pesticide use practices on the environment. International Journal of Biological and Chemical Sciences 9, 723–736. doi:10.4314/ijbcs.v9i2.13.

[2] Amoabeng, B.W., Gurr, G.M., Gitau, C.W., Nicol, H.I., Munyakazi, L., Stevenson, P.C., 2013. Tri-trophic insecticidal effects of african plants against cabbage pests. PLOS ONE 8, e78651. doi:10.1371/annotation/f0351003-b6f8-4249-ace5-bcd84dead916.

[3] Amoabeng, B.W., Gurr, G.M., Gitau, C.W., Stevenson, P.C., 2014. Cost:benefit analysis of botanical insecticide use in cabbage: implications for smallholder farmers in developing countries. Crop Protection 57, 71–76. doi:10.1016/j.cropro.2013.11.019.

[4] Andaloro, J., Rose, K.B., Shelton, A.M., Hoy, C.W., Becker, R.F., 1983. Cabbage growth stages. New York’s Food and Life Sciences Bulletin.

[5] Baidoo, P.K., Adam, J.I., 2012. The Effects of Extracts of Lantana camara (L.) and Azadirachta indica (A. Juss) on the Population
we get $B_y(t_\epsilon) < K_y + \epsilon$. Hence, $B_y(t)$ reaches $K_y + \epsilon$ in finite time $t_y < t_\epsilon$ and will stay below this value afterwards.

Similarly, if we assume that $B_a(0) > K_a + \epsilon$ with $B_y(t) \le K_y + \epsilon$ (which may be valid at time 0 or after the previously described convergence time $t_y$), the second equation of model (3) gives

$$ \dot{B}_a \le \gamma (K_y + \epsilon) - \alpha_a B_a < -\alpha_a \epsilon. $$

And $B_a$ will reach $K_a + \epsilon$ in finite $t_a$ and stay below this value afterwards.

Now, if we suppose that $M(0) > K_M + \epsilon$ with $B_a(t) \le K_a + \epsilon$ (which may be valid at time 0 or after the previously described convergence time $t_a$) the dynamic of $M$ in model (3) gives $\dot{M} \le qh(K_a + \epsilon) - kM < -k\epsilon$. Then $M$ will reach $K_M + \epsilon$ in finite $t_m$ and stay below this value afterwards. Repeating the same approach for $P$, we can conclude at the existence of convergence time $t_p$ at which $P$ reaches $K_P + \epsilon$ and stays below this value afterwards. Finally, we suppose $L(0) > c\psi(K_y + \epsilon)/\mu_L$ with $B_y(t) \le K_y + \epsilon$. Using the third equation of model (3), we have

$$ \dot{L} \le (c\psi(K_y + \epsilon) - \mu_L L) L < -\mu_L \epsilon L. $$

Hence, $L(t) < L(0) \exp(-\mu_L \epsilon t)$. At the time

$$ t_{\epsilon} = \ln(L(0)/(K_L + \epsilon)) / (\mu_L \epsilon), $$

we have $L(t_\epsilon) < K_L + \epsilon$. We deduce that $L$ reaches $K_L + \epsilon$ in finite $t_l < t_\epsilon$ and stay below this value afterwards. We conclude that $\Omega_\epsilon$ is attractive for any $\epsilon > 0$. Hence, $\Omega$ is attractive.

We have just shown that the solution $(B_y(t); B_a(t); L(t); P(t); M(t))$, for any initial satisfying (4) is contained in a compact subset of $\mathbb{R}^5$. Using a classical result of dynamical system, we conclude that the solution $(B_y(t); B_a(t); L(t); P(t); M(t))$ is globally defined in $\mathbb{R}_+$.

## C Proof of Theorem 2

The jacobian matrix of sub-model (6) is given by

$$ \begin{pmatrix} -\alpha_y - 2\mu_y B_y - \psi L & 0 & -\psi B_y & \frac{r_1 k_1}{(1 + \delta_M M)^2} \\ \gamma & -\alpha_a - 2\mu_a B_a & 0 & 0 \\ c\psi L & 0 & c\psi B_y - 2\mu_L L - \phi_L & 0 \\ 0 & qh & 0 & -k \end{pmatrix}. \quad (45) $$

The characteristic polynomial evaluated at the extinction equilibrium $E_0$ is

$$ Q_0(\chi) = (\chi + \phi_L) (\chi^3 + A_2\chi^2 + A_1\chi + A_0), $$

with

$$ \begin{align*} A_2 &= k + \alpha_y + \alpha_a, & A_1 &= \alpha_y \alpha_a + k (\alpha_a + \alpha_y), \\ R_B &= \frac{k_1 \gamma r_1 q h}{k \alpha_y \alpha_a}, & A_0 &= k \alpha_a \alpha_y (1 - R_B). \end{align*} $$

We remark that $-\phi_L$ is an eigenvalue of $Q_0$. The other eigenvalues are solutions of the polynomial function $\chi \mapsto \chi^3 + A_2\chi^2 + A_1\chi + A_0$.

If $R_B > 1$, Descartes rule reveals that $Q_0$ has a unique positive root. Using the Hartman-Grobman theorem in [46], we conclude that $E_0$ is unstable. Assuming $R_B < 1$. We use the well-known Routh-Hurwitz stability criterion for stability by verifying conditions: $A_2, A_1, A_0 > 0$ and $A_2 A_1 - A_0 > 0$.

The positivity of $A_2, A_1$ and $A_0$ are obvious. A direct computation gives

$$ \begin{align*} A_2 A_1 - A_0 &= (k + \alpha_y + \alpha_a) (\alpha_y \alpha_a + k (\alpha_a + \alpha_y)) \\ &\quad - k \alpha_a \alpha_y (1 - R_B) = k^2 (\alpha_y + \alpha_a) \\ &\quad + (\alpha_y + \alpha_a) (\alpha_y \alpha_a + k (\alpha_a + \alpha_y)) + r_1 k_1 q h \gamma > 0. \end{align*} $$

We conclude that $E_0$ is locally asymptotically stable.

We assume that $R_B > 1$. The characteristic polynomial evaluated at the PFE $E_1$ is given by

$$ Q_1(\chi) = (\chi - (c\psi B_{y1} - \phi_L)) (\chi^3 + B_2\chi^2 + B_1\chi + B_0) $$
<figure><img src="image_25.png" /></figure>

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

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

<figure><img src="image_28.png" /><figcaption>**Figure E.3:** Dynamics of state variable $B_y$ of model (3) with: (a) $\alpha = 7 \times 10^{-3}$ day$^{-1}$, (b) $\alpha = 1.35 \times 10^{-2}$ day$^{-1}$ (same value as Fig. 12), (c) $\alpha = 0.01$ day$^{-1}$ (same value as Fig. E.2) and (d) $\alpha = 0.09$ day$^{-1}$ (same value as Fig. E.1). The initial conditions are the same as in Fig. 12. The remaining parameter values are given in Table 1 with $\delta_M = 10^{-5}$ USD$^{-1}$. We have $\alpha' \approx 7.34 \times 10^{-3}$ day$^{-1}$ and $\alpha'' \approx 0.022$ day$^{-1}$.</figcaption></figure>

Dynamics of Plutella xylostella, Brevicoryne brassicae and Hellula undalis on Cabbage. Sustainable Agriculture Research 1. doi:10.5539/sar.v1n2p229.

[6] Baidoo, P.K., Mochiah, M.B., 2016. Comparing the effectiveness of garlic (Allium sativum L.) and hot pepper (Capsicum frutescens L.) in the management of the major pests of cabbage Brassica oleracea (L.). Sustainable Agriculture Research 5. doi:10.5539/sar.v5n2p83.

[7] Barreira, L., Valls, C., 2012. Ordinary differential equations: Qualitative theory. volume 137. American Mathematical Society.

[8] Castillo-Chavez, C., Feng, Z., Huang, W., 2002. On the computation of R0 and its role on global stability. Mathematical Approaches for Emerging and Re-emerging Infectious Diseases: An introduction 125, 31–65.

[9] Charleston, D.S., Kfir, R., Dicke, M., Vet, L.E.M., 2005. Impact of botanical pesticides derived from Melia azedarach and Azadirachta indica on the biology of two parasitoid species of the diamondback moth. Biological Control 33, 131–142. doi:10.1016/j.biocontrol.2005.02.007.

[10] De Oliveira, A.C., Abreu de Siqueira, J.V., De Oliveira, J.V., Da Silva, J.E., Michereff Filho, M., 2011. Resistance of Brazilian diamondback moth populations to insecticides. Scientia Agricola 68, 154–159. doi:10.1590/S0103-90162011000200004.

[11] Van den Driessche, P., Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences 180, 29–48. doi:10.1016/S0025-5564(02)00108-6.

[12] Food and Agriculture Organization (FAO), 2021. World Food and Agriculture - Statistical yearbook 2021. Rome. https://doi.org/10.4060/cb4477en.

[13] Food and Agriculture Organization (FAO), 2022. Agricultural production statistics.2000-2020. FAOSTAT Analytical Brief Series No. 41. Rome. https://www.fao.org/documents/card/en/c/cb9180en.

[14] Furlong, M.J., Wright, D.J., Dosdall, L.M., 2013. Diamondback Moth Ecology and Management: Problems, Progress, and Prospects. Annual Review of Entomology 58, 517–541. doi:10.1146/annurev-ento-120811-153605. pMID: 23020617.

[15] Gautam, M.P., Singh, H., Kumar, S., Kumar, V., Singh, G., Singh, S.N., 2018. Diamondback moth, Plutella xylostella (Linnaeus) (Insecta: Lepidoptera: Plutellidae) a major insect of cabbage in India: A review. Journal of Entomology and Zoology Studies 6, 1394–1399.

[16] Gordon, H.S., 1954. Economic Theory of a Common-Property Resource: the Fishery. Journal of Political Economy 62, 124–142. doi:10.1086/257497.

[17] Guvenc, I., Yildirim, E., 2006. Increasing productivity with intercropping systems in cabbage production. Journal of Sustainable Agriculture 28, 29–44. doi:10.1300/J064v28n04_04.

[18] Huaripata, C., Sánchez, G., et al., 2019. Life Cycle of the Diamondback Moth Plutella xylostella L. (Lepidoptera: Plutelliidae), in Broccoli and Cauliflower under Laboratory Conditions. Peruvian Journal of Agronomy 3, 1–5. doi:10.21704/pja.v3i1.1132.

[19] Isman, M.B., 2006. Botanical insecticides, deterrents, and repellents in modern agriculture and an increasingly regulated world. Annual Review of Entomology 51, 45–66. doi:10.1146/annuev.ento.51.110104.151146.

[20] Jiang, J., Li, X., Wu, X., 2023. The Dynamics of a Bioeconomic Model with Michaelis–Menten Type Prey Harvesting. Bulletin of the Malaysian Mathematical Sciences Society 46, 57. doi:10.1007/s40840-022-01452-4.
where

$$
\begin{align*}
B_2 &= k + \alpha_y + \alpha_a + 2\mu_a B_{a1} + 2\mu_y B_{y1} \\
B_1 &= (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) \\
    &\quad + k (\alpha_a + \alpha_y + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \\
B_0 &= k (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) - \frac{k_1 q h \gamma r_1}{(1 + \delta_M M_1)^2}.
\end{align*}
$$

It is clear that $B_2, B_1 > 0$. Using Equation 11 and the fact that $1 + \delta_M M_1 > 1$, we have

$$
\begin{align*}
B_0 &> \frac{k}{1 + \delta_M M_1} \left( (\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) \right. \\
&\qquad \left. - (\alpha_y (\alpha_a + \mu_a B_{a1}) + \mu (\alpha_a + \mu_a B_{a1})^2 B_{a1}) \right) \\
&> \frac{k}{1 + \delta_M M_1} \left( (\alpha_y + \mu (\alpha_a + \mu_a B_{a1}) B_a) (\alpha_a + \mu_a B_{a1}) \right. \\
&\qquad \left. - (\alpha_y (\alpha_a + \mu_a B_{a1}) + \mu (\alpha_a + \mu_a B_{a1})^2 B_{a1}) \right) = 0.
\end{align*}
$$

Moreover,

$$
\begin{align*}
B_2 B_1 - B_0 &= k^2 (\alpha_a + \alpha_y + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \\
&\quad + (\alpha_y + \alpha_a + 2\mu_a B_{a1} + 2\mu_y B_{y1}) \\
&\quad ((\alpha_y + 2\mu_y B_{y1}) (\alpha_a + 2\mu_a B_{a1}) \\
&\qquad + k (\alpha_a + \alpha_y + 2\mu_a B_{a1} + 2\mu_y B_{y1})) > 0.
\end{align*}
$$

According to Routh-Hurwitz stability criterion, all the eigenvalues of the polynomial function $\chi \mapsto Q_1(\chi)$ have negative real parts. Thus, the local stability of $E_1$ depends only of the sign of $c\psi B_{y1} - \phi_L$. By noticing that $c\psi B_{y1} - \phi_L = \phi_L (\mathcal{R}_0 - 1)$, it comes from Hartman-Grobman theorem that $E_1$ is locally asymptotically stable provided $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

Now supposing $\mathcal{R}_B > 1$ and $\mathcal{R}_0 > 1$. The characteristic polynomial evaluated at $E^*$ is given by

$$
Q_2(\chi) = \chi^4 + C_3\chi^3 + C_2\chi^2 + C_1\chi + C_0 \quad (46)
$$

where

$$
\begin{align*}
J_{11} &= \alpha_y + 2\mu_y B_y^* + \psi L^* \\
J_{22} &= \alpha_a + 2\mu_a B_a^* \\
C_3 &= \mu_L L^* + J_{11} + k + J_{22} \\
C_2 &= c\psi^2 L^* B_y^* + \mu_L L^* J_{11} + k J_{22} \\
    &\quad + (k + J_{22}) (\mu_L L^* + J_{11}) \\
C_1 &= (k + J_{22}) (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) \\
    &\quad + k J_{22} (\mu_L L^* + J_{11}) - \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\
C_0 &= k J_{22} (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) - \mu_L L^* \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2}.
\end{align*}
$$

We want to use again the Routh-Hurwitz stability criterion. It consists in verifying the conditions :
$C_3, C_2, C_1, C_0 > 0$, $C_3C_2 - C_0 > 0$ and $(C_3C_2 - C_1)C_1 - C_3^2C_0 > 0$.

It is clear that $C_3, C_2 > 0$. In addition, using Equation (12)

$$
\begin{align*}
C_1 &= (k + J_{22}) (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) \\
    &\quad + k J_{22} (\mu_L L^* + J_{11}) - \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\
    &> k J_{11} J_{22} - \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\
    &> \frac{k}{1 + \delta_M M^*} \left( \left( \alpha_y + \frac{\mu_y}{\gamma} (\alpha_a + \mu_a B_a^*) B_a^* + \psi L^* \right) \right. \\
    &\qquad \left. (\alpha_a + \mu_a B_a^*) - (\alpha_y + \psi L^*) (\alpha_a + \mu_a B_a^*) \right. \\
    &\qquad \left. - \frac{\mu_y}{\gamma} (\alpha_a + \mu_a B_a^*)^2 B_a^* \right) = 0
\end{align*}
$$
[21] Kim, K.H., Kabir, E., Jahan, S.A., 2017. Exposure to pesticides and the associated human health effects. Science of the total environment 575, 525-535. doi:10.1016/j.scitotenv.2016.09.009.

[22] Kome, G., Enang, R., Yerima, B., 2018. Knowledge and management of soil fertility by farmers in western Cameroon. Geoderma Regional 13, 43–51. doi:10.1016/j.geodrs.2018.02.001.

[23] Kunbhar, S., Rajput, L.B., Gilal, A.A., Channa, G.A., Sahito, J.G.M., 2018. Impact of botanical pesticides against sucking insect pests and their insect predators in brinjal crop. Journal of Entomology and Zoology Studies 6, 83–87.

[24] Kusznierewicz, B., Bartoszek, A., Wolska, L., Drzewiecki, J., Gorin-
stein, S., Namieśnik, J., 2008. Partial characterization of white
cabbages (brassica oleracea var. capitata f. alba) from different regions
by glucosinolates, bioactive compounds, total antioxidant activities
and proteins. LWT-Food Science and Technology 41, 1–9. doi:10.
1016/j.lwt.2007.02.007.

[25] Luenberger, D.G., 1979. Introduction to dynamic systems: theory, models, and applications. volume 1. Wiley.

[26] Machekano, H., Mvumi, B.M., Nyamukondiwa, C., 2017. Dia-
mondback moth, Plutella xylostella (L.) in Southern Africa: research
trends, challenges and insights on sustainable management options.
Sustainability 9, 91. doi:10.3390/su9020091.

[27] Marino, S., Hogue, I.B., Ray, C.J., Kirschner, D.E., 2008. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of Theoretical Biology 254, 178–196. doi:https://doi.org/10.1016/j.jtbi.2008.04.011.

[28] Mason, P.G., Dancau, T., Abram, P.K., Noronha, C., Dixon, P.L., Parsons, C.K., Bahar, M.H., Bennett, A.M.R., Fernández-Triana, J., Brauner, A.M., et al., 2022. The parasitoid complex of diamondback moth, Plutella xylostella (Linnaeus)(Lepidoptera: Plutellidae), in Canada: impact and status. The Canadian Entomologist 154, e12. doi:10.4039/tce.2021.51.

[29] Misra, A.K., Jha, N., Patel, R., 2020. Modeling the effects of insects and insecticides on agricultural crops with NSFD method. Journal of Applied Mathematics and Computing 63, 197–215. doi:10.1007/s12190-019-01314-6.

[30] Misra, A.K., Rai, R.K., Takeuchi, Y., 2018. Modeling the effect of time delay in budget allocation to control an epidemic through awareness. International Journal of Biomathematics 11, 1850027. doi:10.1142/S1793524518500274.

[31] Mkenda, P., Mwanauta, R., Stevenson, P.C., Ndakidemi, P., Mtei, K., Belmain, S.R., 2015. Extracts from Field Margin Weeds Provide Economically Viable and Environmentally Benign Pest Control Compared to Synthetic Pesticides. PLOS ONE 10, e0143530. doi:10.1371/journal.pone.0143530.

[32] Mpumi, N., Machunda, R.S., Mtei, K.M., Ndakidemi, P.A., 2020. Selected insect pests of economic importance to Brassica oleracea, their control strategies and the potential threat to environmental pollution in Africa. Sustainability 12, 3824. doi:10.3390/su12093824.

[33] Mudzingwa, S., Muzemu, S., Chitamba, J., 2013. Pesticidal efficacy of crude aqueous extracts of Tephrosia vogelii L., Allium sativum L. and Solanum incanum L. in controlling aphids (Brevicoryne brassicaceae L.) in rape (Brassica napus L.) n controlling aphids (Brevicoryne brassicaceae L.) in rape (Brassica napus L.). Journal of Research in Agriculture 2, 157–163.

[34] Nazeer, A., Mukhtar, A., Muhammad, S., Hidayat, U., Toheed, I., Khalid, A.A.M., Kiran, S., Rafi, U., Saeed, A., Nibal, A.A.H.A., Hanem, F.K., Muhammad, S., 2021. Botanical Insecticides Are a Non-Toxic Alternative to Conventional Pesticides in the Control of Insects and Pests, in: El-Shafie, H.A.F. (Ed.), Global Decline of Insects. IntechOpen, Rijeka. chapter 6. doi:10.5772/intechopen.100416.

[35] Nofemela, R.S., Kfir, R., 2005. The role of parasitoids in suppressing diamondback moth, Plutella xylostella (L.) (Lepidoptera: Plutellidae), populations on unsprayed cabbage in the North West Province of South Africa. African Entomology 13, 71–83. doi:10.520/EJC32624.

[36] Pavela, R., 2012. Efficacy of three newly developed botanical insecticides based on pongam oil against Plutella xylostella L. larvae. Journal of Biopesticides 5, 62–70.

[37] Prasad, S.K., et al., 1963. Quantitative estimation of damage to crucifers caused by cabbageworm, cabbage looper, diatnonback moth and cabbage aphid. Indian Journal of Entomology 25, 242–259. URL: https://api.semanticscholar.org/CorpusID:83647851.

[38] Rai, R.K., Misra, A.K., Takeuchi, Y., 2019. Modeling the impact of sanitation and awareness on the spread of infectious diseases. Mathematical Biosciences and Engineering 16, 667–700. doi:10.3934/mbe.2019032.

[39] Sarfraz, M., Keddie, A.B., Dosdall, L.M., 2005. Biological control of the diamondback moth, plutella xylostella: A review. Biocontrol Science and Technology 15, 763–789. doi:10.1080/09583150500136956.

[40] Shuai, Z., Van den Driessche, P., 2013. Global stability of infectious disease models using lyapunov functions. SIAM Journal on Applied Mathematics 73, 1513–1532. URL: https://www.jstor.org/stable/24510691.

[41] jardins de la Vallée de la Siagne, L., 2022. Jeunes plants maraîchers bio. Website URL: https://www.associationjvs.com.

[42] Smith, H.L., 1995. Irreducible Cooperative Systems, in: Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems, American mathematical Society, Providence, RI. pp. 55–74.

[43] Smith, V.L., 1969. On Models of Commercial Fishing. Journal of Political Economy 77, 181–198. URL: http://www.jstor.org/stable/1829764.

[44] Sow, G., Diarra, K., Arvanitakis, L., Bordat, D., 2013. The relationship between the diamondback moth, climatic factors, cabbage crops and natural enemies in a tropical area. Folia Horticulturae 25, 3–12. doi:10.2478/fhort-2013-0001.

[45] Steinbrecher, A., Linseisen, J., 2009. Dietary intake of individual glucosinolates in participants of the EPIC-Heidelberg cohort study. Annals of Nutrition and Metabolism 54, 87–96. doi:10.1159/000209266.

[46] Talekar, N.S., Shelton, A.M., 1993. Biology, Ecology, and Management of the Diamondback moth. Annual Review of Entomology 38, 275–301. doi:10.1146/annurev.en.38.010193.001423.

[47] Tekle, B.T., Tesfu, N.B., 2023. Analysis of Head Cabbage Value Chain in Guji Zone, Southern Oromia, Ethiopia. International Journal of Agricultural Economics 8, 14–26. doi:10.11648/j.ijae.20230801.13.

[48] Tetsopgang, S., Fonyuy, F., 2019. Enhancing growth quality and yield of cabbage (Brassica oleracea) while increasing soil pH, chemicals and organic carbon with the application of fines from volcanic pyroclastic materials on a tropical soil in wum, Northwest Cameroon, Africa. Scientific African 6, e00199. doi:10.1016/j.sciaf.2019.e00199.

[49] Timbilla, J.A., Nyarko, K.O., 2004. A survey of cabbage production and constraints in Ghana. Ghana Journal of Agricultural Science 37, 93–101. doi:10.4314/gjas.v37i1.2084.

[50] Waiganjo, M.M., Waturu, C.N., Mureithi, J.M., Muriuki, J., Kamau, J., Munene, R., 2011. Use Of Entomopathogenic Fungi And Neem Bio-Pesticides For Brassica Pest Control And Conservation Of Their Natural Enemies. East African Agricultural and Forestry Journal 1&2.
and

$$C_0 = kJ_{22} (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) - \mu_L L^* \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2} \\ > \mu_L L^* \left( kJ_{11}J_{22} - \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2} \right) > 0.$$

Furthermore,

$$C_3C_2 - C_1 = (k + J_{22}) ((k + J_{22})(\mu_y L^* + J_{11}) + kJ_{22}) \\ + (\mu_y L^* + J_{11}) (c\psi^2 L^* B_y^* + \mu_L L^* J_{11}) \\ + (k_1 + J_{22})(\mu_L L^* + J_{11})) + \frac{qh\gamma r_1 k_1}{(1 + \delta_M M^*)^2} > 0$$

and,

$$(C_3C_2 - C_1) C_1 - C_3^2 C_0 > (k^2 + J_{22}^2)(\mu_L^* + J_{11})(k + J_{22}) \\ + (kJ_{22} + k^2 + J_{22}^2)(\mu_L L^* + J_{11})^2 > 0.$$

We conclude by using Routh-Hurwitz stability criterion that the endemic equilibrium $E^*$, whenever it exists, is always locally asymptotically stable.

# D Proof of Theorem 5

We assume that the solution $(B_y(t); B_a(t); L(t); M(t))$ starts in the domain $\Omega'$ which is described in Theorem 5.

If $B_a = B_{a_1}$, we have $\dot{B}_a = \gamma B_y - \mu_a B_{a_1}^2 - \alpha_a B_{a_1} = \gamma (B_y - B_{y_1}) \le 0$. then $B_a(t) \le B_{a_1}$ for all $t \ge t_0$. Similarly, if $M = M_1$, we obtain $\dot{M} = qhB_a - kM_1 = qh(B_a - B_{a_1})$. Hence, $M(t) \le M_1$ for all $t \ge t_0$. Finally, if $B_y = B_{y_1}$, we obtain

$$\begin{aligned} \frac{dB_y}{dt} &= r_{0y}B_{y_1} - \mu_y B_{y_1}^2 - \psi L B_{y_1} - \gamma B_{y_1} + r_1 \frac{k_1 M}{1 + \delta_M M} \\ &\le r_{0y}B_{y_1} - \mu_y B_{y_1}^2 - \gamma B_y + r_1 \frac{k_1 M}{1 + \delta_M M} \\ &\le r_1 k_1 \left( -\frac{M_1}{1 + \delta_M M_1} + \frac{M}{1 + \delta_M M} \right) \\ &\le r_1 k_1 \frac{M - M_1}{(1 + \delta_M M)(1 + \delta_M M_1)} \\ &\le 0. \end{aligned}$$

Then $B_y(t) \le B_{y_1}$. We conclude that $\Omega'$ is positively invariant.

Before dealing with the stability in $\Omega'$, we first prove the following half-way result.

**Theorem 8.** The pest-free equilibrium $E_1$ is globally asymptotically stable in the sub-region $\mathbb{R}_+^3 \setminus \{0\}$ of hyper-plan $L = 0$ if and only if $\mathcal{R}_0 \le 1 < \mathcal{R}_B$.

*Proof.* This proof is adapted from a result given in [47].

We consider

$$X_1 = g(X_1) \qquad (47)$$

with

$$X_1 = \begin{pmatrix} B_y \\ B_a \\ M \end{pmatrix}, \quad \text{and} \quad g(X_1) = \begin{pmatrix} -\mu_y B_y^2 - \alpha_y B_y + r_1 \frac{k_1 M}{1+\delta_M M} \\ \gamma B_y - \mu_a B_a^2 - \alpha_a B_a \\ qhB_a - kM \end{pmatrix}.$$

System (47) has two equilibria: 0 and $X_1^E = (B_{y_1}; B_{a_1}; M_1)$ on $\mathbb{R}_+^3$ if and only if $1 < \mathcal{R}_B$.

$\mathbb{R}_+^3$ is positively invariant. Furthermore, the Jacobian matrix of $g$ evaluated at $X_1$ is given by

$$\mathcal{J}_g(X_1) = \begin{pmatrix} -\alpha_y - 2\mu_y B_y & 0 & \frac{r_1 k_1}{(1+\delta_M M)^2} \\ \gamma & -\alpha_a - 2\mu_a B_a & 0 \\ 0 & qh & -k \end{pmatrix}. \qquad (48)$$
$\mathcal{J}_g(X_1)$ is a cooperative irreducible matrix. Then, the flow $\phi_t$ of system (47) is strongly monotone in $\mathbb{R}_+^3$. In addition, for every $\epsilon > 0$, the set $Z = [0; z_\epsilon]$ with $z_\epsilon = (K_y + \epsilon; K_a + \epsilon; K_M + \epsilon)^T$, is positively invariant and attractive for the flow $\phi_t$ as $Z$ is the projection of $\Omega_\epsilon$ on this subspace (see (44)).

According to Perron-Frobenius Theorem, the stability modulus $s$ of $\mathcal{J}_g(0)$, is a positive real eigenvalue of $\mathcal{J}_g(0)$ and there is a positive vector $v$ such that $\mathcal{J}_g(0)v = sv$.

Let $x$ be a positive vector of $\mathbb{R}_+^3$. Using attractiveness property of $Z$, we have $y = \phi_{t_1}(x) \in Z$ for large $t_1$.

Then, there is $\eta := \eta(y) > 0$ such that $\phi_t(y) \in [\eta v; z_\epsilon]$ which is a positive invariant order preserving interval and we have $g(z_\epsilon) < 0$ and $g(\eta v) \approx \mathcal{J}_g(0)\eta v = s\eta v > 0$. Using Proposition 3.2.1 in [47], $\phi_t(z_\epsilon)$ and $\phi_t(\eta v)$ converge monotonically both to $X_1^E$.

By monotonicity,

$$ \eta v < \phi_t(\eta v) < \phi_t(y) < \phi_t(z_\epsilon) < z_\epsilon, \quad t > 0. $$

Letting $t \to +\infty$ leads to the $w$-limit set of $y$ is $\{X_1^E\}$. Then $\phi_t(x)$ converges to $X_1^E$ for all positive vector $x$ of $\mathbb{R}_+^3$.

Furthermore, $B_y$ can not stay equal to 0. Indeed, if $B_y = 0$ and $M > 0$ at a time $t_0$, we have $\dot{B}_y = r_1k_1M/(1 + \delta_M M) > 0$. In addition, if $B_y = M = 0$ and $B_a = 0$ at a time $t_0$, one has $\dot{M} = qhB_a > 0$ so that $M$ instantaneously becomes positive and it follows that $\dot{B}_y = r_1k_1M/(1 + \delta_M M) > 0$. One can reproduce the same approach for variables $B_a$ and $M$. Thus, except if $x(t_0) = 0$, $x$ instantaneously enters the positive orthant $\mathbb{R}_+^3$ and is amenable to the previous developments.

Thus, $X_1^E$ is asymptotically stable everywhere in $\mathbb{R}_+^3 \setminus \{0\}$.

For convenience, we write model (6) in the following vector form.

$$ \dot{X}_1 = F(X_1; X_2) \quad \text{and} \quad \dot{X}_2 = G(X_1; X_2) $$

with

$$ X_1 = \begin{pmatrix} B_y \\ B_a \\ M \end{pmatrix}, \quad X_2 = L $$

$$ F(X_1; X_2) = \begin{pmatrix} r_{0y}B_y - \mu_y B_y^2 - \psi L B_y - \gamma B_y + r_1 \frac{k_1 M}{1 + \delta_M M} \\ \gamma B_y + r_{0a} B_a - \mu_a B_a^2 - h B_a \\ qh B_a - k M \end{pmatrix} $$

$$ G(X_1; X_2) = c\psi L B_y - \mu_L L^2 - \phi_L L. $$

In this case, the pest-free equilibrium $E_1$ becomes $(X_1^E; 0)$.

According to Theorem 8, $X_1^E$ is globally asymptotically stable for $\dot{X}_1 = F(X_1; 0)$.

In addition, we remark that

$$ G(X_1; X_2) = AX_2 - \hat{G}(X_1; X_2) $$

where, for all variables in $\Omega'$, one has

$$ A = \frac{\partial G}{\partial X_2}(X_1^E, 0) = \phi_L(\mathcal{R}_0 - 1) $$

$$ \text{and } \hat{G}(X_1; X_2) = c\psi(B_{y_1} - B_y) + \mu_L L^2 \geq 0. $$

Thus, according to [48], the PFE $E_1$ is asymptotically stable in $\Omega'$ if and only if $\mathcal{R}_0 \leq 1 < \mathcal{R}_B$.

We comment that the uniform persistence of sub-model (6) can be obtained in $\Omega'$ when $\mathcal{R}_0 \leq 1 < \mathcal{R}_B$ [49].

## E Other simulations

In this appendix, we simulate sub-model (20) and model (3) for different values of parameter $\alpha$, which represents the investment rate in plant protection. All remaining parameters are given in Table 1.

Figures Fig. 12 and Fig. 13 are similar to Fig. 11 (in which $\alpha = 1.35 \times 10^{-2}$ day$^{-1}$) and depict the dynamics of sub-model (20). Fig. 12 presents a situation where the pest-free equilibrium is always stable and the system does not have an endemic equilibrium. It corresponds to a high investment in plant protection ($\alpha = 0.09$ day$^{-1}$), associated with a cost of 1,530 USD computed according to (38). It generates a seasonal yield of 22.2 tonnes for an income of 15,320 USD, computing according to (19) at the pest-free equilibrium.

Fig. 13 illustrates the situation where there is a unique endemic stable equilibrium and an unstable pest-free equilibrium. It corresponds to a low investment in plant protection ($\alpha = 0.01$ day$^{-1}$), associated with a cost of 79 USD, and it generates a seasonal yield of 10.3 tonnes for an income of 7,111 USD at the endemic equilibrium.
<figure><img src="image_15.png" /><figcaption>Figure 12: Dynamics of state variables (a) $B_y$, (b) $L$, (c) $P$ and (d) $M$ of sub-model (20). The initial conditions and parameter values are the same as in Fig. 11 except for $\alpha = 0.09$ day$^{-1}$ leading to an effective reproductive number $\mathcal{R}_e = 0.19 < \mathcal{R}_e(\alpha^s) \approx 0.2 < 1$. For all initial conditions, the system converges towards the pest-free equilibrium.</figcaption></figure>

The effective reproductive number of original model (3) is given by

$$ \mathcal{R}'_e = \frac{c\psi\alpha_a B_{a1}}{\gamma(\lambda\phi P_1 + \phi_L)(1 + bP_1)}, $$

with $P_1 = a_p\alpha M_1/\phi_P$. The equivalence of the threshold $\alpha^0$ for the original model (3), obtained by solving $\mathcal{R}'_e(\alpha) = 1$, is given by $\alpha^{0'} = \phi_P\Delta/(2a_pM_1\lambda b\phi)$, with $\Delta = -(\lambda\phi + b\phi_L) + \sqrt{(\lambda\phi + b\phi_L)^2 + 4\lambda b\phi (c\psi B_{y1} - \phi_L)}$.

Using parameter values in Table 1, we have $\alpha^{0'} \approx 7.34 \times 10^{-3}$ day$^{-1}$. The lower bound for the existence of endemic equilibria $\alpha^{s'} \approx 0.022$ day$^{-1}$. Fig. 14 shows model (3), for different values of parameter $\alpha$. Comparing these dynamics with those of sub-model (20) depicted in panel (a) of Fig. 11, 12 and 13, we observe that both models have similar dynamics. In fact, we observe a bistability phenomenon for $\alpha \in (7.35 \times 10^{-3}; 0.022)$ day$^{-1}$ (see Fig. 14(b) and Fig. 14(c)). When $\alpha \ge 0.022$ day$^{-1}$, the solutions tend to the pest-free equilibrium of (3) (see Fig. 14(d)), and when $\alpha \le 7.34 \times 10^{-3}$ day$^{-1}$, they converge towards the endemic equilibrium of (3) (see Fig. 14(a)). Furthermore, $\alpha^0 \gg \alpha^{0'}$ means that a pesticide with a double action is more effective and costless than a classical one. Such results encourage vulgarizing the use of botanical pesticides.

## Acknowledgments

This work was carried out within EPITAG (https://team.inria.fr/epitag/), an Inria Associate team part of the LIRIMA, which is supported by INRAE via the TSARA program and by the French Embassy in Cameroon.

## Declaration of competing interest

The authors declare that they have no competing interests that are directly or indirectly related to this work.
<figure><img src="image_16.png" /><figcaption>Figure 13: Dynamics of state variables (a) $B_y$, (b) $L$, (c) $P$ and (d) $M$ of sub-model (20). The initial conditions and parameter values are the same as in Fig. 11 except for $\alpha = 0.01 \text{ day}^{-1}$ leading to an effective reproductive number $R_e = 1.02 > 1$. For all initial conditions, the system converges towards the endemic equilibrium.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 14: Dynamics of state variable $B_y$ of model (3) with: (a) $\alpha = 7 \times 10^{-3}$ day$^{-1}$, (b) $\alpha = 1.35 \times 10^{-2}$ day$^{-1}$ (same value as Fig. 11), (c) $\alpha = 0.01$ day$^{-1}$ (same value as Fig. 13) and (d) $\alpha = 0.09$ day$^{-1}$ (same value as Fig. 12). The initial conditions are the same as in Fig. 11. The remaining parameter values are given in Table 1 with $\delta_M = 10^{-5}$ USD$^{-1}$. We have $\alpha^{0'} \approx 7.34 \times 10^{-3}$ day$^{-1}$ and $\alpha^{s'} \approx 0.022$ day$^{-1}$.</figcaption></figure>