Ratiometric Control of Two Microbial Populations via a Dual Chamber Bioreactor

Maintaining stable coexistence in microbial consortia, particularly when one species grows faster than the other (i.e., the species are non-complementary), poses significant challenges. We introduce a novel control architecture that employs two bioreactors. In this system, the slower-growing species is cultivated separately before being introduced into the main mixing chamber. We analyze the open-loop dynamics of this setup and propose a switching feedback mechanism that controls the dilution rates to ensure robust regulation of population density and composition within the microbial consortium. Validated in silico using parameters identified from real experiments, our approach demonstrates effective and robust maintenance of microbial balance across various strains without requiring genetic modifications.

each population while increasing the overall efficiency of the consortium [2].Over the past decade the translation of these advantages to synthetically engineered biological systems has attracted increasing interest [3] as it is an effective solution for the efficient production of complex chemical compounds [4], [5], [6].
In addition to the design of the phenotype of each population, the realization of functional microbial consortia requires stable, long-term coexistence between the different populations therein.This is a challenging goal to achieve as each population, engineered with different gene regulatory networks, is burdened by a different metabolic load, resulting in heterogeneous growth rates.Left uncontrolled, this imbalance will ultimately lead to the extinction of one or more populations, disrupting the proper functioning of the consortium.Thus, the key issue of regulating the microbial community's composition, known as composition, or ratiometric, control, is essential for creating consortia capable of reliably performing some desired function [7].
A promising scalable approach to co-culturing multiple bacterial populations involves using bioreactors to control both the density and composition of the consortium.Various strategies for managing coexistence within a single bioreactor have been explored.In [8] a deep reinforcement learning algorithm was employed to adjust the growth media composition, ensuring coexistence of two microbial populations dependent on different substrates.Similarly, studies in [9], [10] showed that modifying the dilution rate can support robust coexistence of two complementary strains.These approaches presume that growth conditions can be selectively altered to favor one population over another.However, it is not always feasible to create such conditions, as the growth properties, determined by the cellular resources consumed up to a given time, may remain fixed.In co-cultures, this results in the population that consumes more resources growing slower and potentially facing extinction.In a single bioreactor, this is impossible or cumbersome to avoid even when using a sophisticated feedback control action, because of the competitive exclusion principle [11], which from a dynamical systems standpoint translates into the system being uncontrollable.Therefore, a more general methodology is needed to guarantee long-term coexistence of different microbial species in a bioreactor.
In this letter, we propose an experimental platform endowed with a control strategy for robust regulation of both the density and the composition of a microbial consortium made of two non-complementary strains (one population persistently growing faster than the other).Specifically, we propose to use two communicating bioreactors (turbidostats) (see Fig. 1); one hosting a mixture of both populations, while the other being used as a reservoir where the slower population is grown alone.By using as control inputs the dilution rate D 1 (inlet volumetric flow/volume) associated with fresh media in the first reactor in combination with the dilution rate D 2 , at which the slower population from the reservoir (kept at a reference density therein by means of the dilution rate D 0 ) is added to the mixing chamber, our strategy can guarantee coexistence of the two populations at a desired ratio.
In what follows, we introduce in Section II a mathematical model that describes the growth dynamics of cells within each reactor and analyze its asymptotic behavior.Section III outlines the statement of the control problem, which is addressed in Section IV using an open-loop strategy.A closedloop control law is detailed in Section V.Both strategies are validated in silico in Section VI using a model parametrized with real experimental data from two Chi.Bio turbidostats [12].The proposed approach is versatile, suitable for a broad range of non-complementary strains growing on the same substrate, and does not require genetic interventions.

