Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article Special Issues

Optimal bacterial resource allocation: metabolite production in continuous bioreactors

  • We show novel results addressing the problem of synthesizing a metabolite of interest in continuous bioreactors through resource allocation control. Our approach is based on a coarse-grained self-replicator dynamical model that accounts for microbial culture growth inside the bioreactor, and incorporates a synthetic growth switch that allows to externally modify the RNA polymerase concentration of the bacterial population, thus disrupting the natural process of allocation of available resources in bacteria. Further on, we study its asymptotic behavior using dynamical systems theory, and we provide conditions for the persistence of the bacterial population. We aim to maximize the synthesis of the metabolite of interest during a fixed interval of time in terms of the resource allocation decision. The latter is formulated as an Optimal Control Problem which is then explored using Pontryagin's Maximum Principle. 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 steadystate production objectives in continuous bioreactors: biomass production and metabolite production. For this last purpose, and in addition to the allocation parameter, we control the dilution rate of the bioreactor, and we analyze the results through a numerical approach. The resulting two-dimensional optimization problem is defined in terms of Michaelis-Menten kinetics, and takes into account the constraints for the existence of the equilibrium of interest.

    Citation: Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé. Optimal bacterial resource allocation: metabolite production in continuous bioreactors[J]. Mathematical Biosciences and Engineering, 2020, 17(6): 7074-7100. doi: 10.3934/mbe.2020364

    Related Papers:

    [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
  • We show novel results addressing the problem of synthesizing a metabolite of interest in continuous bioreactors through resource allocation control. Our approach is based on a coarse-grained self-replicator dynamical model that accounts for microbial culture growth inside the bioreactor, and incorporates a synthetic growth switch that allows to externally modify the RNA polymerase concentration of the bacterial population, thus disrupting the natural process of allocation of available resources in bacteria. Further on, we study its asymptotic behavior using dynamical systems theory, and we provide conditions for the persistence of the bacterial population. We aim to maximize the synthesis of the metabolite of interest during a fixed interval of time in terms of the resource allocation decision. The latter is formulated as an Optimal Control Problem which is then explored using Pontryagin's Maximum Principle. 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 steadystate production objectives in continuous bioreactors: biomass production and metabolite production. For this last purpose, and in addition to the allocation parameter, we control the dilution rate of the bioreactor, and we analyze the results through a numerical approach. The resulting two-dimensional optimization problem is defined in terms of Michaelis-Menten kinetics, and takes into account the constraints for the existence of the equilibrium of interest.


    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 Vext [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).

    Figure 1.  Extended coarse-grained self replicator model introduced in [10]. The substrate in the bioreactor S is consumed by the bacterial culture and transformed into precursors P through the action of the metabolic machinery M. Then, precursors are used to make macromolecules of the gene expression machinery R and the metabolic machinery M, with proportions α and 1α respectively; and to synthesize the metabolite X, which is excreted from the cell to the bioreactor. The external control I affects how the precursors P are distributed between both cellular functions M and R.

    The model represents three chemical macroreactions,

    SVMP
    PVRαR+(1α)M
    PVXX

    The first reaction is catalyzed by M and describes the transformation of external substrate S into precursor metabolites P at a rate VM. The second one represents the conversion of precursors into macromolecules R and M and is catalyzed by R, at a rate VR. Finally, the third reaction describes the transformation of precursors P into product X at a rate VX, and catalyzed by M. The natural resource allocation parameter is modeled through the dimensionless parameter α[0,1], that represents the proportion of precursors allocated to the gene expression machinery R, while 1α 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 h1] 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 α 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 α. 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=α). In this work, we limit the analysis to calculating the control input u, without decoupling both individual controls I and α (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=VSinVMVSout,˙P=VMVRVXVPout,˙M=(1u)VRVMout,˙R=uVRVRout,˙X=VXVXout.

    The inflow/outflow rates are defined as

    VSout=DS,VPout=DP,VMout=DM,VRout=DR,VXout=DX,VSin=Fsin,

    where sin [g L1] is the nutrient concentration of the inflow of fresh medium, and D [h1] the dilution rate defined as

    DFVext.

    We define the volume of the cell population V [L] as

    Vβ(M+R), (2.2)

    where β [L g1] 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 V as a simplifying assumption. Then, we express the quantities of the system as concentrations with respect to the volumes,

    pPV,rRV,mMV,sSVext,xXVext, (2.3)

    where p, r and m [g L1] are intracellular concentrations of precursor metabolites, components of the gene expression machinery and metabolic enzymes respectively; and s and x [g L1] the extracellular concentrations of substrate and metabolite. It is worth stressing that intracellular concentrations are defined with respect to the bacterial volume V, while extracellular concentrations with respect to the volume of the bioreactor Vext. Then, using definition (2.2), we obtain that m+r=1/β. We define the growth rate of the bacterial population μ [h1] as the relative variation of cell volume ˙V/V without considering the effect of the outflowing biomass. Replacing with concentrations leads to the system

    S:{˙s=D(sins)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=(1u)vR(p,r)μ(p,r)m,˙x=vX(p,m)VVextDx,˙V=(μ(p,r)D)V, (2.4)

    where vM(s,m), vR(p,r) and vX(p,m) [g L1 h1] are the mass fluxes per unit volume obtained from dividing the rates VM, VR and VX by 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 ˆsin=βsin. Then, dropping all hats yields the following system

    S1:{˙s=D(sins)vM(s,1r)V,˙p=vM(s,1r)vX(p,1r)μ(p,r)(p+1),˙r=(ur)μ(p,r),˙x=vX(p,1r)VDx,˙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 vM(s,m), vR(p,r) and vX(p,m):

    Assumption 1. Functions vM(s,m), vR(p,r) and vX(p,m) meet

    vi(y,z):R+×[0,1]R+

    vi(y,z) continuously differentiable w.r.t. both variables

    vi(0,z)=vi(y,0)=0

    vi() strictly monotonically increasing:

    viy(y,z)>0,(y,z)R>0×(0,1],viz(y,z)>0,(y,z)R>0×(0,1]

    vi() bounded w.r.t. y: limyvi(y,z)=vi,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(0,1), the metabolite synthesis rate vX(p,1r) can be expressed in terms of the macromolecule synthesis rate vR(p,r) (i.e., the growth rate)

    vX(p,1r)=c(r)vR(p,r),

    where c(r):(0,1)R+ is a positive continuously differentiable function.

    As previously described, the reaction vX(p,m) is catalyzed by m, and the reaction vR(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 S1 to

    S1:{˙s=D(sins)vM(s,1r)V,˙p=vM(s,1r)μ(p,r)(p+c(r)+1),˙r=(ur)μ(p,r),˙x=c(r)μ(p,r)VDx,˙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 vM, vR and vX (taking into account that, for these functions, the domain is defined as R+×[0,1β]R+).

    The asymptotic behavior of system S1 describes the “open-loop” operation mode of the continuous bioreactor, where the resource allocation control u(t) is fixed to ˉu(0,1). In the present section we propose a series of mass conservation laws that allow to reduce S1 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 S1 converges to the equilibria of its limiting system. Let us start the analysis of system S1 by defining its invariant region in the following lemma.

    Lemma 1. The set

    Γ={(s,p,r,x,V)R5:sins>0,p0,x0,1r0,V0}

    is positively invariant for the initial value problem.

    Proof. Let us analyze the boundaries of Γ,

    ˙s|s=0=Dsin>0s=0isrepulsive.˙s|s=sin=vM(s,1r)V0s=sinisrepulsive/invariant.˙p|p=0=vM(s,1r)0p=0isrepulsive/invariant.˙r|r=0=0r=0isinvariant.˙r|r=1=(u1)μ(p,1)<0r=1isrepulsive.˙x|x=0=vX(p,m)V0x=0isrepulsive/invariant.˙V|V=0=0V=0isinvariant.

    We will study the initial value problem of system S1 with initial conditions

    sins(0)0,p(0)0,1r(0)>0,x(0)0,V(0)>0, (2.9)

    where two cases have been excluded for being trivial to the analysis: 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 μ(p,0)=0, and therefore it is not possible to self-replicate from that point.

    System S1 can be rewritten as

    {˙φ=D(sinvinvout)+Nvivμμ(p,r),˙V=(μ(p,r)D)V,

    where

    φ[s,p,r,m,x]T is the state vector of concentrations in the bioreactor.

    N is the stoichiometry matrix of the macroreactions.

    vi is the vector of internal synthesis rates.

    vin and vout are the vectors of inflows and outflows respectively, associated to the continuous bioreactor setup.

    vμ is the vector modeling the dilution effect due to variation of the bacterial volume.

    which are defined as

    N[V001110u001u000V],vi[vM(s,1r)vR(p,r)vX(p,1r)],vin[1,0,0,0,0]T,voutdiag(φ)[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

    w1s+(p+m+r)V+x=s+(p+1)V+x,
    w2s+(p+rˉu)V+x.

    The first quantity tends asymptotically to w1=sin as t for every input u(t), as it obeys the dynamical equation

    ˙w1=D(sinw1). (2.10)

    Moreover, when fixing u(t) to ˉu(0,1), the quantity w2 also obeys the same Eq (2.10), meaning that this second quantity also converges to sin, which greatly simplifies the analysis of the asymptotic behavior of the system.

    Lemma 2. The ω-limit set of any solution of system S1 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)=ˉ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 S1 through its limiting system.

    Lemma 2 presents two mass conservation laws that can be used to reduce subsystem S1 by two dimensions. We will first analyze the asymptotic behavior of concentration r: when t, quantities w1=w2=sin, so

    s+(p+1)V+x=s+(p+rˉu)V+xr=ˉu

    meaning that, as t, r will converge to the value ˉu. We can also express x=sins(p+1)V, so that the limiting system of S1 becomes

    S1:{˙s=D(sins)ˉ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),ˉcc(ˉu).

    Details on the convergence of the limiting system S1 to the original one S1 will be addressed later in the article. In next section, we will fully describe the asymptotic behavior of S1.

    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 ˉf(p)>0,ˉf(p)>0,pΓ (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 ˉμ(pw)D:

    - The interior equilibrium Ei exists, is unique and locally stable.

    - The washout equilibrium Ew exists, is unique and locally unstable.

    If ˉμ(pw)<D:

    - The interior equilibrium Ei does not exist.

    - The washout equilibrium Ew exists, is unique and locally stable.

    In the above criterion, the equilibria are defined as follows:

    The interior equilibrium Ei(si,pi,Vi), with

    pi:{pR+:ˉμ(p)=D}, (2.11)
    si:{sR+:ˉvM(s)=ˉf(pi)}, (2.12)
    ViD(sinsi)ˉvM(si). (2.13)

    The washout equilibrium Ew(sin,pw,0), with

    pw:{pR+:ˉ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 pup(t) with dynamical equation

    ˙pupˉvM(sin)ˉf(pup). (2.15)

    In Eq (2.15), pup converges to the value pw satisfying Eq (2.14). Additionally, the vector field at p=pw is always negative (or null when s=sin)

    ˙p|p=pw=ˉvM(s)ˉf(pw)0

    meaning that p=pw is either repulsive or invariant, so a new invariant set ΓΓ can be defined,

    Γ{(s,p,V)R3:sins>0,pwp>0,V0}. (2.16)

    Equilibrium Ei The existence of the interior equilibrium Ei 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 ˉμ(p)=D, ˉμ(p) is strictly monotonically increasing w.r.t. p so, if inequality (2.17) is met, pi should be unique. Similarly, in ˉvM(s)=ˉf(pi), ˉvM(s) is again strictly monotonically increasing so, if inequality (2.18) is met, si 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 Γ, these conditions become

    ˉμ(pw)D, (2.19)

    where the second condition (2.18) is included in the inequality (2.19), as maxˉvM(s)=ˉvM(sin)=ˉf(pw), and the condition ˉf(pw)ˉf(pi) is true if and only if μ(pw)D. The Jacobian matrix is given by

    Ji=[DˉvM(s)Vi0ˉvM(s)ˉvM(s)ˉf(p)00ˉμ(p)Vi0],

    and the characteristic polynomial is

    Pi(λ)=(λ+D+ˉvM(s)Vi)(λ+ˉf(p))λ+ˉvM(s)ˉvM(s)ˉμ(p)Vi (2.20)
    =λ3+λ2(D+ˉvM(s)Vi+ˉf(p))+λ(D+ˉvM(s)Vi)ˉf(p)+ˉvM(s)ˉvM(s)ˉμ(p)Vi.

    Using Definition 2, it can be seen that D is an eigenvalue by replacing in expression (2.20)

    Pi(D)=DˉvM(s)Vi(p+ˉc+1)ˉμ(p)+ˉvM(s)ˉvM(s)ˉμ(p)Vi=ˉvM(s)Viˉμ(p)˙p=0.

    Then, dividing Pi(λ) by λ+D yields

    Pi(λ)1λ+D=λ2+λ(D+ˉvM(s)Vi+(p+ˉc+1)ˉμ(p))>0+ˉvM(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 Ew The washout equilibrium Ew exists for all values of Γ, since the only condition for existence is given by

    maxˉf(p)ˉvM(sin),

    which is always true in Γ since ˉvM(sin)=ˉf(pw) and so the inequality becomes ˉf(pw)ˉf(pw). Again, as ˉf(p) is strictly monotonically increasing, there is a unique solution pw. Its stability is given by

    Jw=[D0ˉcˉvM(s)ˉf(p)000ˉμ(pw)D]

    with characteristic polynomial

    Pw(λ)=(λ+D)(λ+ˉf(p))(λˉμ(pw)+D).

    Eigenvalues λ=D and λ=ˉ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 S1 can be further reduced using another mass conservation law given by the fact that ˉc is constant. This is formalized in the following lemma.

    Lemma 3. The ω-limit set of any solution of the limiting system S1 lies in the hyperplane

    Ω3{(s,p,V)R3:s+(p+ˉc+1)V=sin}

    Proof. The quantity

    w3s+(p+ˉc+1)V (2.21)

    obeys the dynamical equations ˙w3=D(sinw3), which means that the system converges asymptotically to the limit set Ω3.

    Using Lemma 2, we can further reduce system S1 by expressing

    s=sinV(p+ˉc+1), (2.22)

    for the values (p,V) that meet sinV(p+ˉ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.

    Figure 2.  Phase plane of the limiting system S_1'' showing: The case where \bar \mu(p_w) \geq D (left), so that the equilibrium E_i^* exists and attracts all solutions; and the case where \bar \mu(p_w) < D (right) and E_w^* attracts all solutions.

    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.

    Table 1.  Relevant biological system variables.
    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

     | Show Table
    DownLoad: CSV

    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.

    Figure 3.  Results of the numerical simulations using BOCOP [33], showing the optimal control input u and the state variables. The simulation has been done in a time horizon T = 80 , with 5000 time steps and Midpoint discretization method. The Ipopt (interior point nonlinear optimizer) arguments are: max_iter = 1000 , tol = 1.0e-14 , The initial conditions are fixed to: s(0) = 0.1 , p(0) = 0.024 , r(0) = 0.1 , x(0) = 0 , \mathcal{V}(0) = 0.2 .

    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.

    Figure 4.  Results of the numerical simulations using BOCOP, comparing the control u solution of the OCP to the solution of the static optimization problem u_{sp} (both with same initial conditions). The parameters for the numerical simulation match the ones used in Figure 3. The last plot emphasizes the area below the curves x_{sp} and x in OCP, as they are proportional to the total mass of metabolite produced, which is the quantity to be maximized.

    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)
    Figure 5.  Successive derivatives of H_1 obtained with BOCOP. The intervals where the functions vanish are highlighted in light red. All functions vanish along the singular arc except for H_{10001} , highlighted in green, which is negative as required by Kelley condition (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.

    Figure 6.  Numerical results for both static problems. The values for the objective functions J_ \mathcal{V} and J_X are represented through a qualitative colormap. The set \Theta of maximum growth rate delimits the region of existence of the interior equilibrium E_i . Additionally, curves \bar u_{opt}(D) show the optimal allocation \bar u in terms of the dilution rate D .

    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) .

    Figure 7.  Numerical results for both static problems.

    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] R. U. Ibarra, J. S. Edwards, B. O. Palsson, Escherichia coli k-12 undergoes adaptive evolution to achieve in silico predicted optimal growth, Nature, 420 (2002), 186.
    [2] E. Bosdriesz, D. Molenaar, B. Teusink, F. J. Bruggeman, How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization, FEBS J., 282 (2015), 2029-2044. doi: 10.1111/febs.13258
    [3] J. D. Van Elsas, A. V. Semenov, R. Costa, J. T. Trevors, Survival of escherichia coli in the environment: fundamental and public health aspects, ISME J., 5 (2011), 173-183. doi: 10.1038/ismej.2010.80
    [4] N. Giordano, F. Mairet, J.-L. Gouzé, J. Geiselmann, H. De Jong, Dynamical allocation of cellular resources as an optimal control problem: novel insights into microbial growth strategies, PLoS Comput. Biol., 12 (2016), e1004802.
    [5] H. De Jong, S. Casagranda, N. Giordano, E. Cinquemani, D. Ropers, J. Geiselmann, et al., Mathematical modelling of microbes: metabolism, gene expression and growth, J. R. Soc. Interface, 14 (2017), 20170502.
    [6] A. L. Koch, Why can't a cell grow infinitely fast?, Can. J. Microbiol., 34 (1988), 421-426.
    [7] J. Shu, M. Shuler, A mathematical model for the growth of a single cell of E. coli on a glucose/glutamine/ammonium medium, Biotechnol. Bioeng., 33 (1989), 1117-1126. doi: 10.1002/bit.260330907
    [8] D. V. Goeddel, D. G. Kleid, F. Bolivar, H. L. Heyneker, D. G. Yansura, R. Crea, et al., Expression in escherichia coli of chemically synthesized genes for human insulin, Proc. Natl. Acad. Sci., 76 (1979), 106-110.
    [9] L. Huo, J. J. Hug, C. Fu, X. Bian, Y. Zhang, R. Müller, Heterologous expression of bacterial natural product biosynthetic pathways, Nat. Product Rep., 36 (2019), 1412-1436. doi: 10.1039/C8NP00091C
    [10] I. Yegorov, F. Mairet, H. De Jong, J.-L. Gouzé, Optimal control of bacterial growth for the maximization of metabolite production, J. Math. Biol., 78 (2019), 985-1032.
    [11] E. Cinquemani, F. Mairet, I. Yegorov, H. de Jong, J.-L. Gouzé, Optimal control of bacterial growth for metabolite production: the role of timing and costs of control, in 2019 18th European Control Conference (ECC), IEEE, 2019, 2657-2662.
    [12] A. G. Yabo, J.-B. Caillau, J.-L. Gouzé, Singular regimes for the maximization of metabolite production, in 2019 IEEE 58th Conference on Decision and Control (CDC), IEEE, 2019, 31-36.
    [13] J. Izard, C. D. G. Balderas, D. Ropers, S. Lacour, X. Song, Y. Yang, et al., A synthetic growth switch based on controlled expression of rna polymerase, Mol. Syst. Biol., 11 (2015), 840.
    [14] I. Otero-Muras, A. A. Mannan, J. R. Banga, D. A. Oyarzún, Multiobjective optimization of gene circuits for metabolic engineering, IFAC-PapersOnLine, 52 (2019), 13-16.
    [15] D. Molenaar, R. Van Berlo, D. De Ridder, B. Teusink, Shifts in growth strategies reflect tradeoffs in cellular economics, Mol. Syst. Biol., 5 (2009), 323.
    [16] M. T. Wortel, E. Bosdriesz, B. Teusink, F. J. Bruggeman, Evolutionary pressures on microbial metabolic strategies in the chemostat, Sci. Rep., 6 (2016), 29503.
    [17] D. W. Spitzer, Maximization of steady-state bacterial production in a chemostat with ph and substrate control, Biotechnol. Bioeng., 18 (1976), 167-178. doi: 10.1002/bit.260180203
    [18] R. R. Lichtl, M. J. Bazin, D. O. Hall, The biotechnology of hydrogen production by nostoc flagelliforme grown under chemostat conditions, Appl. Microbiol. Biotechnol., 47 (1997), 701-707. doi: 10.1007/s002530050998
    [19] M. C. D'anjou, A. J. Daugulis, A rational approach to improving productivity in recombinant pichia pastoris fermentation, Biotechnol. Bioeng., 72 (2001), 1-11. doi: 10.1002/1097-0290(20010105)72:1<1::AID-BIT1>3.0.CO;2-T
    [20] T. Bayen, F. Mairet, Optimization of the separation of two species in a chemostat, Automatica, 50 (2014), 1243-1248. doi: 10.1016/j.automatica.2014.02.024
    [21] C. Martínez, O. Bernard, F. Mairet, Maximizing microalgae productivity in a light-limited chemostat, IFAC-PapersOnLine, 51 (2018), 735-740. doi: 10.1016/j.ifacol.2018.05.070
    [22] H. R. Thieme, Convergence results and a PoincaréBendixson trichotomy for asymptotically autonomous differential equations, J. Math. Biol., 30 (1992), 755-763.
    [23] L. S. Pontryagin, Mathematical Theory of Optimal Processes, Routledge, 2018.
    [24] H. L. Smith, P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, 1995.
    [25] H. Bremer, P. P. Dennis, Modulation of chemical composition and other parameters of the cell by growth rate, in Escherichia coli and Salmonella: Cellular and Molecular Biology, ASM Press, Washington, (1996), 1553-1569.
    [26] S. Goutelle, M. Maurin, F. Rougier, X. Barbaut, L. Bourguignon, M. Ducher, et al., The hill equation: a review of its capabilities in pharmacological modelling, Fundam. Clin. Pharmacol., 22 (2008), 633-648.
    [27] K. A. Johnson, R. S. Goody, The original Michaelis constant: translation of the 1913 Michaelis-Menten paper, Biochemistry, 50 (2011), 8264-8269. doi: 10.1021/bi201284u
    [28] A. A. Agrachev, Y. Sachkov, Control Theory from the Geometric Viewpoint, Springer Science & Business Media, 2013.
    [29] M. I. Zelikin, V. F. Borisov, Theory of Chattering Control: With Applications to Astronautics, Robotics, Economics, and Engineering, Springer Science & Business Media, 2012.
    [30] A. Van der Schaft, Symmetries in optimal control, SIAM J. Control Optim., 25 (1987), 245-259.
    [31] E. Trélat, E. Zuazua, The turnpike property in finite-dimensional nonlinear optimal control, J. Differ. Equations, 258 (2015), 81-114. doi: 10.1016/j.jde.2014.09.005
    [32] W. Djema, L. Giraldi, S. Maslovskaya, O. Bernard, Turnpike features in optimal selection of species represented by quota models, Submitted.
    [33] I. S. Team Commands, Bocop: an open source toolbox for optimal control, 2017. Available from: http://bocop.org.
    [34] H. Robbins, A generalized Legendre-Clebsch condition for the singular cases of optimal control, IBM J. Res. Dev., 11 (1967), 361-372. doi: 10.1147/rd.114.0361
    [35] M. Scott, C. W. Gunderson, E. M. Mateescu, Z. Zhang, T. Hwa, Interdependence of cell growth and gene expression: origins and consequences, Science, 330 (2010), 1099-1102. doi: 10.1126/science.1192588
    [36] I. Yegorov, F. Mairet, J.-L. Gouzé, Optimal feedback strategies for bacterial growth with degradation, recycling, and effect of temperature, Optim. Control Appl. Methods, 39 (2018), 1084-1109.
  • This article has been cited by:

    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
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4355) PDF downloads(151) Cited by(11)

Figures and Tables

Figures(7)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog