# Stability analysis of a bacterial growth model through computer algebra

Agustín Gabriel Yabo, Mohab Safey El Din, Jean-Baptiste Caillau, Jean-Luc Gouzé

► To cite this version:

Agustín Gabriel Yabo, Mohab Safey El Din, Jean-Baptiste Caillau, Jean-Luc Gouzé. Stability analysis of a bacterial growth model through computer algebra. MathematicS In Action, 2023, 12 (1), pp.175-189. 10.5802/msia.37. hal-03671432v3

HAL Id: hal-03671432

https://hal.science/hal-03671432v3

Submitted on 20 Oct 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
MathematicS

MathS in A.

In Action

AGUSTÍN G. YABO, MOHAB SAFEY EL DIN,
JEAN-BAPTISTE CAILLAU & JEAN-LUC GOUZÉ

Stability analysis of a bacterial growth model through computer algebra

Volume 12 (2023), p. 175-189.

https://doi.org/10.5802/msia.37

© Les auteurs, 2023.

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

Cet article est mis à disposition selon les termes
de la licence CREATIVE COMMONS ATTRIBUTION 4.0.

http://creativecommons.org/licenses/by/4.0/

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

MathematicS In Action *est membre du*
*Centre Mersenne pour l'édition scientifique ouverte*

http://www.centre-mersenne.org/

e-ISSN : 2102-5754
# Stability analysis of a bacterial growth model through computer algebra

AGUSTÍN G. YABO *
MOHAB SAFEY EL DIN **
JEAN-BAPTISTE CAILLAU ***
JEAN-LUC GOUZÉ ‡

* MISTEA, Université Montpellier, INRAE, Institut Agro, Montpellier, France

*E-mail address:* agustin.yabo@inrae.fr

** Sorbonne Université, LIP6 (UMR CNRS 7606), PolSys Team, France

*E-mail address:* mohab.safey@lip6.fr

*** Université Côte d'Azur, CNRS, Inria, LJAD, France

*E-mail address:* jean-baptiste.caillau@univ-cotedazur.fr

‡ Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, Biocore Team, Sophia

Antipolis, France

*E-mail address:* jean-luc.gouze@inria.fr.

## Abstract

We describe microbial growth and production of value-added chemical compounds in a continuous bioreactor through a dynamical system and we study the local stability of the equilibrium of interest by means of the classical Routh-Hurwitz criterion. The mathematical model considers various biological and structural parameters related to the bioprocess (concentration of substrate inflow, constants of the microchemical reactions, steady-state mass fractions of intracellular proteins, etc.) and thus, the stability condition is given in terms of these parameters. This boils down to deciding the consistency of a system of polynomial inequalities over the reals, which is challenging to solve from an analytical perspective, and out of reach even for traditional computational software designed to solve such problems. We show how to adapt classical techniques for solving polynomial systems to cope with this problem within a few minutes by leveraging its structural properties, thus completing the stability analysis of our model. The paper is accompanied by a Maple worksheet available online.

# 1. Introduction

## 1.1. Biological context

The dynamical system analyzed in this paper is based on previous works [31, 32, 37], and it represents a simplified version of a bioprocess used in scientific research and in the chemical and pharmaceutical industries for the production of value-added chemical compounds. The model is twofold. On one side, it considers a bacterial model representing the main cellular functions involved in growth and chemical production: metabolism and production of proteins [9, 35, 36, 38]. On the other hand, it models a continuous bioprocess, a production scheme occurring in a bioreactor that allows steady-state operation for long periods of time, avoiding shutdown for cleaning and maintenance [6]. Thus, through a multi-scale modelling approach, the biological model aims to capture intracellular reactions, extracellular processes, and the interplay between them [33, 34].

<footnote>This work was partially supported by ANR project Maximic (ANR-17-CE40-0024-01), Inria IPL Cosy and Labex SIGNALIFE (ANR-11-LABX-0028-01).</footnote>

<footnote>Keywords: continuous bioreactors, self-replicator model, bacterial growth, local stability.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification: 37N25, 92-08, 13P15, 14Q30, 68W30.</footnote>
The principle behind continuous processing is that the bacterial culture is supplied with a continuous flow of fresh medium rich in nutrients, while waste products and microbial cells are removed at the same volumetric flow rate. This creates a constant volume inside the bioreactor and, provided the appropriate operating conditions are met, a steady-state production regime. The outflow of compounds of interest excreted by bacteria is linear in the volumetric flow, so the process yield can be improved through this operation parameter. However, if too large, it can lead to washout: an undesired condition in which there is no more bacteria in the bioreactor. This usually occurs when the rate of growth of the bacterial culture is not able to “catch-up” with the rate of dilution of the bioreactor. Hence, understanding the asymptotic behavior of the system becomes crucial in successfully operating continuous bioreactors [15].

Additionally, biological models are known to be subject to non-negligible parametric uncertainty [30]. This occurs not only due to the inherent complexity of microorganisms, but also to their genetic variability: generation after generation, bacteria face genetic modifications that can progressively deviate the mathematical model from the real system [10]. Online parameter estimation can compensate for these issues [8], but the operation parameters should be adjusted accordingly. In this context, computing the stability of the equilibria in terms of the internal parameters represents an important advantage.

## 1.2. Reduction to polynomial system solving

In the analyzed dynamical system, the existence of the equilibrium of interest is described by a series of *strict* inequalities in terms of the system parameters; and its local stability is given by the well-known Routh-Hurwitz stability criterion [7], which yields an additional inequality. Deciding (through an algorithmic procedure) the consistency over the real numbers of polynomial systems with real coefficients is a central problem in the area of *effective real algebraic geometry*. Such systems define *basic semi-algebraic sets*. There are two families of algorithms for tackling such a decision problem. The very first one goes back to Hilbert's 17-th problem which asks whether a non-negative multivariate polynomial with real coefficients can always been written as a sum of squares of polynomials with real coefficients. Artin [2] showed that non-negative polynomials with real coefficients are actually sums of squares of rational fractions with real coefficients. A central result of real algebra is the *Positivstellensatz* [13, 27] which exhibits the pattern of a certificate of emptiness for a basic semi-algebraic set. Computing such certificates is proved to be hard. A more popular approach to certify non-negativity is based on reductions to semi-definite programming through the moment method and the so-called Lasserre hierarchy (see e.g. [14, 20]). These approaches are convenient to assess that a given polynomial (or a family of polynomials) is non-negative while in our setting we have strict inequalities. Besides, the certificates that are returned are approximate ones since the resulting semi-definite programs are solved numerically. Lifting them to *exact certificates* is also far from easy [21], especially in our setting (which involve strict inequalities over an unbounded domain).