II. MATHEMATICAL MODEL
We consider two microbial species growing as a continuous culture in a chemostat (labeled as chamber 1 in Fig. 1), with a reservoir (chamber 2 in Fig. 1) in which one species grows separately.The two species are assumed to be independent, that is, they grow on the same substrate but they do not directly condition the growth rate of each other.The mathematical model of the former chemostat with two species and two inputs can be obtained by applying straightforward balance of mass law as in [13], yielding where the variables x 1 , x 2 ∈ R ≥0 and s 1 , s 2 ∈ [0, s in ] denote the concentrations of the biomass of species 1 and 2 and of the substrate in chamber 1 and 2, respectively.The mathematical model of the reservoir (chamber 2) can be written as where x 0 2 ∈ R ≥0 denotes the concentration of the biomass of species 2 in chamber 2.Moreover, μ i (•) is the growth rates of species i (defined below), Y i is the yield coefficient, assumed without loss of generality to be unitary for both species, and s in is the concentration of the substrate in the inlet flows.The control inputs D i are the dilution rates defined as the ratio between the inlet flow rate and the culture volume.Specifically, the control inputs D 1 (t) : R ≥0 → [D 1,min , D 1,max ], with D 1,min ≥ 0, and D 0 (t) : R ≥0 → [D 0,min , D 1,max ], with D 0,min > 0, dilute the concentrations in chambers 1 and 2, respectively, by injecting fresh substrate at concentration s in , while the input D 2 (t) : R ≥0 → [0, D 2,max ] is the dilution rate due to the inlet flow of media from chamber 2 to chamber 1.The dilution rates are assumed to be the same for both species, that is, the culture is well-mixed and mortality and attachment of the bacteria are neglected [10].
The dynamics of systems ( 1)-( 2) is characterized by the growth rate functions μ i (•) : [0, s in ] → R ≥0 , which are typically assumed to be differentiable, strictly increasing, and such that μ i (0) = 0, for i = 1, 2, corresponding to the substrate not having any inhibitory effect at high concentrations [14].The most common analytical growth rate model, that we also consider here, is the so-called Monod law [14] where μ * i > 0 is the maximum growth rate of species i and k i > 0 is a Michaelis-Menten constant.When the two growth functions intersect more than once, it is possible to design some control input D 1 (t) (setting D 2 = 0) such that both species survive [9], [10], [15].On the contrary, when it holds that the growth functions intersect only at 0 and the two species are said to be non-complementary [10], i.e., one species always outgrowing the other.Therefore, in this condition, from the Competitive Exclusion Principle (CEP) [16], it follows that when inflow from the additional reservoir is not available, i.e., D 2 = 0, no control input D 1 (t) exists such that more than one strain survives at steady state.However, in our architecture (Fig. 1) coexistence is still made possible by a controlled injection into the chemostat (chamber 1) of some extra biomass of the slower species taken from the reservoir (chamber 2) by means of the additional control input D 2 (t).

