# Feasibility and optimization results for elimination by mass-trapping in a metapopulation model

Pierre-Alexandre Bliman, Manon de la Tousche, Yves Dumont

► To cite this version:

Pierre-Alexandre Bliman, Manon de la Tousche, Yves Dumont. Feasibility and optimization results for elimination by mass-trapping in a metapopulation model. Applied Mathematical Modelling, 2025, 144, pp.116047. 10.1016/j.apm.2025.116047. hal-04688245

HAL Id: hal-04688245

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

Submitted on 4 Sep 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
Feasibility and optimization results for elimination by mass-trapping in a metapopulation model

Pierre-Alexandre Bliman*¹, Manon de la Tousche†¹,³, and Yves Dumont✉‡²,³,⁴

¹Sorbonne Université, Inria, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions UMR 7598, Equipe MUSCLEES, Paris, France

²CIRAD, Umr AMAP, Pôle de Protection des Plantes, F-97410 Saint Pierre, France

³AMAP, Univ Montpellier, CIRAD, CNRS, INRA, IRD, Montpellier, France

⁴University of Pretoria, Department of Mathematics and Applied Mathematics, Pretoria, South Africa

September 4, 2024

**Abstract**

Having in mind the issue of control of insects vectors or insects pests, we consider in this paper a metapopulation model with patches linearly interconnected, and explore the global effects of the (on purpose) increase of mortality in some of them. Based on previous results by Y. Takeuchi *et al.*, we show that under appropriate conditions, the sign of the stability modulus of the Jacobian of the system at the origin determines the asymptotic behaviour of the solutions. If it is non-positive, then the population becomes extinct in every patch. Conversely, if it is positive, then there exists a unique nonnegative equilibrium, which is positive and globally asymptotically stable.

In the latter case, given a subset of 'controlled' patches where human intervention is allowed, through mass-trapping for instance, we study whether the introduction of additional linear mortality in some of them can result in population elimination in every patch. We characterize this possibility by an algebraic property on the Jacobian at the origin of a so-called residual system. We then assess the minimal globally asymptotically stable equilibrium that may be attained in this way, and when elimination is possible, we study the optimization problem consisting in achieving this task while minimizing

<footnote>*Email: pierre-alexandre.bliman@inria.fr</footnote>

<footnote>†Email: manon.de-la-tousche@inria.fr</footnote>

<footnote>‡Corresponding author, email: yves.dumont@cirad.fr</footnote>
a certain cost function, chosen as a nondecreasing and convex function of the mortality rates added in the controlled patches. We show that such minimization problem admits a global minimizer, which is unique in the relevant cases. An interior point algorithm is proposed to compute the numerical solution.

**Keywords:** Pest; Vector; Metapopulation; Control strategy; Monotone system; Interior-point algorithm; Oriental fruit fly

Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td></td><td>General notations and notions</td><td>5</td></tr><tr><td>2</td><td>Mathematical model and preliminary results</td><td>6</td></tr><tr><td>3</td><td>General properties of a controlled model</td><td>9</td></tr><tr><td>3.1</td><td>Controlled model and assumptions</td><td>9</td></tr><tr><td>3.2</td><td>Study of the control effects</td><td>10</td></tr><tr><td>3.3</td><td>Limit behaviour for large control</td><td>10</td></tr><tr><td>4</td><td>Minimal reachable equilibrium</td><td>12</td></tr><tr><td>4.1</td><td>Irreducible residual graph</td><td>12</td></tr><tr><td>4.2</td><td>Reducible residual graph</td><td>14</td></tr><tr><td>5</td><td>A cost-effective approach for achievable elimination scenarios</td><td>17</td></tr><tr><td>5.1</td><td>Minimization problem and properties</td><td>18</td></tr><tr><td>5.2</td><td>Interior-point algorithm</td><td>18</td></tr><tr><td>6</td><td>Application to <i>Bactrocera dorsalis</i></td><td>21</td></tr><tr><td>6.1</td><td>Three-patch model</td><td>22</td></tr><tr><td>6.1.1</td><td>Intrinsic growth rate and prices effects</td><td>22</td></tr><tr><td>6.1.2</td><td>Exponential convergence rate</td><td>24</td></tr><tr><td>6.1.3</td><td>Diffusion effects</td><td>25</td></tr><tr><td>6.2</td><td>Seven-patch model and the impact of the network configuration</td><td>27</td></tr><tr><td>7</td><td>Conclusion</td><td>32</td></tr><tr><td></td><td>References</td><td>34</td></tr><tr><td></td><td>Appendix – Demonstrations</td><td>38</td></tr><tr><td></td><td>A Fundamental stability result (Thm. 2.2)</td><td>38</td></tr></table>
is widely used, and some methods are specifically designed to target and increase mortality. For instance, gravid traps utilize water infused with decomposing organic material to attract gravid females in search of egg-laying sites, thereafter either retaining or eliminating them [15, 26, 37], and attractive toxic sugar baits also play a significant role [5, 38]. In fruit flies management, Methyl Eugenol, when used together with an insecticide, is frequently used against *Bactrocera dorsalis* to attract and kill male flies [45]. Other mass trapping methods aim to reduce fecundity. For example, ovitraps, particularly used for *Aedes* mosquitoes, are tailored to intercept eggs laid by ovipositing females [30]. As other biological control techniques, the sterile insect technique (SIT) or the incompatible insect technique (IIT) are very promising strategies consisting in massive releases of sterile or incompatible males. After mating with these males, the wild females do not produce viable eggs, thereby contributing to reducing the size of the population [14]. This method has already been implemented successfully in the field [51], and a feasibility project called AttracTIS that combines SIT and other control tools against *Bactrocera dorsalis*, is ongoing in Réunion Island.

These issues deserve deep studies, for which mathematical modelling can be of great help. Mathematical models not only help aggregate and formalize knowledge, but make it possible to identify the processes and parameters important to be studied. These approaches are particularly valuable as an alternative to extensive and potentially impractical simultaneous field studies that are often long, tedious, and costly, or impossible to conduct simultaneously.

In this paper, we delve into modeling mass trapping methods that affect mortality rates. We take into account the spatial aspects, which influence significantly the efficacy of this strategy, since in practice, the surfaces on which mass trapping is applied are limited and heterogeneous. Our focus lies within a metapopulation framework, which consists of a group of populations of the same species occupying distinct homogeneous patches of suitable habitat within a larger landscape. This concept was introduced by Richard Levins in 1969 to describe a model of population dynamics of insect pests in agricultural fields [35]. It has been used mainly by ecologists to show that the connectivity between zones in heterogeneous landscapes promotes species conservation [25]. Subsequently, these metapopulation models were extended to epidemiological issues and other biological processes involving interactions among dispersed populations [1, 27]. Despite biologists' significant interest in metapopulation, such models have been primarily underused in pest management, focusing more on the study of population persistence. For example in [2], the authors proved the existence of a threshold number of source patches necessary to ensure this property.

With the aim of designing efficient control strategies, we consider in the present work the opposite question, exploring the conditions necessary for population elimination. We consider that the control action is limited to a specific number of plots. This limitation might stem e.g. from environmental constraints (gullies, marshes...) that prevent control implementation, reluctance of plot owners to allow any form of intervention, or certain fallow lands remaining untreated. We aim to explore the influence of an untreated land on its neighbours, as well as the impact of treated areas on their surroundings. Specifically, we seek to recognize whether the incapacity to operate in a specific region may detrimentally affect the overall efficacy of the strategy.
The paper is organized as follows. In Section 2, we present a model describing the natural evolution of a population in a patchy environment. We recall and refine Takeuchi's results for such model [36, 47], which state that the elimination or persistence of the population is determined by the sign of the stability modulus of the Jacobian matrix at the origin of the system. In Section 3, we study the feasibility of population elimination by introducing additional mortality terms in a subset of controllable patches. We assess in Section 4 the minimal attainable equilibrium value. When elimination is feasible, we study in Section 5 the optimization problem consisting in achieving this task while minimizing a certain cost function, chosen as a non-decreasing and convex function of the mortality rates added in the controlled patches. We then present an algorithm that computes the numerical solution. In Section 6, we apply the algorithm for a population of *Bactrocera dorsalis* to examine how the model parameters, especially diffusion coefficients, affect the optimal control strategy. For reader's convenience, the proofs of the results are put in Appendix.

## General notations and notions

Let us introduce some general notations and definitions that will be used in this paper. First, the inequalities between vectors are considered in their usual coordinate-wise sense, that is, for any $x = (x_1, ..., x_n)$ and $y = (y_1, ..., y_n)$,

$$ \bullet \ x \le y \iff x_i \le y_i, \ i = 1, ..., n, $$

$$ \bullet \ x < y \iff x \le y, x \ne y, $$

$$ \bullet \ x \ll y \iff x_i < y_i, \ i = 1, ..., n. $$

These definitions are extended to matrices.

**Definition 1.1.** The cardinal of any non-empty subset $C$ of $\{1, ..., n\}$ is denoted $n_C$. For any $x = (x_1, ..., x_n)$, the point of $\mathbb{R}^{nc}$ composed of the $n_C$ components of $x$ with index in $C$ is denoted $x|_C$:

$$ x|_C := (x_i)_{i \in C}. $$

Similarly, for any non-empty subsets $C, D$ of $\{1, ..., n\}$ and any matrix $A \in \mathcal{M}_n(\mathbb{R})$, we define the matrix

$$ A|_{C,D} := (A_{ij})_{i \in C, j \in D}. $$

**Definition 1.2.** We denote

$$ \bar{\mathbb{R}}_+ := \mathbb{R}_+ \cup \{+\infty\} $$

the extended positive real number line.

**Definition 1.3.** [32, p. 31] A square matrix $A \in \mathcal{M}_n(\mathbb{R})$ is said reducible if there exists a permutation matrix $P$ such that $P^TAP$ is block triangular. Otherwise, $A$ is called irreducible.

**Definition 1.4 (Metzler matrix).** A square matrix $A \in \mathcal{M}_n(\mathbb{R})$ is said to be a Metzler matrix if all its non-diagonal entries are non-negative.
Metzler matrices are also called essentially non-negative matrices [13]. Many of their properties are derived by extending properties of non-negative matrices [6].

**Definition 1.5.** The stability modulus of a square matrix $A \in \mathcal{M}_n(\mathbb{R})$ is defined as

$$s(A) := \{\max\{\mathrm{Re}(\lambda)\} : \lambda \text{ is an eigenvalue of } A\}.$$

By convention, when A is a matrix of dimension zero, we put $s(A) = -\infty$.

See e.g. [32, p. 32] for more details. The value $s(A)$ is also called the spectral abscissa [2]. Finally, when A is Metzler and irreducible, $s(A)$ is also referred to as the Perron value or Perron root of A, denoted $r(A)$ (see for example [13]).

**Definition 1.6.** Let $\alpha \ge 0$. The solutions of a dynamical system are said to converge towards the origin at an exponential rate $\alpha$ if $e^{\alpha t}x(t)$ converges to the origin along every trajectory.

**Definition 1.7.** Let $\alpha \in \bar{\mathbb{R}}_+$ The origin of a dynamical system is said $\alpha$-globally asymptotically stable ($\alpha$-GAS) if the following conditions are fulfilled.

(i) The origin of the system is globally asymptotically stable;

(ii) For any $0 < \alpha' < \alpha$, the solutions of the system converge towards the origin at an exponential rate $\alpha'$.

As a consequence of this definition, the origin is said 0-GAS iff it is GAS.

## 2 Mathematical model and preliminary results

We consider the following population model, composed of $n$ patches linearly interconnected and occupied by a single species:

$$\dot{x}_i = x_i g_i(x_i) + \sum_{j=1, j \neq i}^{n} D_{ij} x_j - \sum_{j=1, j \neq i}^{n} D_{ji} x_i, \quad x_i(0) \ge 0, \quad i = 1, \dots, n. \quad (1)$$

For every $i = 1, \dots, n$, $x_i(t)$ is a scalar representing the number of individuals in the $i$-th patch at the time $t$, $\dot{x}_i$ is the derivative of $x_i$ with respect to time, and $g_i$ is the specific growth rate of the population in the $i$-th patch. $D_{ij}$ is a non-negative diffusion coefficient describing the flow of individuals from the patch $j$ to the patch $i$ ($i \neq j$).

We note $x = (x_1, \dots, x_n)$ and we define $F(x)$ as the right hand side of the model (1), i.e.

$$F(x) := (F_1(x), \dots, F_n(x)), \quad F_i(x) := x_i g_i(x_i) + \sum_{j=1, j \neq i}^{n} D_{ij} x_j - \sum_{j=1, j \neq i}^{n} D_{ji} x_i. \quad (2)$$
Let us introduce the matrix $D = (d_{ij})$ where

$$d_{ij} = \begin{cases} D_{ij} & \text{if } j \neq i \\ -\sum_{\substack{k=1 \\ k \neq i}}^{n} D_{ki} & \text{if } i = j. \end{cases} \qquad (3)$$

$D$ is called the connectivity matrix of the system (1), which writes in vector form

$$\dot{x} = (D + \text{diag}(g_i(x_i))) x.$$

Notice that the system (1) is cooperative, and some proofs of the results presented in this paper rely on the theory of monotone cooperative systems. For more details about theoretical properties of such class of models, the reader may refer to [46].

Throughout this article, for any $i = 1, ..., n$, $g_i$ is assumed continuous and locally Lipschitz. We introduce the following hypotheses.

(H1). For any $i = 1, ..., n$, $g_i(0) > 0$.

(H2). For any $i = 1, ..., n$, $g_i$ is decreasing on $\mathbb{R}_+$.

(H3). For any $i = 1, ..., n$, $\lim_{x_i \to +\infty} g_i(x_i) < 0$.

(H4). For any $i = 1, ..., n$, the matrix $D$ is irreducible.

$D$ being irreducible means that the underlying directed graph is strongly connected [6, Theorem 2.7], that is, there always exists a path between two patches. Moreover, it is worth noticing that $D$ is a Metzler matrix.

Assumptions (H1), (H2) and (H3) describe the local behaviour of the population. In absence of migration coupling the patches, (H1) means that when there are very few individuals in the patch, the population in this patch increases. Assumptions (H2) and (H3) express competition effects between individuals from the same patch. The more individuals, the less increase of the population. The latter is indeed bounded, since it decreases when it goes beyond a certain value (the unique zero of $g_i$), called the carrying capacity of the environment. Typical functions satisfying (H1), (H2) and (H3) are the well-known logistic function $g_i(x_i) = r_i(1 - \frac{x_i}{K_i})$, where $K_i > 0$ is the carrying capacity and $r_i > 0$ denotes the intrinsic growth rate; or $g_i(x_i) = b_i e^{-\beta_i x} - m_i$, where $b_i > 0$ and $m_i > 0$ represent respectively the birth rate at low population level and the death rate, such that $b_i - m_i > 0$, and $\beta_i$ takes into account the indirect/direct competition between individuals [8]. For this second model, the carrying capacity is equal to $-\frac{1}{\beta_i} \log \frac{m_i}{b_i}$.

Model (1) has already been the subject of research work, see [16, 17, 18]. In [17], the function $x_i \mapsto x_ig_i(x_i)$ is a logistic function, and a coefficient $\beta$ is set as a factor of the sums in model (1). The author shows that in case of perfect mixing, i.e. when the migration rate $\beta$ tends to infinity, the total population follows a logistic law, with a carrying capacity generally different from the sum of the $n$ carrying capacities of each patch.
In absence of dispersal, model (1) simplifies into the uncoupled system

$$ \dot{x}_i = x_i g_i(x_i), \quad x_i(0) \ge 0, \quad i = 1, \dots, n. \qquad (4) $$

It is easy to prove the following result.

**Theorem 2.1.** Assume (H2) and (H3) hold.

(i) If $g_i(0) \le -\alpha$ for some $\alpha \ge 0$, then the origin is $\alpha$-GAS on $\mathbb{R}_+$ for the system (4).

(ii) If $g_i(0) > 0$, then (4) admits a positive equilibrium point which is GAS on $\mathbb{R}_+^*$.

In this section, we extend this result to model with migration (1), drawing inspiration from the work of Takeuchi in [47] and Lu and Takeuchi in [36], where is studied the similar system

$$ \dot{x}_i = x_i g_i(x_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij}(x_j - x_i), \quad x_i(0) \ge 0, \quad i = 1, \dots, n. \qquad (5) $$

In model (5), the net exchange from the j-th patch to the i-th patch is proportional to the difference $x_j - x_i$ of population densities in patches i and j. The authors showed that the asymptotic behaviour of the solutions of system (5) is determined by the sign of the stability modulus of the Jacobian matrix of the right hand side of the system at the origin.

We now show in Theorem 2.2 that this result also applies to the model (1). Let us define

$$ A := D + \mathrm{diag}(g_i(0)), \qquad (6) $$

which is the Jacobian matrix at the origin of the right hand side of the system (1). The following result classifies the asymptotic behaviour of the solutions of the system (1), according to the value of the stability modulus $s(A)$.

**Theorem 2.2.** Assume Hypothesis (H2) holds.

(i) If $s(A) \le -\alpha$ for some $\alpha \ge 0$, then the origin is $\alpha$-GAS on $\mathbb{R}_+^n$ for the system (1).

(ii) If $s(A) > 0$ and Assumptions (H3) and (H4) also hold, then the system (1) admits a unique non-negative and non-zero equilibrium. Furthermore, the latter is positive and GAS on $\mathbb{R}_+^n \setminus \{0_n\}$.

In particular, Theorem 2.2 gives a simple necessary and sufficient condition for population elimination, namely $s(A) \le 0$.
# 3 General properties of a controlled model

## 3.1 Controlled model and assumptions

The aim of this section is to determine whether the population can be eliminated at a prescribed exponential convergence rate by introduction of linear, diagonal, control terms, and otherwise, to assess the minimal globally asymptotically stable equilibrium that may be attained.

We consider that in the field, the action is limited to a specific number of plots, stemming e.g. from environmental constraints that prevent control implementation or from the reluctance of plot owners to allow any form of intervention. We also assume that the control is realized by introduction of supplementary mortality terms that may be arbitrarily large in some patches, but bounded in others due to practical limitations on achievable mortality rates.

For some $\rho \in \mathbb{R}_+^n$, consider the controlled system

$$ \dot{x}_i = x_i(g_i(x_i) - \rho_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij}x_j - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji}x_i, \quad i = 1, \dots, n. \quad (7) $$

The system (7) is identical to system (1), except that the functions $g_i$ have been replaced by $g_i - \rho_i$, where $\rho_i$ is the additional mortality term induced by the control. Notice that these new functions verify Hypotheses (H2) and (H3) when the functions $g_i$ do. Moreover, the connectivity matrix of the system (7) is equal to $D$, and is thus irreducible when (H4) is verified.

We now introduce a boundedness assumption on the vector $\rho$ in (7).

**Definition 3.1.** Let $\bar{\rho} \in \bar{\mathbb{R}}_+^n$. Define the sets $C_F, C_I$ as follows:

$$ C_F := \{i \in \{1, \dots, n\} : 0 < \bar{\rho}_i < +\infty\}, \quad C_I := \{i \in \{1, \dots, n\} : \bar{\rho}_i = +\infty\} \quad (8a) $$

and

$$ C := C_F \cup C_I. \quad (8b) $$

The sets $C, C_F$ and $C_I$ are called respectively the set of controllable patches, the set of finitely controllable patches and the set of infinitely controllable patches. The complementary set $\bar{C} = \{i \in \{1, \dots, n\} : \bar{\rho}_i = 0\}$ of $C$ is called the set of uncontrollable patches.

**Definition 3.2.** For any $\bar{\rho} \in \bar{\mathbb{R}}_+^n$, we call $\bar{\rho}$-admissible any $\rho \in \mathbb{R}_+^n$ satisfying $\rho \le \bar{\rho}$.

Notice that any $\bar{\rho}$-admissible $\rho$ satisfies $\rho_i = 0$, for every $i \in \bar{C}$.

For any $\bar{\rho}$-admissible $\rho \in \mathbb{R}_+^n$, denote $A(\rho)$ the Jacobian matrix of system (7) at the origin:

$$ A(\rho) := A - \mathrm{diag}(\rho). \quad (9) $$

In particular, $A(0_n) = A$, and $s(A(0_n))$ is simply written $s(A)$.

Theorem 2.2 applies to the system (7), allowing to introduce the following definition.
**Definition 3.3.** Assume (H2), (H3) and (H4) hold. For any $\bar{\rho}$-admissible $\rho$, let $x^*(\rho)$ be the (positive) globally asymptotically stable equilibrium of (7) on $\mathbb{R}_+^n \setminus \{0_n\}$ if $s(A(\rho)) > 0$; and $x^*(\rho) = 0_n$ if $s(A(\rho)) \le 0$.

A key point is that, whatever the value of $s(A(\rho))$, $x^*(\rho)$ is always the equilibrium point that attracts every non-zero trajectory.

## 3.2 Study of the control effects

In this section, we study some properties of the control (Theorem 3.1 and Theorem 3.2). Demonstrations are provided in Appendix B.

**Theorem 3.1.** Assume (H4) holds. The function which associates to any $\bar{\rho}$-admissible $\rho$ the value $s(A(\rho))$ is twice differentiable and decreasing, that is, for any $\bar{\rho}$-admissible $\rho$ and $\rho' \in \mathbb{R}_+^n$,

$$ \rho < \rho' \implies s(A(\rho')) < s(A(\rho)). $$

Moreover, it is convex, and strictly convex if $C \subsetneq \{1, \dots, n\}$.

As long as $x^*(\rho) \gg 0_n$ i.e. as long as $s(A(\rho)) > 0$, the following theorem ensures that the higher the additional mortality terms, the larger the population reduction.

**Theorem 3.2.** Assume (H2), (H3) and (H4) hold. The function which associates to any $\bar{\rho}$-admissible $\rho$ such that $s(A(\rho)) > 0$ the value $x^*(\rho)$ is strictly decreasing in the following sense: for any $\bar{\rho}$-admissible $\rho$, $\rho'$ such that $x^*(\rho') \gg 0_n$,

$$ \rho < \rho' \quad \Rightarrow \quad x^*(\rho') \ll x^*(\rho). $$

Recall that $x^*(\rho)$ is the equilibrium of the system (7) attracting every nonzero trajectory (see Definition 3.3).

Theorem 3.2 implies that introducing mortality in any non-void subset $C \subset \{1, \dots, n\}$ induces population reduction in every patch of the network. In particular, introducing a control in a single patch reduces the positive equilibrium across all patches.

## 3.3 Limit behaviour for large control

A natural issue is now to assess whether, for given $\bar{\rho} \in \mathbb{R}_+^n$ and $\alpha \ge 0$, the origin may be rendered $\alpha$-GAS for some $\bar{\rho}$-admissible controls. As having $s(A(\rho)) \le -\alpha$ for some $\alpha \ge 0$ ensures $\alpha$-global asymptotic stability of the origin for system (7) (see Theorem 2.2), and since the stability modulus decreases with the control (Theorem 3.1), our analysis focuses on evaluating the stability modulus of $A(\rho)$ for large $\bar{\rho}$-admissible $\rho$.

By slight abuse of notation, one denotes $A|_{\bar{c}_I, \bar{c}_I}(\bar{\rho}) = A(\bar{\rho})|_{\bar{c}_I, \bar{c}_I}$. A key technical result is now stated. Its proof in Appendix C.1 is based on singular perturbation theory [34].
**Lemma 3.3.** For $A(\rho)$ defined in (9), one has

$$ \inf\{s(A(\rho)) : \rho \bar{\rho}\text{-admissible}\} = \lim_{\rho \to \bar{\rho}} s(A(\rho)) = s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})), \quad (10) $$

where $C_I$ is defined in (8a). Moreover, the infimum is a minimum if and only if $C_I = \emptyset$.

Building on this lemma, one derives the following theorem. A proof is provided in Appendix C.2.

**Theorem 3.4.** Assume (H2), (H3) and (H4) hold.

* If $C_I$ is not empty, then the following statements are true.

(i) If $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) < 0$, then there exists $\rho' \bar{\rho}$-admissible such that, for any $\rho \bar{\rho}$-admissible satisfying $\rho \ge \rho'$, $x^*(\rho) = 0_n$. Moreover, the origin of (7) is then $\alpha$-GAS for $\alpha := -s(A(\rho))$.

(ii) If $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) \ge 0$, then $x^*(\rho) \gg 0_n$ for any $\bar{\rho}$-admissible $\rho$.

* If $C_I$ is empty, then the following holds.

(iii) If $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) \le 0$, then there exists $\rho' \bar{\rho}$-admissible such that, for any $\rho \bar{\rho}$-admissible satisfying $\rho \ge \rho'$, $x^*(\rho) = 0_n$. Moreover, the origin of (7) is then $\alpha$-GAS for $\alpha := -s(A(\rho))$.

(iv) If $s(A(\bar{\rho})) > 0$, then $x^*(\rho) \gg 0_n$ for any $\bar{\rho}$-admissible $\rho$.

Notice that elimination feasibility does not solely depend on the exchanges between the patches that are either uncontrollable or finitely controllable: the leaks from these patches towards infinitely controllable patches must also be considered, as they appear in the diagonal terms of the matrix $A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})$. As seen below in Section 4, this matrix is the Jacobian matrix in zero of the residual system on the set of patches belonging to $\bar{C}_I$ (defined in (11)).

Theorem 3.4 states in particular that, if every patch is infinitely controllable (that is $C_I = C$), then it is always possible to eliminate the population at any desired exponential convergence rate. As a matter of fact, in this case, $\bar{C}_I = \emptyset$, so that $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) = -\infty$.

Using conditions sufficient to have $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) < 0$ or $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) > 0$ yields the following result, for which a proof is given in Appendix C.3.

**Corollary 3.5.** Assume (H2), (H3) and (H4) hold, and let $C_I \subsetneq \{1, ..., n\}$. Then, the two following statements are true.

(i) If $g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} < 0$ for every $i \in \bar{C}_I$, then there exists $\rho' \bar{\rho}$-admissible such that, for any $\rho \bar{\rho}$-admissible satisfying $\rho \ge \rho'$, $x^*(\rho) = 0_n$. Moreover, the origin of (7) is then $\alpha$-GAS for $\alpha := -s(A(\rho))$.

(ii) If $g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} > 0$ for every $i \in \bar{C}_I$, then $x^*(\rho) \gg 0_n$ for any $\bar{\rho}$-admissible $\rho$.

**Remark 3.1.** When (H2), (H3) or (H4) do not hold, then, adapting ideas from [36, 47], one may show that $s(D + \text{diag}(\sup_{x_i \ge 0} g_i(x_i)) - \text{diag}(\rho_i)) < 0$ is a sufficient condition for the origin of (7) to be $\alpha$-GAS with $\alpha := -s(D + \text{diag}(\sup_{x_i \ge 0} g_i(x_i)) - \text{diag}(\rho_i))$.
# 4 Minimal reachable equilibrium

In this section, we study the minimal globally asymptotically stable equilibrium $x^*(\rho)$ that may be attained for $\bar{\rho}$-admissible $\rho$. It is particularly relevant when $x^*(\rho)$ is positive for any $\bar{\rho}$-admissible $\rho$.

As will be seen shortly, this minimal value depends crucially upon the connectivity matrix<fnref n="1" /> $D|_{\overline{C_I}, \overline{C_I}}$ on the network of patches that are either uncontrollable or finitely controllable. Actually, the minimal value is related to the evolution of the residual system on $\mathbb{R}_+^{n_{\overline{C_I}}}$ defined below in (11). This system is obtained by removing all infinitely controllable patches, and its Jacobian matrix at the origin, called the *residual Jacobian matrix at the origin*, is precisely equal to $A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho})$.

When the underlying graph, called *residual graph*, is irreducible, the same analysis than before holds. On the contrary, when it is reducible, the situation is a little more complex and one has to consider its strongly connected components. These two cases are treated separately in the sequel, respectively in Theorems 4.1 and 4.2.

## 4.1 Irreducible residual graph

When $D|_{\overline{C_I}, \overline{C_I}}$ is irreducible, one has the following Theorem. A proof, as well as a proof of Theorem 4.2, is provided in Appendix D.

**Theorem 4.1.** Assume (H2), (H3) and (H4) hold and $D|_{\overline{C_I}, \overline{C_I}}$ is irreducible. Then the value $x^\infty := \inf\{x^*(\rho) : \rho \ \bar{\rho}\text{-admissible}\}$ fulfils

$$x^\infty = \lim_{\rho \to \bar{\rho}} x^*(\rho).$$

In particular,

$$x^\infty|_{C_I} = 0_{n_{C_I}},$$

and $x^\infty|_{C_I}$ is equal to the equilibrium attracting every non-zero trajectory of the so-called residual system on $\overline{C_I}$, defined as

$$\dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} \right) + \sum_{\substack{j \in \overline{C_I} \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in \overline{C_I} \\ j \neq i}} D_{ji}x_i, \quad i \in \overline{C_I}. \quad (11)$$

Consequently,

* if $s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho})) \le 0$, then

$$x^\infty|_{\overline{C_I}} = 0_{n_{\overline{C_I}}};$$

* if $s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho})) > 0$, then $x^\infty|_{\overline{C_I}}$ is equal to the unique positive equilibrium of the system (11).

<footnote n="1">¹This notation conforms to the convention defined in Definition 1.1.</footnote>
*Proof.* See Appendix D.

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

The value of $x^\infty$ in the infinitely controllable patches is null, as infinite mortality terms introduced in these patches lead to the elimination of the population. The population dynamics within the network is then studied in the patches that are either untreated or finitely controllable, considering that there is no migration from the infinitely controllable patches, due to the elimination of the population in these patches. The migration towards the infinitely controllable patches, as well as the maximum control allowed in $C_F$, produce additional mortality terms in the untreated patches and in the finitely controllable patches.

Recall that when $C_I$ is not empty and $s(A|_{C_I,C_I}(\bar{\rho})) \ge 0$, Theorem 3.4 states that $x^*(\rho) \gg 0_n$ for any $\rho$ $\bar{\rho}$-admissible. Consequently, in this case, the infimum $x^\infty$ is not a minimum. When $C_I$ is empty, $\bar{\rho}$ only takes finite values, and the infimum is a minimum, and the residual system (11) is exactly the system (7) for $\rho = \bar{\rho}$.

An example of irreducible residual graph is shown in Fig. 1. A case is treated numerically
in Section 6.

<figure><img src="image_3.png" /><figcaption>**Fig. 1: Determination of the minimal reachable equilibrium in the case of irreducible residual graph.** The infinitely controllable patches $C_I$, in which the minimal reachable equilibrium is null, are depicted in red. The blue patches are untreated, i.e. belong to $\bar{C}$, and the yellow patches correspond to the finitely controllable patches $C_F$. The red crosses indicate the interruption of migration from the infinitely treated patches to the remaining ones due to population elimination in the former. In the residual graph, the migration towards the infinitely controllable patches are represented by the red arrows. The dashed arrows indicate additional mortality terms corresponding to the maximum control allowed in $C_F$.</figcaption></figure>
## 4.2 Reducible residual graph

We will now consider the case where $D|_{\overline{\mathcal{C}_I, \mathcal{C}_I}}$ is reducible. As in the proof of Theorem 4.1 in Appendix D, determining the minimal globally asymptotically stable equilibrium of system (7) attainable for some $\bar{\rho}$-admissible $\rho$ involves examining the long-term dynamics of the solutions to the residual system (11). However, due to the reducibility of $D|_{\overline{\mathcal{C}_I, \mathcal{C}_I}}$, one cannot apply Theorem 2.2 to (11) for determining the asymptotic behavior of its solutions. Additional work is needed, which uses directed graph theory as central tool [4].

Introducing $\Gamma_{\overline{\mathcal{C}_I}}$, defined as the subgraph whose vertices are the patches belonging to $\overline{\mathcal{C}_I}$, we revert to the irreducible case by decomposing $\Gamma_{\overline{\mathcal{C}_I}}$ into its strongly connected components, which are irreducible subgraphs. Assume $\Gamma_{\overline{\mathcal{C}_I}}$ admits $N$ strongly connected components, where $1 \le N \le n_{\overline{\mathcal{C}_I}}$. The case $N = 1$ corresponds to the situation where $D|_{\overline{\mathcal{C}_I, \overline{\mathcal{C}_I}}}$ is irreducible, and $N = n_{\overline{\mathcal{C}_I}}$ corresponds to the case where $\Gamma_{\overline{\mathcal{C}_I}}$ is composed of $n_{\overline{\mathcal{C}_I}}$ strongly connected subgraphs, each composed of a single patch. It is a classical result that one may arrange the strongly connected components in an acyclic order, that is, into a sequence $G_1, \dots, G_N$, so that there is no arc (no direct path) from $G_j$ to $G_i$ unless possibly if $j < i$ [4, p. 17].

Let us now introduce adequate notations and definitions from graph theory.

**Definition 4.1.** For any subgraph or set of subgraphs $\mathcal{G}$ of $\Gamma_{\overline{\mathcal{C}_I}}$, the set and the number of its vertices are respectively denoted $V_{\mathcal{G}}$ and $n_{V_{\mathcal{G}}}$.

**Definition 4.2.** For any strongly connected component $\mathcal{G}$ of $\Gamma_{\overline{\mathcal{C}_I}}$, the set of its in-neighbouring (resp. out-neighbouring) subgraphs, denoted $\mathcal{G}^-$ (resp. $\mathcal{G}^+$), is the set of strongly connected components from which there exists a direct path leading to $\mathcal{G}$ (resp. that are reachable from $\mathcal{G}$ via a direct path), that is:

$$ \mathcal{H} \in \mathcal{G}^- \iff \text{there exists } i \in V_{\mathcal{H}} \text{ and } j \in V_{\mathcal{H}} \text{ such that } D_{ij} > 0. $$

$$ (\text{resp. } \mathcal{H} \in \mathcal{G}^+ \iff \text{there exists } j \in V_{\mathcal{G}} \text{ and } i \in V_{\mathcal{H}} \text{ such that } D_{ij} > 0.) $$

The in-degree (resp. out-degree) of $\mathcal{G}$ is the number of vertices in $\mathcal{G}^-$ (resp. $\mathcal{G}^+$).

Let us introduce the concept of upstream (resp. downstream) subgraphs. These notions are commonly discussed in literature related to graph theory, algorithms, and network analysis, though it might not always be labeled explicitly as ‘upstream’ (resp. ‘downstream’) subgraph. Instead, related concepts such as ‘reachability’ and ‘subgraph extraction’ are often used [11, 23, 41, 43].

**Definition 4.3.** For any strongly connected subgraph $\mathcal{G}$ of $\Gamma_{\overline{\mathcal{C}_I}}$, define the set of its upstream (resp. downstream) subgraphs, denoted $\mathcal{G}^{up}$ (resp. $\mathcal{G}^{down}$), as the set of strongly connected subgraphs from which there exists a path leading to $\mathcal{G}$ (resp. that are reachable from $\mathcal{G}$ via a path).

With these considerations, one may now introduce the following theorem, which extends Theorem 4.1 to the case of reducible residual graphs.
**Theorem 4.2.** Assume (H2), (H3) and (H4) hold. Then, the value

$$x^{\infty} := \inf\{x^{*}(\rho) : \rho \bar{\rho}\text{-admissible}\}$$

fulfills

$$x^{\infty} = \lim_{\rho \to \bar{\rho}} x^{*}(\rho).$$

In particular,

$$x^{\infty}|_{C_I} = 0_{n_{C_I}},$$

and $x^{\infty}|_{C_I}$ is equal to the equilibrium attracting every positive trajectory of the residual system (11).

Moreover, for any strongly connected component $\mathcal{G}$ of $\Gamma_{\bar{C}_I}$, the following system admits an equilibrium attracting every non-zero trajectory :

$$\dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I \cup V_{G+}} D_{ji} \right) + \sum_{\substack{j \in V_{G} \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in V_{G} \\ j \neq i}} D_{ji}x_i + \sum_{j \in V_{G-}} D_{ij}x_j^{\infty}, \quad i \in V_{G}, \quad (12)$$

and $x^{\infty}|_{V_G}$ is equal to this equilibrium. In particular, if there exists $j \in V_{G-}$ such that $\sum_{j \in V_{G-}} D_{ij}x_j^{\infty} > 0$, $x^{\infty}|_{V_G}$ is positive and GAS on $\mathbb{R}_+^{n_{V_G}}$ for the system (12).

Consequently,

* if $s(A|_{V_G, V_G}(\bar{\rho})) \le 0$ and if one of the following conditions is fulfilled

1. $\mathcal{G}$ is of in-degree zero;

2. $x^{\infty}|_{V_{G-}} = 0_{n_{V_{G-}}}$,

then,

$$x^{\infty}|_{V_G} = 0_{n_{V_G}};$$

* otherwise, $x^{\infty}|_{V_G} \gg 0_{n_{V_G}}$ and $x^{\infty}|_{V_G}$ is equal to the unique positive equilibrium of system (12).

*Proof.* See Appendix D. □

Notice that the first element of the list of the strongly connected components arranged in an acyclic order is always of in-degree zero. The equilibrium values on the strongly connected components of $\Gamma_{\bar{C}_I}$ are determined inductively, following an acyclic order. In particular, if $s(A|_{V_G, V_G}(\bar{\rho})) > 0$ and if one of the conditions 1. or 2. given above is fulfilled, then $x^{\infty}|_{V_G}$ is equal to the unique positive equilibrium of the following system on $V_G$:

$$\dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I \cup V_{G+}} D_{ji} \right) + \sum_{\substack{j \in V_G \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in V_G \\ j \neq i}} D_{ji}x_i, \quad i \in V_G.$$
Moreover, if $x^\infty|_{V_G} \gg 0_{n_{V_G}}$ for some strongly connected component $\mathcal{G}$, then, for any strongly connected component $\mathcal{H} \in \mathcal{G}^{\text{down}}$, the system (12) on $V_\mathcal{H}$ admits a unique equilibrium, which is positive, and $x^\infty|_{V_\mathcal{H}}$ is then equal to this equilibrium.

An example of reducible residual graph is illustrated in Fig. 2, and a numerical case is discussed below in Section 6.

**Remark 4.1.** Theorem 4.1 is a corollary of Theorem 4.2. Indeed, the latter theorem can be applied to the particular case where $D|_{\overline{C_I}, \overline{C_I}}$ is irreducible. In such case, $\Gamma_{\overline{C_I}}$ is composed of a single strongly connected component, which is $\Gamma_{\overline{C_I}}$ itself.
<figure><img src="image_4.png" /><figcaption>**Fig. 2:** Determination of the minimal reachable equilibrium in the case of reducible residual graph. In the underlying graph (left figure), the yellow patches are the patches of $C_F$ and the blue patches are untreated, i.e. belong to $\bar{C}$. The red patches are the infinitely controllable patches $C_I$, in which the minimal reachable equilibrium is null. The red crosses represent the interruption of migration from the infinitely treated patches to the remaining ones, due to population elimination in the former. In the residual system, the migration towards the infinitely controllable patches is represented by the red arrows. The additional mortality terms corresponding to the maximum control allowed in $C_F$ are represented by the dashed arrows. The residual graph is composed of three strongly connected components, $G_1$, $G_2$ and $G_3$ in acyclic order. There is one component of in-degree zero ($G_1$) and $V_{G_1}^{\text{down}} = \{V_{G_2}, V_{G_3}\}$. Moreover, $G_2^- = G_1$, $G_3^- = G_2 = G_1^+$ and $G_2^+ = G_3$. The right figure depicts the various scenarios possible for the equilibrium value on the residual graph</figcaption></figure>

# 5 A cost-effective approach for achievable elimination scenarios

In this section, we assume that Hypotheses (H2), (H3) and (H4) hold, and that the population can be driven to elimination by strong enough action on a subset $C$ of patches (in other terms, some of the conditions for elimination feasibility given in Theorem 3.4 are
fulfilled). We now build a cost-efficient strategy to achieve α-global asymptotic stability of the origin for given α, while minimizing a certain control cost. Here, this cost is taken as a weighted sum of the additional mortality terms ρᵢ introduced in the controllable patches, but any convex function of the latter would be convenient.

## 5.1 Minimization problem and properties

Let us first define the vector $π > 0_n$, where for every $i \in C$, $π_i > 0$ is the relative intervention price<fnref n="2" /> in patch i. For any α such that $0 \le α < -s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho}))$, let us introduce the following problem:

$$ \begin{array}{ll} \text{minimize} & f(\rho) \\ \text{subject to} & \rho \bar{\rho}\text{-admissible} \\ & h(\rho) \le -\alpha, \end{array} \qquad (13) $$

where

$$ f(\rho) := \pi^{\top}\rho, \quad h(\rho) := s(A(\rho)). $$

The function $f(\rho)$ represents the total cost corresponding to the $\bar{\rho}$-admissible control $\rho$. The constraint $h(\rho) \le -\alpha$ ensures α-global asymptotic stability of the origin (see Theorem 2.2).

Note that if $α > -s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho}))$, then the Problem (13) has no solution, since as a consequence of Theorem 3.1 and Lemma 3.3, there is no $\bar{\rho}$-admissible $ρ$ satisfying $h(ρ) \le -α$. The same holds if $α = -s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho}))$, except when $C_I = \emptyset$, in which case the minimization problem has a unique global minimizer, equal to $\bar{\rho}$.

Additionally, it is worth observing that (13) is a problem in $\mathbb{R}_+^{nc}$, since by definition of a $\bar{\rho}$-admissible vector, $ρ_i = 0$ for all $i \in \bar{C}$.

The solvability of problem (13) is addressed in the next result.

**Theorem 5.1.** Assume $s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})) < 0$ and let $α$ satisfy $0 \le α < -s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho}))$. Then Problem (13) admits local minimizers. The latter are also global and their set is convex. Moreover, if $C \subset \{1, ..., n\}$, the minimizer is unique. Last,

* if $s(A) > -α$, then any minimizer $ρ^*$ satisfies $ρ^* > 0_n$ and $s(A(ρ^*)) = -α$;

* if $s(A) \le -α$, the minimizer is unique and equal to $0_n$.

*Proof.* See Appendix E. □

## 5.2 Interior-point algorithm

We propose now an algorithm that solves numerically the Problem (13). It pertains to the class of interior point methods [9, Chapter 11], whose convergence properties have been the subject of many works [50]. The basic idea is to solve the minimization problem by applying Newton's method to a sequence of unconstrained penalized problems with decreasing penalization.

<footnote n="2">²In other words, the price vector fulfils: $π|_C \gg 0_{nc}$.</footnote>
Consider the following barrier problem:

$$ \begin{array}{ll} \text{minimize} & f_{\phi,t}(\rho) \\ \text{subject to} & \rho_i = 0, \ i \in \bar{C}, \end{array} \qquad (14) $$

where

$$ f_{\phi,t}(\rho) := tf(\rho) + \phi(\rho), \quad \phi(\rho) := -\log(-h(\rho) - \alpha) - \sum_{i \in C} \log(\rho_i) - \sum_{i \in C_F} \log(\bar{\rho}_i - \rho_i). $$

The function $\phi$ is a logarithmic barrier function. Its domain is the set of points that satisfy strictly the inequality constraints of (13). The solution of Problem (13) consists in solving a sequence of penalized unconstrained problems of the form (14), obtained while increasing the parameter $t$ at each step, and starting each Newton minimization at the numerical solution of the problem for the previous value of $t$. One must ensure that each solution of the Newton's method for each updated value of $t$ satisfies strictly the inequality constraints of Problem (13), so that the function $\phi$ remains well defined. Note that given the constraints imposed on $\alpha$ in Problem (13), it is always possible to find a $\rho$ that satisfies strictly these inequalities. Moreover, if $C \subset \{1, ..., n\}$, then the function $f_{\phi,t}$ is strictly convex since $h$ is strictly convex (see Theorem 3.1), which implies that for each $t > 0$, the problem (14) admits a unique minimizer.

In order to use Newton's method for the barrier problem (14), one integrates the equality constraints $\rho_i = 0$ for any $i \in \bar{C}$ by reducing the number of effective decision variables. To do so, we identify any $\bar{\rho}$-admissible vector $\rho \in \mathbb{R}_+^n$ with its components in $C$, that is, $\rho := (\rho_i)_{i \in C}$. Similarly, we represent $\pi \in \mathbb{R}_+^n$ using its components in $C$. Consequently, the functions $h = s(A(\cdot))$, $f$, $\phi$ and $f_{\phi,t}$ are now defined in $\mathbb{R}_+^{nc}$.

The proposed interior-point algorithm is given in pseudo-code in Algorithm 1.

**Algorithm 1** Interior-point method

**Require:** $t = t_0 > 0$, $\varepsilon > 0$, $\eta > 0$, $\mu > 1$ and $\rho^0 \gg 0_{nc}$ such that $h(\rho^0) < -\alpha$ and $\rho^0|_{C_F} \ll \bar{\rho}|_{C_F}$.
**while** $\frac{1}{t} > \varepsilon$ **do**
**while** $\|\nabla f_{\phi,t}(\rho^0)\| > \eta$ **do**
$\rho = \rho^0 - d(\rho^0) (\nabla^2 f_{\phi,t}(\rho^0))^{-1} \nabla f_{\phi,t}(\rho^0)$
$\rho^0 \leftarrow \rho$
**end while**
$t \leftarrow \mu t$
**end while**

The expressions $\nabla f_{\phi,t}(\rho^0)$ and $\nabla^2 f_{\phi,t}(\rho^0)$ denote respectively the gradient and the Hessian matrices of the function $f_{\phi,t}$. The parameter $\varepsilon$ is the desired level of accuracy for the approximation (14) of the Problem (13), $\eta$ is the level of accuracy of the solution of the barrier problem on each updated value $t$, and $\mu$ is the factor by which is increased the parameter $t$ at each step.
At the initial step, $\rho^0$ is chosen strictly feasible, in the sense that it satisfies strictly the inequality constraints of Problem (13). Then for each $\rho^0$ computed at a step of the Newton's method, the scalar $d(\rho^0) \in (0, 1]$ is chosen in order to ensure that the next iteration $\rho$ is also strictly feasible. Of course, one may take $d(\rho^0) = 1$ as long as $\rho^0 - (\nabla^2 f_{\phi,t}(\rho^0))^{-1} \nabla f_{\phi,t}(\rho^0)$ satisfies the inequality constraints.

In order to implement this algorithm, it is necessary to compute the gradient and the Hessian of the function $f_{\phi,t}$. For this, one takes advantage of the concept of group inverse of a matrix. The reader is referred to [33] for a thorough presentation of this notion.

**Definition 5.1.** Assume that $B \in \mathcal{M}_n(\mathbb{C})$ is a complex singular matrix, and that its eigenvalue 0 is semisimple, that is, its algebraic and geometric multiplicities coincide. Then, the group inverse of $B$, denoted $B^\#$, is the unique matrix $X \in \mathcal{M}_n(\mathbb{C})$ such that

$$ (i) \ BXB = B; $$

$$ (ii) \ XBX = X; $$

$$ (iii) \ BX = XB. $$

The following technical result shows how to compute the gradient and the Hessian of the function $f_{\phi,t}$. Recall that any $\bar{\rho}$-admissible $\rho$ as well as the vector $\pi$ now belong to $\mathbb{R}_+^{nc}$.

**Lemma 5.2.** Assume with no loss of generality that the control is introduced in the first $n_c$ patches, and among them, the infinite controllable patches are the first $n_{cI}$ patches. Moreover, for any $\rho \in \mathbb{R}_+^{nc}$, let $Q(\rho) := h(\rho)I_n - A(\rho)$, and denote $Q^\#(\rho)$ its group inverse. Then,

$$ \nabla f_{\phi,t}(\rho) = t\pi - \frac{1}{h(\rho) + \alpha} \nabla h(\rho) - \frac{1}{\rho} + \left( 0_{n_{cI}}, \frac{1}{\bar{\rho}|_{c_F} - \rho|_{c_F}} \right)^T, $$

where the gradient $\nabla h(\rho)$ is given by:

$$ \nabla h(\rho) = - \left( I_{n_c} - \text{diag} \left( (Q(\rho)Q^\#(\rho))_{11}, \dots, (Q(\rho)Q^\#(\rho))_{n_c n_c} \right) \right) 1_{n_c}, \quad (15) $$

where $1_{n_c}$ is the vector of $\mathbb{R}^{nc}$ whose components are all equal to 1. Moreover,

$$ \nabla^2 f_{\phi,t}(\rho) = \frac{1}{(h(\rho) + \alpha)^2} \nabla h(\rho) \nabla h(\rho)^\top - \frac{1}{h(\rho) + \alpha} \nabla^2 h(\rho) + \text{diag} \left( \frac{1}{\rho^2} \right) \\ + \text{diag} \left( 0_{n_{cI}}, \frac{1}{(\rho|_{c_F} - \bar{\rho}|_{c_F})^2} \right), $$

where the Hessian matrix $\nabla^2 h(\rho) \in \mathcal{M}_{n_c}(\mathbb{R})$ is such that, for any $i, j = 1, \dots, n_c$,

$$ (\nabla^2 h(\rho))_{ij} = (I_n - Q(\rho)Q^\#(\rho))_{ij} (Q^\#(\rho))_{ji} + (I_n - Q(\rho)Q^\#(\rho))_{ji} (Q^\#(\rho))_{ij}. \quad (16) $$

In Lemma 5.2, we denote by $\frac{1}{\rho}$, $\frac{1}{\bar{\rho}|_{c_F} - \rho|_{c_F}}$, etc. the vectors whose components are the scalar numbers $\frac{1}{\rho_i}$, $i \in \mathcal{C}$, $\frac{1}{\bar{\rho}_i - \rho_i}$, $i \in \mathcal{C}_F$, etc.

To compute the group inverse $Q^\#(\rho)$ of the matrix $Q(\rho) = s(A(\rho))I_n - A(\rho)$, one uses the formula provided in [33, Remark 2.5.3]. This formula is stated therein for matrices written in the form $r(A)I - A$, for $A$ a non-negative irreducible matrix and $r(A)$ its Perron value. One shows easily that it remains valid for an irreducible Metzler matrix $A$.
**Remark 5.1.** In the more general case where (H2), (H3) or (H4) do not hold, one may solve Problem (13) with $h(\rho) := s(A_r)$, $A_r := D + \text{diag}(\sup_{x_i \ge 0} g_i(x_i))$ and $0 < \alpha < -s(A_r|_{\overline{C_I}, \overline{C_I}}(\bar{\rho}))$. Theorem 5.1 still applies, and the interior-point algorithm in Algorithm 1 may be used, with gradient and Hessian as in Lemma 5.2, with now $Q(\rho)$ defined as

$$s(A_r - \text{diag}(\rho_i)) I_n - (A_r - \text{diag}(\rho_i)).$$

Moreover, notice that one does not set $\alpha = 0$, because $h(A_r) = 0$ may not guarantee global asymptotic stability of the origin (see Theorem A.1) in Appendix A.

# 6 Application to *Bactrocera dorsalis*

In this section, we present several numerical simulations to illustrate our results for a population of *Bactrocera dorsalis* evolving across interconnected orchards. This habitat can be represented by a graph, where each patch corresponds to a different orchard. This graph is directed because the propensity to go from an orchard to the other is not necessarily symmetric. The aim is to explore the impact of the different parameters – network structure, diffusion coefficients, intrinsic growth rates, intervention prices and parameter $\alpha$ – on the selection of the most efficient controllable patches and the corresponding optimal control strategy. As indicated in the introduction, we focus on control tools that increase the natural mortality rate, i.e. the use of sticky/deadly traps that attract and kill the fruit flies. We implicitly assume that the control death rate is proportional to the amount of traps that have been set-up in each orchard.

We assume that the network is composed of $n$ patches and we are limited to acting on a maximum of $n_C \le n$ patches. We modify Algorithm 1 so that it solves the minimization Problem (13), for each 'feasible' combination $C$ of $n_C$ patches. This results in a ranked list, arranged in ascending order of cost, of the appropriate combinations of $n_C$ patches along with the corresponding optimal control strategy. Our aim is to study the impact of the parameters of the model on the selection of the most efficient controllable patches. To better highlight this aspect, we suppose that one can add arbitrarily large mortality terms, that is, $C_I = C$ for any combination $C$ of $\{1, ..., n\}$.

For the intrinsic growth rate values, we rely on the data presented in [29, Table 2], [42, Table 3] and [52, Table 3], and consider a combination of numerical values, since orchards are not necessarily mono-specific but a combination of several host fruits. Thus the intrinsic growth rates range, on average, from 0.079 to 0.18. Concerning the carrying capacity of an orchard for *Bactrocera Dorsalis*, it is difficult to find data in literature. In [48], the male population was estimated between 3,300 and 18,000 per hectare. Therefore, we will set the carrying capacity for each patch between 6,000 and 20,000, depending on the size of the patch. Recall that the carrying capacity has no impact on the solution of Problem (13). When the problem is solvable, this quantity only affects the speed at which population converges towards the origin. When it is not solvable, it only affects the value of the minimal reachable positive equilibrium point.
We first consider in Section 6.1 a model with three patches, and explore the impact of the modification of several parameters. We then consider in Section 6.2 a more complex network with seven patches, in an attempt to illustrate the influence of the network configuration.

# 6.1 Three-patch model

We start with a graph containing three patches, each one connected to the others, thus forming a complete directed graph. We explore the effects of changing intrinsic growth rates and prices in Section 6.1.1, of modifying the prescribed convergence rate in Section 6.1.2, and of altering the diffusion coefficients in Section 6.1.3.

## 6.1.1 Intrinsic growth rate and prices effects

In order to examine the effects of the growth rate values and the prices, we first suppose that the diffusion is the same between every pair of patches. This implies in particular that the network is a symmetric directed graph, that is for any $i \neq j$, $D_{ij} = D_{ji}$. As far as we know, there is no data available in literature for migration parameters, but it is likely that ‘mark-release-recapture’ experiments could provide such estimates [31, 44].

For migration terms such that

$$D_{12} = D_{21} = D_{13} = D_{31} = D_{23} = D_{32} = 0.2, \qquad (17)$$

the algorithm outputs, depending on the number $n_c$ of controllable patches, are given in Table 1 for uniform prices across all patches and different intrinsic growth rates. Obviously, the more the controllable patches, the less the optimal cost associated to the optimal combination of patches. Moreover, when the intrinsic growth rate is the same in every patch, and given the rotation invariance of the graph along with identical prices across all patches, the optimal control for each patch is the same. One also notices that for different intrinsic growth rates in the patches and when $n_c = 1$ and $n_c = 2$, the optimal combinations comprise patches with highest growth rates, while the less efficient combination involves those with the lowest growth rates. Additionally, for each combination, the optimal strategy is to introduce a higher control in the patch with the higher intrinsic growth rate.

Moreover, we observe that when $n_c = n = 3$, the optimal control allocated in each patch is equal to the intrinsic growth rate in the patch. Notice that we ran numerical simulations with various symmetrical graphs, and always obtained experimentally that, when $C = \{1, ..., n\}$, then the optimal control in each patch is equal to the intrinsic growth rate.

When the intervention price in the third patch is 7 times more expensive than the others (Table 2), the optimal combination of patches for $n_c = 2$ is $\{1, 2\}$. For the combinations $\{1, 3\}$ and $\{2, 3\}$, almost all the control is allocated respectively to the patch 1 and the patch 2. When $n_c = 3$, all the control is allocated exclusively to patches 1 and 2.
<table><caption>**Table 1:** Optimal combination of patches and the corresponding optimal control and cost for diffusion coefficients in (17), $\alpha = 0$ and $\pi^T = (1, 1, 1)$</caption><thead><tr><th rowspan="3">Intrinsic growth rates in patches<br>{1,2,3}</th><th colspan="9">number n<sub>c</sub> of controllable patches</th></tr><tr><th colspan="3">1</th><th colspan="3">2</th><th colspan="3">3</th></tr><tr><th>Patch</th><th>Control</th><th>Cost</th><th>Patches</th><th>Control</th><th>Cost</th><th>Patches</th><th>Control</th><th>Cost</th></tr></thead><tbody><tr><td rowspan="3">(0.14, 0.14, 0.14)</td><td>1</td><td>1.07</td><td>1.07</td><td>{1,2}</td><td>(0.25,0.25)</td><td>0.50</td><td rowspan="3">{1,2,3}</td><td rowspan="3">(0.14,0.14,0.14)</td><td rowspan="3">0.42</td></tr><tr><td>2</td><td>1.07</td><td>1.07</td><td>{1,3}</td><td>(0.25,0.25)</td><td>0.50</td></tr><tr><td>3</td><td>1.07</td><td>1.07</td><td>{2,3}</td><td>(0.25,0.25)</td><td>0.50</td></tr><tr><td rowspan="3">(0.16, 0.14, 0.16)</td><td>1</td><td>1.37</td><td>1.37</td><td>{1,3}</td><td>(0.27, 0.27)</td><td>0.54</td><td rowspan="3">{1,2,3}</td><td rowspan="3">(0.16, 0.14, 0.16)</td><td rowspan="3">0.46</td></tr><tr><td>3</td><td>1.37</td><td>1.37</td><td>{1,2}</td><td>(0.293, 0.273)</td><td>0.57</td></tr><tr><td>2</td><td>1.74</td><td>1.74</td><td>{2,3}</td><td>(0.273, 0.293)</td><td>0.57</td></tr><tr><td rowspan="3">(0.12, 0.14, 0.16)</td><td>3</td><td>0.91</td><td>0.91</td><td>{2,3}</td><td>(0.226, 0.246)</td><td>0.47</td><td rowspan="3">{1,2,3}</td><td rowspan="3">(0.12, 0.14, 0.16)</td><td rowspan="3">0.42</td></tr><tr><td>2</td><td>1.09</td><td>1.09</td><td>{1,3}</td><td>(0.23, 0.27)</td><td>0.50</td></tr><tr><td>1</td><td>1.33</td><td>1.33</td><td>{1,2}</td><td>(0.253, 0.273)</td><td>0.53</td></tr></tbody></table>

<table><caption>**Table 2:** Optimal combination of patches and the corresponding optimal control and cost for diffusion coefficients in (17), $\alpha = 0$ and $\pi^T = (1, 1, 7)$</caption><thead><tr><th rowspan="3">Intrinsic growth rates in patches<br>{1,2,3}</th><th colspan="9">number n<sub>c</sub> of controllable patches</th></tr><tr><th colspan="3">1</th><th colspan="3">2</th><th colspan="3">3</th></tr><tr><th>Patch</th><th>Control</th><th>Cost</th><th>Patches</th><th>Control</th><th>Cost</th><th>Patches</th><th>Control</th><th>Cost</th></tr></thead><tbody><tr><td rowspan="3">(0.14, 0.14, 0.14)</td><td>1</td><td>1.07</td><td>1.07</td><td>{1,2}</td><td>(0.25,0.25)</td><td>0.50</td><td rowspan="3">{1,2,3}</td><td rowspan="3">(0.25,0.25,0)</td><td rowspan="3">0.50</td></tr><tr><td>2</td><td>1.07</td><td>1.07</td><td>{1,3}</td><td>(0.83,0.03)</td><td>1.02</td></tr><tr><td>3</td><td>1.07</td><td>7.51</td><td>{2,3}</td><td>(0.83,0.03)</td><td>1.02</td></tr><tr><td rowspan="3">(0.16, 0.14, 0.16)</td><td>1</td><td>1.37</td><td>1.37</td><td>{1,2}</td><td>(0.293, 0.273)</td><td>0.57</td><td rowspan="3">{1,2,3}</td><td rowspan="3">(0.293, 0.273, 0)</td><td rowspan="3">0.57</td></tr><tr><td>2</td><td>1.74</td><td>1.74</td><td>{1,3}</td><td>(0.85, 0.05)</td><td>1.18</td></tr><tr><td>3</td><td>1.37</td><td>9.57</td><td>{2,3}</td><td>(0.87, 0.07)</td><td>1.33</td></tr><tr><td rowspan="3">(0.12, 0.14, 0.16)</td><td>2</td><td>1.09</td><td>1.09</td><td>{1,2}</td><td>(0.253, 0.273)</td><td>0.53</td><td rowspan="3">{1,2,3}</td><td rowspan="3">(0.25, 0.27, 0)</td><td rowspan="3">0.53</td></tr><tr><td>1</td><td>1.33</td><td>1.33</td><td>{2,3}</td><td>(0.79, 0.03)</td><td>1.02</td></tr><tr><td>3</td><td>0.91</td><td>6.34</td><td>{1,3}</td><td>(0.81, 0.05)</td><td>1.14</td></tr></tbody></table>
### 6.1.2 Exponential convergence rate

In the previous examples, the constraint $\alpha = 0$ only required the elimination of the population. By setting $\alpha > 0$, we ensure that the elimination occurs at a positive exponential convergence rate, at least equal to $\alpha - \varepsilon$, for any $\varepsilon > 0$. A larger $\alpha$ leads to faster convergence towards zero at the expense of a higher optimal cost. This is illustrated in Fig. 3, for the growth functions

$$g_1(x_1) = 0.14(1 - \frac{x_1}{7,000}), \quad g_2(x_2) = 0.14(1 - \frac{x_2}{9,000}), \quad g_3(x_3) = 0.14(1 - \frac{x_3}{12,000}) \quad (18a)$$

and the initial condition

$$x_1(0) = 8,586.2, \quad x_2(0) = 9,032.5, \quad x_3(0) = 9,501.7. \quad (18b)$$

The latter has been chosen equal to the positive equilibrium point of the uncontrolled system (1), determined numerically for growth functions (18a) and diffusion coefficients (17). For all subsequent numerical simulations throughout the paper, the initial population values are selected to match the positive equilibrium of the respective uncontrolled system.

<figure><img src="image_5.png" /><figcaption>Fig. 3: Optimal cost for the Problem (13) and time for the population to fall below 1% of its initial value, depending on $\alpha$, for $n = 3$, $C = \{1, 2\}$, diffusion coefficients in (17), growth functions in (18a) and initial values in (18b). As $\alpha$ approaches $-s(A|_{\bar{C},\bar{C}})$, the optimal cost tends to infinity and the population drops rapidly below 1% of its initial value.</figcaption></figure>

As $\alpha$ approaches $-s(A|_{\bar{C},\bar{C}})$, the control introduced in the patches within $C$ tends towards infinity, leading to an immediate elimination of the population in those patches. However, in the uncontrolled patches $\bar{C}$, the asymptotic behaviour of the population approaches the asymptotic behaviour of the residual system on $\bar{C}$. Consequently, no matter how large the
control in $C$, the origin in the patches in $\bar{C}$ is at best $-s(A|_{\bar{C},\bar{C}})$-GAS. This behavior is illustrated in Fig. 4 for a value of $\alpha$ near $-s(A|_{\bar{C},\bar{C}})$.

<figure><img src="image_6.png" /><figcaption>Fig. 4: Population trajectories for optimal control $(\rho_1, \rho_2) = (800.2, 800.2)$ introduced at Time = 50 in $C = \{1, 2\}$ associated with $\alpha = -s(A_{\bar{C},\bar{C}}) - 0.0001 = 0.26$. The growth functions are defined in (18a), initial conditions in (18b) and diffusion coefficients in (17).</figcaption></figure>

### 6.1.3 Diffusion effects

Optimal solutions for some diffusion coefficients values are presented in Table 3. When the first orchard is further apart and more difficult to reach, so that for example $D_{13} = D_{12} = 0.001$ $^3$, then it is not possible to eliminate the population controlling only patch 1. As a matter of fact, it is useless to put a control in this patch, even when $n_C \ge 2$. An interpretation is that when the population in patch 1 is reduced to an extremely low level by introducing a sufficient control in this patch, the population dynamics in patches 2 and 3 becomes nearly independent of patch 1, since the migration from and towards this patch is negligible. In these conditions, the control introduced in the patch 1 will have negligible effect on patches 2 and 3. The system for patches 2 and 3 then approaches the following:

$$ \dot{x}_2 = 0.14(1 - \frac{x_2}{9,000})x_2 + 0.5x_3 - 0.5x_2, \quad \dot{x}_3 = 0.14(1 - \frac{x_3}{12,000})x_3 + 0.5x_2 - 0.5x_3. $$

The stability modulus $s(A)$ of this system is equal to $0.14 > 0$, which implies that the population persists. It is then more efficient to act directly on the patches 2 and 3. Moreover, the optimal cost when $n_C = 2$ (resp. $n_C = 2$), equal to 0.28, is 44% (resp. 33%) lower compared to the case where all diffusion parameters are equal to 0.2 (see Table 1). In

<footnote>$^3$We avoid setting $D_{13} = D_{12} = 0$ to maintain the graph's irreducibility.</footnote>
practice, this suggests that integrating mass trapping with measures to prevent migration towards some specific orchards as part of an integrated pest management strategy, would beneficially reduce the overall costs by increasing pest mortality.

<table><caption>**Table 3:** Optimal combination of patches and the corresponding optimal control and total cost for $\alpha = 0$ and $\pi^T = (1, 1, 1)$, and intrinsic growth rates $(r_1, r_2, r_3) = (0.14, 0.14, 0.14)$.</caption><thead><tr><th rowspan="3">Diffusion parameters</th><th colspan="9">number n<sub>C</sub> of controllable patches</th></tr><tr><th colspan="3">1</th><th colspan="3">2</th><th colspan="3">3</th></tr><tr><th>Patch</th><th>Control</th><th>Cost</th><th>Patches</th><th>Control</th><th>Cost</th><th>Patches</th><th>Control</th><th>Cost</th></tr></thead><tbody><tr><td>D<sub>21</sub> = D<sub>23</sub> = D<sub>31</sub> = D<sub>32</sub> = 0.2,<br>D<sub>13</sub> = D<sub>12</sub> = 0.001</td><td>2</td><td>0.61</td><td>0.61</td><td>{2, 3}</td><td>(0.14,0.14)</td><td>0.28</td><td>{1,2,3}</td><td>(0,0.14,0.14)</td><td>0.28</td></tr><tr><td></td><td>3</td><td>0.61</td><td>0.61</td><td>{1,2}</td><td>(0,0.61)</td><td>0.61</td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>{1,3}</td><td>(0,0.61)</td><td>0.61</td><td></td><td></td><td></td></tr><tr><td>D<sub>23</sub> = D<sub>13</sub> = D<sub>21</sub> = 0.08<br>D<sub>32</sub> = D<sub>31</sub> = 0.02, D<sub>12</sub> = 0.12</td><td></td><td></td><td></td><td>{1,2}</td><td>(0.299, 0.259)</td><td>0.558</td><td>{1,2,3}</td><td>(0.178, 0.138, 0.060)</td><td>0.377</td></tr></tbody></table>

For the second set of diffusion parameters, one sees from Table 3 that the population cannot be eliminated if $n_C = 1$. When $n_C = 2$, only one specific combination, patches 1 and 2, guarantees elimination. Therefore, if the control is not allowed in the patch 1 or 2, the control fails. For instance when a control is introduced in patches 2 and 3, the equilibrium point in patch 1 will not fall below a certain threshold despite considerable efforts, as illustrated in Fig. 5 for initial conditions

$$x_1(0) = 9446.4, \quad x_2(0) = 8548.4, \quad x_3(0) = 4762.7. \qquad (19)$$

This is a consequence of Theorem 4.1. Indeed, in this scenario, $\bar{C} = \{1\}$, so the residual system is $\dot{x}_1 = x_1(g_1(x_1) - D_{21} - D_{31})$, with $s(A_{\bar{C},\bar{C}}) = g_1(0) - D_{21} - D_{31} = 0.14 - 0.08 - 0.02 > 0$, which proves the existence of a positive equilibrium threshold $x_1^\infty$ in patch 1. The quantity $x_1^\infty$ is equal to the positive equilibrium of the latter system, i.e. $x_1^\infty = 2,000$ when $g_1(x_1) = 0.14(1 - \frac{x_1}{7,000})$. Moreover, introducing control in all patches results in an optimal cost that is 32% lower than when it is introduced only in patches 1 and 2. This highlights the fact that some patches may play a critical role, the impossibility to control them directly resulting in substantial cost increases or even control failure.
<figure><img src="image_7.png" /><figcaption>Fig. 5: Population trajectories for optimal control $(\rho_1^*, \rho_2^*) = (0.399, 0.259)$ introduced at Time = 50 in patches $\{1, 2\}$ (left), and for large control $(\rho_2, \rho_3) = (1000, 1000)$ in patches $\{2, 3\}$ (right), with $\alpha = 0$, growth functions in (18a), initial conditions in (19) and second set of diffusion coefficients in Table 3. The equilibrium for the right figure is equal to $(2000.56, 0.16, 0.04)$.</figcaption></figure>

## 6.2 Seven-patch model and the impact of the network configuration

Consider now a higher dimension example, involving seven patches. In this section, our aim is to observe the impact of the network spatial structure. To achieve this, we enforce uniform pricing and growth rates across all patches, and only require elimination ($\alpha = 0$). Additionally, we ensure that the diffusion coefficients between patches that are connected are all equal to 0.2. The parameter values are summarized as follows:

$$ \pi^T = (1, 1, 1, 1, 1, 1), \quad g_i(0) = 0.14, \quad i = 1, \dots, 7, \quad \alpha = 0, \quad D_{ij} = 0.2 \text{ if } D_{ij} > 0. \quad (20) $$

As a comparison, when $C = \{1, \dots, 7\}$ and the directed graph is symmetric, then the best strategy consists in introducing mortality rates equal to the intrinsic growth rates in all seven patches, resulting in an optimal cost of $7 \times 0.14 = 0.98$.

We examine in the following three different examples of networks with seven patches (some of these networks being possibly identical, up to vertex permutation). We assume the intervention is limited to three patches (for example because of limited number of traps or because only three orchards owners agree to set-up traps in their fields). Notice that there are 35 possible ways of selecting three patches out of seven.

### Example 1

Let us begin with the example illustrated in Fig. 6, where lines connecting two patches illustrate bidirectional and symmetric migration.
<figure><img src="image_8.png" /><figcaption>**Fig. 6:** A network configuration with seven patches</figcaption></figure>

The algorithm outputs that out of the 35 possible combinations of three patches, 32 can lead to elimination. The optimal combination is patches {2, 3, 7} as illustrated in Fig. 7, with optimal control (0.416, 0.416, 0.416), and optimal cost 1.25. The optimal control is the same in each patch 2, 3 and 7 because we may permute them without changing the structure of the graph, albeit with different numbering.

For the second best combination of three patches, there are three possible configurations, namely {1, 4, 7}, {2, 4, 6} and {3, 4, 5}. This lack of uniqueness also arises because permutations are possible among patches 1, 5, 6, and also among patches 2, 3, 7, without altering the graph structure (albeit with different numbering). For the same reason, $(\rho_1^*, \rho_4^*, \rho_7^*) = (\rho_6^*, \rho_4^*, \rho_2^*) = (\rho_5^*, \rho_4^*, \rho_3^*)$. Additionally, the algorithm indicates that the optimal cost for the second-best combination is 1.53, which is 23% higher than the optimal cost for the best combination of patches.

<figure><img src="image_9.png" /><figcaption>**Fig. 7:** Optimal combination (green) and second-best combinations (orange) of three patches for parameters in (20).</figcaption></figure>
## Example 2

Let us now assume that there is no diffusion between patch 4 and patch 7 (due to landscape configuration, for instance). This tends to lower the number of possible combinations capable of achieving elimination, and the number of possible combinations indeed drops from 32 to 24. Notably, the optimal combination remains unchanged, as illustrated in Fig. 8, but it comes with an optimal cost of 1.28, which is 3% higher than for Example 1. Moreover, the second-best configuration is now unique.

<figure><img src="image_10.png" /><figcaption>**Fig. 8:** Optimal combination (green) and second-best combination (orange) of three patches for parameters in (20). Unlike in Fig. 7, there is no diffusion between patches 4 and 7.</figcaption></figure>

Assume now that there is no diffusion from patch 5 to patch 7 and from patch 7 to patch 6, as illustrated in Fig. 9. The algorithm indicates that there are 23 possible combinations of patches that can lead to elimination with appropriate control. The best combination is {2, 3, 5} with an optimal total cost of 1.11. Notably, the algorithm determines that it is impossible to eliminate the population when $C = \{1, 3, 4\}$.

Let us assess the minimal attainable positive equilibrium point for such $C$. The residual graph is reducible, and composed of three strongly connected components, noted in acyclic order $G_1, G_2$ and $G_3$, where $G_1$ is the graph composed of the unique patch 6, $G_2$ the graph composed of the patch 7, and $G_3$ the graph composed of the patches 2 and 5 and the diffusion between them, as illustrated in Fig. 10a. We apply Theorem 4.2. One has $A|_{V_{G_1}, V_{G_1}} = g_6(0) - D_{36} - D_{76} = -0.26 \le 0$, which implies $x_6^\infty = x^\infty|_{V_{G_1}} = 0$. Similarly, $A|_{V_{G_2}, V_{G_2}} = g_7(0) - D_{57} = -0.06 \le 0$, leading to $x_7^\infty = x^\infty|_{V_{G_2}} = 0$. Next, we examine $A|_{V_{G_3}, V_{G_3}}$ which is given by $A|_{V_{G_3}, V_{G_3}} = \begin{pmatrix} g_2(0) - D_{12} - D_{42} - D_{52} & D_{25} \\ D_{52} & g_5(0) - D_{25} \end{pmatrix} = \begin{pmatrix} -0.46 & 0.2 \\ 0.2 & -0.06 \end{pmatrix}$.

Therefore, $s(A|_{V_{G_3}, V_{G_3}}) = 0.03 > 0$, which implies $(x_2^\infty, x_5^\infty) = x^\infty|_{V_{G_3}} \gg 0_2$. Moreover, still by Theorem 4.2, one has that $(x_2^\infty, x_5^\infty)$ is the unique positive equilibrium of the system

$$ \dot{x}_2 = x_2 (g_2(x_2) - D_{12} - D_{42}) + D_{25}x_5 - D_{52}x_2, \quad \dot{x}_5 = x_5g_5(x_5) + D_{52}x_2 - D_{25}x_5. $$

By applying Newton-Raphson method, one finds numerically that $(x_2^\infty, x_5^\infty) = (469.52, 1097.04)$ for the growth functions

$$ g_i(x_i) = 0.14(1 - \frac{x_i}{K_i}), \quad i = 1, \dots, 7 \qquad (21a) $$
and carrying capacities

$$ \begin{array}{l} K_1 = 7,000, \quad K_2 = 9,000, \quad K_3 = 12,000, \quad K_4 = 16,000, \\ K_5 = 6,000, \quad K_6 = 13,000, \quad K_7 = 19,000. \end{array} \qquad (21b) $$

An illustration when large mortality terms are introduced in the patches {1, 3, 4} is provided in Fig. 10b for the initial condition

$$ \begin{array}{l} x_1(0) = 8,697.3, \quad x_2(0) = 10,0024.5, \quad x_3(0) = 8,846.3, \quad x_4(0) = 10,678.5, \\ x_5(0) = 11,496.7, \quad x_6(0) = 5,535.6, \quad x_7(0) = 8,844.8. \end{array} \qquad (22) $$

<figure><img src="image_11.png" /><figcaption>**Fig. 9:** Network configuration with seven patches. The arrows indicate unidirectional migration.</figcaption></figure>
Residual graph

<figure><img src="image_12.png" /><figcaption>(a) Residual graph when $\mathcal{C} = \mathcal{C}_I = \{1, 3, 4\}$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Population trajectory when $\mathcal{C} = \{1, 3, 4\}$ and $(\rho_1, \rho_3, \rho_4) = (1000, 1000, 1000)$.</figcaption></figure>

**Fig. 10:** Residual graph (a) and population trajectories when large linear mortality terms are introduced at Time = 50 in patches 1,3 and 4 (b) for parameters in (20), (21a) and (21b). The positive equilibrium point for (b) is equal to (0.094, 469.97, 0.0000376, 0.094, 1097.89, 0.00002.89, 0.0000964).

Example 3

Let us now revert to the initial configuration of Example 1, and assume there is no migration between patches 2 and 5, 3 and 4, and 6 and 7, to create a chain structure. This type of structure occurs, for example, when the orchards are positioned side by side. In such cases, it is likely that insects will migrate from one orchard to another sequentially, without skipping any.

The algorithm outputs that 8 combinations of three patches can lead to the population elimination. The optimal combination is depicted in Fig. 11, and the associated optimal total cost increases to 1.78, which is 39% higher than the optimal cost in Example 2 (1.28).
<figure><img src="image_14.png" /><figcaption>**Fig. 11:** Optimal combination (green) and second-best combinations (orange) of three patches for parameters in (20). Unlike in Fig. 7, there is no diffusion between patches 2 and 5, 6 and 7, and 3 and 4.</figcaption></figure>

# 7 Conclusion

Having in mind the issue of control of insects vectors or insects pests, we considered in this paper a metapopulation model with patches linearly interconnected, and explored the global effects of (on purpose) increase of the mortality in some of these patches. The patches are classified in three categories, depending on whether the local value of the mortality rate may be increased or not (controllable or uncontrollable patches); and in the former case, whether the additional mortality rate may be chosen arbitrarily large, or equal at most to a given quantity (infinitely controllable or finitely controllable patches). Typically this setting is supposed to replicate mass trapping campaigns aimed at eliminating the population in an heterogeneous environment where the control application is submitted locally to restrictions.

With the help of results by Y. Takeuchi et al., we showed that the sign of the stability modulus of the Jacobian at the origin determines the asymptotic behaviour of such model. If it is non-positive, then the population becomes extinct in every patch; while if it is positive, then there exists a unique non-negative and non-zero equilibrium, which is positive and globally asymptotically stable.

Based on this preliminary result, the possibility of eliminating the population at a given exponential convergence rate by adequate choice of the mortality rates was investigated. It turns out that elimination feasibility depends upon the value of the stability modulus of the Jacobian matrix at the origin of some *residual system*, obtained by removing the infinitely controllable patches and by adjusting the mortality rates in the finitely controllable patches
at their maximum value allowed. When the stability modulus of the residual system is positive, then the best that can be done is to lower the positive equilibrium to some minimum value. When it is negative, elimination is possible, and we showed how the choice of the mortality rates in the controllable patches may be done in order to achieve elimination while minimizing some control cost.

We then applied these results numerically, using life cycle characteristics of *Bactrocera dorsalis*. We investigated scenarios with different intrinsic growth rates and migration parameters, and highlighted the corresponding effects on the optimal elimination strategies. Our findings underscore the importance of tailoring elimination strategies to seasonal variations, which influence intrinsic growth rates and diffusion coefficients. For example, one could consider the effect of the temperature on the intrinsic growth rate, as done in [20, Table 5].

We also studied the effects induced by modification of the network connectivity. We showed by examples that some patches play a critical role, and that the impossibility to control them directly may lead to significant cost increase, or even control failure. This points out the importance of considering early engagement, thorough explanation, and detailed information dissemination to landowners who are hesitant about intervention. We believe that mathematical modelling and simulations can be of great help to design, size and schedule control strategies, but also to convince reluctant people.

Finally, our work stresses the crucial need of estimation methods for the migration rates in the field in order to develop effective control strategies. Techniques such as Mark-Release-Recapture [44, 31] can be employed to measure the flux between the target areas. While such experiments are lengthy, challenging, and not always successful, they are essential for minimizing the risk of control failure.

This paper may serve as basement for further investigations on vector/pest control tools.

## Acknowledgments

This work is partly supported by the AttracTIS project funded by ECOPHYTO 2021-2022: "Construire avec les outre-mer une agroécologie axée sur la réduction de l'utilisation, des risques et des impacts des produits phytopharmaceutiques". PAB and MdT acknowledge the support of the European Regional Development Fund (ERDF) while visiting 3P-Cirad (Saint-Pierre, La Réunion, France) in March 2024, during the preparation of this work. YD is (partially) supported by the DST/NRF SARChI Chair in Mathematical Models and Methods in Biosciences and Bioengineering at the University of Pretoria (Grant 82770). YD acknowledges the support of the Conseil Régional de la Réunion, the Conseil départemental de la Réunion, the European Regional Development Fund (ERDF), and the Centre de Coopération Internationale en Recherche Agronomique pour le Développement (CIRAD). Last, the financial support of Inria is acknowledged, through the Associate Team MoCoVec 'Modelling and Biological Control of Vector-Borne Diseases: the case of Malaria and Dengue'; and the financial support of the STIC AmSud program, through the project 23-STIC-02 BIO-CIVIP 'Biological control of insects vectors and insects pests'.
References

[1] J. Arino. Diseases in metapopulations. In *Modeling and dynamics of infectious diseases*, pages 64–122. World Scientific, 2009.

[2] J. Arino, N. Bajeux, and S. Kirkland. Number of source patches required for population persistence in a source-sink metapopulation with explicit movement. *Bulletin of Mathematical Biology*, 81:1916–1942, 2019.

[3] R. Baker, G. Gilioli, C. Behring, D. Candiani, A. Gogin, T. Kaluski, M. Kinkar, O. Mosbach-Schulz, F.M. Neri, S. Preti, et al. Bactrocera dorsalis – Pest Report and Datasheet to support ranking of EU candidate priority pests, 2019.

[4] J. Bang-Jensen and G. Z. Gutin. *Digraphs: theory, algorithms and applications*. Springer Science & Business Media, 2008.

[5] J.C. Beier, G.C. Müller, W. Gu, K.L. Arheart, and Y. Schlein. Attractive toxic sugar bait (ATSB) methods decimate populations of Anopheles malaria vectors in arid environments regardless of the local availability of favoured sugar-source blossoms. *Malaria journal*, 11:1–7, 2012.

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

[7] D.I. Bertsekas, A. Nedic, and A. Ozdaglar. *Convex analysis and optimization*, volume 1. Athena Scientific, 2003.

[8] P.-A. Bliman, D. Cardona-Salgado, Y. Dumont, and O. Vasilieva. Implementation of control strategies for sterile insect techniques. *Mathematical biosciences*, 314:43–60, 2019.

[9] S. Boyd and L. Vandenberghe. *Convex optimization*. Cambridge university press, 2004.

[10] A.R. Clarke, K.F. Armstrong, A.E. Carmichael, J.R. Milne, S. Raghu, G.K. Roderick, and D.K. Yeates. Invasive phytophagous pests arising through a recent tropical evolutionary radiation: the Bactrocera dorsalis Complex of Fruit Flies. *Annu. Rev. Entomol.*, 50:293–319, 2005.

[11] T.H. Cormen, C.E. Leiserson, R.L. Rivest, and C. Stein. *Introduction to algorithms*. MIT press, 2022.

[12] D. Cugala, V. Massimiliano, M. Maulid, M. De Meyer, and L. Canhanga. Economic injury level of the Oriental fruit fly, Bactrocera dorsalis (Diptera: Tephritidae), on commercial mango farms in Manica Province, Mozambique. *African entomology*, 28(2):278–289, 2020.
[13] E. Deutsch and M. Neumann. Derivatives of the Perron root at an essentially nonnegative matrix and the group inverse of an *M*-matrix. *Journal of Mathematical Analysis and Applications*, 102(1):1–29, 1984.

[14] V.A. Dyck, J. Hendrichs, and A.S. Robinson. *Sterile Insect Technique: Principles and Practice in Area-Wide Integrated Pest Management*. Taylor & Francis, 2021.

[15] A.E. Eiras, T.S. Buhagiar, and S.A. Ritchie. Development of the gravid Aedes trap for the capture of adult female container-exploiting mosquitoes (Diptera: Culicidae). *Journal of medical entomology*, 51(1):200–209, 2014.

[16] B. Elbetch. Generalized Coupled Source-Sink Models. https://hal.science/hal-04167884/file/hal_version.pdf, July 2023.

[17] B. Elbetch. Generalized logistic equation on networks. *Comptes Rendus. Mathématique*, 361(G5):911–934, 2023.

[18] B. Elbetch, T. Benzekri, D. Massart, and T. Sari. The multi-patch logistic equation. *Discrete and Continuous Dynamical Systems-Series B*, 26(12):6405–6424, 2021.

[19] A. Fereres and B. Raccah. Plant virus transmission by insects. 2015.

[20] K.K.M. Fiaboe, S. Kekeunou, S.N. Nanga, A.F. Kuate, H.E.Z. Tonnang, D. Gnanvos-sou, R. Hanna, et al. Temperature-based phenology model to predict the development, survival, and reproduction of the oriental fruit fly Bactrocera dorsalis. *Journal of Thermal Biology*, 97:102877, 2021.

[21] S. Friedland. Convex spectral functions. *Linear and multilinear algebra*, 9(4):299–316, 1981.

[22] S. A. Gershgorin. Über die abgrenzung der eigenwerte einer matrix. *Izv. akad.Nauk. S.S.S.R.*, (6):749–754, 1931.

[23] A.V. Goldberg, É. Tardos, and R. Tarjan. Network flow algorithm. Technical report, Cornell University Operations Research and Industrial Engineering, 1989.

[24] N.G. Gratz. Critical review of the vector status of Aedes albopictus. *Medical and veterinary entomology*, 18(3):215–227, 2004.

[25] I. Hanski. Metapopulation dynamics. *Nature*, 396(6706):41–49, 1998.

[26] L. Heringer, B.J. Johnson, K. Fikrig, B.A. Oliveira, R.D. Silva, M. Townsend, R. Barrera, Á.E. Eiras, and S.A. Ritchie. Evaluation of Alternative Killing Agents for Aedes aegypti (Diptera: Culicidae) in the Gravid Aedes Trap (GAT). *Journal of Medical Entomology*, 53(4):873–879, 2016.

[27] G. Hess. Disease in metapopulation models: implications for conservation. *Ecology*, 77(5):1617–1632, 1996.
[28] R.A. Horn and C.R. Johnson. *Topics in matrix analysis*. Cambridge university press, 1994.

[29] Y.B.K. Huang and H. Chi. Fitness of Bactrocera dorsalis (Hendel) on seven host plants and an artificial diet. *Turk. Entomol. Derg-Tu*, 38:401–414, 2014.

[30] A. Jaffal, J. Fite, T. Baldet, P. Delaunay, F. Jourdain, R. Mora-Castillo, M.M. Olive, and D. Roiz. Current evidences of the efficacy of mosquito mass-trapping interventions to reduce Aedes aegypti and Aedes albopictus populations and Aedes-borne virus transmission. *PLOS Neglected Tropical Diseases*, 17(3):e0011153, 2023.

[31] E.B. Jang, R.V. Dowell, and N.C. Manoukis. Mark-release-recapture experiments on the effectiveness of methyl eugenol-spinosad male annihilation technique against an invading population of Bactrocera dorsalis. 2017.

[32] E. Kaszkurewicz and A. Bhaya. *Matrix diagonal stability in systems and computation*. Springer Science & Business Media, 2012.

[33] S. J. Kirkland and M. Neumann. *Group inverses of M-matrices and their applications*. CRC Press, 2012.

[34] P. Kokotović, H.K. Khalil, and J. O'reilly. *Singular perturbation methods in control: analysis and design*. SIAM, 1999.

[35] R. Levins. Some demographic and genetic consequences of environmental heterogeneity for biological control. *Bulletin of the ESA*, 15(3):237–240, 1969.

[36] Z. Lu and Y. Takeuchi. Global asymptotic behavior in single-species discrete diffusion systems. *Journal of Mathematical Biology*, 32(1):67–77, 1993.

[37] A.J. Mackay, M. Amador, and R. Barrera. An improved autocidal gravid ovitrap for the control and surveillance of Aedes aegypti. *Parasites & vectors*, 6:1–13, 2013.

[38] G.C. Müller, J.C. Beier, S.F. Traore, M.B. Toure, M.M. Traore, S. Bah, S. Doumbia, and Y. Schlein. Successful field trial of attractive toxic sugar bait (ATSB) plant-spraying methods against malaria vectors in the Anopheles gambiae complex in Mali, West Africa. *Malaria journal*, 9:1–7, 2010.

[39] Y. Nesterov and V.Y. Protasov. Computing closest stable nonnegative matrix. *SIAM Journal on Matrix Analysis and Applications*, 41(1):1–28, 2020.

[40] R.D. Nussbaum. Convexity and log convexity for the spectral radius. *Linear Algebra and its Applications*, 73:59–122, 1986.

[41] S. Parthasarathy, S. Tatikonda, and D. Ucar. A survey of graph mining techniques for biological datasets. *Managing and mining graph data*, pages 547–580, 2010.
[42] W. Pieterse, A. Manrakhan, J.S. Terblanche, and P. Addison. Comparative demography of Bactrocera dorsalis (Hendel) and Ceratitis capitata (Wiedemann) (Diptera: Tephritidae) on deciduous fruit. *Bulletin of entomological research*, 110(2):185-194, 2020.

[43] R. Sedgewick and K. Wayne. *Algorithms*. Addison-Wesley professional, 2011.

[44] T.E. Shelly and J. Edu. Mark-release-recapture of males of Bactrocera cucurbitae and B. dorsalis (Diptera: Tephritidae) in two residential areas of Honolulu. *Journal of Asia-Pacific Entomology*, 13(2):131-137, 2010.

[45] S. Singh and K. Bajaj. Mass Trapping of Fruit Flies Using Methyl Eugenol Based Traps. In *Abiotic and Biotic Stress Management in Plants*, pages 129-153. CRC Press, 2022.

[46] H.L. Smith. *Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems*. Number 41. American Mathematical Soc., 1995.

[47] Y. Takeuchi. Cooperative systems theory and global stability of diffusion models. *Acta Applicandae Mathematica*, 14, 1989.

[48] K.H. Tan and M. Serit. Adult population dynamics of Bactrocera dorsalis (Diptera: Tephritidae) in relation to host phenology and weather in two villages of Penang Island, Malaysia. *Environmental Entomology*, 23(2):267-275, 1994.

[49] L.N. Trefethen and D. Bau. *Numerical linear algebra*. SIAM, 2022.

[50] M. Wright. The interior-point revolution in optimization: history, recent developments, and lasting consequences. *Bulletin of the American mathematical society*, 42(1):39-56, 2005.

[51] X. Zheng, D. Zhang, Y. Li, C. Yang, Y. Wu, X. Liang, Y. Liang, X. Pan, L. Hu, Q. Sun, et al. Incompatible and sterile insect techniques combined eliminate mosquitoes. *Nature*, 572(7767):56-61, 2019.

[52] Y. Zhu, F. Qi, X. Tan, T. Zhang, Z. Teng, Y. Fan, F. Wan, and H. Zhou. Use of Age-Stage, Two-Sex Life Table to Compare the Fitness of Bactrocera dorsalis (Diptera: Tephritidae) on Northern and Southern Host Fruits in China. *Insects*, 13(3), 2022.
# Appendix – Demonstrations

## A Fundamental stability result (Thm. 2.2)

First, let us define the following matrix.

$$A_r := D + \text{diag}(r_i), \qquad (A.23)$$

where for every $i = 1, ..., n$,

$$r_i := \sup_{x_i \ge 0} g_i(x_i).$$

In particular, when (H2) holds, $A_r = A$, and $s(A_r)$ is simply written $s(A)$.

To prove Theorem 2.2, let us introduce the following result which summarizes the asymptotic behaviour of the model (1), depending on Assumptions (H2), (H3) and (H4).

### Theorem A.1.

(i) If $s(A_r) \le -\alpha$ for some $\alpha > 0$, then the origin is $\alpha$-GAS on $\mathbb{R}_+^n$ for system (1). Moreover, if (H2) holds, then this statement is still true for $\alpha = 0$

(ii) Assume (H4) holds. If $s(A) > 0$, then any solution $x$ of the model (1) initiated in a point $x_0 > 0_n$ satisfies

$$\liminf_{t \to +\infty} x(t) \ge 0_n.$$

Furthermore, if (H3) also holds, then the model (1) admits at least one non-negative and non-zero equilibrium $x^*$, which is necessarily positive.

Finally, if in addition to (H3) and (H4), (H2) also holds, then $x^*$ is unique and GAS on $\mathbb{R}_+^n \setminus \{0_n\}$ for the system (1).

*Proof.* The proofs of points (i) for the global asymptotic stability of the origin and (ii) may be obtained by adapting those provided by Takeuchi for the model (5) in [36, 47], under the assumption that (H1), (H2), (H3) and (H4) hold. (In Takeuchi's work, these assumptions are merged.)

Let us now prove the $\alpha$-global asymptotic stability of the origin in point (i) for some $\alpha > 0$. Assume $s(A_r) \le -\alpha < 0$. By definition of $A_r$ in (A.23), $F(x) \le A_r x$. As the linear system $\dot{x} = A_r x$ is cooperative [46], it follows that for any initial condition $x_0 \ge 0_n$, the solution $x(t)$ of the system (1) is upper-bounded by the solution of the linear system, that is,

$$0_n \le x(t) \le e^{A_r t} x_0.$$

On the other hand, there exists a block diagonal matrix $J$, where each block is a Jordan matrix, and an invertible matrix $P$ such that $A_r = PJP^{-1}$. Therefore,

$$\|x(t)\| \le \|e^{A_r t} x_0\| = \|Pe^{Jt}P^{-1}x_0\| \le C\|e^{Jt}\|. \qquad (A.24)$$
If $A_r$ is diagonalizable, then $J$ is the diagonal matrix with the eigenvalues $\lambda_1, ..., \lambda_n$ of $A_r$ on its diagonal. Therefore, $e^{Jt} = \text{diag}(e^{\lambda_1 t}, ..., e^{\lambda_n t})$, so that $\|e^{Jt}\| \le e^{-\alpha t}$. It follows

$$ \|x(t)\| \le Ce^{-\alpha t}. $$

Therefore, for any $\alpha' < \alpha$, one has $e^{\alpha' t} \|x(t)\| \le Ce^{(\alpha' - \alpha)t} \to 0$ as $t \to +\infty$.

If $A_r$ is not diagonalizable, for every eigenvalue $\lambda$ of $A$, let $J_k = \lambda I + N$ be a Jordan block of size $k$, where $N$ is a matrix of all zeros except identities on the first superdiagonal. Since $J$ is a diagonal block matrix, it is enough to bound each block. For a Jordan block $J_k$, $N^n = 0$ for any $n \ge k$, which implies

$$ e^{J_k t} = e^{(\lambda I + N)t} = e^{\lambda t} e^{Nt} = e^{\lambda t} \sum_{n=0}^{k} \frac{N^n}{n!} t^n = e^{(\lambda + \varepsilon)t} \sum_{n=0}^{k} \frac{N^n}{n!} t^n e^{-\varepsilon t}. $$

Therefore, for any $\varepsilon > 0$, $\|e^{J_k t}\| \le C_\varepsilon e^{(-\alpha + \varepsilon)t}$. One deduces from (A.24) that for any $\varepsilon > 0$, there exists $C_\varepsilon > 0$ such that $\|x(t)\| \le C_\varepsilon e^{(-\alpha + \varepsilon)t}$. Moreover, for any $\alpha' < \alpha$, there exists $\varepsilon$ satisfying $0 < \varepsilon < \alpha - \alpha'$. It follows $e^{\alpha' t} \|x(t)\| \le C_\varepsilon e^{(\alpha' - \alpha + \varepsilon)t} \to 0$ as $t \to +\infty$. $\square$

When assumption (H2) holds, Theorem 2.2 derives from Theorem A.1.

# B Control effects properties (Thm. 3.1, 3.2)

This section is dedicated to the proofs of Theorem 3.1 and Theorem 3.2.

## B.1 Proof of Theorem 3.1

Let $C \subset \{1, ..., n\}$. The function: $\rho \in \mathbb{R}_+^n \mapsto s(A(\rho))$ being twice differentiable is a consequence of [13, Lemma 3.1 and Theorem 3.2] extended to irreducible Metzler matrices. Moreover, for any $\rho < \rho'$, one has $A(\rho') < A(\rho)$. This implies that $s(A(\rho')) < s(A(\rho))$ (see [40, Lemma 1.2], also extended to irreducible Metzler matrices).

The convexity of the function $\rho \in \mathbb{R}_+^n \mapsto s(A(\rho))$ is a consequence of Theorem 4.1 in [21], which addresses the (strict) convexity properties of (irreducible) non-negative matrices. Once again, this property can be easily extended to (irreducible) Metzler matrices. When $C \subset \{1, ..., n\}$, for any $\bar{\rho}$-admissible $\rho$ and $\rho' \in \mathbb{R}_+^n$, it is obvious that $\rho - \rho'$ cannot be a scalar multiple of the unit vector, and it also follows from [21, Theorem 4.1] that the function: $\rho \mapsto s(A(\rho))$ is strictly convex.

## B.2 Proof of Theorem 3.2

To prove Theorem 3.2, let us first prove the following result.

**Proposition B.1.** Let $A \in \mathcal{M}_n(\mathbb{R})$ be an irreducible Metzler matrix. Then,

$$ s(A) < 0 \implies A \text{ non singular and } A^{-1} \ll 0_{n \times n}. $$
Proof. We will show the following, statement: for any irreducible matrix A such that $-A$ is Metzler, $s(A) > 0$ implies A non singular and $A^{-1} \gg 0_{n \times n}$.

As a matter of fact, if $-A$ is a Metzler matrix, then $s(A) > 0$ iff A non singular and $A^{-1} > 0_{n \times n}$ (see [28, Theorem 2.5.3]). We make this statement more precise when A is also irreducible.

One has $A = \alpha I - P$ for some irreducible $P > 0_{n \times n}$ and $\alpha > 0$. Without loss of generality, let us suppose that $\alpha = 1$. Then, $A = I - P$, and since $s(A) > 0$, $s(P) < 1$. Thanks to the Perron-Frobenius theorem, $s(P) = \rho(P)$, where $\rho(P)$ denotes the spectral radius of $P$. It follows $\rho(P) < 1$. The Neumann series $I+P+P^2+...$ then converges towards $(I-P)^{-1} = A^{-1}$. Moreover, since $P > 0_{n \times n}$ and is irreducible, one has $I+P+P^2+... \gg 0_{n \times n}$, as a consequence of [6, Theorem 2.1]. Therefore, $A^{-1} \gg 0_{n \times n}$. $\square$

We can now prove Theorem 3.2. For any $i, j = 1, ..., n$, note $E_{ij} \in \mathcal{M}_n(\mathbb{R})$ the matrix whose $(i, j)$-th entry is 1 and whose remaining entries are 0. Since $x^*(\rho)$ is an equilibrium point of system (1), for any $\rho \bar{\rho}$-admissible, one has

$$ \left( D + \mathrm{diag}(g_i(x_i^*(\rho))) - \mathrm{diag}(\rho) \right) x^*(\rho) = 0_n. \qquad (B.25) $$

Let us define

$$ B(\rho) := D + \mathrm{diag}(g_i(x_i^*(\rho))) - \mathrm{diag}(\rho), $$

which is an irreducible Metzler matrix. Since we assume $x^*(\rho) \gg 0_n$, it follows from Perron-Frobenius theorem and (B.25) that $x^*(\rho)$ is a Perron vector of $B(\rho)$. In particular, $s(B(\rho)) = 0$.

On the other hand $\{\rho \in \mathbb{R}^n : x^*(\rho) \gg 0_n\} = \{\rho \in \mathbb{R}^n : s(A(\rho)) > 0\}$ is an open set, by continuity of the function $s \circ A$. By differentiating (B.25) for any $\rho \ge 0_n$ and rearranging sides, one has for any $j \in \mathcal{C}$,

$$ \left[ \mathrm{diag} \left( \frac{\partial x_i^*(\rho)}{\partial \rho_j} g_i'(x_i^*(\rho)) \right) - E_{jj} \right] x^*(\rho) = -B(\rho) \frac{\partial x^*(\rho)}{\partial \rho_j}. $$

Since

$$ \mathrm{diag} \left( \frac{\partial x_i^*(\rho)}{\partial \rho_j} g_i'(x_i^*(\rho)) \right) x^*(\rho) = \mathrm{diag} \left( x_i^*(\rho) g_i'(x_i^*(\rho)) \right) \frac{\partial x^*(\rho)}{\partial \rho_j}, $$

it follows

$$ \left[ B(\rho) + \mathrm{diag} \left( x_i^*(\rho) g_i'(x_i^*(\rho)) \right) \right] \frac{\partial x^*(\rho)}{\partial \rho_j} = E_{jj} x^*(\rho). \qquad (B.26) $$

One has that $B(\rho)+\mathrm{diag}(x_i^*(\rho)g_i'(x_i^*(\rho)))$ is an irreducible Metzler matrix. Since $g_i'(x_i^*(\rho)) < 0$ a.e. du to (H2) and $x_i^*(\rho) > 0$,

$$ B(\rho) + \mathrm{diag} \left( x_i^*(\rho) g_i'(x_i^*(\rho)) \right) < B(\rho) \quad \text{a.e.} $$
It follows by [40, Lemma 1.2] extended to irreducible Metzler matrices that

$$s\left(B(\rho) + \text{diag}\left(x_i^*(\rho)g_i'(x_i^*(\rho))\right)\right) < s(B(\rho)) = 0 \quad \text{a.e.}$$

Ones deduces from Proposition B.1 that $B(\rho) + \text{diag}(x_i^*(\rho)g_i'(x_i^*(\rho)))$ is non singular and

$$\left(B(\rho) + \text{diag}\left(x_i^*(\rho)g_i'(x_i^*(\rho))\right)\right)^{-1} \ll 0_{n \times n} \quad \text{a.e.} \qquad (B.27)$$

By multiplying both sides of (B.26) by $\left(B(\rho) + \text{diag}\left(x_i^*(\rho)g_i'(x_i^*(\rho))\right)\right)^{-1}$, it follows

$$\frac{\partial x^*(\rho)}{\partial \rho_j} = \left(B(\rho) + \text{diag}\left(x_i^*(\rho)g_i'(x_i^*(\rho))\right)\right)^{-1} E_{jj}x^*(\rho). \qquad (B.28)$$

Since $x^*(\rho) \gg 0_n$ and thanks to (B.27), it follows from (B.28)

$$\frac{\partial x^*(\rho)}{\partial \rho_j} \ll 0_n \quad \text{a.e.}$$

This achieves the proof of Theorem 3.2.

# C Limit behaviour for large control (Lemma 3.3, Thm. 3.4, Coro. 3.5)

## C.1 Proof of Lemma 3.3

As $s(A(\rho))$ is strictly decreasing with respect to $\rho$ (see Theorem 3.1), it follows

$$\inf\{s(A(\rho)); \rho \bar{\rho}\text{-admissible}\} = \lim_{\rho \to \bar{\rho}} s(A(\rho)).$$

In particular, the infimum is a minimum if and only if $\bar{\rho}$ does not take infinite values, i.e. if and only if $C_I = \emptyset$.

Denoting the vectors $\bar{\rho}_{C_F}$ and $\mathbf{1}_{C_I}$ as follows:

$$(\bar{\rho}_{C_F})_i := \begin{cases} \bar{\rho}_i & \text{if } i \in C_F \\ 0 & \text{if } i \in \bar{C}_F \end{cases}, \qquad (\mathbf{1}_{C_I})_i := \begin{cases} 1 & \text{if } i \in C_I \\ 0 & \text{if } i \in \bar{C}_I \end{cases}, \qquad (C.29)$$

one has

$$\lim_{\rho \to \bar{\rho}} s(A(\rho)) = \lim_{\lambda \to +\infty} s(A(\lambda \mathbf{1}_{C_I} + \bar{\rho}_{C_F})).$$
### C.1.1 Case 1: $C_I = C = \{1, ..., n\}$

Let us first suppose that $C_I = C = \{1, ..., n\}$. In particular, $\bar{\rho}_{C_F} = 0_n$, which implies $\lim_{\rho \to \bar{\rho}} s(A(\rho)) = \lim_{\lambda \to +\infty} s(A(\lambda \mathbf{1}_n))$. On the other hand, $A(\lambda \mathbf{1}_n) = D + \text{diag}(g_i(0)) - \lambda \text{diag}(\mathbf{1}_n)$. Denoting by $\sigma(A(\lambda \mathbf{1}_n))$ the spectrum of $A(\lambda \mathbf{1}_n)$, one has thanks to the Gershgorin circle theorem [22, 49]

$$ \sigma(A(\lambda \mathbf{1}_n)) \subset \bigcup_{i=1}^{n} \mathcal{B} \left( (A(\lambda \mathbf{1}_n))_{ii}, \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} \right), \qquad (C.30) $$

where $\mathcal{B}$ denotes closed ball. For any $i = 1, ..., n$,

$$ (A(\lambda \mathbf{1}_n))_{ii} + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} = g_i(0) - \lambda \xrightarrow{\lambda \to +\infty} -\infty. $$

It follows from (C.30) that

$$ \lim_{\lambda \to +\infty} s(A(\lambda \mathbf{1}_n)) = -\infty = s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho})). $$

The last equality is a consequence of the convention adopted in Definition 1.5 for $\overline{C_I} = \emptyset$.

### C.1.2 Case 2: $C_I \subsetneq \{1, ..., n\}$

Without loss of generality, let us suppose that the first $n_{C_I}$ patches are infinitely controllable. In any case, it can be achieved by a permutation of components. With such a convention on the ordering of the components, one has

$$ A(\lambda \mathbf{1}_{C_I} + \bar{\rho}_{C_F}) = \begin{pmatrix} A|_{C_I, C_I} - \lambda I_{n_{C_I}} & A|_{C_I, \bar{C}_I} \\ A|_{\bar{C}_I, C_I} & A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho}) \end{pmatrix}, $$

where $I_{n_{C_I}}$ denotes the $n_{C_I} \times n_{C_I}$ identity matrix. Moreover, one has

$$ \lim_{\lambda \to +\infty} s \begin{pmatrix} A|_{C_I, C_I} - \lambda I_{n_{C_I}} & A|_{C_I, \bar{C}_I} \\ A|_{\bar{C}_I, C_I} & A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho}) \end{pmatrix} = s(A|_{\bar{C}_I, \bar{C}_I}(\bar{\rho})). \qquad (C.31) $$

Equality (C.31) is a direct application of the following lemma.

**Lemma C.1.** Let $n \ge 2$ and $m \in \mathbb{N}$ such that $n - m \ge 1$. Let $A_{11} \in \mathcal{M}_m(\mathbb{R})$, $A_{12} \in \mathcal{M}_{m \times (n-m)}(\mathbb{R})$, $A_{21} \in \mathcal{M}_{(n-m) \times m}(\mathbb{R})$ and $A_{22} \in \mathcal{M}_{n-m}(\mathbb{R})$. Then,

$$ \lim_{\varepsilon \to 0^+} s \begin{pmatrix} A_{11} - \frac{1}{\varepsilon} I_m & A_{12} \\ A_{21} & A_{22} \end{pmatrix} = s(A_{22}) $$
*Proof.* Drawing inspiration from singular perturbation theory ([34, pp. 49-55]), we find a change of basis transforming $A(\varepsilon) := \begin{pmatrix} A_{11} - \frac{1}{\varepsilon}I_m & A_{12} \\ A_{21} & A_{22} \end{pmatrix}$ into an upper triangular block matrix $\tilde{A}(\varepsilon)$, whose eigenvalues are thus the combined eigenvalues of the diagonal blocks.

For some $(n-m) \times m$ matrix $L(\varepsilon)$ to be chosen afterwards, let us denote

$$P := \begin{pmatrix} I_m & 0 \\ \varepsilon L(\varepsilon) & I_{n-m} \end{pmatrix}.$$

Straightforward computations show that

$$\tilde{A}(\varepsilon) := PA(\varepsilon)P^{-1} = \begin{pmatrix} A_{11} - \frac{1}{\varepsilon}I_m - \varepsilon A_{12}L(\varepsilon) & A_{12} \\ R(L(\varepsilon), \varepsilon) - L(\varepsilon) & \varepsilon L(\varepsilon)A_{12} + A_{22} \end{pmatrix},$$

where $R(L, \varepsilon) := A_{21} - \varepsilon A_{22}L + \varepsilon LA_{11} - \varepsilon^2 LA_{12}L$. We will now show that, for sufficiently small $\varepsilon^* > 0$, there exists, for any $\varepsilon$ in $[0, \varepsilon^*]$, a fixed-point $L(\varepsilon) \in \mathcal{M}_{(n-m) \times m}(\mathbb{R})$ of $R(L, \varepsilon)$, that is a solution to the equation $R(L, \varepsilon) - L = 0$. For $\varepsilon = 0$, the unique solution of this equation is $L(0) = A_{21}$.

Let $c > \|A_{21}\|$. For any $L$ in the compact ball $\mathcal{B}(0_{(n-m) \times m}, c)$, one has

$$\|R(L, \varepsilon)\| \le \|A_{21}\| + (\|A_{11}\| + \|A_{22}\|)\varepsilon c + \|A_{12}\|\varepsilon^2 c^2.$$

Take $\varepsilon^* > 0$ such that $\|A_{21}\| + (\|A_{11}\| + \|A_{22}\|)\varepsilon^* c + \|A_{12}\|\varepsilon^{*2} c^2 < c$. Then, for any $\varepsilon \in [0, \varepsilon^*]$, the continuous map $L \mapsto R(L, \varepsilon)$ maps $\mathcal{B}(0_{(n-m) \times m}, c)$ into itself. Thus, due to Brouwer's fixed-point theorem, for any $\varepsilon \in [0, \varepsilon^*]$, $R(L, \varepsilon)$ possesses a fixed-point in $\mathcal{B}(0_{(n-m) \times m}, c)$.

Therefore, for any $\varepsilon \in (0, \varepsilon^*]$,

$$\begin{aligned} s(A(\varepsilon)) &= s(P^{-1}\tilde{A}(\varepsilon)P) \\ &= s(\tilde{A}(\varepsilon)) \\ &= \max \left\{ s(A_{11} - \frac{1}{\varepsilon}I_m - \varepsilon L(\varepsilon)A_{12}), s(A_{22} + \varepsilon L(\varepsilon)A_{12}) \right\} \\ &= \max \left\{ s(A_{11} - \varepsilon L(\varepsilon)A_{12}) - \frac{1}{\varepsilon}, s(A_{22} + \varepsilon L(\varepsilon)A_{12}) \right\}. \end{aligned}$$

As the matrix $L(\varepsilon)$ pertain to the bounded set $\mathcal{B}(0_{(n-m) \times m}, c)$, one deduces in conclusion

$$\lim_{\varepsilon \to 0^+} s(A(\varepsilon)) = \max \left\{ \lim_{\varepsilon \to 0^+} \left( s(A_{11}) - \frac{1}{\varepsilon} \right), s(A_{22}) \right\} = s(A_{22}).$$

This achieves the proof of Lemma C.1, and consequently of Lemma 3.3. $\square$

## C.2 Proof of Theorem 3.4

Assume first that $C_I$ is not empty. If $s(A|_{C_I \times C_I}(\bar{\rho})) < 0$, then thanks to (10), there exists, by continuity, $\rho' \bar{\rho}$-admissible such that $s(A(\rho')) \le 0$. Let $\rho \bar{\rho}$-admissible satisfying $\rho \ge \rho'$.
By Theorem 3.1, there exists $\alpha := -s(A(\rho)) \ge -s(A(\rho'))$. One deduces from Theorem 2.2 that the origin is $\alpha$-GAS on $\mathbb{R}_+^n$ for the system (7), which proves the point (i) of Theorem 3.4 when $C_I \ne \emptyset$. If now $s(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho})) \ge 0$, then $\inf\{s(A(\rho)); \rho \bar{\rho}$-admissible\} \ge 0$. As $C_I$ is not empty, the infimum is not a minimum, due to Lemma 3.3, and the point (ii) of Theorem 3.4 when $C_I \ne \emptyset$ also follows from Theorem 2.2.

If now $C_I$ is empty, then the infimum in (10) is a minimum, and $A|_{\overline{C_I},\overline{C_I}}(\bar{\rho}) = A(\bar{\rho})$. Therefore, the origin of (7) is GAS for $\bar{\rho}$-admissible $\rho$ if and only if $s(A(\bar{\rho})) \ge 0$.

### C.3 Proof of Corollary 3.5

Let us first recall that

$$A_{\overline{C_I},\overline{C_I}}(\bar{\rho}) = (D_{ij})_{i \in \overline{C}, j \in \overline{C}} + \text{diag}(g_i(0) - \rho_i, i \in \overline{C}).$$

(i) Let us suppose that $g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} < 0$ for every $i \in \overline{C}$. We show that $s(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho})) < 0$. The Gershgorin circle theorem [22, 49] states that

$$\sigma(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho})) \subset \bigcup_{i \in \overline{C}} \mathcal{B}\left( (A|_{\overline{C_I},\overline{C_I}}(\bar{\rho}))_{i,i}, \sum_{\substack{j \in \overline{C} \\ j \ne i}} D_{ji} \right).$$

Furthermore, for any $i \in \overline{C_I}$,

$$(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho}))_{i,i} = g_i(0) - \bar{\rho}_i - \sum_{\substack{j=1 \\ j \ne i}}^{n} D_{ji},$$

so that

$$(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho}))_{i,i} + \sum_{\substack{j \in \overline{C} \\ j \ne i}} D_{ji} = g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} < 0.$$

Therefore, the Gershgorin circles of $A|_{\overline{C_I},\overline{C_I}}$ stay strictly at the left of the complex axis. This implies that $s(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho})) < 0$. The conclusion then follows from Theorem 3.4.

