
Citation: Wenyuan Xie, Jason Wei Jun Low, Arunmozhiarasi Armugam, Kandiah Jeyaseelan, Yen Wah Tong. Regulation of Aquaporin Z osmotic permeability in ABA tri-block copolymer[J]. AIMS Biophysics, 2015, 2(3): 381-397. doi: 10.3934/biophy.2015.3.381
[1] | Shahab Shamshirband, Javad Hassannataj Joloudari, Sahar Khanjani Shirkharkolaie, Sanaz Mojrian, Fatemeh Rahmani, Seyedakbar Mostafavi, Zulkefli Mansor . Game theory and evolutionary optimization approaches applied to resource allocation problems in computing environments: A survey. Mathematical Biosciences and Engineering, 2021, 18(6): 9190-9232. doi: 10.3934/mbe.2021453 |
[2] | Liang Tian, Fengjun Shang, Chenquan Gan . Optimal control analysis of malware propagation in cloud environments. Mathematical Biosciences and Engineering, 2023, 20(8): 14502-14517. doi: 10.3934/mbe.2023649 |
[3] | K. Renee Fister, Jennifer Hughes Donnelly . Immunotherapy: An Optimal Control Theory Approach. Mathematical Biosciences and Engineering, 2005, 2(3): 499-510. doi: 10.3934/mbe.2005.2.499 |
[4] | Xiaoxiao Yan, Zhong Zhao, Yuanxian Hui, Jingen Yang . Dynamic analysis of a bacterial resistance model with impulsive state feedback control. Mathematical Biosciences and Engineering, 2023, 20(12): 20422-20436. doi: 10.3934/mbe.2023903 |
[5] | Hermann Mena, Lena-Maria Pfurtscheller, Jhoana P. Romero-Leiton . Random perturbations in a mathematical model of bacterial resistance: Analysis and optimal control. Mathematical Biosciences and Engineering, 2020, 17(5): 4477-4499. doi: 10.3934/mbe.2020247 |
[6] | Semu Mitiku Kassa . Three-level global resource allocation model for HIV control: A hierarchical decision system approach. Mathematical Biosciences and Engineering, 2018, 15(1): 255-273. doi: 10.3934/mbe.2018011 |
[7] | Yafei Li, Yuxi Liu . Multi-airport system flight slot optimization method based on absolute fairness. Mathematical Biosciences and Engineering, 2023, 20(10): 17919-17948. doi: 10.3934/mbe.2023797 |
[8] | Jing Jia, Yanfeng Zhao, Zhong Zhao, Bing Liu, Xinyu Song, Yuanxian Hui . Dynamics of a within-host drug resistance model with impulsive state feedback control. Mathematical Biosciences and Engineering, 2023, 20(2): 2219-2231. doi: 10.3934/mbe.2023103 |
[9] | Rinaldo M. Colombo, Mauro Garavello . Stability and optimization in structured population models on graphs. Mathematical Biosciences and Engineering, 2015, 12(2): 311-335. doi: 10.3934/mbe.2015.12.311 |
[10] | Yu Shen, Hecheng Li . A new differential evolution using a bilevel optimization model for solving generalized multi-point dynamic aggregation problems. Mathematical Biosciences and Engineering, 2023, 20(8): 13754-13776. doi: 10.3934/mbe.2023612 |
Microorganisms continuously have to contend with environmental changes in nature, and so they have evolved to accordingly adapt their physiology to cope with this unsteadiness. This is done by reorganizing the gene expression machinery, which is accomplished by dynamically allocating resources to different cellular functions. Microorganisms like bacteria exhibit great genetic variability, which is the main driver of natural selection, a phenomenon that depends on the continuous mutations in living populations. For this reason, they have achieved highly optimized resource management strategies. Such is the case of E. Coli: For specific environmental conditions, they seek to maximize their growth rate, which is a naturally evolved characteristic of this bacterium, and well known in the scientific community [1]. In this regard, optimality of bacterial organisms has been a central hypothesis in several research allocation studies [2].
Nevertheless, most of previous works in the literature consider the resource allocation problem in steady-state, without taking into account the changing conditions of natural environments [3]. This motivated a series of works with a dynamical-oriented approach to the problem [4,5], that aimed to shed light on how bacteria tune their allocation mechanisms when facing changes on the nutrient concentration of the medium. Using Optimal Control theory, they investigate how their internal pathways can be dynamically readjusted so as to maximize their growth rate, obtaining different strategies that can be related to well known natural regulation mechanisms (such as the ppGpp-mediated sensing of the pool of aminoacids). The approach is based on a widely used modeling technique in systems biology: The so-called coarse-grained self-replicator models, used in bacterial growth representations for their simplicity and their capacity to reproduce observed experimental behaviors [6]. Although this kind of single-cell models are somewhat limited when predicting complex phenomena, they can accurately account for bacterial culture growth laws under the right assumptions (e.g., homogeneity of the culture) [7].
From an industry point of view, a natural question triggered by these studies is: How can we divert the natural allocation strategies of bacteria to improve current biotechnological production schemes? Such is the case of the artificial synthesis of metabolites or proteins of interest. The synthesis of such compounds is highly relevant for its wide range of applications: Production of antitumor agents, insulin, antibiotics, immunosuppressive agents and insecticides, among others [8,9]. Motivated by the increasing understanding of the biosynthetic properties of certain microorganisms, research on this area can potentially lead to more efficient and sustainable production schemes. This is the matter addressed in recent work using a strain of Escherichia coli that includes an artificially engineered heterologous pathway for the production of a certain metabolite of interest [10,11,12]. In this approach, a control loop is developed through a bacterial growth switch that allows to externally modify the natural resource allocation decision [13]. The mechanism is implemented by re-engineering the transcriptional control of the expression of RNA polymerase, a key component of the gene expression machinery. This way, it is possible to optimize the productivity of the bioprocess by channeling resources into this new heterologous pathway.
At the same time, the synthesis of these metabolites draws resources from the native pathways used for producing biomass in bacteria, thus leading to an inherent compromise between these two objectives. One possible approach to this trade-off is to model it through different cost functions, thus obtaining multi-objective optimization problems. This is the case of [14], where the authors aim to maximize the production of a metabolite while minimizing the genetic burden caused by pathway expression. In contrast to this method, the work of [10] models the main trade-offs behind the process through a single decision parameter, which considerably reduces the complexity of the optimization problem.
In a similar vein, we address a classic production scheme: The Continuous Stirred-Tank Reactor (CSTR). While resource allocation in bacteria has been vastly studied in constant environments (e.g., under the assumption that there is always enough substrate in the medium), how this goes in continuous bioreactors is not trivial, since a feedback occurs from the physiology of the cell to the environmental conditions given by the interaction bacteria-medium [15]. Examples of self-replicator models in continuous bioreactors can be found in some recent works: In [16], authors use a coarse-grained self-replicator kinetic model of Saccharomyces cerevisiae's in a continuous bioreactor to investigate the trade-off between respiratory and fermentative metabolism, showing that optimal strategies are ‘pure’ metabolic strategies (e.g., either respiration of fermentation, but not respiro-fermentation). Likewise, it is rather classical in the continuous bioreactor scheme to maximize a certain performance measure (i.e., biomass production) in terms of the operational parameters related to the setup, such as dilution rate and/or concentration of the substrate inflow [17,18,19,20,21]. However, incorporating the aforementioned external control—that can disrupt the natural allocation process of the whole culture—provides an extra degree of freedom, which can, in turn, contribute to further improve the classical production scheme.
Based on the presented works, we show novel results addressing the problem of bacterial resource allocation in the CSTR framework. Our approach is based on a coarse-grained self-replicator dynamical model that accounts for the microbial culture growth inside the continuous bioreactor, and incorporates the external allocation control previously described. The novelty of the approach lies in the combination of the resource allocation control scheme, and the capacity to regulate the bacterial growth rate through the dilution rate of the continuous bioreactor. Further on, we study its asymptotic behavior using dynamical systems theory, and we provide conditions for the persistence of the bacterial population. The analysis is carried out by studying the local and global behavior of its limiting system, and relating its convergence to the original model through the theory of asymptotically autonomous systems [22]. Then, we pose the problem of maximizing the synthesis of the metabolite of interest during a fixed interval of time in terms of the resource allocation decision. The latter is expressed as an OCP (Optimal Control Problem) which is then explored using PMP (Pontryagin's Maximum Principle) [23]. We analyze the solution of the problem and propose a sub-optimal control strategy given by a constant allocation decision, which eventually takes the system to the optimal steady-state production regime. On this basis, we study and compare the two most significant steady-state production objectives in CSTRs: Biomass production and metabolite production. For this last purpose, and in addition to the allocation parameter, we control the constant volumetric flow of the bioreactor (or dilution rate), and we analyze the results through a numerical approach. The resulting two-dimensional optimization problem is defined in terms of Michaelis-Menten kinetics with the parameter values of [4], and taking into account the constraints for the existence of the equilibrium of interest.
In this section, we define the coarse-grained self-replicator model including the continuous bioreactor scheme and the allocation parameter. We consider a growing bacterial population in a CSTR bioreactor of constant volume $ \mathcal{V}_\mathrm{ext} $ [L]. The self-replicating system that models the culture is composed of the metabolic machinery $ M $ (transporters, enzymes...) and the gene expression machinery $ R $ (RNA polymerase, ribosomes...), both responsible of the cell growth. The validity of this single-cell model as a representation of a growing bacterial culture depends on a number of simplifying assumptions, one of them being that individual cells share the same macromolecular composition. For the production of the metabolite of interest $ X $, we consider the artificially engineered pathway introduced in [10] (Figure 1).
The model represents three chemical macroreactions,
$ \mathrm{S} \stackrel{\mathrm{V}_{\mathrm{M}}}{\longrightarrow} \mathrm{P} $ |
$ \mathrm{P} \stackrel{\mathrm{V}_{\mathrm{R}}}{\longrightarrow} \alpha \mathrm{R}+(1-\alpha) \mathrm{M} $ |
$ \mathrm{P} \stackrel{\mathrm{V}_{\mathrm{X}}}{\longrightarrow} \mathrm{X} $ |
The first reaction is catalyzed by $ M $ and describes the transformation of external substrate $ S $ into precursor metabolites $ P $ at a rate $ V_M $. The second one represents the conversion of precursors into macromolecules $ R $ and $ M $ and is catalyzed by $ R $, at a rate $ V_R $. Finally, the third reaction describes the transformation of precursors $ P $ into product $ X $ at a rate $ V_X $, and catalyzed by $ M $. The natural resource allocation parameter is modeled through the dimensionless parameter $ \alpha \in [0, 1] $, that represents the proportion of precursors allocated to the gene expression machinery $ R $, while $ 1-\alpha $ indicates that of the metabolic machinery $ M $. The model describes all mass quantities $ S $, $ P $, $ M $, $ R $ and $ X $ in grams, and the rates in grams per hour. A constant volumetric flow rate $ F $ [L h$ ^{-1} $] produces an inflow of fresh medium rich in substrate, and an outflow of bacterial culture and metabolites [24]. As already stated, we include in our scheme the growth switch described in [13] to externally affect the allocation decision $ \alpha $ by varying the inducer concentration in the medium. The mechanism is modeled as
$ u(t)=I(t)α(t),u∈[0,1] $ | (2.1) |
where $ I $ is the external signal, which acts independently from the allocation parameter $ \alpha $. The agammaegated form of Eq (2.1) accounts for the fact that the synthetic switch is capable of adjusting bacterial growth between zero (by setting $ I = 0 $, which yields $ u = 0 $) and the maximal growth rate supported by the medium (given by $ I = 1 $, which yields $ u = \alpha $). In this work, we limit the analysis to calculating the control input $ u $, without decoupling both individual controls $ I $ and $ \alpha $ (refer to Section 5 in [10] for more details on the implementation of this external controller). Then, we write the system of differential equations describing the evolution of the mass of each component,
$ {˙S=VSin−VM−VSout,˙P=VM−VR−VX−VPout,˙M=(1−u)VR−VMout,˙R=uVR−VRout,˙X=VX−VXout. $ |
The inflow/outflow rates are defined as
$ VSout=DS,VPout=DP,VMout=DM,VRout=DR,VXout=DX,VSin=Fsin, $ |
where $ s_{in} $ [g L$ ^{-1} $] is the nutrient concentration of the inflow of fresh medium, and $ D $ [h$ ^{-1} $] the dilution rate defined as
$ D≐FVext. $ |
We define the volume of the cell population $ \mathcal{V} $ [L] as
$ V≐β(M+R), $ | (2.2) |
where $ \beta $ [L g$ ^{-1} $] is the inverse of the cytoplasmic density. The above definition is based on the assumption that the cytoplasmic density of the cells is constant for the whole culture, and it also takes into account the experimental results showing that macromolecules are responsible for most of the biomass in microbial cells [25]. Thus, the mass of precursors $ P $ is excluded from the volume $ \mathcal{V} $ as a simplifying assumption. Then, we express the quantities of the system as concentrations with respect to the volumes,
$ p≐PV,r≐RV,m≐MV,s≐SVext,x≐XVext, $ | (2.3) |
where $ p $, $ r $ and $ m $ [g L$ ^{-1} $] are intracellular concentrations of precursor metabolites, components of the gene expression machinery and metabolic enzymes respectively; and $ s $ and $ x $ [g L$ ^{-1} $] the extracellular concentrations of substrate and metabolite. It is worth stressing that intracellular concentrations are defined with respect to the bacterial volume $ \mathcal{V} $, while extracellular concentrations with respect to the volume of the bioreactor $ \mathcal{V}_\mathrm{ext} $. Then, using definition (2.2), we obtain that $ m + r = 1/\beta $. We define the growth rate of the bacterial population $ \mu $ [h$ ^{-1} $] as the relative variation of cell volume $ \dot{ \mathcal{V}}/ \mathcal{V} $ without considering the effect of the outflowing biomass. Replacing with concentrations leads to the system
$ S:{˙s=D(sin−s)−vM(s,m)VVext,˙p=vM(s,m)−vR(p,r)−vX(p,m)−μ(p,r)p,˙r=uvR(p,r)−μ(p,r)r,˙m=(1−u)vR(p,r)−μ(p,r)m,˙x=vX(p,m)VVext−Dx,˙V=(μ(p,r)−D)V, $ | (2.4) |
where $ v_M(s, m) $, $ v_R(p, r) $ and $ v_X(p, m) $ [g L$ ^{-1} $ h$ ^{-1} $] are the mass fluxes per unit volume obtained from dividing the rates $ V_M $, $ V_R $ and $ V_X $ by $ \mathcal{V} $, function of the concentrations of system $ S $. In this new system, the growth rate becomes
$ μ(p,r)≐˙VV|F=0=˙M+˙RM+R|F=0=βvR(p,r) $ |
showing that the bacterial growth rate is proportional to the macromolecule synthesis rate. We propose a change of variables that simplifies the expressions of the system
$ ˆs=βs,ˆp=βp,ˆr=βr,ˆm=βm,ˆx=βx,ˆV=VVext, $ |
which yields the relation
$ ˆr+ˆm=1, $ | (2.5) |
and we define the non-dimensional synthesis rates
$ ˆvM(ˆs,ˆm)=βvM(s,m),ˆvR(ˆp,ˆr)=βvR(p,r),ˆvX(ˆp,ˆm)=βvX(p,m), $ | (2.6) |
and the non-dimensional substrate concentration $ \hat s_{in} = \beta s_{in} $. Then, dropping all hats yields the following system
$ S1:{˙s=D(sin−s)−vM(s,1−r)V,˙p=vM(s,1−r)−vX(p,1−r)−μ(p,r)(p+1),˙r=(u−r)μ(p,r),˙x=vX(p,1−r)V−Dx,˙V=(μ(p,r)−D)V, $ | (2.7) |
where the dynamical expression of $ m $ has been omitted since it can be computed from $ r $, as shown in Eq (2.5). It can also be seen that both concentrations $ m $ and $ r $ are limited to the interval $ [0, 1] $ due to physical constraints from the relation in Eq (2.5). For this latter, and due to the nature of the reactions involved in the studied problem, we will make some assumptions on the synthesis rates $ v_M(s, m) $, $ v_R(p, r) $ and $ v_X(p, m) $:
Assumption 1. Functions $ v_M(s, m) $, $ v_R(p, r) $ and $ v_X(p, m) $ meet
● $ v_i(y, z): \mathbb{R}^+ \times [0, 1] \rightarrow \mathbb{R}^+ $
● $ v_i(y, z) $ continuously differentiable w.r.t. both variables
● $ v_i(0, z) = v_i(y, 0) = 0 $
● $ v_i(\cdot) $ strictly monotonically increasing:
$ ∂vi∂y(y,z)>0,∀(y,z)∈R>0×(0,1],∂vi∂z(y,z)>0,∀(y,z)∈R>0×(0,1] $ |
● $ v_i(\cdot) $ bounded w.r.t. $ y $: $ \lim_{y\rightarrow \infty}{v_i(y, z)} = {v_i}_{, max}(z) $.
Assumption 1 encompasses all general monotone increasing kinetics models used in biological models, such as Michaelis-Menten or Hill equation-based kinetics [26]. In this first work, we will focus on a particular kind of systems where the synthesis rate related to the metabolite production depends on the growth rate:
Assumption 2. For $ r \in (0, 1) $, the metabolite synthesis rate $ v_X(p, 1-r) $ can be expressed in terms of the macromolecule synthesis rate $ v_R(p, r) $ (i.e., the growth rate)
$ vX(p,1−r)=c(r)vR(p,r), $ |
where $ c(r) : (0, 1) \rightarrow \mathbb{R^+} $ is a positive continuously differentiable function.
As previously described, the reaction $ v_X(p, m) $ is catalyzed by $ m $, and the reaction $ v_R(p, r) $ is catalyzed by $ r $, meaning that the ratio between $ M $ and $ R $ in the microbial culture determines whether the resources are being allocated to the production of biomass or metabolite. This represents the trade-off described in the introduction of this paper, which is here modeled through the function $ c(r) $. The fact that $ c $ does not depend on the concentration of precursors implies that the host cell has the same affinity to synthesize both biomass and metabolite from the precursors, even when the reactions are not expected to consume the precursors in the same proportion. For the particular case of Michaelis-Menten kinetics, this phenomenon is captured by the half-saturation constant [27]. Notably, for fixed values of $ r $, both synthesis rates are simply proportional. This assumption reduces $ S_1 $ to
$ S1:{˙s=D(sin−s)−vM(s,1−r)V,˙p=vM(s,1−r)−μ(p,r)(p+c(r)+1),˙r=(u−r)μ(p,r),˙x=c(r)μ(p,r)V−Dx,˙V=(μ(p,r)−D)V. $ | (2.8) |
We note that both Assumptions 1 and 2 are formulated for the non-dimensional synthesis rates given in definition (2.6), but they also hold for the original functions $ v_M $, $ v_R $ and $ v_X $ (taking into account that, for these functions, the domain is defined as $ \mathbb{R}^+ \times [0, \frac{1}{\beta}] \rightarrow \mathbb{R}^+ $).
The asymptotic behavior of system $ S_1 $ describes the “open-loop” operation mode of the continuous bioreactor, where the resource allocation control $ u(t) $ is fixed to $ \bar u \in (0, 1) $. In the present section we propose a series of mass conservation laws that allow to reduce $ S_1 $ to a 3-dimensional limiting system. Then, we study the local stability of their equilibria, and we show, using the theory of asymptotically autonomous systems, that the full system $ S_1 $ converges to the equilibria of its limiting system. Let us start the analysis of system $ S_1 $ by defining its invariant region in the following lemma.
Lemma 1. The set
$ Γ={(s,p,r,x,V)∈R5:sin≥s>0,p≥0,x≥0,1≥r≥0,V≥0} $ |
is positively invariant for the initial value problem.
Proof. Let us analyze the boundaries of $ \Gamma $,
$ {˙s|s=0=Dsin>0⇒s=0isrepulsive.˙s|s=sin=−vM(s,1−r)V≤0⇒s=sinisrepulsive/invariant.˙p|p=0=vM(s,1−r)≥0⇒p=0isrepulsive/invariant.˙r|r=0=0⇒r=0isinvariant.˙r|r=1=(u−1)μ(p,1)<0⇒r=1isrepulsive.˙x|x=0=vX(p,m)V≥0⇒x=0isrepulsive/invariant.˙V|V=0=0⇒V=0isinvariant.} $ |
We will study the initial value problem of system $ S_1 $ with initial conditions
$ sin≥s(0)≥0,p(0)≥0,1≥r(0)>0,x(0)≥0,V(0)>0, $ | (2.9) |
where two cases have been excluded for being trivial to the analysis: $ \mathcal{V}(0) = 0 $, since it is necessary to have an initial amount of biomass to have bacterial growth; and $ r(0) = 0 $, since an empty gene expression machinery pool implies null growth rate $ \mu(p, 0) = 0 $, and therefore it is not possible to self-replicate from that point.
System $ S_1 $ can be rewritten as
$ {˙φ=D(sinvin−vout)+Nvi−vμμ(p,r),˙V=(μ(p,r)−D)V, $ |
where
● $ \varphi \doteq [s, p, r, m, x]^T $ is the state vector of concentrations in the bioreactor.
● $ N $ is the stoichiometry matrix of the macroreactions.
● $ v_i $ is the vector of internal synthesis rates.
● $ v_{in} $ and $ v_{out} $ are the vectors of inflows and outflows respectively, associated to the continuous bioreactor setup.
● $ v_\mu $ is the vector modeling the dilution effect due to variation of the bacterial volume.
which are defined as
$ N≐[−V001−1−10u001−u000V],vi≐[vM(s,1−r)vR(p,r)vX(p,1−r)],vin≐[1,0,0,0,0]T,vout≐diag(φ)[1,0,0,0,1]T,vμ≐diag(φ)[0,1,1,1,0]T. $ |
By analyzing the left null space of $ N $, it can be seen that there are two mass conservation laws related to the total mass inside the bioreactor.
Definition 1. We define the quantities
$ w1≐s+(p+m+r)V+x=s+(p+1)V+x, $ |
$ w2≐s+(p+rˉu)V+x. $ |
The first quantity tends asymptotically to $ w_1 = s_{in} $ as $ t \rightarrow \infty $ for every input $ u(t) $, as it obeys the dynamical equation
$ ˙w1=D(sin−w1). $ | (2.10) |
Moreover, when fixing $ u(t) $ to $ \bar u \in (0, 1) $, the quantity $ w_2 $ also obeys the same Eq (2.10), meaning that this second quantity also converges to $ s_{in} $, which greatly simplifies the analysis of the asymptotic behavior of the system.
Lemma 2. The $ \omega $-limit set of any solution of system $ S_1 $ given by Eq (2.8) lies in the hyperplane
$ Ω1≐{(s,p,r,x,V)∈R5:s+(p+1)V+x=sin}. $ |
Moreover, under constant input $ u(t) = \bar u $, this is also true for the hyperplane
$ Ω2≐{(s,p,r,x,V)∈R5:s+(p+rˉu)V+x=sin}. $ |
Further on, we will use Lemma 2 to analyze system $ S_1 $ through its limiting system.
Lemma 2 presents two mass conservation laws that can be used to reduce subsystem $ S_1 $ by two dimensions. We will first analyze the asymptotic behavior of concentration $ r $: when $ t \rightarrow \infty $, quantities $ w_1 = w_2 = s_{in} $, so
$ s+(p+1)V+x=s+(p+rˉu)V+x⇒r=ˉu $ |
meaning that, as $ t \rightarrow \infty $, $ r $ will converge to the value $ \bar u $. We can also express $ x = s_{in} -s - (p+1) \mathcal{V} $, so that the limiting system of $ S_1 $ becomes
$ S′1:{˙s=D(sin−s)−ˉvM(s)V,˙p=ˉvM(s)−ˉμ(p)(p+ˉc+1),˙V=(ˉμ(p)−D)V, $ |
where
$ ˉvM(s)≐vM(s,1−ˉu),ˉvR(p)≐vR(p,ˉu),ˉvX(p)≐vX(p,1−ˉu),ˉμ(p)≐μ(p,ˉu),ˉc≐c(ˉu). $ |
Details on the convergence of the limiting system $ S_1' $ to the original one $ S_1 $ will be addressed later in the article. In next section, we will fully describe the asymptotic behavior of $ S_1' $.
In the interest of simplifying the notation, we define the following function.
Definition 2. We define the function
$ ˉf(p)≐ˉvR(p)+ˉvX(p)+ˉμ(p)p=ˉμ(p)(p+ˉc+1), $ |
Function $ f $ meets $ \bar f(p) > 0, \bar f'(p) > 0, \forall p \in \Gamma $ (positive and monotonically increasing).
The main result of local stability study is summarized in Theorem 1.
Theorem 1. The local stability of equilibria is given by the following criterion.
● If $ \bar \mu(p_w) \geq D $:
- The interior equilibrium $ E_i $ exists, is unique and locally stable.
- The washout equilibrium $ E_w $ exists, is unique and locally unstable.
● If $ \bar \mu(p_w) < D $:
- The interior equilibrium $ E_i $ does not exist.
- The washout equilibrium $ E_w $ exists, is unique and locally stable.
In the above criterion, the equilibria are defined as follows:
● The interior equilibrium $ E_i \doteq (s_i, p_i, \mathcal{V}_i) $, with
$ pi:{p∈R+:ˉμ(p)=D}, $ | (2.11) |
$ si:{s∈R+:ˉvM(s)=ˉf(pi)}, $ | (2.12) |
$ Vi≐D(sin−si)ˉvM(si). $ | (2.13) |
● The washout equilibrium $ E_w \doteq (s_{in}, p_w, 0) $, with
$ pw:{p∈R+:ˉf(p)=ˉvM(sin)}. $ | (2.14) |
Proof. First, we will prove the boundedness of $ p $. Using definition (2.14), it is possible to define a constant upper-bound on $ p $ by analyzing the time-varying upper bound $ p_{up}(t) $ with dynamical equation
$ ˙pup≐ˉvM(sin)−ˉf(pup). $ | (2.15) |
In Eq (2.15), $ p_{up} $ converges to the value $ p_w $ satisfying Eq (2.14). Additionally, the vector field at $ p = p_w $ is always negative (or null when $ s = s_{in} $)
$ ˙p|p=pw=ˉvM(s)−ˉf(pw)≤0 $ |
meaning that $ p = p_w $ is either repulsive or invariant, so a new invariant set $ \Gamma' \subset \Gamma $ can be defined,
$ Γ′≐{(s,p,V)∈R3:sin≥s>0,pw≥p>0,V≥0}. $ | (2.16) |
Equilibrium $ E_i $ The existence of the interior equilibrium $ E_i $ is given by the boundedness of the flows,
$ maxˉμ(p)≥D, $ | (2.17) |
$ maxˉvM(s)≥ˉf(pi), $ | (2.18) |
and its uniqueness can be proved through monotonicity arguments: In $ \bar \mu(p) = D $, $ \bar \mu(p) $ is strictly monotonically increasing w.r.t. $ p $ so, if inequality (2.17) is met, $ p_i $ should be unique. Similarly, in $ \bar v_M(s) = \bar f(p_i) $, $ \bar v_M(s) $ is again strictly monotonically increasing so, if inequality (2.18) is met, $ s_i $ should be unique. As is standard in continuous bioreactors, inequality (2.17) implies that the maximal growth rate of the bacterial population should be bigger than the dilution rate $ D $. At the same time, the inequality (2.18) requires the maximal uptake flow to be bigger than the flows responsible for bacterial growth and metabolite production. In $ \Gamma' $, these conditions become
$ ˉμ(pw)≥D, $ | (2.19) |
where the second condition (2.18) is included in the inequality (2.19), as $ max \; \bar v_M(s) = \bar v_M(s_{in}) = \bar f(p_w) $, and the condition $ \bar f(p_w) \geq \bar f(p_i) $ is true if and only if $ \mu(p_w) \geq D $. The Jacobian matrix is given by
$ Ji=[−D−ˉv′M(s)Vi0−ˉvM(s)ˉv′M(s)ˉf′(p)00ˉμ′(p)Vi0], $ |
and the characteristic polynomial is
$ Pi(λ)=(λ+D+ˉv′M(s)Vi)(λ+ˉf′(p))λ+ˉvM(s)ˉv′M(s)ˉμ′(p)Vi $ | (2.20) |
$ =λ3+λ2(D+ˉv′M(s)Vi+ˉf′(p))+λ(D+ˉv′M(s)Vi)ˉf′(p)+ˉvM(s)ˉv′M(s)ˉμ′(p)Vi. $ |
Using Definition 2, it can be seen that $ -D $ is an eigenvalue by replacing in expression (2.20)
$ Pi(−D)=−Dˉv′M(s)Vi(p+ˉc+1)ˉμ′(p)+ˉvM(s)ˉv′M(s)ˉμ′(p)Vi=ˉv′M(s)Viˉμ′(p)˙p=0. $ |
Then, dividing $ P_i(\lambda) $ by $ \lambda + D $ yields
$ Pi(λ)1λ+D=λ2+λ(D+ˉv′M(s)Vi+(p+ˉc+1)ˉμ′(p))⏟>0+ˉv′M(s)Vi(p+ˉc+1)ˉμ′(p)⏟>0, $ |
which, by Routh-Hurwitz criteria, implies that all eigenvalues have negative real part. Then, if it exists, the equilibrium is always stable.
Equilibrium $ E_w $ The washout equilibrium $ E_w $ exists for all values of $ \Gamma' $, since the only condition for existence is given by
$ maxˉf(p)≥ˉvM(sin), $ |
which is always true in $ \Gamma' $ since $ \bar v_M(s_{in}) = \bar f(p_w) $ and so the inequality becomes $ \bar f(p_w) \geq \bar f(p_w) $. Again, as $ \bar f(p) $ is strictly monotonically increasing, there is a unique solution $ p_w $. Its stability is given by
$ Jw=[−D0−ˉcˉv′M(s)−ˉf′(p)000ˉμ(pw)−D] $ |
with characteristic polynomial
$ P_w(\lambda) = (\lambda + D)(\lambda + \bar f'(p))(\lambda - \bar \mu(p_w) + D). $ |
Eigenvalues $ \lambda = -D $ and $ \lambda = - \bar f'(p) $ are always real and negative, and so the stability criterion becomes
$ Ewstableif and only ifˉμ(pw)<D. $ |
In the current section, we show how the limiting system $ S'_1 $ can be further reduced using another mass conservation law given by the fact that $ \bar c $ is constant. This is formalized in the following lemma.
Lemma 3. The $ \omega $-limit set of any solution of the limiting system $ S'_1 $ lies in the hyperplane
$ Ω3≐{(s,p,V)∈R3:s+(p+ˉc+1)V=sin} $ |
Proof. The quantity
$ w3≐s+(p+ˉc+1)V $ | (2.21) |
obeys the dynamical equations $ \dot w_3 = D(s_{in} - w_3) $, which means that the system converges asymptotically to the limit set $ \Omega_3 $.
Using Lemma 2, we can further reduce system $ S_1' $ by expressing
$ s=sin−V(p+ˉc+1), $ | (2.22) |
for the values $ (p, \mathcal{V}) $ that meet $ s_{in} - \mathcal{V}(p + \bar c + 1) > 0 $, and so the new limiting system becomes
$ S1″ $ |
Since $ S''_1 $ is a 2-dimensional continuous system, its global behavior (illustrated in Figure 2) can be studied through the Poincaré-Bendixson trichotomy.
Lemma 4. Every solution of the limiting system $ S_1'' $ with initial conditions (2.9) converges to
● $ E_i^* \doteq (p_i, \mathcal{V}_i) $ if $ \bar \mu(p_w) \geq D $
● $ E_w^* \doteq (p_w, 0) $ if $ \bar \mu(p_w) < D $
Proof. The Poincaré-Bendixson trichotomy ensures that every non-empty compact $ \omega $-limit set of $ S''_1 $ is either a fixed point, a periodic orbit or a cycle of equilibria. Then, by applying Poincaré-Bendixson theorem through the Dulac criterion, we can discard periodic orbits and cycles of equilibria:
$ \begin{align} \frac{\partial}{\partial p} \dot p + \frac{\partial}{\partial \mathcal{V}} \mathcal{\dot V} = &\frac{\partial}{\partial p} \bar v_M\left(s(\cdot)\right) - \bar \mu(p) - \bar \mu'(p) \left(p + \bar c + 1\right) + \bar \mu(p) - D \end{align} $ |
$ \begin{align} = &\bar v_M' \left(s(\cdot)\right) \frac{\partial s(\cdot)}{\partial p} - \bar \mu'(p) \left(p + \bar c + 1\right) -D \lt 0, \end{align} $ | (2.23) |
as $ \partial s(\cdot)/\partial p = - \mathcal{V} $ from Eq (2.22). This ensures that the new limiting system should converge to one of the stable equilibria as $ t \rightarrow \infty $ which, according to Theorem 1, is known to be unique.
Through the theory of asymptotically autonomous systems [22], we can relate the asymptotic behavior of the 2-dimensional limiting system $ S_1'' $ to that of the full 5-dimensional system $ S_1 $. The main result is established in Theorem 2.
Theorem 2. Every solution of system $ S_1 $ with initial conditions (2.9) converges to,
● The extended interior equilibrium $ \hat E_i \doteq (s_i, p_i, \bar u, x_i, \mathcal{V}_i) $ if it exists, with $ x_i \doteq \bar c \mathcal{V}_i $.
● The extended washout equilibrium $ \hat E_w \doteq (s_{in}, p_w, \bar u, 0, 0) $ if $ \hat E_i $ does not exist.
where the condition for the existence of the extended interior equilibrium is $ \bar \mu(p_w) \geq D $.
Proof. We resort to a more particular case of the general theory of asymptotically autonomous systems (introduced in Appendix F of [24]), that requires a certain number of hypotheses to be met related to the limiting system $ S_1'' $, its equilibria and the mass conservation equations defined in Lemma 2:
(A1) The dynamical Eqs (2.10) and (2.21) of quantities $ w_1 $, $ w_2 $ and $ w_3 $ are stable.
(A2) $ S_1'' $ has 2 rest points $ E_i^* $ and $ E_w^* $, which are hyperbolic.
(A3) If $ E_i^* $ exists, $ dim(M^+(E_i^*)) = 2 $ and $ dim(M^+(E_w^*)) = 1 $. If $ E_i^* $ does not exist, $ dim(M^+(E_w^*)) = 2 $.
(A4) $ S_1'' $ has no cycles of rest points, as shown in Eq (2.23).
(A5) $ S_1'' $ has no periodic orbits, as shown in Eq (2.23).
Then, almost all trajectories of the original system $ S_1 $ converge to one of the asymptotically stable equilibria of the limiting system, which is always unique.
In this section, we aim to maximize the production of the metabolite of interest $ X $. We start by modeling the synthesis rates as explicit functions of the concentrations of the system. While this latter is not required for most of the results, it enables the numerical simulations both for the dynamical and for the static optimization problems. In particular, we resort to the Michaelis-Menten kinetics defined in previous works [4] and the same biological constants. Then, we pose the problem of maximizing the production of $ X $ over a fixed period of time, and we approximate the optimal solution with a constant allocation strategy. This finally motivates the static optimization approach, that is solved in terms of the tuple $ (D, \bar u) $.
Michaelis-Menten kinetics are defined as
$ \begin{align} v_M(s,m) &\doteq k_M \, m \frac{s}{K_M + s} = k_M \, (1-r) \frac{s}{K_M + s}, \end{align} $ |
$ \begin{align} v_R(p,r) &\doteq k_R \, r \frac{p}{K_R + p}, \end{align} $ |
$ \begin{align} v_X(p,m) &\doteq k_X \, m \frac{p}{K_X + p} = k_X \, (1-r) \frac{p}{K_X + p}, \end{align} $ |
where the biological constants and their values [4,10] are described in Table 1.
Parameter | Description | Unit | Value |
$ \beta $ | Inverse of the cytoplasmic density | $ L \, g^{-1} $ | $ 0.003 $ |
$ k_M $ | Rate constant of the metabolic macroreaction | $ h^{-1} $ | $ 4.32 $ |
$ K_M $ | Half-saturation constant of the metabolic macroreaction | $ g \, L^{-1} $ | $ 33.33 $ |
$ k_R $ | Rate constant of the protein synthesis reaction | $ h^{-1} $ | $ 3.6 $ |
$ K_R $ | Half-saturation constant of the protein synthesis reaction | $ g \, L^{-1} $ | $ 1 $ |
$ k_X $ | Rate constant of the metabolite synthesis reaction | $ h^{-1} $ | $ 0.5 $ |
$ K_X $ | Half-saturation constant of the metabolite synthesis reaction | $ g \, L^{-1} $ | $ 1 $ |
Assumption 2 is implemented by setting the half-saturation constant $ K_X = K_R $, such that the function $ c(r) $ becomes
$ \begin{align} c(r) = \frac{v_X(p,1-r)}{v_R(p,r)} = \frac{k_X}{k_R} \frac{1 - r}{r}. \end{align} $ |
The substrate concentration of the inflow $ s_{in} $ is set to $ 0.4 $ g/L.
The problem of maximizing the production of the metabolite $ X $ during a fixed interval of time $ T $ is explored in this section. As it is classical in the continuous bioreactor framework, the instantaneous production of metabolite is described by the quantity
$ \begin{align} DX \quad [g \, h^{-1}] \end{align} $ | (3.1) |
which, using definitions (2.3), can be expressed as $ D \mathcal{V}_\mathrm{ext} x $. Then, the total metabolite production over an interval $ T $ amounts to
$ \begin{align} \int_0^T D \, \mathcal{V}_\mathrm{ext} \, x \, \mathrm{d}t \quad [g]. \end{align} $ | (3.2) |
We will be first interested in dynamically adjusting the agammaegated control $ u(t) $ so as to maximize the quantity (3.2) for a fixed dilution rate $ D $. The problem is tackled through an optimal control approach, formulated as
$ \begin{align} (OCP) \, : \, \left\{{\begin{array}{rl} maximize & J_x(u) = D \, \mathcal{V}_\mathrm{ext} \int_0^T x(t) \, \mathrm{d}t \\ subject\; to & dynamics\; of\; S_1, \\ & u(\cdot) \in \mathcal{U}. \end{array}}\right. \end{align} $ |
with $ \mathcal{U} $ the set of admissible controllers, which are Lebesgue measurable real-valued functions defined on the interval $ [0, T] $ and satisfying the constraint $ u(t) \in [0, 1] $. We first note that, since $ D $ and $ \mathcal{V}_\mathrm{ext} $ are constants, the problem reduces to maximizing the integral of the concentration $ x $. Thus, the nature of the OCP allows us to obtain a first result on the existence of the solutions.
Proposition 3. The dynamic maximization problem has at least one solution.
Proof. Since there are no terminal conditions, the set of admissible controls is not empty (any constant control within the prescribed bounds is admissible). For bounded controls ($ u(t) $ belongs to $ [0, 1] $) and for a fixed final time $ T > 0 $, the dynamics (2.8) cannot blow up in finite time so all trajectories remain in a fixed compact. Indeed, all state variables but $ \mathcal{V} $ are bounded (as set (2.16) is invariant), and since
$ \dot{\mathcal{V}} = (\mu(p,r)-D)\mathcal{V} $ |
with $ \mu(p, r) $ bounded, the volume has at worst an exponential rate and is also bounded. As the dynamics is affine in the control, the set of velocities is convex and existence holds by Filippov's theorem [28].
In order to further explore the solution of (OCP), we define the adjoint state vector $ \lambda \doteq (\lambda_s, \lambda_p, \lambda_r, \lambda_x, \lambda_ \mathcal{V}) $ associated to the state vector $ \varphi \doteq (s, p, r, x, \mathcal{V}) $, and we write the Hamiltonian
$ \begin{align} H(\varphi, \delta, u) = \lambda^0x + \langle \lambda, F(\varphi, u) \rangle, \end{align} $ |
where $ F $ represents the right-hand side of system $ S_1 $ given by Eq (2.8). Developing the expression we obtain
$ \begin{align} H = \, &\lambda_s \left( D(s_{in} - s) - v_M(s,1-r) \mathcal{V} \right) + \lambda_p \Big( v_M(s,1-r) - v_X(p,1-r) - \mu(p,r)(p + 1) \Big) \\ &+ \lambda_r \left( u - r \right) \mu(p,r) + \lambda_x \Big( v_X(p,1-r) \mathcal{V} - Dx \Big) + \lambda_ \mathcal{V} (\mu(p,r) - D) \mathcal{V} - \lambda_0 x, \end{align} $ |
which shows that the Hamiltonian is linear in the control $ u $, so it can be rewritten in the affine form
$ \begin{align} H = H_0 + u H_1, \end{align} $ |
where
$ \begin{align} H_0 = \, &\lambda_s \left( D(s_{in} - s) - v_M(s,1-r) \mathcal{V} \right) + \lambda_p \Big( v_M(s,1-r) - v_X(p,1-r) - \mu(p,r)(p + 1) \Big) \end{align} $ |
$ \begin{align} &-r \lambda_r \mu(p,r) + \lambda_x ( v_X(p,1-r) \mathcal{V} - Dx) + \lambda_ \mathcal{V} (\mu(p,r) - D) \mathcal{V} - \lambda_0 x, \\ H_1 = \, &\lambda_r \mu(p,r). \end{align} $ |
In the absence of terminal constraints, there are no abnormal extremals and $ \lambda_0 $ can be set to $ -1 $. Since the constrained optimal control $ u $ should maximize the Hamiltonian, one has
$ \begin{align} u_{OCP}(t) = \left\{ \begin{array}{cl} 0 & if \, H_1 \lt 0, \\ 1 & if \, H_1 \gt 0, \\ u_{s}(t) & if \, H_1 = 0. \end{array}\right. \end{align} $ |
These alternatives account for the possibility of bang and singular arcs, an extremal being in general an arbitrary concatenation of these. Bang controls correspond to pure allocation strategies $ u = 0 $ and $ u = 1 $ (i.e., purely geared towards either the metabolic machinery $ M $ or the gene expression machinery $ R $ respectively), while singular control $ u_s(t) $ occurs when the function $ H_1 $ identically vanishes over some subinterval. Such singular arcs can be further described by successively differentiating the switching function $ H_1 $ until the singular control $ u $ can be explicitly computed. For the present case, this occurs when differentiating four times (“order two” singular arc), showing that the singular regime should be entered and leaved through Fuller phenomenon or “chattering” [29] (an infinite number of switchings between the control bounds $ 0 $ and $ 1 $ over a finite time interval). The theorem hereunder formalizes this statement.
Theorem 4. Any singular arc $ u_s(t) $ of a normal extremal process solution of (OCP) is at least of order two.
Proof. We assume that the process $ (\varphi, u, \lambda) $ that satisfies PMP is a normal extremal, and we set $ \lambda_0 = -1 $. In order to further describe the singular arc, we assume $ H_1 $ vanishes on a whole sub-interval $ \tau = [t_1, t_2] \subset [0, T] $. Then, the switching surface is the set
$ \begin{align} \Sigma = \left\{ \, (\varphi, \lambda) \in \mathbb{R}^{2n} \, | \, H_1 = 0 \, \right\}, \end{align} $ |
where, for this case, $ n = 5 $. The first derivative is
$ \begin{align} \dot H_1 = \dot \lambda_r \mu(p,r) + \lambda_r \Big(\mu_r(p,r) \dot r + \mu_p(p,r) \dot p \Big), \end{align} $ | (3.3) |
where
$ \mu_r(p,r) = \frac{\partial \mu(p,r)}{\partial r}, \quad \mu_p(p,r) = \frac{\partial \mu(p,r)}{\partial p}. $ |
Along the singular arc $ H_1 $ is identically zero, so expression (3.3) also vanishes. In order to compute the successive derivatives of $ H_1 $, we will resort to the Poisson bracket operator [30],
$ \begin{align} \{f, g\} = \sum\limits_{i = 1}^n \left( \frac{\partial f}{\partial \lambda_i}\frac{\partial g}{\partial \varphi_i} - \frac{\partial f}{\partial \varphi_i}\frac{\partial g}{\partial \lambda_i}\right). \end{align} $ |
Applying the latter definition to the derivative of $ H_1 $ we obtain
$ \begin{align} \dot{H_1} & = \frac{\partial H_1}{\partial \varphi} \dot{\varphi} + \frac{\partial H_1}{\partial \lambda} \dot{\lambda} = \sum\limits_{i = 1}^{n}\left( \frac{\partial H}{\partial \lambda_i} \frac{\partial H_1}{\partial \varphi_i} - \frac{\partial H}{\partial \varphi_i} \frac{\partial H_1}{\partial \lambda_i} \right) = \left\{ H,H_1 \right\} = \left\{ H_0 + u H_1 ,H_1 \right\} = \left\{ H_0 ,H_1 \right\} \end{align} $ |
as $ \left\{ u H_1, H_1 \right\} = u \left\{ H_1, H_1 \right\} = 0 $. In this context, we will use the notation $ H_{01} $ to refer to $ \left\{ H_0, H_1 \right\} $, and so forth. The second derivative of $ H_1 $ is
$ \begin{align} \ddot H_1 = \dot H_{01} = \{ H, H_{01} \} = \{ H_0, H_{01} \} + u \{ H_1, H_{01} \} = H_{001} + u H_{101}. \end{align} $ |
Showing that $ H_{101} = 0 $ along the switching surface $ \Sigma $ ensures that the singular arc is at least of order $ 2 $ (otherwise, a singular arc of order 1 could be computed as $ u = -H_{001}/H_{101} $). Since $ H_1 $ depends only on $ \lambda_r $, $ p $ and $ r $, the computation of the Poisson bracket reduces to
$ \begin{align} H_{101} = \{ H_1, H_{01} \} & = \frac{\partial H_1}{\partial \lambda_r} \frac{\partial H_{01}}{\partial r} - \frac{\partial H_1}{\partial r} \frac{\partial H_{01}}{\partial \lambda_r} + \frac{\partial H_1}{\partial \lambda_p} \frac{\partial H_{01}}{\partial p} - \frac{\partial H_1}{\partial p} \frac{\partial H_{01}}{\partial \lambda_p} \\ & = \mu(p,r) \frac{\partial H_{01}}{\partial r} - \lambda_r \mu_r(p,r) \frac{\partial H_{01}}{\partial \lambda_r} - \lambda_r \mu_p(p,r) \frac{\partial H_{01}}{\partial \lambda_p} \end{align} $ | (3.4) |
where
$ \begin{align} \frac{\partial H_{01}}{\partial r} = \dot \lambda_r \mu_r(p,r) + \lambda_r \frac{\partial}{\partial r} \Big(\mu_r(p,r) \dot r + \mu_p(p,r) \dot p \Big). \end{align} $ | (3.5) |
Replacing the expression (3.5) in Eq (3.4) yields
$ \begin{align} H_{101} = \mu(p,r) \left[ \dot \lambda_r \mu_r(p,r) + \lambda_r \frac{\partial}{\partial r} \Big(\mu_r(p,r) \dot r + \mu_p(p,r) \dot p \Big) \right] - \lambda_r \mu_r(p,r) \frac{\partial H_{01}}{\partial \lambda_r} - \lambda_r \mu_p(p,r) \frac{\partial H_{01}}{\partial \lambda_p} \end{align} $ |
which is also equal to $ 0 $ along the switching surface since every term is multiplied either by $ \lambda_r $ or $ \dot \lambda_r $, both identically zero along a singular arc (see $ H_1 $ expression), showing that the singular arc is at least of order 2.
It is noteworthy that the above proof holds for all general flows $ v_M(s, 1-r) $, $ v_X(p, 1-r) $ and $ v_R(p, r) $ considered in Assumption 1, with no need to use the defined Michaelis-Menten kinetics. Numerical results shown in Figure 3 confirm that, when the conditions for the existence of the interior equilibrium $ E_i $ are met, the optimal control strategy consists on a series of bang-bang arcs, and a singular arc that is entered and left through the chattering phenomenon. While, as already mentioned, it is compulsory to enter and leave order two singular arcs through chattering, a more precise description of the structure of the extremal would require a deeper analysis (for instance to prove that there is only one singular arc). Moreover, when the simulation time is long enough, we observe that the singular control $ u_s(t) $ converges to a constant value, eventually taking the system to steady-state. {This is related to the so-called turnpike phenomenon [31] that relates the singular arc of the dynamic optimization problem with the solution of the static one. See, e.g., [32] for a preliminary analysis on a similar case.
Additionally, we observe that the steady-state approached by the system during the singular regime maximizes the integrand of the cost function, which is in fact the instantaneous production of metabolite described in expression (3.1) (otherwise, it would not be the optimal strategy). Therefore, the static optimal control to which the dynamic one converges along the singular arc is the solution of the optimization problem
$ \begin{align} (SP) \, : \, \left\{{\begin{array}{rll} maximize && J_{x}(\bar u) \doteq x \\ subject \;to && dynamics\; of\; S_1, \\ && \dot \varphi = 0, \\ && 0 \leq \bar u \leq 1. \end{array}}\right. \end{align} $ |
The dynamical optimal solution acts as the gold-standard in terms of what can be achieved through control techniques. However, implementing such strategy is, in practice, unfeasible. Consequently, it is possible to design a static suboptimal strategy consisting of a constant allocation $ u_{sp} $ that takes the system to the same steady-state which maximizes the integrand of the cost function. Numerical simulations of such strategy are shown in Figure 4, where it can be seen that the area below curves $ x_{sp} $ and $ x $ in OCP only differ slightly ($ < 1\% $ for this particular simulation). The constant control $ u_{sp} $ basically disregards the initial and final bang arcs, as well as the chattering phenomena, which constitutes a small fraction of the complete time horizon $ T $. Moreover, this fraction gets smaller as $ T $ becomes large, which is a typical feature of the turnpike effect for a specific class of optimal control problems [31]. This way, the difference between strategies becomes marginal in long-term production schemes.
Additionally, we provide a numerical computation of the successive derivatives of the switching function $ H_1 $ in Figure 5. It can be seen that the fourth derivative $ H_{10001} \neq 0 $ over the interval where $ H_1 = 0 $, which shows that the singular arc is of order $ 2 $. Furthermore, we verify that the Kelley (or generalized Legendre-Clebsch [34]) condition
$ \begin{align} (-1)^k \frac{\partial}{\partial u} \left[ \frac{d^{2k}}{dt^{2k}} \left( \frac{\partial H}{\partial u} \right) \right] \lt 0, \quad \forall t \in \mathcal{I} \end{align} $ | (3.6) |
is met along the singular arc, which in this case is equivalent to $ H_{10001} < 0, \, \forall t \in \mathcal{I} $. A check of this condition, necessary for optimality, is also shown in Figure 5. Although there is no available sufficient condition to test local optimality of extremals with Fuller arcs, verifying the Legendre-Clebsch condition along the singular arc only ensures that we do not compute a too crude local minimizer. Besides, numerically only a small finite number of bang arcs are retrieved by the optimizer for the chattering parts before and after the singular arc. This is usually sufficient to give a very good approximation of the solution.
We have shown that a constant allocation decision $ \bar u $ represents a simplified alternative to the optimal control solution, a strategy composed of bang arcs, a time-varying singular arc and the chattering artifact. In this section, we will further explore the static optimization problem by adding a second degree of freedom to the problem: the dilution rate $ D $. In addition to that, we investigate two objectives: the production of biomass $ \mathcal{V} $ and of metabolite $ X $. The static biomass maximization problem (SP$ _{ \mathcal{V}} $) can be written as
$ \begin{align} (SP_{ \mathcal{V}}) \, : \, \left\{ {\begin{array}{rll} maximize && J_{ \mathcal{V}}(\bar u, D) \doteq D \mathcal{V} \\ subject \;to && dynamics \;of\; S_1, \\ && \dot \varphi = 0, \\ && 0 \leq \bar u \leq 1. \end{array}}\right. \end{align} $ |
Analogously, the product maximization problem can be defined as
$ \begin{align} (SP_X) \, : \, \left\{{\begin{array}{rll} maximize && J_{X}(\bar u, D) \doteq DX = D \mathcal{V}_\mathrm{ext} x\\ subject \;to && dynamics\; of\; S_1, \\ && \dot \varphi = 0 ,\\ && 0 \leq \bar u \leq 1. \end{array}}\right. \end{align} $ |
Since we look for the steady-states that maximize each objective, the washout equilibrium $ E_w $ can be excluded from the analysis since, as shown in Theorem 2, the equilibrium corresponds to the steady-state values $ \mathcal{V} = 0 $ and $ x = 0 $. Therefore, the static problems are reduced to finding the equilibria $ E_i $ in terms of the pair $ (D, \bar u) $ that maximize each objective function. Moreover, it can be shown that the optimal solution cannot belong to the boundary of the equilibrium $ E_i $
$ \begin{align} \Theta \doteq \left\{ (\bar u, D) \in \mathbb{R}^2 \, : \, {\bar \mu(p_w) = D} \right\}. \end{align} $ |
Indeed, using the definitions (2.11), (2.12) and (2.13), it can be seen that on $ \Theta $ there is no bacterial population, as $ s_i = s_{in} $, which means that $ \mathcal{V}_i = x_i = 0 $. We formalize the latter reasoning in the following proposition.
Proposition 5. Every solution $ (D, \bar u) $ of the static optimization problems (SP$ _ \mathcal{V} $) and (SP$ _X $) is in the region of existence of $ E_i $ given by the condition $ \bar \mu(p_w) > D $.
When considering the same self-replicator scheme under constant environmental conditions (which could describe fed-batch cultivation) the solution for both static problems corresponds to the steady-state with maximal growth rate [10]. Proposition 5 shows that this is not the case in continuous bioreactors, where maximal growth rate is attained at the boundary set of existence of the equilibrium $ E_i $. Two interesting particular cases within $ \Theta $ are the pure static allocation strategies $ \bar u = 0 $ and $ \bar u = 1 $: A pure metabolic strategy $ \bar u = 0 $ will lead to a bacterial population with no RNA polymerase (i.e., $ r = 0 $), which will eventually stop the production of biomass (as $ v_R(p, 0) = 0 $), leading to washout in the bioreactor. Analogously, allocating all resources to the gene expression machinery will finally empty the metabolic machinery $ m $, halting the absorption of substrate from the environment (as $ v_M(p, 0) = 0 $) and depleting the bacteria from precursors, which will also lead to washout. From this analysis, we can conclude that any optimal steady-state allocation $ \bar u $ should belong to $ (0, 1) $.
We recall that, in continuous bioreactors, the growth rate $ \bar \mu(p_i) $ is fixed by the dilution rate $ D $, as shown in Eq (2.11). Figure 6 illustrates a numerical analysis of the static problems. It is interesting to notice in both subfigures 6a and 6b that the model accounts for the classical quasi-linear relation between the maximal growth rate (that lies on the boundary $ \Theta $), and the control $ \bar u $ which regulates the RNA/protein mass ratio of the bacterial population [35]. The latter is a phenomenon which has been first observed experimentally, and later used to develop dynamical self-replicator models for natural and biotechnological purposes [4,10,11]. However, in the present case, the growth rate is fixed through $ D $, so it is not a result of the nutrient quality in the environment, which enables the multivariate approach.
Biomass maximization objective Results for this problem are shown in Figure 6a. For this case, the curve $ \bar u_{opt} $ remains over $ 0.6 $ for all values of the growth rate, tending to $ \bar u_{opt} = 1 $ as the growth rate goes to $ 0 $. This shows that, in order to maximize $ \mathcal{V} $, the allocation strategy should prioritize the synthesis of macromolecules of the gene expression machinery $ R $ over the metabolic machinery $ M $, which catalyzes the production of biomass and not the synthesis of metabolites. Moreover, the maximum $ D \mathcal{V}_i $ is attained through a fairly high growth rate ($ \approx 85\% $ of the maximal growth rate).
Product maximization objective Results for this case are shown in Figure 6b. We can see that, in opposition to the first case, the optimal solution requires allocating as much precursors as possible into the metabolic machinery $ M $, no matter the value of the dilution rate. In other words, the allocation control $ \bar u $ should be as low as possible within the region of existence of the equilibrium $ E_i $ (but not in $ \Theta $). This result is consistent with the fact that the metabolic machinery $ M $ catalyzes the synthesis of metabolite $ X $. It is noteworthy that the optimal point $ D X_i $ is accomplished through a rather lower dilution rate $ D $ in comparison to the biomass production case, resulting in a continuous production at a rate of about $ 35 \% $ of the maximal growth rate. The latter might appear counter-intuitive, as it is well established in the literature that high dilution rates in continuous bioreactors imply high production rates. This characteristic can be attributed to the compromise between allocating resources to the metabolic machinery $ M $ and increasing the dilution rate, linked to the maximum value of the function $ c(r) $. In other words, the lower the dilution rate of operation, the wider the interval of existence of the equilibrium $ (u_{\min}, u_{\max}) $, which enables the possibility of further promoting the synthesis of compounds of the metabolic machinery (by artificially lowering the allocation parameter $ \bar u $) without going to washout. Figure 7a illustrates the resulting synthesis rates for each solution. In (SP$ _X $), the synthesis rate of the metabolite $ v_X $ is increased at the expense of reducing both precursor and biomass synthesis rates $ v_M $ and $ v_R $, which is in large part due to the reduction of the dilution rate $ D $. This difference in the flux distribution impacts directly on the mass quantities inside the bioreactor at steady-state—depicted in Figure 7b—for each problem: for the (SP$ _X $) case, there is a reduction in the biomass $ \mathcal{V} $ of only $ 20\% $ w.r.t. the (SP$ _ \mathcal{V} $) case. However, we can see an increase on the amount of metabolite $ X $ of about $ 5 $ times that of the (SP$ _ \mathcal{V} $) case. As already said, the allocation parameter $ \bar u $ becomes quite lower in the metabolite synthesis objective. In turn, this inhibits the synthesis of macromolecules, which translates into a decrease of the bacterial population's growth rate. To compensate this effect, the steady-state pool of precursors $ P $ required for producing the metabolite $ X $ becomes considerably bigger than that of the biomass production objective (about $ 10 $ times), therefore increasing the rate of fabrication of biomass $ v_R(p, r) $.
The objective in this paper was to synthesize a certain metabolite of interest by re-adjusting the natural allocation of resources in bacteria. This is done by drawing resources from the native pathways of the cell originally used for producing biomass, and allocating them into the production of the compound $ X $. We addressed the problem in the continuous bioreactor framework, which allows a steady-state production regime. Based on previous dynamical systems approaches [4,10,12], we proposed a coarse-grained self-replicator model capable of accounting for well-studied bacterial growth laws in a simplified way, and we studied its asymptotic behavior. We tackled two different production objectives: Biomass maximization $ \mathcal{V} $ and metabolite maximization $ X $, and we compared results in order to better understand the potential control strategies required to that effect. We rely on novel bio-engineering techniques capable of delivering groundbreaking control schemes: A synthetic growth switch that allows to control the transcription of RNA polymerase through an inducible promoter. In addition to this regulation mechanism, we include the dilution rate as a control input, which yields a multi-variable optimization problem. We concluded by showing very contrasting results, but in accordance with our previous understanding of microbial resource allocation. The biomass-oriented strategy involves an almost maximal dilution rate, and prioritizes investing resources into the gene expression machinery. Conversely, producing the compound of interest requires a rather low value of the dilution rate, which allows an allocation strategy more geared towards the synthesis of components of the metabolic machinery (i.e., a lower value of $ \bar u $). The latter shows there exist a compromise between augmenting the dilution rate and artificially boosting the production of enzymes and transporters.
The capacity to interfere with the natural bacterial behaviors at molecular level represents a promising tool to improve biotechnological processes. We are interested in further exploring the details of the implementation of such control schemes, by considering different models of the external growth switch. A next step in this regard would involve taking into account the nature of the external signal $ I $ and its physical constraints, in order to adapt our current results towards a more implementable control loop. Additionally, our approach is based on very simplified representations of bacterial growth. Such biological processes usually involve numerous reactions and variables. In our case, we clustered all macromolecules into only 2 different classes, and we purposely omitted a number of known phenomena in bacteria, such as cell division, protein degradation and the influence of temperature. However, some of these effects have been proven to affect only marginally the results regarding the allocation problem [36]. Thus, in our case, a simplified representation becomes useful to emphasize the effects of optimally dealing with the internal resource distribution of bacteria in industrial frameworks. Eventually, such optimal strategies could provide guidance for developing online feedback solutions based on real-time measuring of the process.
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). We acknowledge the support of the FMJH Program PGMO and the support to this program from EDF-THALES-ORANGE. We also thank Hidde de Jong (INRIA Grenoble - Rhône-Alpes) for discussions on the manuscript.
All authors declare no conflicts of interest in this paper.
[1] |
Kumar M, Grzelakowski M, Zilles J, et al. (2007) Highly permeable polymeri membranes based on the incorporation of the functional water channel protein Aquaporin Z. PNAS 104: 20719-20724. doi: 10.1073/pnas.0708762104
![]() |
[2] | Service RF (2006) Desalination freshens up. Science 313:1088-1090. |
[3] | Tal A (2006) Seeking sustainability: Israel's evolving water management strategy. Science 313: 1081-1084. |
[4] |
Discher BM, Won YY, Ege DS, et al. (1999) Polymersomes: tough vesicles made from diblock copolymers. Science 284:1143-1146. doi: 10.1126/science.284.5417.1143
![]() |
[5] |
Calamita G, Kempf B, Rudd KE, et al. (1997) The aquaporin-Z water channel gene of Escherichia coli: Structure, organization and phylogeny. Biology of the Cell 89: 321-329. doi: 10.1111/j.1768-322X.1997.tb01029.x
![]() |
[6] |
Calamita G, Bishai WR, Preston GM, et al. (1995) Molecular cloning and characterization of AqpZ, a water channel from Escherichia coli. J Biol Chem 270: 29063-29066. doi: 10.1074/jbc.270.49.29063
![]() |
[7] | Scheuring S, Ringler P, Borgnia M, et al. (1999) High resultion AFM topographs of the Escherichia coli water channel aquaporin Z. EMBO J 18: 4981-4987. |
[8] | Gorin MB, Yancey SB, Cline J, et al. (1984) The major intrinsic protein (MIP) of the bovine lens fiber membrane: Characterization and structure based on cDNA cloning. Cell 39: 49-59. |
[9] |
Ishibashi K, Kuwahara M, Gu Y, et al. (1997) Cloning and functional expression of a new water channel abundantly expressed in the testis permeable to water, glycerol and urea. J Biol Chem 272: 20782-20786. doi: 10.1074/jbc.272.33.20782
![]() |
[10] | Soupene E, King N, Lee H, et al. (2001) Aquaporin Z of Escherichia coli: Reassessment of Its Regulation and Physiological Role. J Bacter 184: 4304-4307. |
[11] |
Calamita G, Kempf B, Bonhivers B, et al. (1998) Regulation of the Escherichia coli water channel gene AqpZ. Proc Natl Acad Sci U S A 95: 3627-3631. doi: 10.1073/pnas.95.7.3627
![]() |
[12] | Borgnia MJ, Kozono D, Calamita G, et al. (1999) Funcation Reconstitution and Characterization of AqpZ, the E. coli Water Channel Protein. J Mol Biol 291: 1169-1179. |
[13] |
Kozono D, Yasui M, King LS, et al. (2002). Aquaporin water channels: atomic structure molecular dynamics meet clinical medicine. J Clin Inves 109: 1395-1399. doi: 10.1172/JCI0215851
![]() |
[14] |
Nemeth-Cahalan KL, Hall JE (2000) pH and Calcium Regulate the Water Permeability of Aquaporin 0. J Biol Chem 275: 6777-6782. doi: 10.1074/jbc.275.10.6777
![]() |
[15] | Cahalan K, Kalman K, Hall JE (2004) Molecular Basis of pH and Ca2+ Regulation of Aquaporin Water Permeability. J Gen Physiol 123: 573-580 |
[16] | Zhou W, Jones SW (1996) The effects of external pH on calcium channel currents in bullfrog sympathetic neurons. Biophys J 70: 1326-1334 |
[17] | Gonen T, Walz T (2006) The structure of aquaporins. Q Rev Biophys 39: 361-396. |
[18] |
Chaumont F, Moshelion F, Daniels MJ (2005) Regulation of plant aquaporin activity. Biol Cell 97: 749-764. doi: 10.1042/BC20040133
![]() |
[19] |
Tong J, Canty JT, Briggs MM, et al. (2013) The water permeability of lens aquaporin-0 depends on its lipid bilayer environment. Exp Eye Res 113: 32-40. doi: 10.1016/j.exer.2013.04.022
![]() |
[20] |
Andersen OS, Bruno MJ, Sun H, et al. (2007) Single-molecule methods for monitoring changes in bilayer elastic properties. Meth Mol Biol 400: 543-570 doi: 10.1007/978-1-59745-519-0_37
![]() |
[21] |
Hong H, Tamm LK (2004) Elastic coupling of integral membrane protein stability to lipid bilayer forces. Proc Natl Acad Sci U S A 101: 4065-4070. doi: 10.1073/pnas.0400358101
![]() |
[22] |
Nyholm TK, Ozdirekcan S, Killian JA (2007) How protein transmembrane segments sense the lipid environment. Biochemistry 46: 1457-1465. doi: 10.1021/bi061941c
![]() |
[23] |
Phillips R, Ursell T, Wiggins P, et al. (2009) Emerging roles for lipids in shaping membrane-protein function. Nature 459: 379-385. doi: 10.1038/nature08147
![]() |
[24] | Yuan C, O'Connell RJ, Jacob RF, et al. (2007) Regulation of the gating of BKCa channel by lipid bilayer thickness. J Biol Chem 282: 7276-7286. |
[25] |
Dumas F, Tocanne JF, Leblanc G, et al. (2000) Consequences of hydrophobic mismatch between lipids and melibiose permease on melibiose transport. Biochem 39: 4846-4854. doi: 10.1021/bi992634s
![]() |
[26] |
Perozo E, Kloda A, Cortes DM, et al. (2002) Physical principles underlying the transduction of bilayer deformation forces during mechano senditive channel gating. Nat Struct Biol 9: 696-703. doi: 10.1038/nsb827
![]() |
[27] |
Xie W, He F, Wang B, et al. (2013) An aquaporin-based vesicle-embedded polymeric membrane for low energy water filtration. J Mater Chem A 1: 7592-7600. doi: 10.1039/c3ta10731k
![]() |
[28] | Wang H, Chung TS, Tong YW, et al. (2011) Preparation and characterization of pore-suspending biomimetic membranes embedded with Aquaporin Z on carboxylated polyethylene glycol polymer cushion. Soft Matter 7: 7274-7280. |
[29] | Wang H, Chung TS, Tong YW, et al. (2012) Highly permeable and selective pore-spanning biomimetic membrane embedded with aquaporin Z. Small 8: 1185-1190, 1125. |
[30] | Duong PHH, Chung TS, Jeyaseelan K, et al. (2012) Planar biomimetic aquaporin-incorporated triblock copolymer membranes on porous alumina supports for nanofiltration. J Membr Sci 409: 34-43. |
[31] | Zhong PS, Chung TS, Jeyaseelan K, et al. (2012) Aquaporin-embedded biomimetic membranes for nanofiltration. J Membr Sci 407: 27-33. |
[32] | Savage DF, Egea PF, Colmenares YR, et al. (2013) Architecture and selectivity in aquaporins: 2.5A X-Ray Structure of Aquaporin Z. PLoS Biol 1: 334-340. |
[33] |
Nielsen CH (2009) Biomimetic membranes for sensor and separation applications. Bioanal Chem 395: 697-718. doi: 10.1007/s00216-009-2960-0
![]() |
[34] | Discher DE, Eisenberg A(2002) Polymer vesicles. Science 297: 967-973. |
[35] |
Ahmed F, Photos PJ, Discher DE (2006) Polymersomes as viral capsid mimics. Drug Develop Res 67: 4-14. doi: 10.1002/ddr.20062
![]() |
[36] |
Lewis BA, Engelman DM (1983) Lipid Bilayer Thickness Varies Linearly with Acyl Chain Length in Fluid Phosphatidylcholine Vesicles. J Mol Biol 166: 211-217. doi: 10.1016/S0022-2836(83)80007-2
![]() |
[37] |
Dave PC, Tiburu EK, Damodaran K, et al. (2004) Investigating Structural Changes in the Lipid Bilayer upon Insertion of the Transmembrane Domain of the Membrane-Bound Protein Phospholamban Utilizing 31P and 2H Solid-State NMR Spectroscopy. Biophys J 86: 1564-1573. doi: 10.1016/S0006-3495(04)74224-1
![]() |
[38] | Marsh D (2008) Energetics of Hydrophobic Matching in Lipid-Protein Interactions. Biophys J 94: 3996-4013. |
[39] |
Xu Q, Kim M, David Ho KW, et al. (2008) Membrane Hydrocarbon Thickness Modulates the Dynamics of a Membrane Transport Protein. Biophys J 95: 2849-2858. doi: 10.1529/biophysj.108.133629
![]() |
[40] |
He F, Tong YW (2014) A mechanistic study on amphiphilic block co-polymer poly(butadiene-b-(ethylene oxide)) vesicles reveals the water permeation mechanism through a polymeric bilayer. RSC Adv 4: 15304-15313. doi: 10.1039/c3ra48063a
![]() |
[41] |
Yang B, Verkman AS (1997) Water and Glycerol Permeabilities of Aquaporins 1-5 and MIP Determined Quantitatively by Expression of Epitope-tagged Constructs in Xenopus Oocytes. J Biol Chem 272: 16140-16146. doi: 10.1074/jbc.272.26.16140
![]() |
[42] |
Mehdizadeh H, Dickson JM, Eriksson PK (1989) Temperature effects on the performance of thin-film composite, aromatic polyamide membranes. Ind Eng Chem Res 28: 814-824. doi: 10.1021/ie00090a025
![]() |
1. | Agustín G. Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé, Hidde de Jong, Francis Mairet, Dynamical Analysis and Optimization of a Generalized Resource Allocation Model of Microbial Growth, 2022, 21, 1536-0040, 137, 10.1137/21M141097X | |
2. | Jean-Baptiste Caillau, Walid Djema, Jean-Luc Gouzé, Sofya Maslovskaya, Jean-Baptiste Pomet, Turnpike Property in Optimal Microbial Metabolite Production, 2022, 194, 0022-3239, 375, 10.1007/s10957-022-02023-0 | |
3. | Agustin Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouze, 2021, Hierarchical MPC applied to bacterial resource allocation and metabolite synthesis, 978-1-6654-3659-5, 667, 10.1109/CDC45484.2021.9683025 | |
4. | Nicolas Augier, Agustín G. Yabo, Time‐optimal control of piecewise affine bistable gene‐regulatory networks, 2022, 1049-8923, 10.1002/rnc.6012 | |
5. | Agustin Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouze, 2022, Optimal allocation of bacterial resources in fed-batch reactors, 978-3-9071-4407-7, 1466, 10.23919/ECC55457.2022.9838346 | |
6. | Haihui Cheng, Xinzhu Meng, Multistability and Hopf bifurcation analysis for a three-strategy evolutionary game with environmental feedback and delay, 2023, 620, 03784371, 128766, 10.1016/j.physa.2023.128766 | |
7. | Agustín G. Yabo, Mohab Safey El Din, Jean-Baptiste Caillau, Jean-Luc Gouzé, Stability analysis of a bacterial growth model through computer algebra, 2023, 12, 2102-5754, 175, 10.5802/msia.37 | |
8. | Agustín G. Yabo, Optimal control strategies in a generic class of bacterial growth models with multiple substrates, 2025, 171, 00051098, 111881, 10.1016/j.automatica.2024.111881 | |
9. | Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé, Optimal Bacterial Resource Allocation Strategies in Batch Processing, 2024, 84, 0036-1399, S567, 10.1137/22M1506328 | |
10. | Agustín G. Yabo, Predicting microbial cell composition and diauxic growth as optimal control strategies, 2023, 56, 24058963, 6217, 10.1016/j.ifacol.2023.10.745 | |
11. | Haowen Gong, Huijun Xiang, Yifei Wang, Huaijin Gao, Xinzhu Meng, Strategy evolution of a novel cooperative game model induced by reward feedback and a time delay, 2024, 9, 2473-6988, 33161, 10.3934/math.20241583 |
Parameter | Description | Unit | Value |
$ \beta $ | Inverse of the cytoplasmic density | $ L \, g^{-1} $ | $ 0.003 $ |
$ k_M $ | Rate constant of the metabolic macroreaction | $ h^{-1} $ | $ 4.32 $ |
$ K_M $ | Half-saturation constant of the metabolic macroreaction | $ g \, L^{-1} $ | $ 33.33 $ |
$ k_R $ | Rate constant of the protein synthesis reaction | $ h^{-1} $ | $ 3.6 $ |
$ K_R $ | Half-saturation constant of the protein synthesis reaction | $ g \, L^{-1} $ | $ 1 $ |
$ k_X $ | Rate constant of the metabolite synthesis reaction | $ h^{-1} $ | $ 0.5 $ |
$ K_X $ | Half-saturation constant of the metabolite synthesis reaction | $ g \, L^{-1} $ | $ 1 $ |
Parameter | Description | Unit | Value |
$ \beta $ | Inverse of the cytoplasmic density | $ L \, g^{-1} $ | $ 0.003 $ |
$ k_M $ | Rate constant of the metabolic macroreaction | $ h^{-1} $ | $ 4.32 $ |
$ K_M $ | Half-saturation constant of the metabolic macroreaction | $ g \, L^{-1} $ | $ 33.33 $ |
$ k_R $ | Rate constant of the protein synthesis reaction | $ h^{-1} $ | $ 3.6 $ |
$ K_R $ | Half-saturation constant of the protein synthesis reaction | $ g \, L^{-1} $ | $ 1 $ |
$ k_X $ | Rate constant of the metabolite synthesis reaction | $ h^{-1} $ | $ 0.5 $ |
$ K_X $ | Half-saturation constant of the metabolite synthesis reaction | $ g \, L^{-1} $ | $ 1 $ |