A. Reduced System
Under the assumption that the sum of the control inputs D 0 (t) and D 2 (t) is persistently exciting, that is, such that Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
are attracted to the invariant subspace S 0 := {(x 0 2 , s 2 ) : x 0 2 + s 2 = s in }.This can be demonstrated using invariant set theory, as illustrated in [17], similarly to the approach used in [15, Proposition 1] (see also [18] for further details).Thus, if the initial condition x 0 2 (0) belongs to S 0 , the dynamics of (2) on S 0 can be reduced to Likewise, under the assumption that the sum of D 1 (t) and D 2 (t) is persistently exciting, it can be proved using similar arguments that the solutions of system (1) are attracted to the invariant subspace S := {(x 1 , x 2 , s 1 ) : Thus, solutions rooted therein can be obtained from the reduced order model

B. Stability Analysis of the Reservoir Model
Depending on the value of (D 0 + D 2 ), system (5) can have one or two equilibria.One equilibrium is always at 0, while the other is x0 The point x0 2 is admissible and stable for 0 < (D 0 + D 2 ) < μ 2 (s in ).It becomes non-admissible, exchanging its stability with the point at the origin, when (D 0 +D 2 ) = μ 2 (s in ), through a transcritical bifurcation.

C. Stability Analysis of the Chemostat Model
As shown below, system (6) has a richer dynamics depending on the values of D 1 and D 2 . 1) 12 < μ 1 (s in ), there exists a stable equilibrium point in S r such that both species coexist.This point is located at the intersection between the x 1 -nullcline (ẋ 1 = 0) and the x 2nullcline (ẋ 2 = 0), obtained by solving the system of equations Moreover, there exists another equilibrium point, that is unstable, lying on the x 2 -axis (x 1 = 0) at the intersection with (9).When D 1 + D 2 is above some critical value μ c 12 , the coexistence point collides with the unstable point on the x 2 -axis, exchanging stability via a transcritical bifurcation, and ceases to be admissible.At this point only species 2 survives, because the total dilution rate in the bioreactor is too high for species 1 to survive.However, note that the complete extinction of species 1 is not possible in practice, as experimentally some cells will not be flushed out.Due to the nonlinearity of ( 9), the analytical expression of μ c 12 is cumbersome and is omitted here for the sake of brevity.Finally, the orange dashed lines indicate the safety bounds for each species (i.e., x 1,min , x 2,min ).b) Graphical representation of the closedloop control strategy.S r is divided in four regions (R 1 , R 2 , R 3 , R 4 ) by the switching surfaces 1 and 2 (blue lines).The control actions D 1 and D 2 steer the trajectory of the system (blue arrows) onto x d (red star).c) Block diagram of the closed-loop control architecture.The switching controller SC regulates the concentrations of the biomasses x 1 and x 2 in chamber 1 to the desired levels x 1,d and x 2,d , respectively.The control input D 0 , regulating x 0 2 in chamber 2 to the desired level x0 2 , is computed by a PI controller with a feed-forward disturbance compensation.
2) D 1 ≈ 0, D 2 > 0: The dynamics is similar to the previous case, that is, there exists a stable equilibrium point associated with coexistence whose position depends only on D 2 .However, when D 2 > μ 1 (s in ), that is, the dilution rate from the reservoir is higher than the maximum growth rate of species 1, the faster species is flushed out from the bioreactor (x 1 ≈ 0) and as D 2 increases the concentration of species 2 reaches asymptotically the same value as in the reservoir (i.e., x2 ≈ x 0 2 ). 3) D 1 > 0, D 2 = 0: Model (6) becomes the same model described in [9], and as D 1 varies from 0 to D 1,max > μ 1 (s in ) we retrieve the same dynamics depicted in the phase portraits shown in Fig. 2 from [9] (namely cases I, IV, V, VI).
4) D 1 = 0, D 2 = 0: All solutions converge to the attractive equilibrium set E 0 := {(x 1 , x 2 ) ∈ S r :x 2 = −x 1 + s in } corresponding to the biomass being in starvation and therefore it is not an admissible working condition in continuous culture (see case I in [9, Fig. 2]).

III. CONTROL PROBLEM
Given the two-chambers bioreactor architecture depicted in Fig. 1, with two non-complementary microbial species, whose mathematical model and dynamics have been described in Section II, our control objective is to control the dilution rates such that at steady state the two species coexist and their concentrations are regulated at a desired ratio.That is, we want to solve the following problem.
Control problem: Design some control inputs D 0 (t), D 1 (t), and D 2 (t) for system (5)-( 6), under condition (4), such that: 1) coexistence of the two species is guaranteed for all time, that is, x i (t) ≥ x i,min , ∀t > 0, i = 1, 2, where x i,min > 0 is some safety value to avoid extinction of species i; 2) the total biomass is regulated to a desired value to guarantee efficient utilization of the resources, that is, where OD d is the desired optical density in chamber 1, assumed for the sake of simplicity to be an exact proxy of the total biomass therein; 3) the ratio x 2 /x 1 between the concentration of the two species in chamber 1 is robustly regulated at steady state to a desired value r d , that is, Notice that the above control problem corresponds to requiring that all solutions of system ( 5)-( 6), starting from any initial conditions in S 0 ∪ S r , converge to the point of intersection between the two curves defined in ( 10) and ( 11) (Fig. 2.a), that is, to the point In the following, we assume that the concentrations x 1 , x 2 and x 0 2 are either directly or indirectly measurable, for example via fluorescent reporters or state observers.

IV. OPEN-LOOP CONTROL
A simple open-loop controller that solves the problem described above can be designed by exploiting the dynamics of the system described in Section II, that is, by choosing constant values D0 , D1 and D2 so as to make the nullclines ( 8) and ( 9) intersect at the desired point (12).Specifically, setting D12 := D1 + D2 by using ( 8) and (10) we obtain Moreover, by using ( 9) and ( 12) we get where the steady-state value x0 2 in the reservoir can be reached by setting the control input D 0 to D0 = μ 2 (s in − x0 2 ) − D2 , obtained from setting the first term in (5) to zero.By definition, we also have that D1 = D12 − D2 .Further details about the derivation of ( 13)-( 14) are reported in [18].Although this feed-forward control action is simple and solves the control problem in nominal conditions, as we will see later in Section VI, it neither guarantees robustness to perturbation nor fast convergence.The coexistence requirement, as defined in Section III, is fulfilled by the open-loop controller provided that the initial condition (x 1 (0), x 2 (0)) is chosen beneath the intersection between the nullcline (8) and the constraint x 1 ≥ x 1,min .

V. CLOSED-LOOP CONTROL
In biological applications robustness is as important as precision in the regulation and open-loop strategies cannot often reject external perturbations and uncertainties in the model parameters.Therefore, we propose next a closed-loop switching controller that can globally and robustly guarantee regulation of the desired density and composition.
The controller is designed by defining two switching surfaces, 1 := {(x 1 , x 2 ) ∈ S r : x 1 = x 1,d } and 2 := {(x 1 , x 2 ) ∈ S r : x 2 = x 2,d }, which divide the domain S r into four regions (see Fig. 2.b), R 1 where x 1 < x 1,d and x 2 > x 2,d , R 2 where x 1 > x 1,d and x 2 > x 2,d , R 3 where x 1 < x 1,d and x 2 < x 2,d , and R 4 where x 1 > x 1,d and x 2 < x 2,d .We denote by ij the boundary between regions R i and R j .
Theorem 1: Choosing x0 2 > x 2,d as set-point for the reservoir, and the switching control inputs as so that the bioreactor is in flush-out mode, • D 1 = D 2 = 0 when x ∈ R 3 , so that way the biomass is left to grow freely, (12), is attractive for all solutions of system (6) in S r .
Proof: This can be proved by showing that the intersection of the two switching surfaces, 1 ∩ 2 is attracting.Let n 1 = [1 0] T , n 2 = [0 1] T be the normal vectors to 1 and 2 , respectively, and f 1 (x), f 2 (x), f 3 (x), and f 4 (x) the vector fields of system (6) in each of the four regions determined by our choice of the switching control inputs, that is, The boundary 13 is a stable sliding region because n T 2 f 1 (x) < 0 and n T 2 f 3 (x) > 0, since (μ 2 (s 1 ) − D + 1 ) < (μ 2 (s 1 ) − μ 2 (s in )) and μ 2 (s 1 ) > 0, respectively.Moreover, solutions sliding on 13 are attracted towards d .This is because the intersection between the nullcline (8) and 2 lies above x d , and μ s (s 1 ) > 0, respectively.Similarly, 34 is a stable sliding region because n T . Moreover, solutions on 34 are attracted towards Finally, the boundary 12 has always a crossing region, in which solutions from R 2 cross the boundary because n T 1 f 2 (x) < 0 while n T 1 f 1 (x) changes sign at the intersection of 12 with the nullcline (8).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
However, all crossing solutions will then converge onto 13 since ẋ2 < 0, and therefore they will move towards x d , as proved above.This concludes the proof.
Although as stated earlier, complete extinction of species 1 is not possible, some trajectories starting within R 1 ⊂ S r may still violate the coexistence bounds defined in Section III.The set of initial conditions (x 1 (0), x 2 (0)) ∈ R 1 yielding trajectories that temporarily violate those bounds can be removed by choosing D + 1 = μ 2 (s in )+ , with > 0 arbitrarily small, so that the nullcline (8) intersects at the uppermost point the line x 1 = x 1,min .Therefore, initial conditions must be chosen so that The controller implemented in the reservoir (chamber 2) to regulate x 0 2 is a PI controller with a disturbance compensation.Specifically, the dilution rate D 0 in ( 5) is chosen as D 0 = D PI − D 2 , where D PI is an action computed using the same proportional-integral controller presented in [12] for the regulation of the turbidity.This controller ensures that the error e PI (t) = x0 2 − x 0 2 (t) converges to zero, so that the density of species 2 in the reservoir x 0 2 is regulated to the desired value x0 2 .