We then focus on another family of algorithms which are “root finding” methods that can handle systems involving strict inequalities. These procedures will compute points in the solution set of the polynomial system under consideration whenever such points exist. When the solution set is empty, they just return an empty list (hence, without a witness of emptiness). To ensure exactness, these algorithms make use of computer algebra methods which manipulate symbolically the input polynomials with an exact arithmetic. Within this framework, there are two families of algorithms. The very first one, initiated by Collins [5], computes a partition of the semi-algebraic set defined by the input into pieces which are homeomorphic to $[0, 1]^i$, for $i$ ranging from 0 to the dimension of the ambient space, and which are arranged cylindrically through repeated projections on coordinate subspaces. The computational cost of this approach is doubly exponential in the dimension of the ambient space and polynomial in the number of input polynomials and their maximum degree. A more modern approach, introduced in [11] and
refined in [3] allows one to compute sample points per connected components of basic semi-
algebraic sets in time which is singly exponential in the number of variables and polynomial
in the number of input polynomials and their maximum degree. We refer the reader to [4] for
the foundations of such approaches and to [23] for more modern algorithms based on the crit-
ical point method and to the Maple package RAGLIB [22] which implements them. The two
main functions which are provided by this Maple package RAGLIB are HasRealSolutions and
PointsPerComponents. Both take as input polynomial systems of equations and strict inequal-
ities with coefficients in Q. They respectively decide whether the solution set is empty (when
it is not a sample point assessing the non-emptiness is provided) and compute sample points
per connected components. Roughly speaking, these functions “reduce” the problem of solving
polynomial systems of equations and inequalities over the reals to the one of solving families of
polynomial systems with finitely many complex solutions. These systems will have real solutions
if and only if the initial one has. They can then be rewritten in a triangular form, which is
the internal exact representation of the solutions, from which certified approximations of the
real solutions can be extracted using univariate real root isolation and backward substitution.
The polynomial system we need to solve is out of reach for the softwares implementing the
aforementioned methods as none of them was able to solve our problem within a full month of
computation.

1.3. Contributions

We investigate some structural properties of the polynomial system to solve and in particular
the degree pattern of our polynomial constraints w.r.t. the involved variables. By leveraging this
degree pattern, we simplify significantly the first steps of the Cylindrical Algebraic Decomposi-
tion algorithm and show that on the system we need to solve, the classical complexity growth one
faces does not occur. We push forward this investigation and show how to exploit even better the
degree patterns of intermediate data to find even stronger simplifications. In the end, after those
simplifications, we obtain a polynomial system with less variables which can be solved rather
easily (within a few seconds) by the RAGLIB Maple package. These simplifications are based
on a careful geometric analysis of the solution set to the system we need to solve. This analy-
sis enables us to find more compact ways to describe their projection on coordinate subspaces
thanks to the degree patterns we already mentioned. In Section 2, we recall the self-replicator
model of the bioreactor under consideration and show how using Routh-Hurwitz criterion leads
to polynomial system solving issues. In Section 3, we describe how we take advantage of the
degree patterns in our problem to solve our polynomial system.

2. Self-replicator model for bacterial growth

2.1. Model definition

We consider a self-replicator mathematical model representing bacterial growth [37], which is
described in this paper in a simplified manner. A more detailed description of the dynamical
system can be found in [32], with a thorough analysis of the biological principles and assumptions
behind its derivation. In the model, the bacterial culture in the bioreactor has constant volume
and is subject to an inflow of fresh medium rich in substrate. The medium is supplied at constant
volumetric flow rate and with a certain substrate concentration. The bioreactor is also subject
to an outflow of same volumetric flow rate of bacterial culture, containing substrate, microbial
cells and metabolites of interest. In the continuous bioreactor framework, this is represented
through a parameter *D* called *dilution rate*, which is defined as the ratio of in-flow/out-flow
rate to culture volume. The quantities in the bioreactor are expressed as fractions of the culture
volume:
* *s*: the volume fraction of substrate, which is used by bacteria to grow and to produce the compounds of interest.

* *V*: the volume fraction of the bacterial population.

* *x*: the volume fraction of the compound of interest.

The composition of the bacterial culture is described at the cell level as in [38], and thus the quantities are mass fractions of the total bacterial mass:

* *p*: mass fraction of precursor metabolites, used to produce biomass and the compound of interest.

* *r*: mass fraction of proteins of the gene expression machinery, responsible for the production of biomass.

* *m*: mass fraction of proteins of the metabolic machinery, responsible for the uptake of substrate from the medium, and the production of compounds of interest.

As, in nature, proteins account for most of the biomass of bacterial cells, the cellular mass is assumed to be described by *m* and *r*, such that *m* + *r* = 1. However, if required for biological purposes, the quantity *m* + *r* could be easily rescaled to be equal to a constant smaller than 1, as done in [36]. Thus, the assumption does not imply that *p* occupies a negligible fraction of the cell. The system can be externally controlled through two essential parameters: the above-defined dilution rate *D*, which is subject to a positivity constraint, and an allocation control *u*, satisfying *u* ∈ [0, 1]. The role of the allocation control is to decide whether the precursors are being used to produce proteins *r* or *m*. In nature, the allocation of resources is governed by natural mechanisms optimized by evolution. In biotechnological processes, this parameter can be controlled through biosynthetic methods able to affect the expression of RNA polymerase [12, 17]. As this paper focuses on the steady-state behavior of the system, *u* is a constant parameter.

<figure><img src="image_4.png" /><figcaption>FIGURE 2.1. Extended coarse-grained self replicator model [38].</figcaption></figure>
Following [32, 37], the system of differential equations that describes the evolution of the state variables is given by

$$ egin{cases} 
\dot{s} = D(1-s) - v_M(s,m)\mathcal{V}, \\
\dot{p} = v_M(s,m) - v_X(p,m) - v_R(p,r)(p+1), \\
\dot{r} = (u-r)v_R(p,r), \\
\dot{m} = (1-u-m)v_R(p,r), \\
\dot{x} = v_X(p,m)\mathcal{V} - Dx, \\
\dot{\mathcal{V}} = (v_R(p,r) - D)\mathcal{V},
\end{cases} 
\quad (\text{S}_{\text{original}}) $$

where the functions $v_M$, $v_R$ and $v_X$ correspond to the relative synthesis rates of precursor metabolites, macromolecules and compounds of interest, respectively. Based on [9], we define the synthesis rates as linear in the concentrations $m$ and $r$ [24], and using Michaelis–Menten kinetics [19] for the dependence on the component used for each reaction:

$$ egin{aligned} v_M(s, m) &\doteq k_2 m \frac{s}{K_2 + s}, \\ v_R(p, r) &\doteq r \frac{p}{K + p}, \\ v_X(p, m) &\doteq k_1 m \frac{p}{K_1 + p}. \end{aligned} \quad (2.1) $$

As the system is analyzed in a steady-state regime, the control is fixed to a constant parameter $u(t) = u^*$. Thus, the set of parameters is defined as $\theta \doteq (K, k_1, K_1, k_2, K_2, u^*, D)$, which is subject to the constraints

$$ K > 0, \quad k_1 > 0, \quad K_1 > 0, \quad k_2 > 0, \quad K_2 > 0, \quad 1 \ge u^* \ge 0, \quad D > 0. \quad (2.2) $$

