Time-Optimal Adaptation in Metabolic Network Models

Analysis of metabolic models using constraint-based optimization has emerged as an important computational technique to elucidate and eventually predict cellular metabolism and growth. In this work, we introduce time-optimal adaptation (TOA), a new constraint-based modeling approach that allows us to evaluate the fastest possible adaptation to a pre-defined cellular state while fulfilling a given set of dynamic and static constraints. TOA falls into the mathematical problem class of time-optimal control problems, and, in its general form, can be broadly applied and thereby extends most existing constraint-based modeling frameworks. Specifically, we introduce a general mathematical framework that captures many existing constraint-based methods and define TOA within this framework. We then exemplify TOA using a coarse-grained self-replicator model and demonstrate that TOA allows us to explain several well-known experimental phenomena that are difficult to explore using existing constraint-based analysis methods. We show that TOA predicts accumulation of storage compounds in constant environments, as well as overshoot uptake metabolism after periods of nutrient scarcity. TOA shows that organisms with internal temporal degrees of freedom, such as storage, can in most environments outperform organisms with a static intracellular composition. Furthermore, TOA reveals that organisms adapted to better growth conditions than present in the environment (“optimists”) typically outperform organisms adapted to poorer growth conditions (“pessimists”).

1 Introduction

Over the past decades, various modeling frameworks have been proposed to understand the organization and functioning of cellular metabolism and growth. Among the most popular approaches are constraint-based methods, in particular flux balance analysis (FBA) (Orth et al., 2010). Constraint-based methods typically make use of optimality principles that are motivated by evolutionary arguments. That is, instead of requiring a detailed mechanistic understanding of the underlying regulatory machinery, properties of cellular metabolism, such as exchange fluxes or biomass accumulation, are predicted based on the assumption that metabolism has evolved according to certain evolutionary optimality principles.

More recently, constraint-based methods have been extended to quantitatively account for the synthesis costs of the biological macromolecules that are required for cellular metabolism and growth, giving rise to resource balance analysis (RBA) (Goelzer et al., 2011) and integrated reconstructions of Metabolism and macromolecular Expression (ME) (Lerman et al., 2012). While the initial approaches were restricted to time-invariant environments and subject to steady-state conditions, various dynamic extensions have also been proposed, such as dynamic FBA (dFBA) Mahadevan et al. (2002), dynamic enzyme-cost FBA (deFBA) (Waldherr et al., 2015), conditional FBA (cFBA) (Rügen et al., 2015; Reimers et al., 2017), dynamic RBA (dRBA) (Jeanne et al., 2018), dynamic ME (Yang et al., 2019), and regulatory dynamic enzyme-cost FBA (r-deFBA) (Liu and Bockmayr, 2020). These dynamic frameworks are computationally more expensive and allow predicting time courses over a given time interval, such that the variables fulfil a given (linear) optimality principle. Typically, within these frameworks, the time intervals over which the solutions are considered are predefined.

In this work, we extend these existing approaches and propose time-adaptation (TOA) as a new constraint-based modeling framework that allows us to evaluate the fastest possible adaptation to a pre-defined cellular state while fulfilling a given set of dynamic and static constraints. If the underlying dynamics of the biological system are governed by ordinary differential equations (ODEs) subject to algebraic constraints such as positivity, that is, so-called differential-algebraic equations (DAEs), time-optimal adaptation falls into the mathematical problem class of time-optimal control problems, which are optimal control problems where the time-interval is part of the objective (Hermes and Lasalle, 1969). In its general form, TOA can be applied in a very broad sense and thereby extends most of the existing constraint-based modeling frameworks.

Our approach allows us to compute feasible time courses to simulate or predict adaptations of cellular metabolism to environmental shifts. Potential applications include an analysis of cellular doubling, i.e., to analyze the optimal metabolic trajectory that results in a doubling of all cellular components in the shortest time, as well as an analysis of the temporal adaptation to changing nutrient availability.

We exemplify TOA using a coarse-grained self-replicator model (Molenaar et al., 2009; Giordano et al., 2016; Yegorov et al., 2018; Yabo et al., 2022) and demonstrate that TOA allows us to explain several known experimental phenomena that are difficult to investigate using existing static or dynamic constraint-based analysis methods. In particular, we demonstrate that TOA can explain the accumulation of storage compounds also in time-invariant environments–a counterintuitive fact that cannot be predicted using RBA and related methods. Likewise, we demonstrate that “luxury uptake” of nutrients, i.e., the fact that microorganisms may take up more of a limiting resource than strictly required for steady-state growth, can be explained by TOA and does not necessarily require competition within a microbial community. Furthermore, our analysis shows that organisms with internal temporal degrees of freedom, such as storage, can in most environments outperform organisms with a static intracellular composition. Finally, TOA shows that in constant (or slowly changing) environments, organisms adapted to better growth conditions (“optimists”) outperfom organisms adapted to poorer growth conditions (“pessimists”) when placed in the same environment.

The manuscript is organized as follows: Within Sections 2.1 and 2.2 we introduce notation and define a general constraint-based framework to describe cellular metabolism and growth. This framework captures most current examples of dynamic constraint-based modeling, in particular dynamic FBA (Mahadevan et al., 2002), dynamic enzyme-cost FBA (Waldherr et al., 2015) and conditional FBA (Reimers et al., 2017). In Section 2.3, we formally introduce time-optimal adaptation (TOA) and discuss two relevant applications in Section 2.4: cell doubling in minimal time, as well as transition after a nutrient shift. The latter is formulated as a two-objective optimization problem (in the sense of Pareto) that considers a minimal time for the transition versus a total increase in biomass. In Sections 2.5–2.7, we discuss numerical aspects, variability analysis, and implementation, respectively.

Readers not interested in the mathematical details may skip most of Materials and Methods and focus on Results. In Sections 3.1 and 3.2, we describe the coarse-grained self-replicator model and its properties using RBA. In Section 3.3, we then apply TOA to describe cell doubling in minimal time in a constant environment. In Section 3.4, we discuss the role of “expectation”, i.e., the consequences of being mis-adapted to a given environment. In Section 3.5, we apply TOA to simulate the metabolic response after a nutrient shift. In the final Sections 4 and 5, we discuss the biological implications of our results, and provide conclusions.

2 Materials and Methods

2.1 Introduction and Notation

The dynamic simulation of metabolic networks by means of a fully parameterized ODE/DAE model is an ideal scenario that, in most cases, cannot be met due to the inherent incompleteness and uncertainty of the description and the involved parameters. Constraint-based modeling (Bordbar et al., 2014) has therefore become an important paradigm for the computational description of cellular metabolism and growth. The general idea can be framed as follows: instead of making use of a fully mechanistic description of biochemical dependencies by means of reaction rate equations, the system is characterized by a set of constraints/inclusions, typically defined by (in-)equalities that constrain the dynamics over a time interval [t₀, t_end] of interest.

Before capturing our approach in mathematical terms in Section 2.2, we introduce some notation, see also Supplementary Appendix S1. The function $\mathit{y}:[t_0,t_{\text{end}}]\to {\mathbb{R}}^{n_{\mathit{y}}}$ is used to describe the cellular dynamics by the total amounts y(t) of intracellular compounds at time t (typically measured in number of molecules, mol), with $\dot{\mathit{y}}(t)=\frac{\mathrm{d}}{\mathrm{d}t}\mathit{y}(t)$ denoting the time-derivative. For simplicity, we focus on the dynamics of intracellular compounds only, extracellular compounds (e.g., nutrient or waste product concentrations) are not included in y. Our framework, however, can be readily adapted to include the dynamics of extracellular compounds (see the Supplementary Appendix S2.3 for details). Furthermore, our description is based on the assumption of a well-stirred metabolism, i.e., the spatial distribution of compounds is not considered.

We distinguish the total amounts of molecules y(t) from their concentrations c(t), defined by

$\mathit{c}\left(t\right)=\frac{\mathit{y}\left(t\right)}{bio\left(t\right)},$(2.1)

where the term bio(t)≔w^⊤⋅y(t) denotes the total biomass of the system (measured in Gram cellular dry mass). The vector $\mathit{w}\in {\mathbb{R}}^{n_{\mathit{y}}}$ denoting the molar masses of the entities of y (measured in gram cellular dry mass per mol).

The time evolution of the state vector y(t) can be described by means of ordinary differential equations.

$\dot{\mathit{y}}\left(t\right)=\mathbf{S}\cdot \mathit{v}\left(t\right),$(2.2a)

where $\mathbf{S}\in {\mathbb{R}}^{n_{\mathit{y}}\times n_{\mathit{v}}}$ denotes the stoichiometric matrix and $\mathit{v}:[t_0,t_{\text{end}}]\to {\mathbb{R}}^{n_{\mathit{v}}}$ the flux rates of the reactions. The flux rates v(t) may in general also depend on the environment the cells are exposed to. Typically, and specifically for large networks, the stoichiometric matrix $\mathbf{S}=\left(\begin{array}{c}\hfill {\mathbf{S}}_{\mathit{y}}\hfill \\ \hfill {\mathbf{S}}_{\mathit{x}}\hfill \end{array}\right)$ is split up such that “fast” and “slow” intracellular compounds, usually metabolites resp. macromolecules, are described separately and (2.2a) is replaced by.

$\begin{array}{ll}\hfill \dot{\mathit{y}}\left(t\right)& ={\mathbf{S}}_{\mathit{y}}\cdot \mathit{v}\left(t\right),\hfill \\ \hfill \mathbf{0}& ={\mathbf{S}}_{\mathit{x}}\cdot \mathit{v}\left(t\right),\hfill \end{array}$(2.2b)

where the fast compounds, corresponding to the rows of S_x, are subject to a quasi steady-state approximation (QSSA) (Segel and Slemrod, 1989). In this case, for simplicity of notation, the fast components will be removed from the vector y(t). We note that the splitting into “slow” and “fast” compounds is not a necessary step and its validity has to be verified in any particular application.

2.2 Constraint-Based Modeling

To capture the broad range of simulation frameworks that time-optimal adaptation is able to cover, we abstractly denote the constraints defining the specific constraint-based description of a cell via.

$\text{for}\left(\mathrm{a}\mathrm{l}\mathrm{m}\mathrm{o}\mathrm{s}\mathrm{t}\right)\mathrm{a}\mathrm{l}\mathrm{l}t:\left(\dot{\mathit{y}}\left(t\right),\mathit{y}\left(t\right),\mathit{u}\left(t\right)\right)\in \mathcal{A}\left(t\right),$(2.3a)