(ii) Let us suppose that $g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} > 0$ for every $i \in \overline{C}$. We will show that $s(A|_{\overline{C_I},\overline{C_I}}(\bar{\rho})) > 0$. Let us suppose without loss of generality that the patches in $\overline{C_I}$ are the first $n_{\overline{C_I}}$ patches. One has,

$$\begin{align*}
\mathbf{1}_{n_{\overline{C_I}}}^{\mathrm{T}} A|_{\overline{C_I},\overline{C_I}}(\bar{\rho}) &= (g_1(0) - \bar{\rho}_1 - \sum_{j \in C_I} D_{j1}, \dots, g_{n_{\overline{C}}}(0) - \bar{\rho}_{n_{\overline{C}}}) - \sum_{j \in C_I} D_{jn_{\overline{C}}} \\
&\ge \min_{i \in \overline{C_I}} \{g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji}\} \mathbf{1}_{n_{\overline{C_I}}}^{\mathrm{T}}.
\end{align*}$$
Since $\min_{i \in \overline{C}_I} \{g_i(0) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji}\} > 0$ and since $\mathbf{1}_{n_{\overline{C}_I}} > 0_{n_{\overline{C}_I}}$, we can apply [39, Lemma 8],
which can be easily extended to Metzler matrices. Therefore,