Using $m + r = 1$, we can express $m = 1 - r$, and thus remove $m$ from the dynamical model. Additionally, by analyzing the dynamics of the quantities $s + (p+1)\mathcal{V} + x$ and $s + (p+r/u^*)\mathcal{V} + x$, we can see that they both tend asymptotically to 1 when $t \to \infty^1$. Thus, the $\omega$-limit set (i.e. the set of points that can be limit of subtrajectories of the system) is characterized as follows.

**Proposition 2.1.** The $\omega$-limit set of any solution of system (S$_{original}$) lies in the hyperplanes

$$ \Omega_1 \doteq \{(s, p, r, x, \mathcal{V}) \in \mathbb{R}^5 : s + (p+1)\mathcal{V} + x = 1\}, \quad (2.3) $$

$$ \Omega_2 \doteq \{(s, p, r, x, \mathcal{V}) \in \mathbb{R}^5 : s + \left(p + \frac{r}{u^*}\right) \mathcal{V} + x = 1\}. \quad (2.4) $$

The latter implies that, for every trajectory, $r$ converges to $u^*$ asymptotically, and $x$ can be expressed in terms of the remaining states as

$$ x = 1 - s - (p+1)\mathcal{V}. \quad (2.5) $$

The case $u^* = 0$ is excluded from the study for its triviality, as it is rather simple to analyze. Indeed, for this case, the hyperplane $s + (p + m/(1 - u^*))\mathcal{V} + x = 1$ can be used for a similar study. Thus, the asymptotic behavior of the original system (S$_{original}$) can be studied through its limiting system given by

$$ egin{cases} 
\dot{s} = D(1-s) - v_M(s, 1-u^*)\mathcal{V}, \\
\dot{p} = v_M(s, 1-u^*) - v_X(p, 1-u^*) - (p+1)v_R(p, u^*), \\
\dot{\mathcal{V}} = (v_R(p, u^*) - D)\mathcal{V}. 
\end{cases} 
\quad (S) $$

<footnote n="1">¹While in the present non-dimensional formulation these quantities have no physical meaning, they are based on their equivalent on the original model formulation, which represents total mass in the bioreactor.</footnote>
In order to ensure that the original system ($S_{original}$) converges to the equilibria of the limiting system (S), an additional analysis is required. The reader can refer to [32] for an application of the theory of asymptotically autonomous systems [29].

## 2.2. Stability of the equilibrium of interest

In (S), there is a single equilibrium of interest ($s^*, p^*, V^*$), as the other existing equilibrium (usually referred to as the washout equilibrium) is characterized by having $V^* = 0$ (and thus $s^* = 1$), which does not allow for bacterial growth and metabolite production. By solving $\dot{s} = \dot{p} = \dot{V} = 0$, we obtain the steady-state values in terms of the system parameters $\theta$.

**Proposition 2.2.** *The equilibrium of interest is*

$$ egin{cases} s^*(\theta) = \frac{K_2\gamma(\theta)}{k_2(1-u^*) - \gamma(\theta)}, \\ p^*(\theta) = \frac{DK}{u^* - D}, \\ V^*(\theta) = \frac{D(1-s^*(\theta))(K_2 + s^*(\theta))}{k_2(1-u^*)s^*(\theta)}, \end{cases} \qquad (2.6) $$

where

$$ \gamma(\theta) \doteq k_1(1-u^*) \frac{p^*(\theta)}{K_1 + p^*(\theta)} + (p^*(\theta) + 1)u^* \frac{p^*(\theta)}{K + p^*(\theta)}. \qquad (2.7) $$

As the variables represent concentrations and biomass, they are non-negative quantities. Thus, additional inequalities are required for the existence of the equilibrium, which can be obtained by enforcing $(s^*, p^*, V^*) > 0$:

$$ 1 > u^* > D, \quad 1 > s^*(\theta) > 0. \qquad (2.8) $$

**Definition 2.3.** A set of parameters $\theta$ is called admissible if it satisfies (2.2) and (2.8).

The steady state is thus determined by 7 variables. Preliminary numerical explorations on a wide range of admissible parameters suggest that the system enjoys strong stability properties; see Figure 2.2. Not having been able to prove its global stability despite some structural properties of the dynamics, we focus on local stability of the previous equilibrium. To this end, we compute the eigenvalues over $\mathbb{C}$ of the Jacobian matrix and check whether their real part is negative. For this case, the Jacobian matrix is defined as

$$ \begin{bmatrix} -D - \frac{\partial}{\partial s^*} v_M(s^*, 1-u^*) V^* & 0 & -v_M(s^*, 1-u^*) \\ \frac{\partial}{\partial s^*} v_M(s^*, 1-u^*) & \beta(\theta) & 0 \\ 0 & \frac{\partial}{\partial p^*} v_R(p^*, u^*) V^* & v_R(p^*, u^*) - D \end{bmatrix} \qquad (2.9) $$

with

$$ \beta(\theta) = \frac{\partial}{\partial p^*} v_X(p^*, 1-u^*) - v_R(p^*, u^*) - (p^* + 1) \frac{\partial}{\partial p^*} v_R(p^*, u^*) \qquad (2.10) $$

We define the functions

$$ egin{align} \phi_1(\theta) &= \frac{V^* k_2^2 K K_2 s^* (1 - u^*)^2 u^*}{(K + p^*)^3 (K_2 + s^*)^2}, & \phi_2(\theta) &= \frac{K (p^* + 1) u^*}{(K + p^*)^2}, \tag{2.11} \\ \phi_3(\theta) &= \frac{k_1 K_1 (1 - u^*)}{(K_1 + p^*)^2}, & \phi_4(\theta) &= \frac{V^* k_2 K_2 (1 - u^*)}{(K_2 + s^*)^2}, \nonumber \end{align} $$
<figure><img src="image_5.png" /><figcaption>FIGURE 2.2. Numerical trajectories of (S) obtained for fixed initial conditions and random set of admissible parameters $\theta$ (satisfying so the set of inequalities (2.2) and (2.8)). In every case, the system converges to the equilibrium of interest, which is characterized by the persistence of the bacterial population $V^* > 0$.</figcaption></figure>

that satisfy $\phi_i(\theta) > 0$ for $i \in \{1, \dots, 4\}$. Computing the characteristic polynomial evaluated at the equilibrium point given by (2.6) yields

$$ 
\begin{aligned}
P(\lambda) &= \phi_1 + (\phi_2 + \phi_3 + \lambda + D)(\phi_4 + D + \lambda)\lambda \\
&= \lambda^3 + (2D + \phi_2 + \phi_3 + \phi_4)\lambda^2 + (D^2 + D\phi_2 + D\phi_3 + D\phi_4 + \phi_2\phi_4 + \phi_3\phi_4)\lambda + \phi_1.
\end{aligned} 
\quad (2.12) $$

The Routh-Hurwitz criterion for degree-three polynomials states that the roots of a degree-three polynomial $\lambda^3 + \alpha_2\lambda^2 + \alpha_1\lambda + \alpha_0$ belong to the open left complex plane if and only if

$$ \alpha_2 > 0, \quad \alpha_0 > 0 \quad \text{and} \quad \alpha_1\alpha_2 > \alpha_0. $$

From (2.12), since every function $\phi_i$ is positive, the stability criterion reads

$$ (2D + \phi_2 + \phi_3 + \phi_4)(D^2 + D\phi_2 + D\phi_3 + D\phi_4 + \phi_2\phi_4 + \phi_3\phi_4) - \phi_1 > 0. \quad (2.13) $$

In order to prove the latter, we can prove that no state-parameter combination satisfies the negation of the condition,

$$ -(2D + \phi_2 + \phi_3 + \phi_4)(D^2 + D\phi_2 + D\phi_3 + D\phi_4 + \phi_2\phi_4 + \phi_3\phi_4) + \phi_1 \geq 0. \quad (\text{POL}) $$

The stability of the steady state of interest is given by the validity of the latter inequality, which depends on the value of $\theta$.

Routh-Hurwitz criterion reduces the study of the local stability of an equilibrium to an algebraic problem. In some cases, depending on the complexity of the dynamical system (such as its dimension and non-linearity), such a problem may be solved analytically. The difficulty in our case stems from the large number of non-zero coefficients of the left-hand side of (POL) in the standard monomial basis. The tailored computer algebra methods described in the rest of the paper allows us to obtain the following stability property of the family of equilibria of interest:

**Theorem 2.4.** On a dense subset of the set of admissible parameters $\theta = (K, k_1, K_1, k_2, K_2, u^*, D)$, the equilibrium $(s^*(\theta), p^*(\theta), V^*(\theta))$ is locally stable.

The next section is devoted to the proof of this result. The Maple worksheet used to prove the result is available online<fnref n="2" /> as a companion notebook to the current paper. It makes an intensive use of the RAGLIB package.

<footnote n="2">²See https://ct.gitlabpages.inria.fr/gallery</footnote>
### 3. Computer algebra analysis of stability

We start the proof of Theorem 2.4 with a preparatory change of variables. As is clear from (2.6) and (2.7), in the range of admissible parameters one can express $K$ in terms of $(p^*, u^*, D)$, and $k_1$ in terms of $(\gamma, u^*, p^*, K_1, K)$ while $\gamma$ itself is obtained as a function of $(s^*, K_2, k_2, u^*)$. Accordingly, the following is true:

**Lemma 3.1.** There is a birational change of variables

(i) from $\theta = (K, k_1, K_1, k_2, K_2, u^*, D)$ to $(p^*, s^*, K_1, k_2, K_2, u^*, D)$,

(ii) then from $(p^*, s^*, K_1, k_2, K_2, u^*, D)$ to $\theta' := (s^*, p^*, K'_1, k_2, K'_2, u^*, D)$ where
$$K'_1 := K_1 + p^*, \quad K'_2 := K_2 + s^*.$$

The resulting system of inequalities becomes

$$K > 0, \quad k_1 > 0, \quad K'_1 - p^* > 0, \quad k_2 > 0, \quad K'_2 - s^* > 0, \quad 1 > u^* > D > 0, \quad 1 > s^* > 0, \quad C < 0, \quad (3.1)$$

where $C$ is a polynomial of total degree 13 involving the above 7 variables, whose number of non-zero coefficients in the standard monomial basis is 164, that is associated with the numerator of (POL) left-hand side. (We refer to the companion notebook for the detailed expression of $C$.) Note that we restrict the analysis to strict inequalities. Clearly, it suffices to prove the proposition below to obtain Theorem 2.4.

**Proposition 3.2.** The semi-algebraic set defined by (3.1) is empty.