where the set $\mathcal{A}(t)\subseteq {\mathbb{R}}^{n_{\mathit{y}}}\times {\mathbb{R}}^{n_{\mathit{y}}}\times {\mathbb{R}}^{n_{\mathit{u}}}$ is typically defined through (in-)equalities such as steady-state assumptions and/or positivity requirements. The particular form of the set $\mathcal{A}(t)$ usually depends on the chosen modeling framework and its granularity. For the present work, we model the influence of the external conditions via the explicit time-dependence of $\mathcal{A}(t)$. The vector-valued function $\mathit{u}=\mathit{u}(t):[t_0,t_{\text{end}}]\to {\mathbb{R}}^{n_{\mathit{u}}}$ signifies the degrees-of-freedom of the cell, i.e., quantities that are not uniquely determined from the current state of the cell and its environment. In the context of control theory, u(t) defines the controls; on the biochemical level, it can for example stand for flux rates v(t) but also for parameters within the model.

The formal statement (2.3a) is usually not enough to sufficiently constrain the solutions, because the feasible region is too large to obtain biochemical insight. To get biochemically meaningful results, (2.3a) is therefore often accompanied by boundary conditions and an optimality principle, i.e., a global objective function f to be optimized:

${\mathit{\varphi }}_{\text{bndry}}\left(\mathit{y}\left(t_0\right),\mathit{y}\left(t_{\text{end}}\right),\mathit{u}\left(t_0\right),\mathit{u}\left(t_{\text{end}}\right)\right)\le \mathit{0},$(2.3b)

$\underset{\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right)}{\min }f\left(\mathit{y},\mathit{u}\right)$(2.3c)

The boundary conditions (2.3b) are defined by means of inequalities to allow for more generality of this description. Usually, the boundary conditions will only contain initial values, provided by equality constraints, i.e., two inequalities. In some cases, optimality principles are already incorporated into the constraint set $\mathcal{A}(t)$, see the following examples.

In the context of optimal control-based methods with ODE/DAE constraints, the flux rates at any fixed point in time cannot (mathematically) be determined as they enter the problem as control variables (Gerdts, 2011). This is why (2.3a) technically can only be enforced for almost all times. Numerically or with respect to the biochemical reasoning, however, this has no further implications. In the following, we illustrate how (2.3) provides an abstract framework to describe established examples of constraint-based modeling.

Example 2.1. (Dynamic FBA, dFBA). Dynamic (or iterative) flux balance analysis (Varma and Palsson, 1994; Mahadevan et al., 2002), although one of the most commonly used dynamic frameworks within constraint-based modeling, is not consistently defined in the literature. Here, we refer to the formulation in (Höffner et al., 2016), see also (Höffner et al., 2012), for the characterization of dynamic FBA as a “dynamical system with a linear program embedded.”The control quantities u(t) can in this case be directly identified with the flux rates in the metabolic network model, i.e., v(t) = u(t). The overall dynamics are governed by (2.2a), positivity requirements on y(t) and flux bounds $\mathit{l}\mathit{b},\mathit{u}\mathit{b}\in {\mathbb{R}}^{n_{\mathit{v}}}$, which might be dependent on the time t:

$\begin{array}{lll}\hfill \dot{\mathit{y}}\left(t\right)& ={\mathbf{S}}_{\mathit{y}}\cdot \mathit{v}\left(t\right),\hfill & \hfill \left(\text{dynamics, often just biomass}\right)\\ \hfill \mathbf{0}& ={\mathbf{S}}_{\mathit{x}}\cdot \mathit{v}\left(t\right),\hfill & \hfill \left(\text{quasi steady}-\text{state}\right)\\ \hfill \mathbf{0}& \le \mathit{y}\left(t\right),\hfill & \hfill \left(\text{positivity}\right)\\ \hfill \mathit{l}\mathit{b}\left(t\right)& \le \mathit{v}\left(t\right)\le \mathit{u}\mathit{b}\left(t\right),\hfill & \hfill \left(\text{flux bounds, dependent on environment}\right)\end{array}$

with given initial conditions

$\begin{array}{ll}\hfill \mathit{y}\left(t_0\right)& ={\mathit{y}}_0\in {\mathbb{R}}^{n_{\mathit{y}}}.\hfill \end{array}$

The flux rates are determined through optimization of a linear functional (often the flux through the biomass reaction, assembled in a vector ${\mathit{w}}_{\text{obj}}\in {\mathbb{R}}^{n_{\mathit{v}}}$)

$\begin{array}{ll}\hfill \mathit{v}\left(t\right)& \in \underset{\mathit{v}}{\mathrm{argmin}}{\mathit{w}}_{\text{obj}}^{\top }\cdot \mathit{v}.\hfill \end{array}$

The quantities in (2.3) can be identified as:

$\begin{array}{ll}\hfill \mathcal{A}\left(t\right)& =\left\{\left(\dot{\mathit{y}},\mathit{y},\mathit{u}\right):\begin{array}{cc}\hfill \dot{\mathit{y}}& \hfill ={\mathbf{S}}_{\mathit{y}}\cdot \mathit{u},\hfill \\ \hfill \mathbf{0}& \hfill ={\mathbf{S}}_{\mathit{x}}\cdot \mathit{u},\hfill \\ \hfill \mathbf{0}& \le \mathit{y},\hfill \\ \hfill \mathit{l}\mathit{b}\left(t\right)& \hfill \le \mathit{u}\le \mathit{u}\mathit{b}\left(t\right),\hfill \\ \hfill \mathit{u}& \in \hfill & {\mathrm{argmin}}_{\mathit{v}}{\mathit{w}}_{\text{obj}}^{\top }\cdot \mathit{v}\hfill \end{array}\right\},\hfill \\ \hfill {\mathit{\varphi }}_{\text{bndry}}& =\left(\begin{array}{c}\hfill \mathit{y}\left(t_0\right)-{\mathit{y}}_0\hfill \\ \hfill {\mathit{y}}_0-\mathit{y}\left(t_0\right)\hfill \end{array}\right),\hfill \end{array}$

while typically no additional (global) objective function is present. Note that the defining condition on the fluxes $\mathit{u}\in \underset{\mathit{v}}{\mathrm{argmin}}{\mathit{w}}_{\text{obj}}^{\top }\cdot \mathit{v}$ is an inclusion, such that the solutions to dynamic FBA problems are, in general, not unique. To remedy this, flux variability analysis (FVA) (Mahadevan and Schilling, 2003) was introduced as a computational tool to explore the range of possible solutions of the static sub-problems.

Example 2.2. (Dynamic enzyme-cost FBA, deFBA). Dynamic enzyme-cost FBA (Waldherr et al., 2015) is a dynamic extension of FBA that takes into account the temporal development and function of the enzymes. This is modeled by a system of linear inequalities

${\mathbf{H}}_{\mathit{y}}\left(t\right)\cdot \mathit{y}\left(t\right)+{\mathbf{H}}_{\mathit{v}}\left(t\right)\cdot \mathit{v}\left(t\right)\le \mathit{h}\left(t\right),$(2.4)

with

${\mathbf{H}}_{\mathit{y}}\in {\mathbb{R}}^{n_{\mathit{h}}\times n_{\mathit{y}}},{\mathbf{H}}_{\mathit{v}}\in {\mathbb{R}}^{n_{\mathit{h}}\times n_{\mathit{v}}}.$(2.5)

The model is usually formulated as an initial-value problem

$\mathit{y}\left(t_0\right)={\mathit{y}}_0\in {\mathbb{R}}^{n_{\mathit{y}}}.$

Similar to FBA, deFBA assumes that a certain objective function is to be optimized. Since the framework entails a fully dynamic model over the whole time range of interest, the objective function contains “global” information, expressed as an optimal control objective of Boltza-type (Gerdts, 2011),

$\underset{\mathit{y},\mathit{u}}{\min }{\int }_{t_0}^{t\mathrm{e}\mathrm{n}\mathrm{d}}{\mathit{\varphi }}_{\mathit{y}}^{\top }\left(t\right)\cdot \mathit{y}\left(t\right)dt+{\mathit{\varphi }}_{\text{end}}^{\top }\cdot \mathit{y}\left(t_{\text{end}}\right),\text{\hspace{0.17em}with\hspace{0.17em}}{\mathit{\varphi }}_{\mathit{y}},{\mathit{\varphi }}_{\text{end}}\in {\mathbb{R}}^{n_{\mathit{y}}}.$

Like in dFBA, the control variables in deFBA can be identified with the flux rates and the description in terms of (2.3) is given by

$\begin{array}{ll}\hfill \mathit{u}\left(t\right)& =\mathit{v}\left(t\right),\hfill \\ \hfill \mathcal{A}\left(t\right)& =\left\{\left(\dot{\mathit{y}},\mathit{y},\mathit{u}\right):\begin{array}{cc}\hfill \dot{\mathit{y}}& \hfill ={\mathbf{S}}_{\mathit{y}}\cdot \mathit{u},\hfill \\ \hfill \mathbf{0}& \hfill ={\mathbf{S}}_{\mathit{x}}\cdot \mathit{u},\hfill \\ \hfill \mathbf{0}& \hfill \le \mathit{y}\hfill \\ \hfill \mathit{l}\mathit{b}\left(t\right)& \hfill \le \mathit{u}\le \mathit{u}\mathit{b}\left(t\right),\hfill \\ \hfill \mathit{h}\left(t\right)& \hfill \ge {\mathbf{H}}_{\mathit{y}}\left(t\right)\cdot \mathit{y}+{\mathbf{H}}_{\mathit{u}}\left(t\right)\cdot \mathit{u}\hfill \end{array}\right\}\hfill \\ \hfill {\mathit{\varphi }}_{\text{bndry}}& =\left(\begin{array}{c}\hfill \mathit{y}\left(t_0\right)-{\mathit{y}}_0\hfill \\ \hfill {\mathit{y}}_0-\mathit{y}\left(t_0\right)\hfill \end{array}\right),\hfill \\ \hfill f\left(\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right)\right)& ={\int }_{t_0}^{t_{\text{end}}}{\mathit{\varphi }}_{\mathit{y}}^{\top }\left(t\right)\cdot \mathit{y}\left(t\right)dt+{\mathit{\varphi }}_{\text{end}}^{\top }\cdot \mathit{y}\left(t_{\text{end}}\right).\hfill \end{array}$

Example 2.3. (Conditional FBA, cFBA). This framework (Rügen et al., 2015; Reimers et al., 2017) is again a dynamic extension of resource balance analysis (RBA) (Goelzer et al., 2011). Like in deFBA, enzymatic constraints (potentially alongside further constraints, e.g., on the cell’s density) are included via (2.4). The boundary values in cFBA, however, are defined through a periodicity condition that accounts for the growth of the cell:

$\mathit{c}\left(t_0\right)=\frac{1}{bio\left(t_0\right)}\cdot \mathit{y}\left(t_0\right)=\frac{1}{bio\left(t_{\text{end}}\right)}\cdot \mathit{y}\left(t_{\text{end}}\right)=\mathit{c}\left(t_{\text{end}}\right).$(2.6)

Instead of using the biomass production on all time points, the objective in cFBA is the total growth of the cell until t_end. In terms of (2.3), cFBA can be summarized as

$\begin{array}{ll}\hfill \mathit{u}\left(t\right)& =\left(\begin{array}{c}\hfill {\mathit{u}}_1\hfill \\ \hfill {\mathit{u}}_{2:}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill \alpha \hfill \\ \hfill \mathit{v}\left(t\right)\hfill \end{array}\right),\alpha \in \mathbb{R},\hfill \\ \hfill \mathcal{A}\left(t\right)& =\left\{\left(\dot{\mathit{y}},\mathit{y},\mathit{u}\right):\begin{array}{cc}\hfill \dot{\mathit{y}}& \hfill ={\mathbf{S}}_{\mathit{y}}\cdot {\mathit{u}}_{2:},\hfill \\ \hfill \mathbf{0}& \hfill ={\mathbf{S}}_{\mathit{x}}\cdot {\mathit{u}}_{2:},\hfill \\ \hfill \mathbf{0}& \hfill \le \mathit{y},\hfill \\ \hfill \mathit{l}\mathit{b}\left(t\right)& \hfill \le {\mathit{u}}_{2:}\le \mathit{u}\mathit{b}\left(t\right),\hfill \\ \hfill \mathit{h}\left(t\right)& \hfill \ge {\mathbf{H}}_{\mathit{y}}\left(t\right)\cdot \mathit{y}+{\mathbf{H}}_{\mathit{v}}\left(t\right)\cdot {\mathit{u}}_{2:}\hfill \end{array}\right\},\hfill \\ \hfill {\mathit{\varphi }}_{\text{bndry}}& =\left(\begin{array}{c}\hfill {\mathit{u}}_1\cdot \mathit{y}\left(t_0\right)-\mathit{y}\left(t_{\text{end}}\right)\hfill \\ \hfill \mathit{y}\left(t_{\text{end}}\right)-{\mathit{u}}_1\cdot \mathit{y}\left(t_0\right)\hfill \end{array}\right),\hfill \\ \hfill f\left(\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right)\right)& =-{\mathit{u}}_1,\hfill \end{array}$

where u₁ refers to the first component of the vector u and u_2: to the vector of the remaining entries. If no constraints on the cell density are included in (2.4), the inequalities defining cFBA are often scale-invariant in the sense that for each solution y(t) and each number β ≥ 0, the function β ⋅y(t) is also a solution. To exclude trivial solutions, the boundary conditions are therefore often extended such that the biomass at t₀ is equal to one. Note that cFBA is inherently nonlinear as the products u₁⋅y in the boundary value constraints contribute quadratically in the unknowns y and u. Like in RBA, the numerical solution of cFBA problems therefore comprises a series of linear programs that have to be solved after a discretization of the dynamics by means of, for example, a collocation scheme.

Example 2.4. (Iterative RBA, (Liu, 2020), see also dynamic ME (Yang et al., 2019)). Just as dynamic FBA can be seen as a dynamic extension of classical FBA by iteratively applying the algorithm with constraints following the external conditions, resource balance analysis (RBA, see Goelzer et al. (2011)) can also be applied consecutively. In doing that, the limit case of infinitesimally short sub-intervals leads to a fully dynamic framework. Numerically, this limiting process is skipped and one only solves RBA problems on a series of short—but finite—time intervals. Note that, as cFBA, RBA uses periodicity conditions like (2.6) which implies that, in constant external conditions, only one RBA problem needs to be solved. The full solution in this case is given by an exponential curve for y(t). Note that there are fewer degrees-of-freedom for the cell when compared to deFBA or cFBA, as the fixed concentration values for the metabolites in the case of iterative RBA also block internal dynamics of the metabolic network.In the notation of the constraint-based framework (2.3), iterative RBA can be written as

$\begin{array}{ll}\hfill \mathit{u}\left(t\right)& =\left(\begin{array}{c}\hfill \alpha \hfill \\ \hfill {\mathit{v}}_0\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\mathit{u}}_1\hfill \\ \hfill {\mathit{u}}_{2:}\hfill \end{array}\right),\hfill \\ \hfill \mathcal{A}\left(t\right)& =\left\{\left(\dot{\mathit{y}},\mathit{y},\mathit{u}\right):\begin{array}{cc}\hfill \dot{\mathit{y}}& ={\mathbf{S}}_{\mathit{y}}\cdot \mathit{v}\hfill \\ \hfill \mathbf{0}& ={\mathbf{S}}_{\mathit{x}}\cdot \mathit{v}\hfill \\ \hfill \mathbf{0}& \le \hfill & \mathit{y}\hfill \\ \hfill \mathit{l}\mathit{b}\left(t\right)& \le \mathit{v}\le \mathit{u}\mathit{b}\left(t\right)\hfill \\ \hfill \mathit{h}\left(t\right)& \hfill \ge {\mathbf{H}}_{\mathit{y}}\left(t\right)\cdot \mathit{y}+{\mathbf{H}}_{\mathit{v}}\left(t\right)\cdot \mathit{v}\hfill \\ \hfill \mathit{v}& ={\mathit{u}}_{2:}\cdot \exp \left(\lambda \cdot \left(t-t_0\right)\right)\hfill \\ \hfill \lambda & =\ln \left({\mathit{u}}_1\right)/\left(t_{\text{end}}-t_0\right)\hfill \end{array}\right\},\hfill \\ \hfill {\mathit{\varphi }}_{\text{bndry}}& =\left(\begin{array}{c}\hfill \mathit{y}\left(t_0\right)-{\mathit{y}}_0\hfill \\ \hfill {\mathit{y}}_0-\mathit{y}\left(t_0\right)\hfill \end{array}\right),\hfill \\ \hfill f\left(\mathit{y}\left(t\right),\mathit{u}\right)& =-{\mathit{u}}_1.\hfill \end{array}$

Note that the control variables u are not time-dependent, i.e., they enter the model as control parameters rather than functions that need to be optimized in the sense of optimal control.

2.3 Time-Optimal Adaptation: Definition and Forms

Previous frameworks for constraint-based optimization did not explicitly include the time interval as part of the optimization objective. In the following, we introduce Time-Optimal Adaptation (TOA) as a framework to analyze transition between different cellular states in the shortest possible time. TOA is motivated by the assumption that under certain environmental conditions, cells may have evolved to reach target amounts y^goal in the shortest possible time, starting from initial amounts y^init. This transition might either take place in a variable environment, encoded by a time-dependent set $\mathcal{A}(t)$, or in a constant environment. Likewise, the target and initial amounts may either have to fulfill additional optimality criteria, or may correspond to pre-defined or experimentally measured states. Mathematically, we capture such a strategy in the following way.

Time-Optimal Adaptation

Given an initial/current amount of molecules ${\mathit{y}}^{\text{init}}\in {\mathbb{R}}^{n_{\mathit{y}}}$ and a target amount ${\mathit{y}}^{\text{goal}}\in {\mathbb{R}}^{n_{\mathit{y}}}$, the optimization objective is to transition from the former to the latter as quickly as possible.

$\underset{\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right),T>t_0}{\min }T$(2.7a)

$\begin{array}{ll}\hfill \text{s.\hspace{0.17em}t.\hspace{0.17em}}& \begin{array}{cc}\hfill \mathit{y}\left(t_0\right)& ={\mathit{y}}^{\text{init}},\hfill \\ \hfill \mathit{y}\left(T\right)& ={\mathit{y}}^{\text{goal}}\hfill \end{array}\hfill \end{array}$(2.7b)

$\text{and}\left(\dot{\mathit{y}}\left(t\right),\mathit{y}\left(t\right),\mathit{u}\left(t\right)\right)\in \mathcal{A}\left(t\right)\text{for}\left(\text{almost}\right)\text{all}t\in \left[t_0,T\right],see(2.3a)$(2.7c)

The constraints (2.7c) and (2.7b) can be framed within the abstract constraint-based framework (2.3) by including y^init and y^goal using

${\mathit{\varphi }}_{\text{bndry}}\left(\mathit{y}\left(t_0\right),\mathit{y}\left(t_{\text{end}}\right),\mathit{u}\left(t_0\right),\mathit{u}\left(t_{\text{end}}\right)\right)=\left(\begin{array}{c}\hfill \mathit{y}\left(t_0\right)-{\mathit{y}}^{\text{init}}\hfill \\ \hfill {\mathit{y}}^{\text{init}}-\mathit{y}\left(t_0\right)\hfill \\ \hfill \mathit{y}\left(t_{\text{end}}\right)-{\mathit{y}}^{\text{goal}}\hfill \\ \hfill {\mathit{y}}^{\text{goal}}-\mathit{y}\left(t_{\text{end}}\right)\hfill \end{array}\right),$

whereas the global objective function, cf. f in (2.3), does not explicitly contain any of the variables y or u. Instead, the general framework of constraint-based modeling (2.3) is extended through time-optimal adaptation by using the end point of the time interval of interest itself as the optimization objective function. In contrast to the frameworks with non-time-dependent objective function as defined in (2.3c), TOA provides solutions (y(t), u(t)) only on the time interval [t₀, T] instead of (arbitrary) [t₀, t_end].

Remark 2.5. Within this work, we assume that the target amounts y^goal are accessible. Specifically, we assume that a time t_end ≥ T exists such that all values within the optimization problem defining TOA are well-defined. We note that the accessibility of the target state is a classical problem in time-optimal control, and accessibility is a prerequisite for applying TOA. In practice, the target state will often be defined by means of an RBA solution and we conjecture that these target states will be accessible.

