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 y : [ t 0 , t end ] → R n 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 y ̇ ( t ) = d d t 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 c t = y t b i o t , (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 w ∈ R n 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. y ̇ t = S ⋅ v t , (2.2a) where S ∈ R n y × n v denotes the stoichiometric matrix and v : [ t 0 , t end ] → R n 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 S = S y S x is split up such that “fast” and “slow” intracellular compounds, usually metabolites resp. macromolecules, are described separately and (2.2a) is replaced by. y ̇ t = S y ⋅ v t , 0 = S x ⋅ v t , (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.

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. for a l m o s t a l l t : y ̇ t , y t , u t ∈ A t , (2.3a) where the set A ( t ) ⊆ R n y × R n y × R n u is typically defined through (in-)equalities such as steady-state assumptions and/or positivity requirements. The particular form of the set 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 A ( t ) . The vector-valued function u = u ( t ) : [ t 0 , t end ] → R n 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: ϕ bndry y t 0 , y t end , u t 0 , u t end ≤ 0 , (2.3b) min y ⋅ , u ⋅ f y , u (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 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 l b , u b ∈ R n v , which might be dependent on the time t: y ̇ t = S y ⋅ v t , dynamics, often just biomass 0 = S x ⋅ v t , quasi steady - state 0 ≤ y t , positivity l b t ≤ v t ≤ u b t , flux bounds, dependent on environment with given initial conditions y t 0 = y 0 ∈ R n y . The flux rates are determined through optimization of a linear functional (often the flux through the biomass reaction, assembled in a vector w obj ∈ R n v ) v t ∈ argmin v w obj ⊤ ⋅ v .

The quantities in (2.3) can be identified as: A t = y ̇ , y , u : y ̇ = S y ⋅ u , 0 = S x ⋅ u , 0 ≤ y , l b t ≤ u ≤ u b t , u ∈ argmin v w obj ⊤ ⋅ v , ϕ bndry = y t 0 − y 0 y 0 − y t 0 , while typically no additional (global) objective function is present. Note that the defining condition on the fluxes u ∈ argmin v w obj ⊤ ⋅ 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 H y t ⋅ y t + H v t ⋅ v t ≤ h t , (2.4) with H y ∈ R n h × n y , H v ∈ R n h × n v . (2.5) The model is usually formulated as an initial-value problem y t 0 = y 0 ∈ R n 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), min y , u ∫ t 0 t e n d ϕ y ⊤ t ⋅ y t d t + ϕ end ⊤ ⋅ y t end ,  with  ϕ y , ϕ end ∈ R n 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 u t = v t , A t = y ̇ , y , u : y ̇ = S y ⋅ u , 0 = S x ⋅ u , 0 ≤ y l b t ≤ u ≤ u b t , h t ≥ H y t ⋅ y + H u t ⋅ u ϕ bndry = y t 0 − y 0 y 0 − y t 0 , f y ⋅ , u ⋅ = ∫ t 0 t end ϕ y ⊤ t ⋅ y t d t + ϕ end ⊤ ⋅ y t end .

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: c t 0 = 1 b i o t 0 ⋅ y t 0 = 1 b i o t end ⋅ y t end = c t end . (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 u t = u 1 u 2 : = α v t , α ∈ R , A t = y ̇ , y , u : y ̇ = S y ⋅ u 2 : , 0 = S x ⋅ u 2 : , 0 ≤ y , l b t ≤ u 2 : ≤ u b t , h t ≥ H y t ⋅ y + H v t ⋅ u 2 : , ϕ bndry = u 1 ⋅ y t 0 − y t end y t end − u 1 ⋅ y t 0 , f y ⋅ , u ⋅ = − u 1 , 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 u t = α v 0 = u 1 u 2 : , A t = y ̇ , y , u : y ̇ = S y ⋅ v 0 = S x ⋅ v 0 ≤ y l b t ≤ v ≤ u b t h t ≥ H y t ⋅ y + H v t ⋅ v v = u 2 : ⋅ exp λ ⋅ t − t 0 λ = ln u 1 / t end − t 0 , ϕ bndry = y t 0 − y 0 y 0 − y t 0 , f y t , u = − u 1 . 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.

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

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 y goal = 2 ⋅ y 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), y t 0 = y init , c T = 1 b i o T ⋅ y T = 1 w ⊤ ⋅ y goal ⋅ y goal = c goal . (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: min y ⋅ , u ⋅ , T > 0 T − α s.t.  y 0 = y init y T = α ⋅ y goal and   y ̇ t , y t , u t ∈ A t for almost all t ∈ 0 , T , (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 b i o t end ≥ b i o t 0  or  w ⊤ ⋅ y t end ≥ w ⊤ ⋅ y t 0 can be added, illustrating again the power of constraint-based modeling. Whenever necessary, this was done for the in silico experiments in Section 3.

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

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 y ( t ̃ ) for all time instances 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 : min u ⋅ , y ⋅ ± ∫ t 0 T y i t 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.

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.

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.

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.

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

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.

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

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.

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 m e a n M = 1 T ċ ∫ 0 T M t w ⊤ ⋅ y t d t , (3.1) where y (t) is a solution obtained by relative TOA-VA, cf. Section 2.6.

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.

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.

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.

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.

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 * ≈ 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 * ≈ 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.