MathematicS In Action

We describe microbial growth and production of value-added chemical compounds in a continuous biore-actor 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.


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 ocurring 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].
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.

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 semialgebraic 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 critical 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 inequalities 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.

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 Decomposition 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 analysis 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.

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.
Following [32,37], the system of differential equations that describes the evolution of the state variables is given by 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: 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 θ (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)V +x and s+(p+r/u * )V +x, we can see that they both tend asymptotically to 1 when t → ∞1 .Thus, the ω-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 ω-limit set of any solution of system (S original ) lies in the hyperplanes The latter implies that, for every trajectory, r converges to u * asymptotically, and x can be expressed in terms of the remaining states as (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 * ))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 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].

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 ṡ = ṗ = V = 0, we obtain the steady-state values in terms of the system parameters θ.

Proposition 2.2. The equilibrium of interest is
where 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: 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 C of the Jacobian matrix and check whether their real part is negative.For this case, the Jacobian matrix is defined as We define the functions Numerical trajectories of (S) obtained for fixed initial conditions and random set of admissible parameters θ (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.
that satisfy ϕ i (θ) > 0 for i ∈ {1, . . ., 4}.Computing the characteristic polynomial evaluated at the equilibrium point given by (2.6) yields (2.12) The Routh-Hurwitz criterion for degree-three polynomials states that the roots of a degree-three polynomial λ 3 + α 2 λ2 + α 1 λ + α 0 belong to the open left complex plane if and only if From (2.12), since every function ϕ i is positive, the stability criterion reads In order to prove the latter, we can prove that no state-parameter combination satisfies the negation of the condition, The stability of the steady state of interest is given by the validity of the latter inequality, which depends on the value of θ.
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
The next section is devoted to the proof of this result.The Maple worksheet used to prove the result is available online 2 as a companion notebook to the current paper.It makes an intensive use of the RAGlib package.

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 (γ, u * , p * , K 1 , K) while γ 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 The resulting system of inequalities becomes where C is a polynomial of total degree 13 involving the above 7 variables, whose number of nonzero 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. 3The 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 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 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 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 i f and x j g for 0 ≤ i ≤ q − 1 and 0 ≤ j ≤ 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 ∂f ∂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 Q[x 1 , . . ., x n ], one can see them as polynomials in (as does the leading coefficient of f w.r.t.x n which we denote by 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 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 , . . ., f s be polynomials in Q[x 1 , . . ., x n ] and let S ⊂ R n be the semi-algebraic set defined by S is open for the Euclidean topology, its projection on the (x 1 , . . ., x n−1 ) coordinate space is open too.Let S ′ be the semi-algebraic set defined by By [28], there exist connected components C 1 , . . ., C ℓ of S ′ such that C 1 ∪ • • • ∪ C ℓ is semialgebraic and dense in the projection of S on the (x 1 , . . ., x n−1 ) coordinate space.We call the OpenCAD, the procedure which, given the polynomial system f 1 > 0, . . ., 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 α 1 , . . ., α r be those sample points, solve the univariate system of polynomial inequalities f t+1 (α i , x n ) > 0, . . ., f s (α i , x n ) > 0, for each 1 ≤ i ≤ 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 ≤ i ≤ s − 1.We let C be a connected component of S ′ that is contained in the projection of S on the (x 1 , . . ., x n−1 ) coordinate space.In other words, for any α ∈ C ⊂ R n−1 , there exists ϑ ∈ R such that f i (α, ϑ) > 0 for t + 1 ≤ i ≤ s.By definition of S ′ , the discriminant of f s does not vanish over C. We then consider the two possible cases: and is then positive (because as explained above, there exists ϑ ∈ R such that f s (α, ϑ) > 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 α ∈ C, f s (α, x n ) has exactly two real roots with multiplicity one and there exists ϑ such that f s (α, ϑ) 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 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 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 α ∈ C which does not cancel lc(f s , x n ) lies in the projection of S on the (x 1 , . . ., x n−1 ) coordinate space.
Proof.Consider the semi-algebraic set defined by f 1 > 0, . . ., f s−1 > 0. Running the procedure OpenCAD on the input [f 1 , . . ., 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 ) ̸ = 0 and resultant 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) ξ 1 , . . ., ξ k from C to R such that the cylinder C × R is the disjoint union of cells defined as • either the graph of one of the maps ξ i for 1 ≤ i ≤ k, • or a band of some cylinder for 0 ≤ i ≤ k with ξ 0 = −∞ and ξ k+1 = +∞ by convention, over which the polynomials f 1 , . . ., f s−1 are sign invariant (because the system output by the call to the procedure OpenCAD on the input [f 1 , . . ., 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 α ∈ C, ξ 1 (α), . . ., ξ k (α) is the ordered sequence of roots of the polynomials f t+1 (α, x n ), . . ., f s−1 (α, x n ) (note that these are univariate).Hence, the semi-algebraic maps ξ i 's are defined as maps which send α ∈ C to some root of some f j (α, x n ) (those roots being maintained ordered when α ranges over C).
Assume now that C has a non-empty intersection with the projection of S on the (x 1 , . . ., 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 α be in this intersection such that lc(f s , x n ) does not vanish at α (such a point α exists because C is open for the Euclidean topology as well as the projection of S and the solution set to lc(f s , x n ) ̸ = 0).In other words, there exists ϑ such that f 1 (α, ϑ) > 0, . . ., f s (α, ϑ) > 0. Let B i be the band which contains (α, ϑ).Consider now α ′ ∈ C and assume that lc(f s , x n )(α ′ ) ̸ = 0. We need to prove that there exists ϑ ′ ∈ R such that (α ′ , ϑ ′ ) ∈ S. Note that, by the sign invariance property over the band B i of f 1 , . . ., f s−1 , we already know that there exists ϑ ′ such that f j (α ′ , ϑ ′ ) > 0 for 1 ≤ j ≤ s − 1.Hence, what is missing is a control on the sign of f s , precisely that there exists ϑ ′ ∈ R such that f s (α ′ , ϑ ′ ) > 0. The next lemma is a first step towards this.We postpone its proof to next section.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 , . . ., 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 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.

Proof of auxiliary lemmas
Proof of Lemma 3.5.Note that for any v ∈ [0, 1], and any ξ i (γ(v)) < ϑ v < ξ i+1 (γ(v)), the following holds: 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 ∈ [0, 1], f s (γ(v), x n ) has no root R (because we have assumed that lc(f s , x n ) does not vanish at γ(v)).Since, by construction, f s (γ(0), x n ) is positive over R, we deduce that for any v ∈ [0, 1], f s (γ(v), x n ) is positive over R (else, by continuity of γ, there would exist some v ′ such that f s (γ(v ′ ), x n ) has some real root while lc(f s , x n ) does not vanish at γ(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 ∈ [0, 1], there exists ϑ v ∈ B i such that f s (γ(v), ϑ v ) > 0. We start with v = 0. We already know by construction that there exists ξ i (γ(0)) < ϑ < ξ i+1 (γ(0)) such that f s (γ(0), ϑ) > 0. Assume by contradiction that there exists v ′ ∈ [0, 1] such that for any ξ By continuity of γ, ξ i and ξ 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 ξ i 's we deduce that f s (γ(v ′ min ), ξ i (γ(v ′ min ))) = 0 or f s (γ(v ′ min ), ξ i+1 (γ(v ′ min ))) = 0. Recall that there exists some t + 1 ≤ j ≤ s − 1 such that ξ i maps α ∈ C to some root of f j .In other words, the gcd of f s (γ(v ′ min ), x n ) and f j (γ(v ′ min ), x n ) has degree ≥ 1.This implies that for some j, the resultant polynomial resultant(f j , f s , x n ) vanishes at γ(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 ) ̸ = 0.This ends the proof of Lemma 3.5.□ Proof of Lemma 3.6.Take w ∈ ]0, 1[.Up to scaling [0, w] to [0, 1], one applies Lemma 3.5 to deduce that γ(w) lies in the projection of S on the (x 1 , . . ., x n−1 ) coordinate subspace.Since γ(w) lies in the projection of S on the (x 1 , . . ., x n−1 ) coordinate subspace, there exists Note that ϑ w can be chosen in the interval ]ξ i (γ(w)), ξ i+1 (γ(w))[ where ξ i is one of the "root" functions introduced in the previous section.By continuity of γ, the equality γ(1) = lim w→1 (γ(w)) holds.We are going to prove that γ(1) lies in the projection of S on the (x 1 , . . ., x n−1 ) coordinate subspace.By assumption, lc(f s , x n ) vanishes at γ (1).Note that, by definition of the discriminant of a quadratic polynomial, this implies that discriminant(f s , x n ) is non-negative at γ (1).Also, discriminant(f s , x n ) is sign invariant over γ([0, 1]) (since the inequation constraint discriminant(f s , x n ) ̸ = 0 is part of the polynomial system (3.3)).Hence, we deduce that discriminant(f s , x n ) is positive over γ([0, 1]).Recall that we established that, for all w ∈ ]0, 1[, γ(w) lies in the projection of S on the (x 1 , . . ., x n−1 ) coordinate space.Hence, for w ∈ ]0, 1[, the univariate polynomial f s (γ(w), x n ) has two real roots (let us denote them by ρ 1 (γ(w)) and ρ 2 (γ(w))) and the locus where it is positive meets the interval ]ξ i (γ(w)), ξ i+1 (γ(w))[.Assume by contradiction that γ(1) does not lie in the projection of S on the (x 1 , . . ., x n−1 ) coordinate space.

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.
Figure2.2.Numerical trajectories of (S) obtained for fixed initial conditions and random set of admissible parameters θ (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.