Remark 2.6. Within this work, we use the term “adaptation” in a control-theoretic sense. That is, the term refers to changes in the intracellular amounts or concentrations in response to the environmental conditions, respecting the given constraints. In an evolutionary context, such changes are typically considered as “acclimation”.

Remark 2.7. We do not require the constraint set $\mathcal{A}(t)$ in (2.3a) to have any specific form. This means that time-optimal adaptation can be defined irrespective of the concrete modeling paradigm underneath the simulation. Practically, even discrete time/state systems fit well within TOA. To be concise, however, we concentrate in the following on frameworks closely related to deFBA and cFBA. In Example 2.8, we therefore introduce TOA also in a simplified setting that directly builds upon d(e)FBA, cf. (Waldherr et al., 2015; Höffner et al., 2016). From the viewpoint of the general framework (2.3), this is a special case of deFBA with a modified objective function.

Example 2.8. (TOA as an extension of deFBA). Assume that there is no distinction between “fast” and “slow” components within the metabolic network. In this case, the dynamics of its molecular amounts can be described purely by ordinary differential equations $\dot{\mathit{y}}(t)=\mathbf{S}\cdot \mathit{v}(t)$. As for classical flux balance analysis, the fluxes are constrained by upper and lower bounds lb, ub that might depend on the possibly changing environment, i.e., ub(t) ≤ v(t) ≤ ub(t). If y contains compounds with enzymatic function, the flux rates (or weighted sums thereof) may additionally be constrained by (weighted sums of) components of y. Such bounds can be collected into a single set of linear inequalities by introducing suitable matrices/vectors H_y(t), H_v(t), h(t), i.e.,

${\mathbf{H}}_{\mathit{y}}\left(t\right)\cdot \mathit{y}\left(t\right)+{\mathbf{H}}_{\mathit{v}}\left(t\right)\cdot \mathit{v}\left(t\right)\le \mathit{h}\left(t\right),$

see (Waldherr et al., 2015, Section 2.3) for a detailed description. To account for y being total amounts, y is constrained to positive values, i.e., y(t) ≥0. As outlined above, TOA requires fixed initial and terminal values for the molecular amounts, mathematically captured by

$\mathit{y}\left(t_0\right)={\mathit{y}}^{\text{init}},\mathit{y}\left(T\right)={\mathit{y}}^{\text{goal}}.$

In summary, TOA can be aggregated in this simplified case to the following constrained optimization problem

$\begin{array}{cc}\hfill \underset{\mathit{y}\left(\cdot \right),\mathit{v}\left(\cdot \right),T>t_0}{\min }& T\hfill \\ \hfill \text{s.t.}& \left.\begin{array}{cc}\hfill \dot{\mathit{y}}\left(t\right)& \hfill =\mathbf{S}\cdot \mathit{v}\left(t\right)\hfill \\ \hfill \mathit{l}\mathit{b}\left(t\right)& \hfill \le \mathit{v}\left(t\right)\le \mathit{u}\mathit{b}\left(t\right)\hfill \\ \hfill \mathit{h}\left(t\right)& \hfill \ge {\mathbf{H}}_{\mathit{y}}\left(t\right)\cdot \mathit{y}\left(t\right)+{\mathbf{H}}_{\mathit{v}}\left(t\right)\cdot \mathit{v}\left(t\right)\hfill \\ \hfill \mathit{y}\left(t\right)& \hfill \ge \mathbf{0}\hfill \end{array}\right\}\hfill \\ \hfill & \left.\begin{array}{cc}\hfill \dot{\mathit{y}}\left(t_0\right)& \hfill ={\mathit{y}}^{\text{init}}\hfill \\ \hfill \dot{\mathit{y}}\left(T\right)& \hfill ={\mathit{y}}^{\text{goal}}\hfill \end{array}\right\}\hfill \end{array}$

The notation “${\min }_{\mathit{y}(\cdot ),\mathit{v}(\cdot ),T>t_0}$” can be understood in the sense of optimal control, i.e., one is searching for the optimal objective value among all (differentiable) functions y(t), t ∈ [t₀, T], and (measurable) functions v(t), t ∈ [t₀, T]. The framework identifies possible time courses for the fluxes v(t) and amounts y(t) such that (i) stoichiometry, (ii) flux bounds, and (iii) enzyme activities are included in the model and such that the transition from one given amount to another is as fast as possible.

2.4 Applications and Case Studies

Next we introduce two particularly relevant applications of TOA.

Application 2.1. (Cell Doubling). A first natural application of TOA is cell doubling, where the objective is to double all cellular components in minimal time, such that

${\mathit{y}}^{\text{goal}}=2\cdot {\mathit{y}}^{\text{init}}.$

The resulting trajectory thus can be interpreted as one cell cycle. Neither the initial, nor the target amount have to be optimal with regard to other objectives. Within the TOA framework, cell doubling can be considered either in a constant environment, or with time-dependent external conditions. We note that applications of constraint-based optimization of metabolism typically do not distinguish between solutions for a single cell vs. solutions for a homogeneous population of cells. Similarly, the time courses for cell doubling predicted by TOA can either be interpreted for a single cell or a homogeneous, synchronized population of cells. If cells are not synchronized, that is each cell within the population is at a different time point with respect to its cell cycle, we have to average over the population or, equivalently, over a full cell cycle, to obtain in silico measurements of a population.

Application 2.2. (Transitions after a nutrient shift). A second important application of TOA is to consider a sudden change in the external conditions, i.e., from a given constant nutrient availability for t < 0 to a different one for t ≥ 0. In this scenario, TOA can be utilized to predict the transition of the intracellular amounts y^init to new target amounts y^goal. The new target amounts might either be optimal with respect to the new environmental conditions (as defined by RBA), or be provided otherwise (for example by experimental observations). In both scenarios, the target amounts are typically defined in terms of concentrations instead of absolute amounts. Hence, we must also formulate the boundary conditions in terms of c(t),

$\begin{array}{cc}\hfill \mathit{y}\left(t_0\right)& ={\mathit{y}}^{\text{init}},\hfill \\ \hfill \mathit{c}\left(T\right)& =\frac{1}{bio\left(T\right)}\cdot \mathit{y}\left(T\right)=\frac{1}{{\mathit{w}}^{\top }\cdot {\mathit{y}}^{\text{goal}}}\cdot {\mathit{y}}^{\text{goal}}={\mathit{c}}^{\text{goal}}.\hfill \end{array}$(2.8)

As shown in Supplementary Appendix S3, it is possible to rearrange conditions (2.8) such that a linear equality system in the unknowns (y(t₀), y(T)) is obtained. Therefore, the concentration-based definition has no immediate drawbacks regarding the numerical solution.We must further consider that an as-quick-as-possible transition from one intracellular concentration to another does not incorporate the overall (i.e., biomass-) growth of the cell and thus might not represent an evolutionarily plausible strategy. Rather, the transition to new external conditions involves a balance between fast transition to a (better adapted) novel state and the requirement to increase (or not decrease) the total biomass of the cell. To obtain a general framework, we therefore propose a two-objective optimization problem:

$\begin{array}{cc}\hfill \underset{\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right),T>0}{\min }& \left(\begin{array}{c}\hfill T\hfill \\ \hfill -\alpha \hfill \end{array}\right)\hfill \\ \hfill \mathit{\text{s.t.\hspace{0.17em}}}\mathit{y}\left(0\right)& ={\mathit{y}}^{\text{init}}\hfill \\ \hfill \mathit{y}\left(T\right)& =\alpha \cdot {\mathit{y}}^{\text{goal}}\hfill \\ \hfill \mathrm{and}\left(\dot{\mathit{y}}\left(t\right),\mathit{y}\left(t\right),\mathit{u}\left(t\right)\right)& \in \mathcal{A}\left(t\right)\text{for}\left(\text{almost}\right)\text{all}t\in \left[0,T\right],\hfill \end{array}$(2.9)

where y^init denotes a normalized vector of intracellular amounts which describe the cells for the environment t < 0. The “normalization” here can, for example, be understood as w^⊤⋅y^init = 1. Accordingly, y^goal denotes a normalized vector for the environmental conditions after the nutrient shift. “Minimality” in (2.9) is to be understood in the sense of Pareto: a triple (y(t), u(t), T) is optimal if T cannot be decreased without decreasing α such that y(T) = α ⋅y^goal, and vice-versa if α cannot be increased without also increasing the end time T. The set of all optimal solutions of (2.9) describes the different compromises between fast adaptation and continued growth.Remark 2.9. Note that the boundary conditions (2.8) do not entail any direct condition concerning bio(t_end) = w^⊤⋅y(t_end). If the metabolic network allows for a quick degradation of compounds, it might be optimal (in the sense of TOA) to shrink (in terms of absolute biomass) before actually adapting to the new concentrations, or even to completely disintegrate all metabolic compounds to zero. Such a behavior would be in line with the description of time-optimal time courses as induced by (2.8). To remedy this, a linear inequality

$\begin{array}{cc}\hfill bio\left(t_{\text{end}}\right)& \ge bio\left(t_0\right)\hfill \\ \hfill \text{\hspace{0.17em}or\hspace{0.17em}}{\mathit{w}}^{\top }\cdot \mathit{y}\left(t_{\text{end}}\right)& \ge {\mathit{w}}^{\top }\cdot \mathit{y}\left(t_0\right)\hfill \end{array}$

can be added, illustrating again the power of constraint-based modeling. Whenever necessary, this was done for the in silico experiments in Section 3.

2.5 Numerical Solution

The optimization problem (2.7) of TOA contains a general condition on the dynamics of y in the form of (2.7c). To design an algorithm able to cope with this generality, we assume that a numerical method is available that can simulate this dynamic behavior subject to boundary conditions on a given fixed time interval [t₁, t₂] ⊆ [t₀, T] and/or to determine whether such a solution exists. Provided this condition (and tacitly assuming that the relevant feasible end time points T lie in a connected set), the minimal value for T can be found using any one-dimensional root finding algorithm. For its simplicity and guaranteed convergence, we propose to use the following bisection method for the determination of T in TOA:

For the initial time interval [t_min, t_max], one needs to assume that (2.7b) and (2.7c) define an infeasible problem on [t₀, t_min], while the corresponding problem on [t₀, t_max] is feasible. The quick convergence of the bisection method entails that an already very good initial guess is not crucial for an efficient implementation, as long as the simulation task is not too computationally expensive.

If there is legitimate doubt about the result, the algorithm can be re-started with another initial interval or one can change to a more fine-grained sampling for the evaluation of feasible and infeasible points. The numerical results in Section 3 were preceded by an exhaustive scan of end time points, which indicated that the set of feasible end time points do indeed form a single interval (i.e., a connected set) in all shown examples.

Remark 2.10. The bisection method was chosen here for several reasons over more “classical” methods in time-optimal control: firstly, the “simplicity” aspect of the bisection method does not only refer to it being easily applicable for various extensions of the framework (like time- and/or state-discrete systems, or a framework that incorporates heterogeneity within a community or in space) but also to the implementation. Many existing toolboxes include interfaces to (MI-)LP solvers. Algorithms for dynamic simulations are moreover often highly optimized, such that checking for feasibility over a given time range can be more efficient than implementing a new interface to an optimal control library.

Secondly, the inherently linear structure of problems like deFBA should be preserved. For time-dependent constraint sets $\mathcal{A}(t)$ this is only possible if the time variable t is treated as the independent variable in the optimal control algorithm. In existing optimal control libraries like BOCOP (Bonnans et al., 2017), time-optimal control problems are often transferred to optimal control problems on a unit interval by introducing an artificial independent variable. If the time-dependency of some of the constraints is non-linear, this translates to the optimization problems that need to be solved within the optimal control routine.

We note, however, that in the non-linear case the application of “classical algorithms” for time-optimal control problems like shooting-methods, or those based on the Pontryagin principle might generally outperform the bisection approach taken here.

Remark 2.11. For the solution of the Pareto problem (2.9) it is not necessary to implement algorithms for maximizing α, i.e., optimizing the second objective. Instead, one can continue using the algorithm for time-optimal adaptation while simultaneously fixing feasible values of α. With respect to the definition of Pareto-optimality, this means that for any feasible value of one objective, the other one is optimized, corresponding to the so-called ϵ-constraint method in multi-objective optimization, cf. (Ehrgott, 2000).

2.6 Time-Optimal Adaptation Variability Analysis

Minimizing T need not suffice to uniquely determine the time courses in y. If this is the case, the variability over time can be captured by enumerating possible time series once the optimal end time point was found. We will refer to this procedure as TOA-Variability Analysis (TOA-VA). In contrast to static flux variability analysis (FVA), there are several ways to define what such an enumeration means. One way would be to determine the maximum and minimum possible value for all components of y and separately at each time point. This, however, would not only lead to time-consuming computations, but would also be difficult to interpret: a numerical solution that is constructed via putting together maxima or minima $\mathit{y}(\tilde{t})$ for all time instances $\tilde{t}$ does not have to fulfill the dynamics defined by the original model. Here, we understand TOA-VA as the minimization and maximization of the integral over all components of y, i.e., for all i = 1, 2, …, n_y:

$\underset{\mathit{u}\left(\cdot \right),\mathit{y}\left(\cdot \right)}{\min }\pm {\int }_{t_0}^T{y}_i\left(t\right)\mathrm{d}t$(2.10)

subject to the dynamic and/or boundary constraints in the original problem. Note that the overall time courses might still not be uniquely defined from (2.10).

To explore the variability of the time courses for the concentrations c(t), we use the following variant of TOA-VA, which is called relative TOA-VA:

(i) Compute the optimal end time point T of time-optimal adaptation.

(ii) For all i = 1, 2, … , n_y: Use TOA-VA as in (2.10) to obtain a minimal value I_min,i and a maximal value I_max,i for the integral of y_i over [t₀, T].

(iii) Calculate for all i = 1, 2, … , n_y the (maximal and minimal) concentrations c_i(t) as given by (2.1) where y(t) is calculated from

$\begin{array}{lll}\hfill \underset{\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right)}{\max }{\int }_{t_0}^{T}bio\left(t\right)\mathrm{d}t& \hfill & \hfill \text{s.t.\hspace{0.17em}}{\int }_{t_0}^T{y}_i\left(t\right)\mathrm{d}t\le I_{\text{min},i},\\ \hfill \underset{\mathit{y}\left(\cdot \right),\mathit{u}\left(\cdot \right)}{\min }{\int }_{t_0}^{T}bio\left(t\right)\mathrm{d}t& \hfill & \hfill \text{s.t.\hspace{0.17em}}{\int }_{t_0}^T{y}_i\left(t\right)\mathrm{d}t\ge I_{\text{max},i}\end{array}$

(again subject to the original constraints of the problem).

There is still no guarantee that the solutions to these problems are unique. However, since the concentrations are defined as the ratio of total amounts to the biomass, the above definition is reasonable as one is maximized whilst minimizing the other. Note, that this definition implies that the weighted sum of all (maximal or minimal) concentrations no longer needs to add to the total biomass.

2.7 Implementation

The calculations for all experiments in Section 3 were done in Python 3.8.1 on a laptop computer. Scripts that reproduce the numerical experiments below are available on GitHub, https://github.com/MarkusKoebis/StaticTOA_py The numerical solutions were determined from a complete parameterization (using the trapezoidal rule) of the compounds and fluxes over the entire time range of interest using n = 100 steps on an equidistant grid. This leads to a sparse LP problem which was solved using gurobipy on Gurobi 9.0.1 solver (Gurobi, 2021) with standard settings (concerning problem formulation and tolerances). Most experiments were repeated (for verification) with tight error tolerances without notable differences. For time-optimal adaptation, no objective vector for the LPs is necessary, so we used the null vector 0. For (relative) TOA-VA, the integrals in the objective or constraints were approximated using the same time grid and also the trapezoidal rule.

3 Results

3.1 A Coarse-Grained Self-Replicator Model

We illustrate TOA by means of a coarse-grained self-replicator model (Molenaar et al., 2009; Giordano et al., 2016; Yabo et al., 2022). The model, cf. Figure 1, consists of three compounds: M (intracellular metabolic precursor), Tr (transporter), and R (ribosome), as well as five biochemical reactions, and the external nutrient N. The uptake of the external nutrient N is catalyzed by the transporter Tr and depends on the availability of N via a Michaelis-Menten rate equation. Depending on the application, the concentration of the external nutrient N may either be constant or vary over time. The synthesis of the catalytic macromolecules Tr and R is limited by the ribosome amount. Within the model, macromolecules can be disassembled into the precursor M. For energetic consistency, however, disassembly results in fewer precursor molecules than required for synthesis, reflecting the energy expenditure of protein synthesis and thereby avoiding futile cycles. We note that within the model, no compound is subject to the quasi steady-state assumption, and the metabolic precursor M can accumulate over time. Hence M also serves as a storage compound. All constraints of the model can be formulated in terms of linear inequalities. A detailed definition is provided in Supplementary Appendix S2.1.

FIGURE 1

FIGURE 1. A schematic illustration of the coarse-grained self-replicator model; solid lines represent biochemical reactions between the nodes (biochemical compounds), dashed dark-blue lines indicate that a reaction is catalyzed by the respective compound. Abbreviations: N, external nutrient; M, metabolic precursor/storage; Tr, transporter; R: ribosome; v_N, nutrient uptake reaction; v_R, ribosome production reaction; vd_R, ribosome degradation reaction; v_Tr, transporter production reactions; vd_Tr, transporter degradation reaction.

3.2 Constant Environments and RBA

Before the dynamic behavior of the model is studied by means of TOA, we summarize the steady-state properties of the model in a constant environment using Resource Balance Analysis (RBA). RBA provides a method to calculate the steady-state amounts of the cell that maximize the growth rate under constant external conditions, i.e., for a constant external nutrient concentration. In the following, extracellular nutrient is measured relative to the Michaelis constant K_M of the uptake reaction, with N/K_M as a dimensionless parameter.

Figure 2A shows the maximal growth rate λ as a function of the relative nutrient availability. The growth rate follows a Monod equation with a maximum λ^max ≈ 0.435 h⁻¹ and an effective (dimensionless) affinity constant K_A ≈ 0.347, corresponding to the value of the relative nutrient availability N/K_M at which the cell grows at half the maximal growth rate λ^max.

FIGURE 2

FIGURE 2. (A) Maximal growth rate λ as a function of the (relative) extracellular nutrient availability as predicted by RBA. (B) Cellular amounts of intracellular compounds as functions of relative nutrient availability. Extracellular nutrient is measured relative to the Michaelis constant K_M of the uptake reaction.

Figure 2B shows the total amounts of the three intracellular components M, Tr, and R as a function of the (relative) nutrient availability. The amounts were scaled such that the total biomass always equals one unit (e.g., 1 g cellular dry mass). As expected, when maximizing the growth rate, the level of the precursor/storage component M is always zero. This reflects the fact that the precursor M has no catalytic activity, and any non-zero amount of M would consume resources that otherwise could be allocated to transport or protein translation.

The amounts of the other intracellular components Tr and R follow the well-known growth laws of microbiology (Scott and Hwa, 2011). The concentrations are a function of the growth rate, and hence the external nutrient availability, the well-known linear relationship is shown in Supplementary Appendix S4. With increasing nutrient availability, the relative amount of transporter decreases, whereas the relative amount of ribosome increases.

3.3 TOA in Constant Environments

Our first case study using TOA is to consider the doubling of a microbial cell in minimal time. We assume that the self-replicator model in Figure 1 has pre-described initial amounts y(t₀) = y₀ which simultaneously identify the pre-defined initial state y^init. The objective is to double all cellular components as fast as possible, cf. Application 2.1. The environment is assumed to be constant with a relative (external) nutrient availability N/K_M = 1. The initial (and final) amounts are not assumed to be optimal for the given environment. Instead, y(t₀) is obtained by solving an RBA problem corresponding to N/K_M = 2.0. In other words, the cell is assumed to be adapted to a higher nutrient level than is present in the current environment. In the following, we will refer to such cells as “optimists”.

Figure 3 shows the time course of intracellular components for one cell doubling. The predicted time-optimal amounts of metabolic compounds are shown as solid lines (red, blue, and yellow), the total biomass is shown in green. The dashed lines correspond to a solution obtained by iterative RBA (cf. Example 2.4), which corresponds to exponential growth of all cellular components with no further internal degrees of freedom. Figure 3B shows the respective flux rates over the simulated time range. Solid lines again indicate the solution of TOA, while dashed lines (exponential curves) correspond to the solution found with iterative RBA.