$$s(A_{\overline{C}_I, \overline{C}_I}^{\top}(\bar{\rho})) \geq \min_{i \in \overline{C}_I} \{g_i(0) - \bar{\rho}_i - \sum_{j \in C_i} D_{ji}\} > 0.$$

As $s(A|_{\overline{C}_I, \overline{C}_I}(\bar{\rho})) = s(A_{\overline{C}_I, \overline{C}_I}^{\top}(\bar{\rho}))$, it follows $s(A_{\overline{C}_I, \overline{C}_I}) > 0$. The conclusion then follows from Theorem 3.4.

## D Minimal globally asymptotically stable equilibrium (Thm. 4.1 and Thm. 4.2)

First note that if $s(A|_{\overline{C}_I, \overline{C}_I}(\bar{\rho})) < 0$, or $s(A|_{\overline{C}_I, \overline{C}_I}(\bar{\rho})) \le 0$ and $C_I = \emptyset$, it directly follows from Theorem 3.4 that

$$x^{\infty} = \min\{x^{*}(\rho) : \rho \ \bar{\rho}\text{-admissible}\} = 0_{n}.$$

More generally, as $x^*(\rho)$ is a non-decreasing function of $\rho$ and using similar arguments as for the proof of Lemma 3.3 in Appendix C.1, one has:

$$x^{\infty} = \lim_{\rho \to \bar{\rho}} x^{*}(\rho) = \lim_{\lambda \to +\infty} x^{*}(\lambda \mathbf{1}_{C_I} + \bar{\rho}_{C_F}),$$

for $\bar{\rho}_{C_F}$ and $\mathbf{1}_{C_I}$ respectively defined in (C.29).

For any $\lambda \ge 0$, $x^*(\lambda \mathbf{1}_C + \bar{\rho}_{C_F})$ belongs to the set

$$E := \{x \in \mathbb{R}_+^n : x \le x^*(0_n)\}. \qquad (D.32)$$

Note that this set is independent of $\lambda$ and is compact.

On the other hand, for any $\lambda > 0$, $x^*(\lambda \mathbf{1}_{C_I} + \bar{\rho}_{C_F})$ is the largest non-negative solution of the following system

$$\left\{ \begin{array}{ll} -x_i + \frac{1}{\lambda} \left( x_i g_i(x_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij} x_j - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} x_i \right) = 0 & i \in C_I \\ x_i (g_i(x_i) - \bar{\rho}_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij} x_j - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} x_i = 0 & i \in \overline{C}_I. \end{array} \right.$$

Due to the fact that $x^*(\lambda \mathbf{1}_C + \bar{\rho}_{C_F})$ belongs to the bounded set $E$, the limit equations when $\lambda \to +\infty$ are given by

$$\left\{ \begin{array}{ll} x_i = 0 & i \in C_I \\ x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} \right) + \sum_{\substack{j \in C_I \\ j \neq i}} D_{ij} x_j - \sum_{\substack{j \in C_I \\ j \neq i}} D_{ji} x_i = 0 & i \in \overline{C}_I. \end{array} \right. \qquad (D.33)$$
Notice that any solution $x$ of this system satisfies

$$x|_{c_I} = 0_{nc_I} \qquad (D.34)$$

and $x|_{\bar{c}_I}$ is an equilibrium point of the residual system (11).

## D.1 Irreducible residual system

We now assume that $D|_{\bar{c}_I, \bar{c}_I}$ is irreducible, so that we can apply Theorem 2.2 to the system (11). Since $A|_{\bar{c}_I, \bar{c}_I}(\bar{\rho})$ is the Jacobian in zero of the system (11), it follows that

* If $s(A|_{\bar{c}_I, \bar{c}_I}(\bar{\rho})) \le 0$, then the origin of the system (11) is GAS on $\mathbb{R}_+^{n_{\bar{c}}}$.

* If $s(A|_{\bar{c}_I, \bar{c}_I}(\bar{\rho})) > 0$, then the system (11) admits a positive equilibrium point, which is globally asymptotically stable on $\mathbb{R}_+^{n_{\bar{c}}} \setminus \{0_{n_{\bar{c}}}\}$.

In any case, the globally asymptotically stable equilibrium point belongs to the set $E$ defined in (D.32) and is the largest solution of the system (D.33). Moreover, since $E$ is compact, the equilibrium $x^*(\lambda \mathbf{1}_{c_I} + \bar{\rho}_{c_F})$ has at least a limit point in $E$ when $\lambda$ goes to infinity. Therefore,

* If $s(A|_{\bar{c}_I, \bar{c}_I}(\bar{\rho})) \le 0$, then

$$x^\infty|_{\bar{c}_I} = \lim_{\lambda \to +\infty} x^*(\lambda \mathbf{1}_{c_I} + \bar{\rho}_{c_F})|_{\bar{c}_I} = 0_{n_{\bar{c}_I}}.$$

* If $s(A|_{\bar{c}_I, \bar{c}_I}(\bar{\rho})) > 0$, then

$$x^\infty|_{\bar{c}_I} = \lim_{\lambda \to +\infty} x^*(\lambda \mathbf{1}_{c_I} + \bar{\rho}_{c_F})|_{\bar{c}_I}$$

is equal to the positive equilibrium of the residual system (11).

Moreover, in both cases, it follows from (D.34)

$$x^\infty|_{c_I} = \lim_{\lambda \to +\infty} x^*(\lambda \mathbf{1}_{c_I} + \bar{\rho}_{c_F})|_{c_I} = 0_{nc_I}.$$

This ends the proof of Theorem 4.1.

## D.2 Reducible residual system

To prove Theorem 4.2, we have to study the equilibrium points of the residual system (11) on $\bar{c}_I$, which may not be irreducible. For reasons that will be explained later, in Section D.2.2, this requires to analyze an auxiliary model to the model (1), with a constant inflow. This model is studied in Section D.2.1, where the existence of a positive GAS equilibrium is established (Theorem D.1). To prove this result, we first demonstrate the existence of a positive GAS equilibrium for the one-patch model (Lemma D.2 in Section D.2.1.1). We then derive a persistence result for the multi-patch model (Lemma D.3 in Section D.2.1.2), which allows to prove Theorem D.1 in Section D.2.1.3. The proof of Theorem 4.2, which relies on Theorem D.1 is then given in Section D.2.2.
### D.2.1 Auxiliary model with constant inflow

Let us study the following auxiliary model, which corresponds to the model (1) with the addition of a non-negative constant vector $C = (C_1, ..., C_n)$.

$$ \dot{x}_i = x_i g_i(x_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij} x_j - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} x_i + C_i, \quad x_i(0) \ge 0, \quad i = 1, ..., n. \quad (D.35) $$

Let us define $F^C(x)$ as the right hand side of the model (D.35), that is,

$$ F^C(x) := (F_1^{C_1}(x), ..., F_n^{C_n}(x)), \quad F_i^{C_i}(x) := x_i g_i(x_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij} x_j - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} x_i + C_i. $$

Notice that $F^C(x) = F(x) + C$, where $F$, defined in (2), is the vector form of the right-hand side of the system (1). When $C = 0_n$, the system (D.35) is identical to the system (1). When $C > 0_n$, the aim of this subsection is to prove the following theorem.

**Theorem D.1.** Assume (H2), (H3) and (H4) hold. If $C > 0_n$, then the system (D.35) admits a unique equilibrium point $x^*$, which is positive and GAS on $\mathbb{R}_+^n$.

To prove this theorem, one first studies the one-patch model.

#### D.2.1.1 One-patch model

The one-patch model writes

$$ \dot{x}_1 = x_1 g_1(x_1) + C_1, \quad x_1(0) \ge 0. \quad (D.36) $$

When $C_1 > 0$, the following lemma gives the asymptotic behaviour of the solutions of the model (D.36).

**Lemma D.2.** Assume $C_1 > 0$. Then any solution $x_1$ of (D.36) initiated in $x_1(0) \ge 0$ satisfies

$$ \liminf_{t \to +\infty} x_1(t) > 0. \quad (D.37) $$

If in addition (H3) holds, then the system (D.36) admits a positive equilibrium point. Moreover, if (H2) also holds, then this equilibrium point is GAS on $\mathbb{R}_+$.

*Proof.* Since $C_1 > 0$, one has $F_1(0) = C_1 > 0$. By continuity of $F_1$, the exists $\bar{a} > 0$ such that

$$ F(a) > 0 \quad \text{for any } a \in [0, \bar{a}). $$

Moreover, for any $a \in [0, \bar{a})$, the set $\{x_1 \ge a\}$ is positively invariant.

Let $x_1(0) > 0$. There exists $a \in (0, \bar{a})$ such that $x_1(0) \ge a$. Since $\{x_1 \ge a\}$ is positively invariant, it follows $x_1(t) \ge a$ for any $t \ge 0$. This implies (D.37).
If $x_1(0) = 0$, then since $F(x_1(0)) = F(0) > 0$, there exists $\eta$ such that for any $t \in [0, \eta)$, $x_1(t) > 0$. As above, one shows that (D.37) is veritifed.

Let us now assume that (H3) holds. The function $x_1 \mapsto x_1g_1(x_1)$ is null when $x_1 = 0$, and converges towards $-\infty$ as $x_1 \to +\infty$. Therefore, by continuity, there exists $x_1^* > 0$ satisfying $x_1^*g_1(x_1^*) = -C_1 < 0$. This proves the existence of a positive equilibrium point for the system (D.36). Moreover, since $x_1^* > 0$, we have in particular that $g_1(x_1^*) < 0$. Let us assume that (H2) holds and let $x_1 > x_1^*$. Then,

$$ 
\begin{aligned}
F_1^{C_1}(x_1) &= x_1g_1(x_1) + C_1 \\
&< x_1g_1(x_1^*) + C_1 \quad \text{since } g_1 \text{ is decreasing} \\
&< x_1^*g_1(x_1^*) + C_1 \quad \text{since } x_1 > x_1^* > 0 \text{ and } g_1(x_1^*) < 0 \\
&= 0.
\end{aligned}
$$

Therefore,

$$ F_1^{C_1}(x_1) < 0, \quad \text{for any } x_1 > x_1^*. \qquad (D.38) $$

Following the same methodology for $x_1 < x_1^*$, one finds

$$ F_1^{C_1}(x_1) > 0, \quad \text{for any } x_1 < x_1^*. \qquad (D.39) $$

One easily deduces from (D.38) and (D.39) that $x_1^*$ is globally asymptotically stable on $\mathbb{R}_+$.

### D.2.1.2 Multi-patch model

Drawing inspiration from [36, 47], one now studies the general case where $n \ge 2$.

**Lemma D.3.** If (H3) holds, then all the solutions of (D.35) are uniformly bounded.

If (H4) holds and $C > 0_n$, then any solution $x$ of (D.35) initiated in a point $x(0) \ge 0_n$ satisfies

$$ \liminf_{t \to +\infty} x(t) \gg 0_n. \qquad (D.40) $$

If (H3) and (H4) both hold, and if $\lim_{x_i \to +\infty} g_i(x_i) < g_i(0)$ for every $i = 1, \dots, n$, then the system (D.35) admits a positive equilibrium point.

*Proof.*

#### (i) Boundedness

Assume (H3) holds and let us compute the sum of the $\dot{x}_i$.

$$ \sum_{i=1}^{n} \dot{x}_i = \sum_{i=1}^{n} x_i g_i(x_i) + \sum_{i=1}^{n} C_i. $$

Thanks to (H3), $\lim_{x_i \to +\infty} g_i(x_i) < 0$, so that $x_i g_i(x_i) + C_i < 0$ for $x_i$ large enough. Therefore, $\sum_{i=1}^{n} \dot{x}_i < 0$ for sufficiently large $x_i$, $i = 1, \dots, n$.
(ii) Persistence

Assume (H4) holds and let $x(0) \ge 0_n$.

**Step 1:** Let us prove that for every $i = 1, ..., n$ such that $C_i > 0$, $\liminf_{t \to +\infty} x_i(t) > 0$ for any $x_i(0) \ge 0$. Let $i \in \{1, ..., n\}$ such that $C_i > 0$, then

$$ \dot{x}_i = x_ig_i(x_i) + \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ij}x_j - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji}x_i + C_i \geq x_i \left( g_i(x_i) - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} \right) + C_i. $$

Let us define the system

$$ \dot{x}_i = x_i \left( g_i(x_i) - \sum_{\substack{j=1 \\ j \neq i}}^{n} D_{ji} \right) + C_i. \quad (D.41) $$