VI. IN SILICO VALIDATION
We validated the designed controllers in silico to assess their performance, stability, and robustness.Specifically, we simulated the experimental platform shown in Fig. 1 implemented using two Chi.Bio reactors [12] and included the technological constraints imposed by the experimental platform.Each reactor measures the density of the biomasses at every minute by measuring the turbidity of the liquid (i.e., the optical density) and updates the duty-cycle of the peristaltic pumps that dilute the cultures with fresh media.The dilution rates are saturated to 0.065 min −1 for safety reasons.In all control experiments we set the desired reference signals to OD d = 0.7 and r d = 1.The parameters of the simulated set-up were parametrized from real data as explained below.

A. Model Parametrization
For the identification of the parameters in system (6), we conducted open-loop characterization experiments using the commercially available reactor Chi.Bio [12].Specifically, we characterized the growth dynamics of two E. coli strains embedding two different implementations of a genetic toggleswitch presented in [19] and [20], respectively.In each experiment, using a single Chi.Bio reactor, we grew a single strain in LB media supplemented with either 50 μg/mL kanamycin or 25μg/mL Chloramphenicol, and 1 mM isopropyl β-D-1-thiogalactopyranoside (IPTG), at 37 • C, shaking at half the maximum speed available on the platform.During these experiments, the culture was diluted using the maximum available dilution rate until the optical density dropped below 0.3.Then, cells underwent a recovery phase where the dilution rate was kept close to zero for 60 minutes.Finally, the cells' growth was excited changing the dilution rates every 15 to 30 minutes for 1 hour (Fig. 3).
Using Genetic Algorithm (GA) optimization (as implemented in [21]), we estimated the growth parameters for each strain to be: μ * 1 = 0.026, μ * 2 = 0.013, k 1 = 20.32,k 2 = 31.25,and s in = 88.75.The ranges for the μ * i were set according to typical values in literature [20], while the other parameters' ranges were obtained heuristically.The intervals were set as follows: [0.001, 0.05] for μ * 1 and μ * 2 , [0.01, 100] for k 1 and k 2 , [1,100] for s in to start the GA optimization.The total number of iterations was 200, and the initial population of the GA algorithm was set to be 150 individuals.In Fig. 3 the model's predictions are depicted alongside the experimental data, demonstrating the model ability to capture quantitatively the growth dynamics of the microbial species under varying conditions.The data used for validation were different from those employed for the parameters identification.The quality of the predictions was quantified using the Root Mean Squared Error (RMSE), which proved to be low for the estimation of both populations' trajectories (i.e., RMSE = 3.18 for species 1 and RMSE = 1.73 for species 2) and comparable with the RMSE of the identification data set (i.e., RMSE = 1.86 for species 1 and RMSE = 1.09 for species 2).

B. Control Validation
Firstly, we simulated model ( 1)-( 2) under the action of the open-loop controller, described in Section IV (with dilution rates equal to D1 = 0.016, D2 = 0.0051, D0 = 0.0044), including also the technological constraints described above, observing that the desired density and configuration were reached asymptotically (see Fig. 4.a).However, the convergence time was excessively long, namely 97 • 10 3 min and 116 • 10 3 min to reach the desired optical density and ratio, respectively, making this strategy unfeasible in practice for in vivo implementation.
Instead, as shown in Fig. 4.b the switching control strategy developed in Section V was able to guarantee a much faster settling time (at 5%), namely, 40.5 min and 52.6 min for the optical density and the ratio between the two species, respectively, at the cost of a relatively small residual error at steady state.We quantified the relative residual errors by evaluating the average deviation from the desired values in the last 60 minutes of the experiments, that is, ēr = |(r − r d )/r d |, where r = 1/60 • t f t f −60 (x 2 (τ )/x 1 (τ ))dτ , with t f = 150 min, and similarly for ēOD , obtaining 0.026 and 0.0004, respectively.
Finally, we assessed the robustness of the proposed switched controller against variations in the parameters of the growth functions μ i (s), i = 1, 2. Specifically, the values of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.).The results of this analysis, reported in Fig. 5, highlight that, even when a consistent mismatch is present between the identified parameters and their true value (CV = 30%), the average steady-state residual errors remain below 0.05 for the ratio and 0.01 for the optical density.