FIGURE 3

FIGURE 3. Cell cycle of an “optimistic” cell; (A) amounts and biomass as a function of time, (B) flux rates as a function of time; solid lines indicate the solution of TOA, dashed lines indicate iterative RBA (exponential growth) with the same “optimistic” initial values.

Using TOA, the time for one cellular doubling is T = 2.17 h. In contrast, the solution based on iterative RBA results in a slightly longer doubling time of T = 2.34  h, showing that internal degrees of freedom shorten the calculated division time. The time course of y(t) over one cell doubling can be subdivided into four time intervals (marked as I-IV in Figure 3A). At the beginning (marked as interval “I”), cell growth is limited by the lack of transporter Tr due to the “optimistic” initial configuration of the cell. Hence, ribosome R is actively disassembled into precursor M to increase the synthesis of Tr. In interval “II”, the cell is perfectly adapted to the given nutrient environment and grows exponentially, before the re-adaptation to the target composition y^goal = 2 y^init begins in interval “III”. Within interval “III”, the cell still has an overabundance of Tr, which allows it to accumulate the precursor M. In the final interval “IV”, synthesis of transporter Tr ceases and all resources are devoted to the synthesis of the ribosome R, until the target amounts y^goal are reached.

The biological plausibility of these time courses is discussed in Section 4. Here we only summarize the following results: Given the initial amounts y^init, cell doubling using TOA in time-invariant environments gives rise to complex intracellular dynamics different from solutions obtained by iterative RBA. Importantly, these solutions involve a transient accumulation of the precursor M as a storage compound–a phenomenon not observed with iterative RBA. The minimal division time predicted by time-optimal adaptation is shorter than division times obtained by iterative RBA.

So far, we considered a particular initial amount y^init such that the cell was adapted to a higher nutrient availability than actually present in the environment (“optimist”). To obtain a broader view, we evaluated cell doubling using TOA in different time-invariant environments with initial (and final) amounts adapted to different external nutrient availability. The results are shown in Figure 4. Solid lines correspond to intracellular amounts using TOA, dashed lines correspond to a solution obtained with iterative RBA (exponential growth without internal degrees of freedom). Shaded areas correspond to variability in the sense of TOA-VA (cf. Section 2.6), i.e., possible solutions that equally satisfy all constraints and the optimality criterion. In this case, the solid lines display a “nominal” solution, i.e., one that was provided by the algorithm before an additional variability analysis (we note that since the numerical solution is based on a feasibility problem, the LP solver has no incentive to favor a smooth solution to any other).

FIGURE 4

FIGURE 4. Time course solutions of time-optimal adaptation and a cell doubling experiment under different constant external nutrient conditions; solid lines: TOA, shaded areas: TOA-VA, dashed lines: iterative RBA (simulated until cell doubling was achieved); upper row: pessimistically adapted, middle row: perfectly adapted (recovery of iterative RBA), bottom row: optimistically adapted for constant relative nutrient availability of N/K_M = 0.5 (left column), N/K_M = 1.0 (middle column), and N/K_M = 5.0 (right column).

Columns in Figure 4 correspond to different relative nutrient availability levels: the first column to a nutrient availability N(t)/K_M ≡ 0.5; the second column to N(t)/K_M ≡ 1.0, and the third to N(t)/K_M ≡ 5.0. The rows in Figure 4 correspond to different “expectations” of the cells, that is, which external nutrient availability the initial (and final) amounts are adapted to. Specifically, the first row corresponds to “pessimists”. That is, cells adapted to a nutrient availability below the one present in the environment, while retaining the objective to double all cellular components in minimal time. The second row corresponds to cells perfectly adapted to the environmental nutrient availability. The final row corresponds to “optimists”, i.e., cells adapted to a higher nutrient availability than present in the environment.

The latter scenario corresponds to the example already shown in Figure 3. We again observe an initial increase in the transporter synthesis, followed by a delayed onset of ribosome synthesis. Importantly, in each case, we can see a transient accumulation of storage M(t) that is absent in solutions obtained by iterative RBA. In the case of perfectly adapted cells (middle row), solutions obtained by TOA are equivalent to solutions obtained by iterative RBA. For “pessimistic” cells (top row), we again observe complex time courses. In particular, cells adapted to lower nutrient levels than present in the environment exhibit an overabundance of transporter. Hence, we observe an initial rapid uptake of nutrient and transient accumulation of the precursor M. In the initial interval, resources are primarily allocated to the synthesis of ribosomes. Only in the later interval, the transporter is synthesized to the required amounts (even at the expense of ribosomes that may be disassembled into precursors). The transient accumulation of precursor M exhibits considerable variability and the solutions of TOA are no longer unique.

A detailed discussion about the biological plausibility of these time courses is again relegated to Section 4. Here we only note that, despite the simplicity of the model, the solutions exhibit a wide variety of qualitatively different complex temporal behaviors, including the transient accumulation of the precursor M.

3.4 The Role of Expectation: Optimists vs. Pessimists

We further investigate two key observations obtained in the previous experiments: the transient accumulation of precursor M as a storage compound, as well as the impact of the initial cellular state on the predicted doubling time.

Firstly, Figure 5 shows the average storage concentration predicted for a population of cells adapted to different nutrient availabilities (N/K_M ∈ (0.2, 2.0), x-axis) in an environment with an actual relative nutrient availability N/K_M ≡ 1.0. To calculate the average storage concentration predicted by TOA for a population of cells, we assume that the (in silico) measurements are taken from a heterogeneous population of unsynchronized cells that are (equidistributed) at various stages of their cell cycle. To take this non-uniform age distribution into account, the population average was computed, cf. (Powell, 1956), as

$\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\left(M\right)=\frac{1}{T}ċ{\int }_0^T\frac{M\left(t\right)}{{\mathit{w}}^{\top }\cdot \mathit{y}\left(t\right)}\mathrm{d}t,$(3.1)

where y(t) is a solution obtained by relative TOA-VA, cf. Section 2.6.

FIGURE 5

FIGURE 5. Influence of optimistic and pessimistic goal states in cell doubling: Main plot: mean relative storage accumulation, see (3.1), as a function of nutrient adaptation level. Bottom row: Three selected time courses, cf. Figure 4, for nutrient adaptation levels N/K_M of 0.2, 1.1, and 2.0. For N/K_M < 1, the quantity mean(M) is no longer unique such that a shaded area indicates the possible range, as TOA-VA also predicts a range of possible solutions (shaded area in the bottom left plot).

As shown in Figure 5, we observe (the possibility of) a nonzero average storage concentration for all cellular states that are not perfectly adapted to the respective environment. For optimistic cells adapted to a higher nutrient availability than present in the environment, the average storage concentration increases slightly with the distance to the perfectly adapted state. The effect is more pronounced for pessimistic cells adapted to a lower nutrient availability than present in the environment. In this case, the solutions of TOA are not unique and the range of average storage is indicated as a shaded area. For “pessimist” cells, the large average storage is due to a high abundance of transporter molecules, which implies that uptake and accumulation of precursor is not restricted.

Secondly, Figure 6 shows the predicted growth rate for cells adapted to a different relative nutrient availability (N/K_M ∈ (0.2, 2.0)) than present in the environment (N/K_M ≡ 1.0). The straight line indicates the growth rate of cells that are perfectly adapted, resulting in a maximal growth rate λ = λ_env ≈ 0.32  h⁻¹. The maximal growth rates for cells adapted to a different environment (misadaptation) are shown as a solid green line for solutions obtained with TOA and as a purple line for solutions obtained with iterative RBA.

FIGURE 6

FIGURE 6. Growth rate of differently adapted cells as predicted by cell-doubling experiments using TOA and iterative RBA in an environment with relative nutrient availability N/K_M = 1.0; λ_env ≈ 0.32  h⁻¹ denotes the maximal growth rate as predicted by RBA.

We observe that misadaptation always results in a reduced growth rate, as compared to a perfectly adapted cell. However, solutions obtained by TOA always outperform solutions obtained by iterative RBA, demonstrating that internal degrees of freedom and transient accumulation of storage shorten the predicted doubling time. Furthermore, the decrease in growth rate is more pronounced for “pessimistic” adaptation, that is, for cells that are adapted to a lower nutrient level than present in the environment. In contrast, “optimistic” adaptation, that is, cells are adapted to a higher levels than present in the environment, together with TOA results in growth rates close to perfectly adapted cells–indicating that “optimistic” adaptation carries a lower evolutionary cost than “pessimistic” adaptation.

3.5 Time-Optimal Adaptation at a Nutrient Shift

As our second application, we consider a nutrient shift, i.e., a sudden change in the external conditions from a given constant nutrient availability for t < 0 to a different one for t ≥ 0. TOA is utilized to predict the time-optimal transition of a cell perfectly adapted to the initial state at t < 0 to a state perfectly adapted to maximize growth in the new environment for t ≥ 0. As noted in Section 2.4, the target state for the new environment is typically defined in terms of concentrations rather than absolute amounts, because it is unknown whether or how much the cells are able to grow during adaptation.

Figure 7 shows the resulting time courses for the coarse-grained self-replicator model used in the previous sections. Shown are time-optimal shifts from a low nutrient availability to a higher nutrient availability (left column in Figure 7), as well as time-optimal shifts from a high nutrient availability to a lower nutrient availability (right column in Figure 7). Non-unique solutions are again displayed as shaded areas indicating the maximum and minimum range in which solutions can be found (TOA-FVA, see Section 2.6). We observe that the time-optimal transition from lower to higher nutrient availability again entails a transient accumulation of storage.

FIGURE 7

FIGURE 7. Adaptation to a single nutrient jump (shown as a dashed green line), left column: adaptations from poorer to richer medium, right column: adaptation to scarcer environment; shaded areas: solutions in the sense of TOA-VA. We note that for t < 0, TOA makes no assumptions about y(t).

As detailed in Section 2.4, time-optimal adaptation alone is not sufficient as an evolutionary principle to explain cellular adaptation after a nutrient shift. Rather, we consider a two-objective optimization (in the sense of Pareto) with the conflicting objectives of a fastest possible adaptation to the new state vs. a maximal increase in total cellular biomass.