Let $x_i$ and $\tilde{x}_i$ be respectively the solutions of system (D.35) and (D.41) initiated in $x_i(0) \ge 0$. By comparison,

$$ x_i(t) \ge \tilde{x}_i(t) \quad \text{for any } t \ge 0, $$

Therefore,

$$ \liminf_{t \to +\infty} x_i(t) \ge \liminf_{t \to +\infty} \tilde{x}_i(t) > 0, $$

the last inequality being a consequence of Lemma D.2.

**Step 2:** Let us prove that if there exists $i$ and $j$ in $\{1, ..., n\}$ such that $i \neq j$, $\liminf_{t \to +\infty} x_j(t) > 0$ and $D_{ij} > 0$, then $\liminf_{t \to +\infty} x_i(t) > 0$.

Since $\liminf_{t \to +\infty} x_j(t) > 0$, there exists $\varepsilon > 0$ and $T \ge 0$ such that

$$ x_j(t) \ge \varepsilon \quad \text{for any } t \ge T. $$

Then, for any $t \ge T$,

$$ \begin{aligned} \dot{x}_i &= x_ig_i(x_i) + \sum_{\substack{k=1 \\ k \neq i}}^{n} D_{ik}x_k - \sum_{\substack{k=1 \\ k \neq i}}^{n} D_{ki}x_i + C_i \\ &\ge x_ig_i(x_i) + D_{ij}x_j - \sum_{\substack{k=1 \\ k \neq i}}^{n} D_{ki}x_i + C_i \\ &\ge x_ig_i(x_i) + D_{ij}\varepsilon - \sum_{\substack{k=1 \\ k \neq i}}^{n} D_{ki}x_i + C_i \\ &= x_i \left( g_i(x_i) - \sum_{\substack{k=1 \\ k \neq i}}^{n} D_{ki} \right) + D_{ij}\varepsilon + C_i, \end{aligned} $$
with $D_{ij}\epsilon + C_i > 0$. The proof of the persistence in the i-th patch is then the same as in Step 1.

**Step 3:** Let us prove that $\liminf_{t \to +\infty} x(t) \gg 0_n$. Since $C > 0_n$, there exists $i \in \{1, ..., n\}$ such that $C_i > 0$. Moreover, for any $k \in \{1, ..., n\}$ such that $k \neq i$, there exists by irreducibility of $D$ a path between patches $k$ and $i$, that is, there exists $m \in \{2, ..., n\}$ and $a_2, ..., a_{m-1}$ all different such that

$$D_{ka_{m-1}}...D_{a_3a_2}D_{a_2i} > 0.$$

It follows by induction and from step 1 and step 2 that (D.40) holds for any $x(0) \ge 0_n$.

### (iii) Positive equilibrium point

Let us now assume that (H3) and (H4) both hold, and that $\lim_{x_i \to +\infty} g_i(x_i) < g_i(0)$ for every $i = 1, ..., n$. Since $A$ is irreducible thanks to (H4), one may apply the Perron-Frobenius theorem [6]. It states that there exists a vector $v \gg 0_n$ such that

$$Av = s(A)v. \qquad (D.42)$$

On the other hand, with some computations, we find that for any scalar $a > 0$,

$$F_i^{C_i}(av) = av_i(g_i(av_i) - g_i(0)) + a(Av)_i + C_i. \qquad (D.43)$$

Let us first assume $s(A) > 0$. Thanks to (D.42) and since $v \gg 0_n$, we have $Av \gg 0_n$. For every $i = 1, ..., n$, one has thanks to (D.43) that

$$\left. \frac{\partial F_i^{C_i}(av)}{\partial a} \right|_{a=0} = (Av)_i > 0.$$

and

$$F_i^{C_i}(av)|_{a=0} = C_i \ge 0.$$

Therefore, there exists a scalar $\underline{a} > 0$ such that

$$F^C(av) \gg 0_n \text{ for any } 0 < a < \underline{a}.$$

Let $x$ be the solution of the system (D.35) initiated in $av$ for $0 < a < \underline{a}$. Since $F^C(av) \gg 0_n$ and since the system (D.35) is cooperative, $x$ is increasing. As it is bounded thanks to Lemma D.3, it converges towards a positive equilibrium. Therefore, if $s(A) > 0$, there exists a positive equilibrium point.

Let us now assume $s(A) \le 0$. Thanks to (D.42) and since $v \gg 0_n$, one has $Av \le 0_n$. Since $\lim_{x_i \to +\infty} g_i(x_i) < g_i(0)$, it follows as a consequence of (D.43) that, for every $i = 1, ..., n$,

$$F_i^{C_i}(av) < 0 \text{ for a large enough.} \qquad (D.44)$$

Therefore, there exists a scalar $\bar{a} > 0$ such that

$$F^C(av) \ll 0_n \text{ for any } a > \bar{a}.$$

Let $x$ be the solution of the system (D.35) initiated in $av$ for $a > \bar{a}$. Since $F^C(av) \ll 0_n$, $x$ is decreasing along this trajectory. By boundedness and since (D.40) holds, it converges towards a positive equilibrium point. This proves the existence of a positive equilibrium point if $s(A) \ge 0$ and achieves the proof of Lemma D.3. $\square$
### D.2.1.3 Proof of Theorem D.1

If (H2) holds, it implies that $\lim_{x_i \to +\infty} g_i(x_i) < g_i(0)$ for every $i = 1, \dots, n$. If (H3) and (H4) also hold, it follows from Lemma D.3 that the system (D.35) admits a positive equilibrium point.

Let $x^*$ be a positive equilibrium$^4$. With some computations, one finds the identity

$$F_i^{C_i}(ax) = ax_i(g_i(ax_i) - g_i(x_i)) + aF_i^{C_i}(x) + (1-a)C_i,$$

valid for any $x \in \mathbb{R}_+^n$ and any scalar $a \ge 0$. In particular, for any scalar $a > 0$,

$$F_i^{C_i}(ax^*) = ax_i^*(g_i(ax_i^*) - g_i(x_i^*)) + (1-a)C_i.$$

Since $g_i$ is decreasing and $C_i \ge 0$, this leads to

$$\begin{cases} F_i^{C_i}(ax^*) > 0 & \text{for any } 0 < a < 1 \\ F_i^{C_i}(ax^*) < 0 & \text{for any } a > 1. \end{cases}$$

By [36, Lemma 2], $x^*$ is then globally asymptotically stable on $\mathbb{R}_+^n \setminus \{0_n\}$. Moreover, since $F^C(0_n) = C > 0_n$, the solution of (D.35) initiated in the origin is increasing. By boundedness of the solutions, it converges towards the unique equilibrium point $x^*$. It implies that $x^*$ is GAS on $\mathbb{R}_+^n$.

### D.2.2 Proof of Theorem 4.2

First, as in the proof of Theorem 4.1 in Appendix D.1, one has

$$x^\infty|_{C_I} = 0_{nC_I}$$

and $x^\infty|_{\overline{C_I}}$ is equal to the equilibrium point that attracts every positive trajectory of the residual system (11).

$$\dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I} D_{ji} \right) + \sum_{\substack{j \in \overline{C_I} \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in \overline{C_I} \\ j \neq i}} D_{ji}x_i, \quad i \in \overline{C_I}.$$

To prove the existence of this equilibrium point, one splits $\Gamma_{\overline{C_I}}$ into strongly connected components and order them such that they follow an acyclic ordering. Therefore, for any strongly connected subgraph $\mathcal{G}$ and any of its upstream subgraphs $\mathcal{G}^-$, one has $D_{ji} = 0$ for every $j \in V_{\mathcal{G}^-}$ and $i \in V_{\mathcal{G}}$. Similarly, for any subgraph $\mathcal{G}^+$ in its downstream subgraphs, one has $D_{ij} = 0$ for every $j \in V_{\mathcal{G}^+}$ and $i \in V_{\mathcal{G}}$. It follows that

$$\dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{\substack{j \in C_I \cup V_{\mathcal{G}^+} \\ j \neq i}} D_{ji} \right) + \sum_{\substack{j \in V_{\mathcal{G}^-} \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in V_{\mathcal{G}^-} \\ j \neq i}} D_{ji}x_i + \sum_{\substack{j \in V_{\mathcal{G}^-} \\ j \neq i}} D_{ij}x_j, \quad i \in V_{\mathcal{G}}.$$

<footnote>$^4$Notice that by (D.40), any equilibrium of system (D.35) is necessarily positive when (H4) holds.</footnote>
To find the value of $x^\infty|_{V_G}$, we replace $x_j$ by $x_j^\infty$ for every $j \in V_G^-$. This results in the following system:

$$ \dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I \cup V_G^+} D_{ji} \right) + \sum_{\substack{j \in V_G \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in V_G \\ j \neq i}} D_{ji}x_i + \sum_{j \in V_G^-} D_{ij}x_j^\infty, \quad i \in V_G. \qquad (D.45) $$

The Jacobian of this system evaluated in zero, is equal to $A|_{V_G, V_G}(\bar{\rho})$. Moreover, since its connectivity matrix is irreducible, we can apply Theorem 2.2 and Theorem D.1. In particular, since $x^\infty|_{V_G}$ is the equilibrium attracting every positive trajectory of (D.45), one has either $x^\infty|_{V_G} = 0_{n_{V_G}}$ or $x^\infty|_{V_G} \gg 0_{n_{V_G}}$.

Assume first that either $V_G^-$ is of in-degree zero, or $x^\infty|_{V_G^-} = 0_{n_{V_G^-}}$. In such a case, the system (D.45) is equal to

$$ \dot{x}_i = x_i \left( g_i(x_i) - \bar{\rho}_i - \sum_{j \in C_I \cup V_G^+} D_{ji} \right) + \sum_{\substack{j \in V_G \\ j \neq i}} D_{ij}x_j - \sum_{\substack{j \in V_G \\ j \neq i}} D_{ji}x_i, \quad i \in V_G. \qquad (D.46) $$

System (D.46) is in the form of system (1), and applying Theorem 2.2 yields:

* if $s(A|_{V_G, V_G}(\bar{\rho})) > 0$, then the system (D.46) admits a positive equilibrium point, which is GAS on $\mathbb{R}_+^{n_{V_G}} \setminus \{0_{n_{V_G}}\}$. Then $x^\infty|_{V_G}$ is equal to this positive equilibrium point.

* if $s(A|_{V_G, V_G}(\bar{\rho})) \le 0$, then the origin of (D.46) is GAS on $\mathbb{R}_+^{n_{V_G}}$, yielding $x^\infty|_{V_G} = 0_{n_{V_G}}$.

Assume now that $V_G^- \neq \emptyset$ and $x^\infty|_{V_G^-} > 0_{n_{V_G^-}}$. In such case, there exists at least one strongly connected component $\mathcal{H}$ in $\mathcal{G}^-$ such that $x^\infty|_{V_{\mathcal{H}}} \gg 0_{n_{V_{\mathcal{H}}}}$. Therefore, there exists $i \in V_G$ such that

$$ \sum_{j \in V_G^-} D_{ij}x_j^\infty > 0. $$

Then, by applying Theorem D.1, the system (D.45) admits a unique equilibrium, which is positive and GAS on $\mathbb{R}_+^{n_{V_G}}$, and $x^\infty|_{V_G} \gg 0_{n_{V_G}}$ is equal to this equilibrium.

# E Properties of the minimization problem (Thm. 5.1 and Lemma 5.2)

## E.1 Proof of Theorem 5.1

The set of the $\bar{\rho}$-admissible vectors $\rho$ is convex, as well as the function $f(\rho)$. Furthermore, Theorem 3.1 states that $\rho \mapsto h(\rho)$ is convex. It follows that Problem (13) is a convex problem.
Now assume $s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho})) < 0$ and let $\alpha$ satisfy $0 \le \alpha < -s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho}))$. The set of the controls $\rho$ satisfying the constraints of Problem (13) is a closed set, and is non empty by assumption on $\alpha$ and $s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho}))$. Indeed, by Lemma 3.3, $\inf\{h(\rho); \rho \text{-admissible}\} = s(A|_{\overline{C_I}, \overline{C_I}}(\bar{\rho})) < -\alpha$, so that by continuity of the function $h$, there exists a $\bar{\rho}$-admissible control $\rho$ satisfying $h(\rho) \le -\alpha$. Since $f(\rho) \to +\infty$ as $\|\rho\| \to +\infty$, the Problem (13) admits at least one local minimizer $\rho^*$ (see [7]), which is necessarily global by convexity. The set of minimizers is convex, by convexity of the problem.

Let $\rho^*$ be a minimizer of Problem (13). If $s(A) = s(A(0_n)) > -\alpha$, then $\rho^* > 0_n$ since $\rho^*$ is non-negative. Moreover, assume by contradiction that $h(\rho^*) < -\alpha$. By continuity of $h$, there exists $\varepsilon \in (0, 1)$ such that $h((1 - \varepsilon)\rho^*) < -\alpha$. Then, $(1 - \varepsilon)\rho^*$ satisfies the constraints of Problem (13). As $\rho^* > 0_n$, it follows $f(\rho^*) > 0$. Therefore, $f((1 - \varepsilon)\rho^*) = (1 - \varepsilon)f(\rho^*) < f(\rho^*)$. This contradicts the fact that $\rho^*$ is a global minimizer of Problem (13). Therefore, $h(\rho^*) = \alpha$.

If $s(A) \le -\alpha$, then $0_n$ is the unique minimizer. Indeed, assume there exists another minimizer $\rho^{**} > 0_n$. In such case, one has $f(\rho^{**}) > 0_n = f(0_n)$, which is in contradiction with the fact that $\rho^{**}$ is a global minimizer.

Assume now that $C \subsetneq \{1, ..., n\}$ and let us show the uniqueness of the minimizer. Denote $f(\rho^*) := f^*$. When $s(A) \le -\alpha$, one has $\rho^* = 0_n$ and it is the unique minimizer, as seen above. Now assume that $s(A) > -\alpha$. Since $\rho^* > 0_n$, it follows that $f^* > 0$. Let us suppose by contradiction that Problem (13) admits another global minimizer $\rho^{**} \ne \rho^*$. Since $C \ne \{1, ..., n\}$, $h(\rho)$ is strictly convex by Theorem 3.1. Then,

$$h\left(\frac{\rho^* + \rho^{**}}{2}\right) < \frac{1}{2}h(\rho^*) + \frac{1}{2}h(\rho^{**}) \le -\alpha.$$

Consequently, by continuity of $r$, there exists a scalar $\varepsilon \in (0, 1)$ such that

$$h\left((1 - \varepsilon)\frac{\rho^* + \rho^{**}}{2}\right) < -\alpha,$$

with $(1 - \varepsilon)\frac{\rho^* + \rho^{**}}{2}$ satisfying the constraints of Problem (13). By linearity of $f$ and since $f(\rho^*) = f(\rho^{**}) = f^* > 0$,

$$f((1 - \varepsilon)\frac{\rho^* + \rho^{**}}{2}) = (1 - \varepsilon)f^* < f^*.$$

It is in contradiction with $\rho^*$ being a global minimizer of Problem (13). This achieves the proof of Theorem 5.1.

## E.2 Proof of Lemma 5.2

Let us compute the gradient (with respect to the variable $\rho$) of the function $f_{\phi,t} = tf + \phi$. Recall that we assume with no loss of generality that the control is introduced in the first $n_C$ patches, and among them, the infinite controllable patches are the first $n_{C_I}$ patches. In any case, this may be achieved by a permutation of components. Moreover, we identify in this
section any $\bar{\rho}$-admissible $\rho$ as well as the vector $\pi$ with their components in $\mathcal{C}$. Consequently, these vectors belong to $\mathbb{R}_+^{nc}$, and the functions $h = s(A(\cdot))$, $f, \phi$ and $f_{\phi,t}$ are now defined in $\mathbb{R}_+^{nc}$.

By composition,

$$ \nabla\phi(\rho) = -\frac{1}{h(\rho) + \alpha}\nabla h(\rho) - \frac{1}{\rho} + \left(0_{n_{c_I}}, \frac{1}{\bar{\rho}|_{c_F} - \rho|_{c_F}}\right)^T. $$

Since $\nabla f(\rho) = \pi$, one deduces that

$$ \nabla f_{\phi,t}(\rho) = t\pi - \frac{1}{h(\rho) + \alpha}\nabla h(\rho) - \frac{1}{\rho} + \left(0_{n_{c_I}}, \frac{1}{\bar{\rho}|_{c_F} - \rho|_{c_F}}\right)^T. $$

Moreover, with some computations, one finds that for any $i = 1, \dots, n_c$,

$$ (\nabla h(\rho))_i = -\frac{\partial s(A(\rho))}{\partial \nu_{ii}}, $$

where $\frac{\partial s(A(\rho))}{\partial \nu_{ii}}$ is the derivative of $s(A(\rho))$ with respect to the $(i, i)$-th component of $A(\rho)$. Let $Q(\rho) := h(\rho)I - A(\rho)$ and note $Q^{\#}(\rho)$ its group inverse. Equation (15) is then deduced from the following identity, extracted from [13, Lemma 3.1]:

$$ \frac{\partial s(A(\rho))}{\partial \nu_{ii}} = (I - Q(\rho)Q^{\#}(\rho))_{ii}. $$

Let us now compute the Hessian matrix of the function $f_{\phi,t}$. By composition,

$$ \nabla^2\phi(\rho) = \frac{1}{(h(\rho) + \alpha)^2}\nabla h(\rho)\nabla h(\rho)^\mathrm{T} - \frac{1}{h(\rho) + \alpha}\nabla^2h(\rho) + \mathrm{diag}\left(\frac{1}{\rho^2}\right) \\ + \mathrm{diag}\left(0_{n_{c_I}}, \frac{1}{(\rho|_{c_F} - \bar{\rho}|_{c_F})^2}\right). $$

Since $\nabla^2 f(\rho) = 0_{n_c \times n_c}$, the conclusion follows for $\nabla^2 f_{\phi,t}(x)$. Let us now compute $\nabla^2 h(\rho)$. For any $i, j = 1, \dots, n_c$,

$$ (\nabla^2 h(\rho))_{ij} = \frac{\partial^2 s}{\partial \nu_{jj} \partial \nu_{ii}}(A(\rho)). $$

Thanks to [13, Theorem 3.2], Equality (16) then follows. This achieves the proof of Lemma 5.2.