Feeding HasRealSolutions with system (3.1) without further simplification does not allow to obtain a solution in reasonable time.<fnref n="3" /> The crucial observation to obtain a computationally tractable proof of Proposition 3.2 is related to the degree pattern of the input system. Indeed, observe that the system (3.1) involves 11 polynomial inequalities in the polynomial ring $\mathbb{Q}[s^*, p^*, K'_1, k_2, K'_2, u^*, D]$ but only two of them have positive degree in $K'_1$. In other words, amongst these inequalities, 9 lie actually in $\mathbb{Q}[s^*, p^*, k_2, K'_2, u^*, D]$. This gives rise to the following idea: one may reduce our initial problem (deciding the consistency of the system (3.1)) to the study of some semi-algebraic set defined by polynomial inequalities in $\mathbb{Q}[s^*, p^*, k_2, K'_2, u^*, D]$. Hence, we are in the process of eliminating the variable $K'_1$ from (3.1), that is computing polynomial inequalities that would define the projection of the solution set to (3.1) on the $(s^*, p^*, k_2, K'_2, u^*, D)$ coordinate space. A key tool for algebraic elimination is the resultant of two univariate polynomials, which we now recall before dwelving into the proof of Proposition 3.2.

Let $R$ be a ring and $f$ and $g$ be two polynomials in $R[x]$ of respective degree $p$ and $q$. The resultant associated to $(f, g)$ is the determinant of the Sylvester matrix which is the one obtained by stacking on each row the coefficients of the polynomials $x^if$ and $x^jg$ for $0 \le i \le q-1$ and $0 \le j \le p-1$ (see e.g. [4, Chapter 4]). This resultant is 0 if and only if $f$ and $g$ have a gcd of positive degree. The discriminant of a polynomial $f$ is the resultant of $f$ and its derivative $\frac{\partial f}{\partial x}$ divided by an appropriate power of the leading coefficient of $f$ (see e.g. [4, Chapter 4]). When dealing with polynomials $f$ and $g$ in $\mathbb{Q}[x_1, \dots, x_n]$, one can see them as polynomials in $\mathbb{Q}[x_1, \dots, x_{n-1}][x_n]$. The resultant (resp. discriminant) associated to $(f, g)$ (resp. $f$) is then denoted by resultant($f, g, x_n$) (resp. discriminant($f, x_n$)). Note that they lie in $\mathbb{Q}[x_1, \dots, x_{n-1}]$ (as does the leading coefficient of $f$ w.r.t. $x_n$ which we denote by $\text{lc}(f, x_n)$). These are the tools we use to describe the projection on a subspace of coordinates of a semi-algebraic set. These tools eliminating one variable from a finite sequence of polynomials are the ones used by the Cylindrical

<footnote n="3">³Computation not terminated after one month on an Intel(R) Xeon(R) Gold 6244 CPU @ 3.60GHz with 1.5 To of RAM.</footnote>
Algebraic Decomposition recursively algorithm and its so-called “open” Cylindrical Algebraic Decomposition which applies to systems of strict inequalities [28]. We will use properties of resultants and discriminants further in our approach.

Let now $f_1, \dots, f_s$ be polynomials in $\mathbb{Q}[x_1, \dots, x_n]$ and let $S \subset \mathbb{R}^n$ be the semi-algebraic set defined by

$$f_1 > 0, \dots, f_s > 0 \qquad (3.2)$$

with $f_i \in \mathbb{Q}[x_1, \dots, x_{n-1}]$ for $1 \le i \le t$ and $\deg(f_j, x_n) > 0$ for all $t+1 \le j \le s$. Note that since $S$ is open for the Euclidean topology, its projection on the $(x_1, \dots, x_{n-1})$ coordinate space is open too. Let $S'$ be the semi-algebraic set defined by

$$\begin{array}{ll} f_i > 0 & \text{for } 1 \le i \le t, \text{ and} \\ \text{resultant}(f_i, f_j, x_n) \neq 0, \text{ discriminant}(f_i, x_n) \neq 0 & \text{for } t+1 \le i \le j \le s \\ \text{lc}(f_i, x_n) \neq 0 & \text{for } t+1 \le i \le s. \end{array} \qquad (3.3)$$

By [28], there exist connected components $C_1, \dots, C_\ell$ of $S'$ such that $C_1 \cup \dots \cup C_\ell$ is semi-algebraic and dense in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space. We call the OpenCAD, the procedure which, given the polynomial system $f_1 > 0, \dots, f_s > 0$ and the variable $x_n$ as inputs, returns the polynomial system (3.3). Hence, in order to decide the consistency over the reals of the system of inequalities (3.2), it suffices to

(a) compute from it the system of polynomial inequations and inequalities (3.3) using the procedure OpenCAD;

(b) compute at least one sample point per connected component in the semi-algebraic set $S'$ defined by (3.3) (note that we may get several points per connected component). This is what the command PointsPerComponents from the RAGLIB package is designed for;

(c) letting $\alpha_1, \dots, \alpha_r$ be those sample points, solve the univariate system of polynomial inequalities

$$f_{t+1}(\alpha_i, x_n) > 0, \dots, f_s(\alpha_i, x_n) > 0,$$

for each $1 \le i \le r$. If the system is not consistent over the reals then one can conclude that $S$ is empty (in other words, the system (3.2) is not consistent over the reals).

Launching this computation on system (3.1) is still not computationally tractable. To further reduce computations, we rely on the following property of the system, which can easily be checked (e.g., on the MAPLE worksheet attached with the paper).

**Lemma 3.3.** In system (3.1), only two polynomials have positive degree in $K'_1$: one has degree 2 in $K'_1$, and the other one has degree 1 in $K'_1$.

We assume from now on that $\deg(f_s, x_n) = 2$ (so $f_s$ has real roots if and only if its discriminant is non-negative) and $\deg(f_i, x_n) = 1$ for $t+1 \le i \le s-1$. We let $C$ be a connected component of $S'$ that is contained in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space. In other words, for any $\alpha \in C \subset \mathbb{R}^{n-1}$, there exists $\vartheta \in \mathbb{R}$ such that $f_i(\alpha, \vartheta) > 0$ for $t+1 \le i \le s$. By definition of $S'$, the discriminant of $f_s$ does not vanish over $C$. We then consider the two possible cases:

(i) when discriminant($f_s, x_n$) is negative over $C$, for any $\alpha \in C$, $f_s(\alpha, x_n)$ has no real root and is then positive (because as explained above, there exists $\vartheta \in \mathbb{R}$ such that $f_s(\alpha, \vartheta) > 0$ and since $f_s$ has no real root, our positivity statement holds);

(ii) when discriminant($f_s, x_n$) is positive over $C$, for any $\alpha \in C$, $f_s(\alpha, x_n)$ has exactly two real roots with multiplicity one and there exists $\vartheta$ such that $f_s(\alpha, \vartheta)$ is positive.
Note that the leading coefficient of $f_s$ did not play any role in the above observations. This leads us to consider the semi-algebraic set $S''$ which is defined by

$$ \begin{array}{ll} f_i > 0 & \text{for } 1 \le i \le t, \text{ and} \\ \text{resultant}(f_i, f_j, x_n) \ne 0, \text{ discriminant}(f_i, x_n) \ne 0 & \text{for } t+1 \le i \le j \le s \\ \text{lc}(f_i, x_n) \ne 0 & \text{for } t+1 \le i \le s-1. \end{array} \quad (3.4) $$

Note also that the only difference is that we do not consider the leading coefficient of $f_s$ anymore (the last inequalities in (3.4) range from $i = t+1$ to $s-1$). The connected components of $S''$ are unions of some connected components of $S'$ and some connected semi-algebraic subsets of the vanishing set of $\text{lc}(f_s, x_n)$ (over which none of the other polynomials used to define $S''$ vanishes). Moreover, by definition of $S''$, over any connected component $C$ of $S''$, the discriminant polynomial discriminant($f_s, x_n$) is sign invariant.

**Proposition 3.4.** Assume that the semi-algebraic set $S$ is non-empty. Then, there exists a connected component $C$ of $S''$ such that any $\alpha \in C$ which does not cancel $\text{lc}(f_s, x_n)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space.

*Proof.* Consider the semi-algebraic set defined by $f_1 > 0, \dots, f_{s-1} > 0$. Running the procedure OpenCAD on the input $[f_1, \dots, f_{s-1}], x_n$ builds a polynomial system which is the same as system (3.3) without the constraints discriminant($f_s, x_n$) $\ne 0$ and resultant($f_i, f_s, x_n$) $\ne 0$ for $1 \le i \le s-1$. Let $C$ be a connected component of $S''$. By [28, Lemma 3.6 and the discussion below on the proof of Algo. 3.5], there exist finitely many continuous semi-algebraic maps (i.e. maps whose graphs are semi-algebraic sets) $\xi_1, \dots, \xi_k$ from $C$ to $\mathbb{R}$ such that the cylinder $C \times \mathbb{R}$ is the disjoint union of cells defined as

* either the graph of one of the maps $\xi_i$ for $1 \le i \le k$,

* or a band of some cylinder

$$ B_i = \{(\alpha, \vartheta) \mid \alpha \in C \text{ and } \xi_i(\alpha) < \vartheta < \xi_{i+1}(\alpha)\} \\ \text{for } 0 \le i \le k \text{ with } \xi_0 = -\infty \text{ and } \xi_{k+1} = +\infty \text{ by convention,} $$

over which the polynomials $f_1, \dots, f_{s-1}$ are sign invariant (because the system output by the call to the procedure OpenCAD on the input $[f_1, \dots, f_{s-1}], x_n$ is contained in (3.3), and then defines a dense open semi-algebraic subset of $S''$).

Note that, following the argumentation of the correctness of [28, Algorithm 3.5], for $\alpha \in C$, $\xi_1(\alpha), \dots, \xi_k(\alpha)$ is the ordered sequence of roots of the polynomials $f_{t+1}(\alpha, x_n), \dots, f_{s-1}(\alpha, x_n)$ (note that these are univariate). Hence, the semi-algebraic maps $\xi_i$'s are defined as maps which send $\alpha \in C$ to some root of some $f_j(\alpha, x_n)$ (those roots being maintained ordered when $\alpha$ ranges over $C$).

Assume now that $C$ has a non-empty intersection with the projection of $S$ on the $(x_1, \dots, x_{n-1})$ subspace of coordinates (note that if $S$ is non-empty, there exists such a connected component $C$ of $S''$ since $S$ and $S''$ are both open for the Euclidean topology). Let $\alpha$ be in this intersection such that $\text{lc}(f_s, x_n)$ does not vanish at $\alpha$ (such a point $\alpha$ exists because $C$ is open for the Euclidean topology as well as the projection of $S$ and the solution set to $\text{lc}(f_s, x_n) \ne 0$). In other words, there exists $\vartheta$ such that

$$ f_1(\alpha, \vartheta) > 0, \dots, f_s(\alpha, \vartheta) > 0. $$

Let $B_i$ be the band which contains $(\alpha, \vartheta)$. Consider now $\alpha' \in C$ and assume that $\text{lc}(f_s, x_n)(\alpha') \ne 0$. We need to prove that there exists $\vartheta' \in \mathbb{R}$ such that $(\alpha', \vartheta') \in S$. Note that, by the sign invariance property over the band $B_i$ of $f_1, \dots, f_{s-1}$, we already know that there exists $\vartheta'$ such that $f_j(\alpha', \vartheta') > 0$ for $1 \le j \le s-1$. Hence, what is missing is a control on the sign of $f_s$, precisely that there exists $\vartheta' \in \mathbb{R}$ such that $f_s(\alpha', \vartheta') > 0$. The next lemma is a first step towards this. We postpone its proof to next section.
**Lemma 3.5.** Assume that there exists a continuous semi-algebraic path $\gamma : [0, 1] \to C$ such that $\gamma(0) = \alpha$ and $\gamma(1) = \alpha'$ such that for any $v \in [0, 1]$, $\mathrm{lc}(f_s, x_n)$ does not vanish at $\gamma(v)$, with $\alpha$ in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space. Then, for any $v \in [0, 1]$, $\gamma(v)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$-space.

Since $\alpha$ and $\alpha'$ both lie in $C$ which is connected, they can be connected with semi-algebraic paths in $C$. Assume now that any continuous semi-algebraic path in $C$ linking $\alpha$ to $\alpha'$ meets the vanishing set of $\mathrm{lc}(f_s, x_n)$. Note that since the vanishing set of $\mathrm{lc}(f_s, x_n)$ has dimension at most $n-1$, one can, without loss of generality, assume that the intersection of the continuous semi-algebraic path under consideration, linking $\alpha$ to $\alpha'$, with the vanishing set of $\mathrm{lc}(f_s, x_n)$ is finite. Hence, there exists $\{v_1, \dots, v_N\} \subset [0, 1]$ such that

$$ \{v_1, \dots, v_N\} = \{v \in [0, 1] \mid \mathrm{lc}(f_s, x_n)(\gamma(t)) = 0\} $$

and we define $v_0 = 0$ and $v_{N+1} = 1$. We let $\gamma_j$ be the semi-algebraic continuous path $v \in [v_j, v_{j+1}] \to \gamma(v)$ for $0 \le j \le N$. Note that $\gamma_j$ can be extended by continuity to $v_j$ and $v_{j+1}$ by just taking $\gamma_j(v_j) = \gamma(v_j)$ and $\gamma_j(v_{j+1}) = v_{j+1}$ (since $\gamma$ is continuous). By a slight abuse of notation, we still denote by $\gamma_j$ that semi-algebraic map which is extended to $[v_j, v_{j+1}]$. Note that $\gamma_0(v_0) = \alpha$ which, by assumption, lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$-space. We use the following lemma (also proved in the next section):

**Lemma 3.6.** Let $\gamma : [0, 1] \to C$ be a semi-algebraic continuous map such that

* for all $v \in [0, 1]$, $\mathrm{lc}(f_s, x_n)$ does not vanish at $\gamma(v)$ but does at $\gamma(1)$;

* $\gamma(0)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space.

Then, for all $v \in [0, 1]$, $\gamma(v)$ also lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space.

We apply inductively Lemma 3.6 to the semi-algebraic continuous paths $\gamma_j$ (up to translating and scaling $[v_j, v_{j+1}]$ to $[0, 1]$) as follows. First note that for $j = 0$, $\gamma_0$ satisfies the assumptions of Lemma 3.6 (again, scaling $[0, v_1]$ to $[0, 1]$) since $\gamma_0(0) = \gamma(0) = \alpha$ which lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$-space by assumption. Hence, Lemma 3.6 allows us to deduce that for all $v \in [0, v_1]$, $\gamma(v)$ lies in the aforementioned projection of $S$. This implies that $\gamma_1$ satisfies the assumptions of Lemma 3.6.

Assume now that for $j < N$, $\gamma_j$ satisfies the assumptions of Lemma 3.6 (again, by translation and scaling, $[v_j, v_{j+1}]$ is sent to $[0, 1]$. Following the same reasoning as above, we deduce that $\gamma_{j+1}$ satisfies the assumptions of Lemma 3.6. We deduce that $\gamma_N$ does satisfy the assumptions of Lemma 3.6 which then implies that $\gamma_N(1) = \gamma(1) = \alpha'$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$-space, as requested. This ends the proof of Proposition 3.4. $\square$

Denote by ModifiedOpenCAD the procedure which takes as input $f_1, \dots, f_s$ as above and returns the system of polynomial constraints as in (3.3). By Proposition 3.4, it suffices to compute sample points $\alpha_1, \dots, \alpha_r$ per connected components of the semi-algebraic set defined by (3.3) at which $\mathrm{lc}(f_s, x_n)$ does not vanish (recall that doing so is the purpose of the PointsPerComponents command of RAGLIB) and decide if the univariate system of polynomial inequalities

$$ f_{t+1}(\alpha_i, x_n) > 0, \dots, f_s(\alpha_i, x_n) > 0 $$

has real solutions (which is easily done with real root isolation algorithms) for each $1 \le i \le r$. Running this computation does not provide any result in a reasonable time, even when one takes advantage of special simplifications which we make explicit in the next section. However, note that one can use the same strategy in one dimension less, and apply the procedure OpenCAD to (3.3) as input which provides a new system of polynomial inequalities/inequations. Computing sample points per connected components for the semi-algebraic set defined in $\mathbb{R}^{n-2}$ by this
new system is done within a minute on a standard laptop with PointsPerComponents (when exploiting some simplifications which we make explicit in the next section). Next, one can “lift” those points by instantiating $x_1, \dots, x_{n-2}$ in (3.3) and solve the obtained univariate polynomial systems of constraints for $x_{n-1}$. This will provide sample points per connected components of the semi-algebraic set defined by (3.3), which can be lifted to $\mathbb{R}^n$ as sketched above and finally solve our decision problem. These latter steps take a few second. All in all, these computations prove that the semi-algebraic set defined by (3.1) is empty, which allows to conclude the proof of Proposition 3.2. Theorem 2.4 is proved.

## 4. Proof of auxiliary lemmas

*Proof of Lemma 3.5.* Note that for any $v \in [0, 1]$, and any $\xi_i(\gamma(v)) < \vartheta_v < \xi_{i+1}(\gamma(v))$, the following holds:

$$f_1(\gamma(v), \vartheta_v) > 0, \dots, f_{s-1}(\gamma(v), \vartheta_v) > 0.$$

Recall also that, by definition of $S''$, the discriminant discriminant($f_s, x_n$) is sign invariant over $C$. If it is negative, then for any $v \in [0, 1]$, $f_s(\gamma(v), x_n)$ has no root in $\mathbb{R}$ (because we have assumed that $\mathrm{lc}(f_s, x_n)$ does not vanish at $\gamma(v)$). Since, by construction, $f_s(\gamma(0), x_n)$ is positive over $\mathbb{R}$, we deduce that for any $v \in [0, 1]$, $f_s(\gamma(v), x_n)$ is positive over $\mathbb{R}$ (else, by continuity of $\gamma$, there would exist some $v'$ such that $f_s(\gamma(v'), x_n)$ has some real root while $\mathrm{lc}(f_s, x_n)$ does not vanish at $\gamma(v')$, which would contradict our assumption of the negativity of discriminant($f_s, x_n$) over $C$). Assume now that discriminant($f_s, x_n$) is positive over $C$. As above, we need to prove that for any $v \in [0, 1]$, there exists $\vartheta_v \in B_i$ such that $f_s(\gamma(v), \vartheta_v) > 0$. We start with $v = 0$. We already know by construction that there exists $\xi_i(\gamma(0)) < \vartheta < \xi_{i+1}(\gamma(0))$ such that $f_s(\gamma(0), \vartheta) > 0$. Assume by contradiction that there exists $v' \in [0, 1]$ such that for any $\xi_i(\gamma(v')) < \vartheta' < \xi_{i+1}(\gamma(v'))$, it holds that $f_s(\gamma(v'), \vartheta') \le 0$. By continuity of $\gamma$, $\xi_i$ and $\xi_{i+1}$, the set of such reals $v'$ is closed in $[0, 1]$. We let $v'_{\min}$ be the smallest element of this set. By continuity of the $\xi_i$'s we deduce that $f_s(\gamma(v'_{\min}), \xi_i(\gamma(v'_{\min}))) = 0$ or $f_s(\gamma(v'_{\min}), \xi_{i+1}(\gamma(v'_{\min}))) = 0$. Recall that there exists some $t+1 \le j \le s-1$ such that $\xi_i$ maps $\alpha \in C$ to some root of $f_j$. In other words, the $\mathrm{gcd}$ of $f_s(\gamma(v'_{\min}), x_n)$ and $f_j(\gamma(v'_{\min}), x_n)$ has degree $\ge 1$. This implies that for some $j$, the resultant polynomial resultant($f_j, f_s, x_n$) vanishes at $\gamma(v'_{\min})$. This is a contradiction since $C$ is a connected component of $S''$, defined by (3.3) which contains the constraint resultant($f_j, f_s, x_n$) $\ne 0$. This ends the proof of Lemma 3.5. $\square$

*Proof of Lemma 3.6.* Take $w \in ]0, 1[$. Up to scaling $[0, w]$ to $[0, 1]$, one applies Lemma 3.5 to deduce that $\gamma(w)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate subspace.

Since $\gamma(w)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate subspace, there exists $\vartheta_w \in \mathbb{R}$ such that $(\gamma(w), \vartheta_w) \in S$

$$f_1(\gamma(w), \vartheta_w) > 0, \dots, f_{s-1}(\gamma(w), \vartheta_w) > 0, f_s(\gamma(w), \vartheta_w) > 0$$

Note that $\vartheta_w$ can be chosen in the interval $]\xi_i(\gamma(w)), \xi_{i+1}(\gamma(w))[$ where $\xi_i$ is one of the “root” functions introduced in the previous section. By continuity of $\gamma$, the equality $\gamma(1) = \lim_{w \to 1}(\gamma(w))$ holds. We are going to prove that $\gamma(1)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate subspace. By assumption, $\mathrm{lc}(f_s, x_n)$ vanishes at $\gamma(1)$. Note that, by definition of the discriminant of a quadratic polynomial, this implies that discriminant($f_s, x_n$) is non-negative at $\gamma(1)$. Also, discriminant($f_s, x_n$) is sign invariant over $\gamma([0, 1])$ (since the inequality constraint discriminant($f_s, x_n$) $\ne 0$ is part of the polynomial system (3.3)). Hence, we deduce that discriminant($f_s, x_n$) is positive over $\gamma([0, 1])$. Recall that we established that, for all $w \in ]0, 1[$, $\gamma(w)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space. Hence, for $w \in ]0, 1[$, the univariate polynomial $f_s(\gamma(w), x_n)$ has two real roots (let us denote them by $\rho_1(\gamma(w))$ and $\rho_2(\gamma(w)))$ and the locus where it is positive meets the interval $]\xi_i(\gamma(w)), \xi_{i+1}(\gamma(w))[$. Assume by contradiction that $\gamma(1)$ does not lie in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space.
Then, by continuity of $\gamma, \xi_i, \xi_{i+1}$ and the two roots $\rho_1$ and $\rho_2$ of $f_s(\gamma(\cdot), x_n)$ we deduce that one of the two roots $\rho_1(\gamma(1))$, or $\rho_2(\gamma(1))$ coincides with either $\xi_i(\gamma(1))$ or $\xi_{i+1}(\gamma(1))$. By definition of the $\xi_i$'s, this implies that the resultant polynomial resultant($f_r, f_s, x_n$) (for $t+1 \le r \le s-1$) vanishes at $\gamma(1)$. This is a contradiction since this resultant polynomial cannot vanish over $\gamma([0, 1])$ (by definition of $S''$). Hence we deduce that $\gamma(1)$ lies in the projection of $S$ on the $(x_1, \dots, x_{n-1})$ coordinate space as requested. $\square$

## 5. Conclusion

System (S) has a specific structure that is rather common in models of mathematical biology. As is clear by inspecting the Jacobian matrix (2.9), the system has a negative feedback loop. Local and global stability properties of such systems have been widely studied (see, e.g., [1, 16]), taking advantage of monotone systems properties [25]. A raw picture in dimension three is that such a system either oscillates or is globally stable around the unique positive equilibrium ([16]). We tried to apply these techniques to our system, but up to now did not manage to conclude. So we believe that the local stability approach (instead of a global analysis) that we have developed in this paper is a meaningful step towards a more general result.

The relevance of the study of systems describing bacterial growth in the context of continuous bioreactors is a subject of research with a long history [18, 26] and a wide range of applications. In particular, the biosynthesis of chemical compounds with the aid of bacterial cells is central in numerous applications such as the production of insulin, antitumor agents, antibiotics and insecticides, among others. Numerous questions regarding the study of the dynamical behavior of such systems remain unanswered. In this context, this paper represents an alternative multi-disciplinary approach to the standard methods existing in the research community of mathematical biology.

## References

[1] David Angeli and Eduardo D Sontag. Monotone control systems. *IEEE Trans. Autom. Control*, 48(10):1684-1698, 2003.

[2] Emil Artin. Über die Zerlegung definiter Funktionen in Quadrate. *Abh. Math. Semin. Univ. Hamb.*, 5(1):100-115, 1927.

[3] Saugata Basu, Richard Pollack, and Marie-Françoise Roy. A new algorithm to find a point in every cell defined by a family of polynomials. In *Quantifier elimination and cylindrical algebraic decomposition*, Texts and Monographs in Symbolic Computation. Springer, 1998.

[4] Saugata Basu, Richard Pollack, and Marie-Françoise Roy. *Algorithms in Real Algebraic Geometry*, volume 10 of *Algorithms and Computation in Mathematics*. Springer, 2006.

[5] George E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decomposition. In *ATFL 2nd GI Conf. Kaiserslautern*, volume 33 of *Lecture Notes in Computer Science*, pages 134-183, 1975.

[6] Charles L. Cooney and Konstantin B. Konstantinov. White paper on continuous bioprocessing. In *International Symposium on Continuous Manufacturing of Pharmaceuticals*, 2014.

[7] Edmund X. DeJesus and Charles Kaufman. Routh-Hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinary differential equations. *Phys. Rev. A*, 35(12):5288, 1987.

[8] Attila Gábor and Julio R. Banga. Robust and efficient parameter estimation in dynamic models of biological systems. *BMC Syst. Biol.*, 9(1):1-25, 2015.

[9] Nils Giordano, Francis Mairet, Jean-Luc Gouzé, Johannes Geiselmann, and Hidde de Jong. Dynamical Allocation of Cellular Resources as an Optimal Control Problem: Novel Insights into Microbial Growth Strategies. *PLoS Comput. Biol.*, 2016.
[10] Antoine Giraud, Miroslav Radman, Ivan Matic, and François Taddei. The rise and fall of mutator bacteria. *Curr. Opin. Microbiol.*, 4(5):582–585, 2001.

[11] Dimitri Grigoriev and Nikolai Vorobjov. Solving systems of polynomial inequalities in subexponential time. *J. Symb. Comput.*, 5(1-2):37–64, 1988.

[12] Jérôme Izard, Cindy D. C. Gomez Balderas, Delphine Ropers, Stephan Lacour, Xiaohu Song, Yifan Yang, Ariel B. Lindner, Johannes Geiselmann, and Hidde de Jong. A synthetic growth switch based on controlled expression of RNA polymerase. *Mol. Syst. Biol.*, 11(11):840, 2015.

[13] Jean-Louis Krivine. Anneaux préordonnés. *J. Anal. Math.*, 12:307–326, 1964.

[14] Jean-Bernard Lasserre. Global Optimization with Polynomials and the Problem of Moments. *SIAM J. Optim.*, 11(3):796–817, 2001.

[15] Bingtuan Li. Global asymptotic behavior of the chemostat: general response functions and different removal rates. *SIAM J. Appl. Math.*, 59(2):411–422, 1998.

[16] John Mallet-Paret and Hal L. Smith. The Poincaré-Bendixson theorem for monotone cyclic feedback systems. *J. Dyn. Differ. Equations*, 2(4):367–421, 1990.

[17] Andreas Milias-Argeitis, Marc Rullan, Stephanie K. Aoki, Peter Buchmann, and Mustafa Khammash. Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth. *Nat. Commun.*, 7(1):1–11, 2016.

[18] Aaron Novick and Leo Szilard. Description of the chemostat. *Science*, 112(2920):715–716, 1950.

[19] Bernhard O. Palsson and Edwin N. Lightfoot. Mathematical modelling of dynamics and control in metabolic networks. I. On Michaelis-Menten kinetics. *J. Theor. Biol.*, 111(2):273–302, 1984.

[20] Pablo A. Parrilo. *Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization*. PhD thesis, California Institute of Technology, Pasadena, California, 2000.

[21] Helfried Peyrl and Pablo A. Parrilo. Computing sum of squares decompositions with rational coefficients. *Theor. Comput. Sci.*, 409(2):269–281, 2008.

[22] Mohab Safey El Din. RAGlib – A library for real solving polynomial systems of equations and inequalities. https://www-polsys.lip6.fr/~safey/RAGlib/.

[23] Mohab Safey El Din and Éric Schost. Polar varieties and computation of one point in each connected component of a smooth real algebraic set. In *ISSAC'03*, pages 224–231. ACM Press, 2003.

[24] Matthew Scott, Stefan Klumpp, Eduard M. Mateescu, and Terence Hwa. Emergence of robust growth laws from optimal regulation of ribosome synthesis. *Mol. Syst. Biol.*, 10(8):747, 2014.

[25] Hal L. Smith. *Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems*, volume 41 of *Mathematical Surveys and Monographs*. American Mathematical Society, 2008.

[26] Hal L. Smith and Paul Waltman. *The theory of the chemostat: dynamics of microbial competition*, volume 13 of *Cambridge Studies in Mathematical Biology*. Cambridge University Press, 1995.

[27] Gilbert Stengle. A nullstellensatz and a positivstellensatz in semialgebraic geometry. *Math. Ann.*, 207:87–97, 1974.

[28] Adam Strzeboński. Solving systems of strict polynomial inequalities. *J. Symb. Comput.*, 29(3):471–480, 2000.

[29] Horst R. Thieme. Convergence results and a Poincaré-Bendixson trichotomy for asymptotically autonomous differential equations. *J. Math. Biol.*, 30(7):755–763, 1992.

[30] Darren J. Wilkinson. Bayesian methods in bioinformatics and computational systems biology. *Brief. Bioinfor.*, 8(2):109–116, 2007.

[31] Agustín Gabriel Yabo. *Optimal resource allocation in bacterial growth: theoretical study and applications to metabolite production*. PhD thesis, Université Côte d'Azur, Nice, France, 2021.
[32] Agustín Gabriel Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Optimal bacterial resource allocation: metabolite production in continuous bioreactors. *Math. Biosci. Eng.*, 17(6):7074–7100, 2020.

[33] Agustín Gabriel Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Hierarchical MPC applied to bacterial resource allocation and metabolite synthesis. In *2021 IEEE 60th Conference on Decision and Control (CDC)*. IEEE, 2021.

[34] Agustín Gabriel Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Optimal allocation of bacterial resources in fed-batch reactors. In *2022 European Control Conference (ECC)*. IEEE, 2022.

[35] Agustín Gabriel Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Optimal bacterial resource allocation strategies in batch processing. preprint, https://hal.inria.fr/hal-03710681, 2022.

[36] Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé, Hidde de Jong, and Francis Mairet. Dynamical analysis and optimization of a generalized resource allocation model of microbial growth. *SIAM J. Appl. Dyn. Syst.*, 21(1):137–165, 2022.

[37] Agustín Gabriel Yabo and Jean-Luc Gouzé. Optimizing bacterial resource allocation: metabolite production in continuous bioreactors. *IFAC-PapersOnLine*, 53(2):16753–16758, 2020.

[38] Ivan Yegorov, Francis Mairet, Hidde de Jong, and Jean-Luc Gouzé. Optimal control of bacterial growth for the maximization of metabolite production. *J. Math. Biol.*, 2018.