Figure 8 (main panel) shows the resulting Pareto fronts for different transitions in terms of the minimal time T* for adaptation vs. the maximal increase in cellular biomass given by the factor α, cf. (2.9). Panels A–D in Figure 8 show selected time courses of intracellular amounts at different positions of the Pareto front. In the subplots A and B, the shaded areas indicate that the cell is perfectly adapted to the environment in the sense of RBA, i.e., from the start of the shaded areas, the cell exhibits balanced exponential growth at the maximal growth rate and no further internal dynamics take place. The absence of internal dynamics explains that, for larger values of α or T*, the lines in the main plot become asymptotically parallel.

FIGURE 8

FIGURE 8. Two-objective optimization of adaption time T* and total biomass growth factor α for time-optimal adaptation at a nutrient shift. Main plot: Pareto fronts for three different initial adaptations (measured in N/K_M) of 0.2, 1.0 and 10. Subplots (A–D): Time courses at different points on the Pareto fronts, cf. Figure 7. The shaded areas in subplots A and B indicate time intervals where the cell is perfectly adapted (exponential growth).

In the absence of a nutrient shift (i.e., the transition N/K_M: 1 → 1, blue line in Figure 8), the minimal time for adaptation is T* = 0 with a growth factor α = 1, in this case the relationship between transition times T* > 0 and increase in biomass is consistent with exponential growth (note the logarithmic scale on the y-axis).

For a nutrient shift from high to low nutrient availability (N/K_M: 10 → 1, black line) the minimal transition time is $T^{*}=T_1^{*}\approx$ 0.44 h. Figures 8A–C show two representative transitions on the Pareto front with panel A corresponding to a scenario that prioritizes an increase in biomass (factor α) over the transition time T*, and panel C corresponding to a scenario that prioritizes a minimal transition time over the accumulation of biomass.

For a nutrient shift from low to higher nutrient availability (N/K_M: 0.2 → 1, green line) the minimal transition time is $T^{*}=T_2^{*}\approx$ 1.14 h. Figures 8B,D show two representative transitions on the Pareto front with panel B corresponding to a scenario that prioritizes an increase in biomass (factor α) over the transition time T*, and panel D corresponding to a scenario that prioritizes a minimal transition time over the accumulation of biomass. In either case, the optimal transition involves a transient accumulation of the storage compound M.

Consistent with results in the previous section, Figure 8 also shows that “optimistic” adaptation carries a lower evolutionary cost than “pessimistic” adaptation. A cell adapted to high nutrient availability exhibits only a slightly reduced biomass increase when transitioning into a low nutrient environment, as compared to a cell already adapted to this environment. In contrast, a cell adapted to a lower nutrient environment exhibits a more pronounced reduction in accumulated biomass when transitioning into higher nutrient availability, as compared to either a cell that is already adapted to the higher nutrient availability, or likewise as compared to a cell that was previously adapted to even higher nutrient availability.

4 Discussion

In this work we introduced TOA, a novel approach to simulate and predict time-optimal adaptation of microbial metabolism and growth. While time-optimal modeling has been considered before, see, among others, (Klipp et al., 2002) (temporal gene expression), (Pavlov and Ehrenberg, 2013) (fast proteome adaptation to environmental change), (Waldherr et al., 2015) (maximize survival time under nutrient depletion), (Basan et al., 2020) (minimization of lag/response-time), and (Djema et al., 2020) (bio-reactor applications), our work builds upon the recent advances in dynamic constraint-based modeling, such as dFBA, deFBA and cFBA, cf. Section 2.2. TOA is versatile and extends most approaches currently employed in constraint-based modeling of microbial metabolism and growth.

In particular, while the analysis of balanced steady-state growth dominates current experimental and computational studies, in most natural environments microbes have to continuously adapt to perturbations and changes in nutrient availability. TOA allows us to study such transitions in the context of established constraint-based models of microbial metabolism. Similar to other constraint-based methods, the solutions obtained from TOA are not based on mechanistic understanding of the regulatory system that governs the respective transition, but are derived from the assumption that, under certain conditions, a time-optimal transition may be evolutionary beneficial. We emphasize that an application of TOA does not necessarily imply that a time-optimal transition is the only or most important evolutionary objective. Rather, and again similar to other optimality-based methods, the solutions of TOA provide a computational “gold standard”, (Giordano et al., 2016), to which experimentally observed behavior can be compared.

Within this work, we exemplified the use of TOA by considering two prototypical applications: the doubling of a cell in a constant environment (cf. Application 2.1), as well as the time-optimal adaptation to a nutrient shift (cf. Application 2.2). Following previous works (Molenaar et al., 2009; Giordano et al., 2016; Yabo et al., 2022), the application of TOA was illustrated using a coarse-grained self-replicator model. The results illustrate the utility of TOA to generate and explore biological hypotheses.

The premise underlying the in silico experiments of our first application, cell doubling in a constant environment, was that microbial cells are not necessarily precisely adapted to the given environment, but may nonetheless have evolved a regulatory scheme that allows them to double their intracellular composition in minimal time. Based on this premise, the application of TOA gives rise to several predictions, we observe 1) complex intracellular dynamics different from solutions obtained by iterative RBA, 2) that transient accumulation of storage compounds reduces the predicted doubling time, and 3) that (mis-)adaptation to a higher nutrient availability than actually present in the environment carries a lower evolutionary cost than (mis-)adaptation to a lower nutrient availability.

Due to the simplicity of the coarse-grained model, we do not necessarily expect the specific time courses obtained for the model to be exact predictions of biological reality. In particular, we acknowledge that the coarse-grained model lacks further intracellular constraints that affect progress through the cell cycle (for example, checkpoints and a detailed representation of DNA replication and segregation) that also impact metabolic processes. Nonetheless, we are confident that the results reveal several insights that reflect biological reality. Specifically, the role of storage compounds in cellular metabolism is difficult to explore using existing constraint-based models. Here, the application of TOA demonstrates that, beyond the role of storage in diurnal oscillations, cf. (Rügen et al., 2015; Reimers et al., 2017) and as a safeguard for periods of nutrient scarcity, storage may play an important role even under constant environmental conditions. As shown with TOA, intracellular dynamics and transient accumulation of nutrients may contribute to a reduction of doubling time. Indeed, and different from typical steady-state solutions of current constraint-based methods, cells do exhibit coordinated metabolic dynamics over a cell cycle (Papagiannakis et al., 2017).

The application of TOA was further exemplified by simulations of time-optimal cellular adaptation to a nutrient shift. Similar to the results obtained for constant environments, TOA demonstrates that transient accumulation of storage can reduce the time required for adaptation–a finding supported by experimental evidence that storage compounds, such as glycogen, indeed provide short-term benefits in changing environments (Sekar et al., 2020).

In particular, the rapid uptake and storage of nutrients following an upshift in nutrient supply (as shown in Figure 7, left column) is reminiscent of “luxury uptake” or “over-compensation”. The latter phenomenon is well known (Powell et al., 2009) and occurs when cells are starved and re-exposed to a limiting nutrient, such as phosphate. “Luxury uptake” and “over-compensation” after starvation can be exploited, for example, for nutrient removal from wastewater (Powell et al., 2009). Our analysis shows that such “over-compensation” or “overshoot” phenomena are readily explained using principles of (optimal) cellular resource allocation, and do not necessarily require explanations that invoke competition between individuals to rationalize rapid nutrient uptake after starvation.

We conjecture that, while the specific trajectories of the cellular response to environmental shifts might be different under specific conditions, for example, due to additional constraints not present in the model, many of the principles revealed by TOA remain valid in more elaborate models of cellular growth transitions–and thereby provide an important reference to identify optimal vs. suboptimal behavior. Indeed, it was previously shown that growth transition kinetics of E. coli are indeed suboptimal under the studied conditions (Erickson et al., 2017)–a finding that could only be obtained by comparison to an optimal reference solution. As shown in this work, TOA can also be readily incorporated into a multi-objective framework (in the sense of Pareto) that allows us to incorporate additional objectives.

Finally, the results of TOA demonstrate that the costs of mis-adaptation to an environment are not symmetric, neither for cell doubling in a constant environment (Figure 6), nor for adaptation after a nutrient shift (Figure 8). In either case, a cell that is adapted to a higher level of (extracellular) nutrient than available in the environment (“optimist”) has only a minor disadvantage compared to an already perfectly adapted cell. Vice versa, however, cells that are adapted to a lower level of (extracellular) nutrient than available in the environment (“pessimist”) have a pronounced disadvantage compared to a perfectly adapted cell. This asymmetry indicates that adaptation to a low nutrient environment is only advantageous if the low nutrient state persists for an extended period of time. This asymmetry is supported by experimental evidence. For example, it has been suggested that some microorganisms, such as Lactococcus lactis, preserve a large overcapacity of ribosomes and glycolytic enzymes to be ready to rapidly respond and grow when conditions improve (Goel et al., 2015), and thereby implement an “optimistic strategy”.

5 Conclusions and Outlook

Constraint-based optimization plays an important role to elucidate and eventually predict cellular behavior. As an extension of previous modeling frameworks, we introduced time-optimal adaptation. TOA is motivated by the assumption that under certain conditions it is evolutionary favorable to adapt to a new cellular state in minimal time. In its general form, TOA can be applied in a very broad sense and thereby extends most of the existing constraint-based modeling frameworks.

As shown in this work, TOA allowed us to obtain insight into several biological phenomena, such as the accumulation of storage in constant environments and “overshoot” accumulation of nutrients after starvation, which cannot be readily explained using existing methods–thereby demonstrating the utility of TOA for future analysis.

While the examples discussed within this work focused on constant environments and simple nutrient shifts, TOA can also be applied in time-dependent environments and can be readily extended to include further constraints. Likewise, as shown in this work, TOA can be included within multi-objective optimization in the sense of Pareto.

Possible further extensions include “t-max adaptation”, i.e., to maximize, for example, survival time under nutrient starvation, as well as more general constraints on the target state (for example, to attain a minimal amount of a specific intermediate in minimal time, while the amounts other cellular components are not specified).

We are therefore confident that TOA and its possible extensions are a valuable contribution in the context of constraint-based modeling with manifold applications beyond the examples discussed in this work.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

All authors contributed equally to the conceptualization of the method, the writing and editing of the manuscript. MK conducted the numerical experiments. RS provided interpretations and discussion.

Funding

The work of MK was carried out during the tenure of an ERCIM “Alain Bensoussan” Fellowship of the author at the Norwegian University of Science and Technology. The work of RS is funded by the grant STE 2062/2-1 of the German Research Foundation (DFG). We acknowledge support by the German Research Foundation (DFG) and the Open Access Publication Fund of Humboldt-Universität zu Berlin.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.866676/full#supplementary-material