VII. CONCLUSION
We introduced a dual-reactor architecture to steer the overall density and composition of a two-strain consortium.We analyzed its open-loop dynamics and then formulated two tailored control strategies.These strategies were validated in silico with model parameters identified experimentally.Our approach enables long-term coexistence within a consortium of non-complementary cell populations without the integration of additional synthetic circuits into the cells.This capability is crucial for supporting the simultaneous growth of different populations in a consortium, potentially across several different species such as bacteria and yeasts.We are currently implementing a prototype of the proposed architecture shown in Fig. 1, featuring two Chi.Bio reactors.Ongoing in vivo experiments aim to further validate the effectiveness of our control strategies, with results to be detailed elsewere.

Fig. 2 .
Fig.2.a) Graphical representation of the control objective.The green and blue lines represent the sets where x 1 + x 2 = OD d and x 2 = r d x 1 , respectively.The black dashed line delimits the invariant set S r .Finally, the orange dashed lines indicate the safety bounds for each species (i.e., x 1,min , x 2,min ).b) Graphical representation of the closedloop control strategy.S r is divided in four regions (R 1 , R 2 , R 3 , R 4 ) by the switching surfaces 1 and 2 (blue lines).The control actions D 1 and D 2 steer the trajectory of the system (blue arrows) onto x d (red star).c) Block diagram of the closed-loop control architecture.The switching controller SC regulates the concentrations of the biomasses x 1 and x 2 in chamber 1 to the desired levels x 1,d and x 2,d , respectively.The control input D 0 , regulating x 0 2 in chamber 2 to the desired level x0 2 , is computed by a PI controller with a feed-forward disturbance compensation.

Fig. 3 .
Fig. 3. Comparison between model predictions and the experimental data obtained from open-loop characterization experiments conducted using the Chi.Bio reactor.Two different E. coli strains implementing genetic toggle-switches from [19] (panel a) and [20] (panel b) were used in the experiments.The top panels depict the temporal evolution of the biomass in the in vivo experiment (orange) and in the numerical simulation (gray).The bottom panels show the inputs applied to the bioreactors.The experimental data are collected using a single reactor (i.e., setting D 0 = D 2 = 0).

Fig. 4 .
Fig. 4. In silico validation of the open-loop (a) and closed-loop (b) control strategies.The top panels show the evolution in time of x 2 /x 1 and x 1 + x 2 , depicted as gray solid lines, with their reference values (dashed green lines) and, for those in (b), the average trajectory over a time window of 30 minutes (blue solid line).The reference values are set as OD d = 0.7 and r d = 1.The bottom panels show the inputs applied to the first chamber.The control parameters of the closed-loop controller are set to D + 1 = 0.013, D − 1 = D − 2 = 0.065.

Fig. 5 .
Fig. 5. Robustness analysis of the closed-loop controller.Mean (bar height) and standard deviation (whisker) of ēr (panel a) and ēOD (panel b), for increasing value of the coefficient of variation CV .Each bar is obtained by averaging 100 simulations each computed with different values of parameters (μ * 1 , μ * 2 , k 1 , k 2 ).parameters(μ * 1 , μ * 2 , k 1 , k 2 )were randomly drawn from normal distributions, each centered on their nominal value, say η, with a standard deviation equal to σ = CV • η, where CV is the coefficient of variation.For each value of CV, taken in the interval [0.05, 0.3], we evaluated the mean and standard deviation of the residual errors ēOD and ēr , over 100 simulations, each computed with different values of the parameters(μ * 1 , μ * 2 , k 1 , k 2 ).The results of this analysis, reported in Fig.5, highlight that, even when a consistent mismatch is present between the identified parameters and their true value (CV = 30%), the average steady-state residual errors remain below 0.05 for the ratio and 0.01 for the optical density.
Schematic representation of the two-reactor architecture.The slower species (species 2, in pink) is grown separately in chamber 2 and added at rate D 2 in chamber 1, where it is mixed with species 1 (in green).Species 2 is regulated at the desired concentration x0 2 in chamber 2 by means of the dilution rate D 0 .The dilution rate D 1 is used in combination with D 2 to keep the two species at the desired concentrations in chamber 1.
c 2024 The Authors.This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/Authorizedlicensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 1.