References

Basan, M., Honda, T., Christodoulou, D., Hörl, M., Chang, Y.-F., Leoncini, E., et al. (2020). A Universal Trade-Off between Growth and Lag in Fluctuating Environments. Nature 584, 470–474. doi:10.1038/s41586-020-2505-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnans, F., Martinon, P., and Grélard, V. (2017). BOCOP: An Open Source Toolbox for Optimal Control.

Google Scholar

Bordbar, A., Monk, J. M., King, Z. A., and Palsson, B. O. (2014). Constraint-based Models Predict Metabolic and Associated Cellular Functions. Nat. Rev. Genet. 15, 107–120. doi:10.1038/nrg3643

PubMed Abstract | CrossRef Full Text | Google Scholar

Djema, W., Bernard, O., and Giraldi, L. (2020). Separating Two Species of Microalgae in Photobioreactors in Minimal Time. J. Process Control 87, 120–129. doi:10.1016/j.jprocont.2020.01.003

CrossRef Full Text | Google Scholar

Ehrgott, M. (2000). Multicriteria Optimization. Springer Berlin Heidelberg. doi:10.1007/978-3-662-22199-0

CrossRef Full Text | Google Scholar

Erickson, D. W., Schink, S. J., Patsalo, V., Williamson, J. R., Gerland, U., and Hwa, T. (2017). A Global Resource Allocation Strategy Governs Growth Transition Kinetics of Escherichia coli. Nature 551, 119–123. doi:10.1038/nature24299

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerdts, M. (2011). Optimal Control of ODEs and DAEs. Berlin, Boston: De Gruyter.

Google Scholar

Giordano, N., Mairet, F., Gouzé, J.-L., Geiselmann, J., and de Jong, H. (2016). Dynamical Allocation of Cellular Resources as an Optimal Control Problem: Novel Insights into Microbial Growth Strategies. PLoS Comput. Biol. 12 (28), e1004802. doi:10.1371/journal.pcbi.1004802

PubMed Abstract | CrossRef Full Text | Google Scholar

Goel, A., Eckhardt, T. H., Puri, P., de Jong, A., Branco dos Santos, F., Giera, M., et al. (2015). Protein Costs Do Not Explain Evolution of Metabolic Strategies and Regulation of Ribosomal Content: Does Protein Investment Explain an Anaerobic Bacterial Crabtree Effect? Mol. Microbiol. 97, 77–92. doi:10.1111/mmi.13012

PubMed Abstract | CrossRef Full Text | Google Scholar

Goelzer, A., Fromion, V., and Scorletti, G. (2011). Cell Design in Bacteria as a Convex Optimization Problem. Automatica 47, 1210–1218. doi:10.1016/j.automatica.2011.02.038

CrossRef Full Text | Google Scholar

Gurobi (2021). Gurobi Optimizer Reference Manual. Tech. Rep. Beaverton, OR: Gurobi Optimization LLC.

Google Scholar

Hermes, H., and Lasalle, J. P. (1969). Functional Analysis and Time Optimal Control. New York: Academic Press.

Google Scholar

Höffner, K., Harwood, S. M., and Barton, P. I. (2012). A Reliable Simulator for Dynamic Flux Balance Analysis. Biotechnol. Bioeng. 110, 792–802. doi:10.1002/bit.24748

PubMed Abstract | CrossRef Full Text | Google Scholar

Höffner, K., Khan, K. A., and Barton, P. I. (2016). Generalized Derivatives of Dynamic Systems with a Linear Program Embedded. Automatica 63, 198–208. doi:10.1016/j.automatica.2015.10.026

CrossRef Full Text | Google Scholar

Jeanne, G., Goelzer, A., Tebbani, S., Dumur, D., and Fromion, V. (2018). Dynamical Resource Allocation Models for Bioreactor Optimization. IFAC-PapersOnLine 51, 20–23. doi:10.1016/j.ifacol.2018.09.020

CrossRef Full Text | Google Scholar

Klipp, E., Heinrich, R., and Holzhütter, H.-G. (2002). Prediction of Temporal Gene Expression. Eur. J. Biochem. 269, 5406–5413. doi:10.1046/j.1432-1033.2002.03223.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerman, J. A., Hyduke, D. R., Latif, H., Portnoy, V. A., Lewis, N. E., Orth, J. D., et al. (2012). In Silico method for Modelling Metabolism and Gene Product Expression at Genome Scale. Nat. Commun. 3. doi:10.1038/ncomms1928

CrossRef Full Text | Google Scholar

Liu, L., and Bockmayr, A. (2020). Regulatory Dynamic Enzyme-Cost Flux Balance Analysis: A Unifying Framework for Constraint-Based Modeling. J. Theor. Biol. 501, 110317. doi:10.1016/j.jtbi.2020.110317

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L. (2020). Unifying Metabolic Networks, Regulatory Constraints, and Resource Allocation. Berlin, Germany: Freie Universität Berlin. Ph.D. thesis. doi:10.17169/REFUBIUM-27061

CrossRef Full Text | Google Scholar

Mahadevan, R., Edwards, J. S., and Doyle, F. J. (2002). Dynamic Flux Balance Analysis of Diauxic Growth in Escherichia coli. Biophysical J. 83, 1331–1340. doi:10.1016/s0006-3495(02)73903-9

CrossRef Full Text | Google Scholar

Mahadevan, R., and Schilling, C. H. (2003). The Effects of Alternate Optimal Solutions in Constraint-Based Genome-Scale Metabolic Models. Metab. Eng. 5, 264–276. doi:10.1016/j.ymben.2003.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Molenaar, D., van Berlo, R., de Ridder, D., and Teusink, B. (2009). Shifts in Growth Strategies Reflect Tradeoffs in Cellular Economics. Mol. Syst. Biol. 5, 323. doi:10.1038/msb.2009.82

PubMed Abstract | CrossRef Full Text | Google Scholar

Orth, J. D., Thiele, I., and Palsson, B. Ø. (2010). What Is Flux Balance Analysis? Nat. Biotechnol. 28, 245–248. doi:10.1038/nbt.1614

PubMed Abstract | CrossRef Full Text | Google Scholar

Papagiannakis, A., Niebel, B., Wit, E. C., and Heinemann, M. (2017). Autonomous Metabolic Oscillations Robustly Gate the Early and Late Cell Cycle. Mol. Cell 65, 285–295. doi:10.1016/j.molcel.2016.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavlov, M. Y., and Ehrenberg, M. (2013). Optimal Control of Gene Expression for Fast Proteome Adaptation to Environmental Change. Proc. Natl. Acad. Sci. U.S.A. 110, 20527–20532. doi:10.1073/pnas.1309356110

PubMed Abstract | CrossRef Full Text | Google Scholar

Powell, E. O. (1956). Growth Rate and Generation Time of Bacteria, with Special Reference to Continuous Culture. J. General Microbiol. 15, 492–511. doi:10.1099/00221287-15-3-492

PubMed Abstract | CrossRef Full Text | Google Scholar

Powell, N., Shilton, A., Chisti, Y., and Pratt, S. (2009). Towards a Luxury Uptake Process via Microalgae - Defining the Polyphosphate Dynamics. Water Res. 43, 4207–4213. doi:10.1016/j.watres.2009.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimers, A. M., Knoop, H., Bockmayr, A., and Steuer, R. (2017). Cellular Trade-Offs and Optimal Resource Allocation during Cyanobacterial Diurnal Growth. Proc. Natl. Acad. Sci. U. S. A. 114, E6457–E6465. doi:10.1073/pnas.1617508114

PubMed Abstract | CrossRef Full Text | Google Scholar

Rügen, M., Bockmayr, A., and Steuer, R. (2015). Elucidating Temporal Resource Allocation and Diurnal Dynamics in Phototrophic Metabolism Using Conditional FBA. Sci. Rep. 5 (16), 15247. doi:10.1038/srep15247

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, M., and Hwa, T. (2011). Bacterial Growth Laws and Their Applications. Curr. Opin. Biotechnol. 22, 559–565. doi:10.1016/j.copbio.2011.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Segel, L. A., and Slemrod, M. (1989). The Quasi-Steady-State Assumption: A Case Study in Perturbation. SIAM Rev. 31, 446–477. doi:10.1137/1031091

CrossRef Full Text | Google Scholar

Sekar, K., Linker, S. M., Nguyen, J., Grünhagen, A., Stocker, R., and Sauer, U. (2020). Bacterial Glycogen Provides Short-Term Benefits in Changing Environments. Appl. Environ. Microbiol. 86. doi:10.1128/aem.00049-20

PubMed Abstract | CrossRef Full Text | Google Scholar

Varma, A., and Palsson, B. O. (1994). Stoichiometric Flux Balance Models Quantitatively Predict Growth and Metabolic By-Product Secretion in Wild-type Escherichia coli W3110. Appl. Environ. Microbiol. 60, 3724–3731. doi:10.1128/aem.60.10.3724-3731.1994

PubMed Abstract | CrossRef Full Text | Google Scholar

Waldherr, S., Oyarzún, D. A., and Bockmayr, A. (2015). Dynamic Optimization of Metabolic Networks Coupled with Gene Expression. J. Theor. Biol. 365, 469–485. doi:10.1016/j.jtbi.2014.10.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Yabo, A. G., Caillau, J.-B., Gouzé, J.-L., de Jong, H., and Mairet, F. (2022). Dynamical Analysis and Optimization of a Generalized Resource Allocation Model of Microbial Growth. SIAM J. Appl. Dyn. Syst. 21, 137–165. doi:10.1137/21m141097x

CrossRef Full Text | Google Scholar

Yang, L., Ebrahim, A., Lloyd, C. J., Saunders, M. A., and Palsson, B. O. (2019). DynamicME: Dynamic Simulation and Refinement of Integrated Models of Metabolism and Protein Expression. BMC Syst. Biol. 13, 2. doi:10.1186/s12918-018-0675-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Yegorov, I., Mairet, F., and Gouzé, J.-L. (2018). Optimal Feedback Strategies for Bacterial Growth with Degradation, Recycling, and Effect of Temperature. Optim. Control Appl. Meth 39, 1084–1109. doi:10.1002/oca.2398

CrossRef Full Text | Google Scholar