5. differentiability: The optimal value of the intracellular optimization problem may not be continuously differentiable, which poses an obstacle when dFBA is used for optimal control or parameter estimation.

6. non-unique solutions: The solution of the intracellular optimization problem is usually not unique which can make fluxes unrealistically “jump” between different optimal solutions.

**Static Optimization Approach (SOA)** divides the total time horizon of the dFBA simulation into several smaller time intervals. The optimization problem is solved to obtain the flux distribution at the beginning of each time interval, and then the ODEs are integrated over the time interval with this fixed flux distribution. The dynamics calculated from this time step are used to constrain the optimization problem solved at the beginning of the next time interval, and the process is repeated until the end of the simulation time is reached. SOA can be implemented easily with the use of an Euler scheme for integrating the system and a suitable existing LP solver for solving the FBA at each time step. SOA is also implemented in the constraint-based reconstruction and analysis (COBRA) toolbox for MATLAB [338] which can perform dFBA simulations. Since its implementation is relatively simple, SOA has been widely used in studies for the diauxic [329, 331], aerobic and fermentative [328, 322] growth of *E. coli*, for *S. cerevisiae* fermentations [333, 326, 327], as well as for the metabolism underlying plant growth [334]. Many of these applications include larger-scale or genome-scale metabolic networks, due to the scalability of SOA. However, the main drawback is that SOA is inefficient and can become computationally expensive because it has to solve the optimization problem at each time step. This can be challenging for most dFBA problems which are stiff and require small time steps to ensure accuracy, convergence, and stability of the solution.

**Dynamic Optimization Approach (DOA)** follows closely the general dynamic optimization framework described in Section 11.2: an objective function that depends on the dynamic states of the system over the complete time horizon of interest is formulated, and the dynamics (11.26)–(11.27) and algebraic constraints (11.28)–(11.29) are added as optimization constraints. In other words, DOA discretizes the total time horizon of the dFBA simulation, and then transforms the dynamic optimization problem into a non-linear programming (NLP) problem, which is solved once by simultaneously optimizing over the entire time of the simulation. In this way, DOA obtains the time profiles of fluxes and metabolite concentrations in the system, and allows the formulation of a dynamic objective function, which could provide useful information about the design of genetically modified metabolic networks or the maximization of bioprocess productivity. Because of this characteristic, DOA is often used in dynamic metabolic engineering, parameter estimation and optimal control applications. For example, DOA has been used for simulating the diauxic growth of *E. coli* [329, 339]. On the downside, even though the optimization problem does not need to be repeatedly solved like in SOA, the single NLP of DOA can become easily intractable, as its dimension increases with the fineness of time discretization. Additionally, DOA has been mainly limited to small-scale metabolic networks, since it cannot be easily applied to genome-scale metabolic networks due to the large number of variables and constraints that are introduced in the NLP as the size of the network increases.

**Direct Approach (DA)** has been formulated more recently than SOA and DOA, and directly includes the LP solver for the FBA in the right-hand side evaluator function of the ODEs. In this way, it can take advantage of existing ODE integrators with adaptive step size and error control that can reduce the number of integration steps and provide better solution accuracy compared to SOA. DA has been implemented in the ORCA toolbox [340], which complements the constraint-based reconstruction and analysis (COBRA) toolbox for MATLAB [338]. Furthermore, DA has been used for studying the (an)aerobic growth of wild-type and engineered *E. coli* strains [330], the aerobic growth of *Corynebacterium glutamicum* on glucose and xylose in biorefinery simulations [332], as well as the aerobic and anaerobic growth of wild type and engineered *S. cerevisiae* strains [324, 325, 326]. Some of these applications involve dynamic metabolic engineering for product maximization, and many of them include genome-scale metabolic networks, since DA is relatively easily scalable like SOA.

However, DA requires the LP to be resolved at least once, every time the right-hand side of the ODEs is evaluated [341]. This can make DA computationally demanding, especially for larger metabolic networks. Another major challenge is that when evaluating the right-hand side of the ODEs close to the boundary of feasibility, the LP can become infeasible and make the dFBA simulation fail. The LP can become infeasible
Jun S, Si F, Pugatch R, Scott M. Fundamental principles in bacterial physiology-history, recent progress, and the future with focus on cell size control: a review. Rep Prog Phys. 2018 May;81(5):056601. doi: 10.1088/1361-6633/aaa628. Epub 2018 Jan 9. PMID: 29313526; PMCID: PMC5897229.

Serbanescu D, Ojkic N, Banerjee S. Cellular resource allocation strategies for cell size and shape control in bacteria. FEBS J. 2021 Oct 19. doi: 10.1111/febs.16234. Epub ahead of print. PMID: 34665933.

Amir A, Männik J, Woldringh CL, Zaritsky A. Editorial: The Bacterial Cell: Coupling between Growth, Nucleoid Replication, Cell Division, and Shape Volume 2. Front Microbiol. 2019 Sep 4;10:2056. doi: 10.3389/fmicb.2019.02056. PMID: 31551989; PMCID: PMC6737747.

Kleckner NE, Chatzi K, White MA, Fisher JK, Stouf M. Coordination of Growth, Chromosome Replication/Segregation, and Cell Division in E. coli. Front Microbiol. 2018 Jul 9;9:1469. doi: 10.3389/fmicb.2018.01469. PMID: 30038602; PMCID: PMC6046412.
$\Sigma$ Math box 11.B DFBAlab

In order to address some of the computational challenges of dFBA, Höffner, Harwood, and Barton proposed a simulator for dFBA, which was initially coded in FORTRAN [343], but gained popularity when implemented in MATLAB with the name Dynamic Flux Balance Analysis laboratory (DFBAlab) [337], and more recently in Python [344]. It is noted that the DFBAlab is compatible with the COBRA toolbox [338]. Based on this dFBA simulator, it is not necessary to resolve the LP each time the right-hand side of the ODEs is evaluated and consequently, the solution process becomes faster. This is possible because the FBA solution at an initial time can be used to compute future optimal solutions by detecting changes in the active set or by computing the optimal basis of the FBA solution [345]. Unfortunately, such formulations need to continuously monitor the active set of the LP, which increases with the size of the metabolic network, or need to choose a basis for the optimal solution that is most likely to remain optimal as the simulation proceeds [345]. The latter is challenging since the optimal basis can be non-unique even for a unique optimal solution. Nevertheless, DFBAlab manages to reduce the number of times that the LP is resolved, and also avoids obtaining infeasible LPs and numerical failure by using the LP feasibility problem and the Karush-Kuhn-Tucker (KKT) optimality conditions of the FBA problem (see below). In addition, the differentiability problem could be solved with the help of non-smooth analysis which provides optimality conditions in terms of sub-gradients or generalized gradients, for convex and non-convex functions respectively. Furthermore, to tackle the issue of primal multiplicity of the FBA problem, DFBAlab performs lexicographic optimization.

either because it is really infeasible and the simulation should be terminated, or because the ODE integrator becomes unable to evaluate the right-hand side of the ODEs and the simulation is discontinued, or erroneous death phase messages are being displayed. The latter can happen as dFBA simulations involve discrete events that correspond to switches in the active set of the LP solution. More specifically, different bases for the optimal solution of the LP can emerge at each time step. Moreover, at the points of change of the active set, the dFBA model is not differentiable, since the optimal value of the LP as a function of the right-hand side of the constraints is not continuously differentiable. This is a problem because the first and second derivatives of the model must be computed when dFBA is used for optimal control or parameter estimation applications. Finally, another drawback emerges due to the primal multiplicity of the LP. As it is well-known, FBA is formulated as an underdetermined problem and therefore, the LP does not have a unique solution [342]. Non-unique optimal reaction fluxes can lead different ODE integrators to different results.

$\Sigma$ Math box 11.C Lexicographic Optimization

Lexicographic or hierarchical optimization involves the solution of a series of LPs with auxiliary objectives ranked in priority order. The use of auxiliary objectives reduces the feasible space and leads to a unique optimal solution, while the auxiliary objectives can have specific biological meanings and can be selected based on prior knowledge about the organism [297]. Some of the most-used auxiliary objectives are based on the assumption that evolution leads to the exclusion of inefficient pathways so that cells can biosynthesize the smaller possible number of enzymes. Examples of such auxiliary objectives include the minimization of enzyme cost [346], the minimization of the total reaction flux [347], and the minimization of the number of active reactions [348]. However, it has been shown that such objectives may not be suitable for some engineered cells. In general, it is not trivial to find a series of auxiliary objectives that are consistent with experimental data, assure uniqueness and preserve continuity of the optimal solution. In some cases, even when all auxiliary objectives have been used, hierarchical optimization cannot ensure the uniqueness of the dFBA solution. Apart from the use of auxiliary objectives, auxiliary rules or auxiliary parameters have also been proposed to address the primal multiplicity of FBA. For example, geometric methods have been proposed to identify a unique distribution of reaction fluxes for FBA [349], even though there is no biological evidence to justify such methods.

## 11.6 Concluding remarks

One approach for understanding the response of microorganisms to changes in their environment is to assume that this response has been optimized by evolution. That is, the regulatory mechanisms control-
<figure><img src="image_111.png" /><figcaption>Figure 11.4: Example result of a dFBA simulation for the E. coli core model [350] performed with DFBAlab (A) Illustration of dynamic metabolites and relevant exchange fluxes. Extracellular metabolites use mass concentration (g L⁻¹), exchange fluxes are in molar amount per dry biomass and time (mmol gDW⁻¹ h⁻¹). (B) Concentration-dependent constraints applied to the exchange fluxes during the dFBA simulation. (C) Differential equation model for biomass concentration X and metabolite mass concentrations. Growth rate μ and exchange fluxes v* (in red) are optimal values from the underlying FBA model. mO₂, mG, mE, mA are molar masses of O₂, glucose, ethanol, and acetate, respectively. kLa = 8.5 h⁻¹ is the volumetric mass transfer coefficient for oxygen. (D) Simulation results for concentrations of biomass and metabolites. We can observe four growth phases: aerobic growth on glucose with production of acetate, anaerobic growth on glucose with production of ethanol and acetate, aerobic growth on acetate, and a stationary phase. (E) Simulation result for oxygen concentration in the liquid medium. Oxygen is depleted in the second growth phase due to mass transfer limitations, but replenishes at the start of the third phase. (F) Penalty function time course. Increases of the penalty function indicate periods where the underlying FBA model is infeasible. Here, this occurs when the ATP maintenance constraint cannot be satisfied due to a lack of substrates, and happens in this simulation during a brief period where the switch from glucose to acetate as a substrate takes place, since oxygen needs to be replenished first, and in the stationary phase.</figcaption></figure>

ling the response optimize an objective, or a trade-off between competing objectives, subject to a variety of physical and biochemical constraints. This approach gives rise to dynamic optimization problems (11.1)-(11.2) that can be solved by techniques from optimal control theory. Three examples of such problems were considered in this chapter: dynamic optimization of enzyme expression in metabolic pathways, dynamic optimization of coarse-grained models of cellular growth, and dynamic flux balance analysis. This does not exhaust the range of possible problems that can be considered. One example is the combination of the resource allocation perspective with dynamic flux balance analysis [339].

Some of the predictions obtained by means of dynamic optimization seem to be supported by available experimental data, such as the time-ordering of enzyme expression in a linear pathway. Other predictions cannot currently be tested or may not be consistent with the available experimental data, such as the dynamical ribosomal protein synthesis pattern. The contradiction between a predicted optimal response and the observations is interesting, because it indicates that some of the assumptions underlying the problem need to be revised. The model may not account for all relevant processes taking place in the cell, important constraints may have been ignored, or the objective may not capture the actual processes taking place.
Chapter 13

Cell behavior in the face of uncertainty

David Lacoste, Olivier Rivoire, and David Tourigny

Chapter highlights

* Organisms that grow and survive in uncertain environments may need to change their physiological state as the environment changes.

* When the environment is uncertain, one strategy known as bet-hedging is to make these changes randomly and independently of the environment, to ensure that at least part of the population is well adapted.

* Organisms that collect information from their environment may also use this information to modulate their changes of physiological states.

* We review these different strategies and point out parallels with the theory of optimal financial investments.

To a large extent, the content of this textbook prior to the current chapter has dealt with models of microorganisms under the implicit assumption that the dynamics of both environmental factors and intracellular components are deterministic, and that behavior is optimized uniformly across cells in a population. On longer time scales however, natural selection also acts on populations and these populations may encounter environments that fluctuate across both time and space. Under these conditions, natural selection may not favor a homogeneous deterministic cellular response across the population, but rather select for a certain level of population diversity and heterogeneity, including behaviors arising from mechanisms that are fundamentally stochastic. Stochasticity is inherent to intracellular processes such as gene expression and signal transduction due to the small number of molecules that they involve. It is often referred to as “noise”, but this terminology can be misleading because stochasticity may also fulfill an essential role in cellular function and survival, for example during growth in uncertain environmental conditions. The purpose of this chapter is to highlight this role, introduce the mathematical models necessary for understanding it, and draw a new economic analogy with problems of investment in finance.

Before expanding upon the role that uncertainty plays in shaping cellular behavior, we briefly point out some general limitations of deterministic models based on optimal regulation of behavior in time as described in Chapter 10. In that chapter, it was assumed that microorganisms have evolved, under selective pressures exerted by the environment, to optimize a specific objective criterion or combination of objective criteria that were shared by all cells of a population. This assumption was then incorporated into an optimal control framework to explain how cellular behavior (e.g., enzyme expression) is optimally regulated in time depending on deterministic interactions between a microbial population and its environment. In particular, we consider optimal control strategies across a prescribed time window. Defining optimality in such case assumes the organism has perfect information on how the environment will change (including in response to actions taken) over time. In an uncertain environment, this information is simply not available. An alternative is instantaneous optimization of growth rate at each time point but this is a shortsighted strategy that excludes any partial
The solution of a dynamic optimization problem may be different from the concatenation of the solutions of repeated static optimization problems defined over short, consecutive time intervals making up the time horizon. For instance, the optimal pattern of resource allocation over a time horizon may involve the accumulation of a reserve of unused resources that, while being wasteful in the short run, is beneficial when the whole time interval of interest is considered. One example was given in the introduction of this chapter, concerning the expression of maltose enzymes in the presence of lactose [293]. Another example is the accumulation of glycogen in cyanobacteria during daylight, providing the energetic resources for maintenance metabolism in the night time [290].

The analysis of the growth of microorganisms using dynamic optimization critically depends on the choice of an appropriate objective function. In the examples discussed above, optimization was performed with respect to a single objective, e.g., the maximal accumulation of biomass over a time interval or the minimal time to deplete a given amount of substrate. It is plausible, however, that microorganisms have evolved under the necessity to satisfy several objectives simultaneously. This can be taken into account by formulating a weighted sum of the different objective functions, such as the simultaneous minimization of throughput time and investment in enzymes in Eq. (11.16). Another approach is to generalize the optimization problem to a multi-objective optimization problem, with sets of Pareto optimal solutions, each providing a trade-off between mutually conflicting objectives. One example of such a multi-objective optimization problem is given by a generalization of the dynamic optimization of enzyme expression in metabolic pathways in Section 11.3, with the double objective of minimizing the time to consume a given amount of substrate and minimizing the concentration of (possibly toxic) intermediate metabolites [351].

Instead of making *a-priori* assumptions about the objectives presumably optimized by microorganisms, one could try to infer the latter from the experimental data. This inverse optimization approach leverages the large amounts of time-course data on the dynamic response of microorganisms to environmental perturbations that have accumulated in the past decade. Inverse optimization requires the solution of complex inverse optimal control problems that have been little explored until now [296].

## Recommended readings

### General references introducing the problem of dynamic optimization in biology:

* Sutherland, The best solution, *Nature*, 2005 [352]. Influential statement of the use of optimization as an explanatory and predictive principle in biology.

* Banga, Optimization in computational systems biology, *BMC Systems Biology*, 2008 [353]. General review of optimization approaches in systems biology.

* Heinrich, Schuster, and Holzhütter, Mathematical analysis of enzymic reaction systems using optimization principles, *European Journal of Biochemistry*, 1991 [291]. Classical review on the use of dynamic optimization in the analysis of biochemical systems.

* Ewald, Bartl, and Kaleta. Deciphering the regulation of metabolism with dynamic optimization: an overview of recent advances. *Biochemical Society Transactions*, 2017 [304]. Recent review on the use of dynamic optimization in the analysis of metabolism.

### Methodological foundations of dynamic optimization and optimal control:

* Gerdts, *Optimal Control of ODEs and DAEs*, De Gruyter, 2011 [294]. Introduction to optimal control theory.

* Kirk, *Optimal Control Theory: An Introduction*, Courier Corporation, 2004 [295].

### Dynamic optimization of enzyme expression in metabolic pathways:

* Klipp, Heinrich, and Holzhütter. Prediction of temporal gene expression: metabolic optimization by re-distribution of enzyme activities. *European Journal of Biochemistry*, 2002 [301]. Classical paper on dynamic optimization of enzyme expression in unbranched pathways.

* Oyarzún, Ingalls, Middleton, and Kalamatianos. Sequential activation of metabolic pathways: a dynamic optimization approach. *Bulletin of Mathematical Biology*, 2009 [305]. Recent example revisiting the problem of dynamic optimization of enzyme expression.

### Dynamic optimization of resource allocation in coarse-grained models of cellular growth:
information on future environmental states that the organism may have acquired over the course of evolution. Such deterministic models may be suitable for deterministically changing environments, but cannot account for stochastic behaviors that may be advantageous to population growth in uncertain environments.

In this chapter, it will be shown how principles of optimality can be formulated to study the behavior of organisms growing under uncertainty. Unlike the deterministic setting however, optimality will instead need to be defined in terms of probabilities and expected returns. Analogous to the general unification of deterministic models for cellular behavior using an optimal control theory framework, models including uncertainty are unified by the subject of *stochastic optimal control*. Beyond biology, this subject has wide-reaching applications to engineering but the most relevant analogy is with finance where stochastic strategies of portfolios diversification mirror stochastic strategies of cellular diversification. This will add a new economic analogy to the economic analogies of previous chapters.

## 13.1. Strategies to cope with uncertainty: a financial analogy

We will use the topic of bacterial persistence as a recurring example throughout this chapter (Figure 13.1). When a clonal population of bacteria is exposed to an antibiotic, not all cells within the population are killed – a small sub-population, although genetically identical to the rest, may nevertheless be in a distinct phenotypic state that is growth-dormant and resistant to treatment (Figure 13.1A). While the peers of this dormant sub-population previously grew well in the absence of antibiotic, upon exposure to treatment these growing cells are killed, and only the dormant cells (the persisters) remain alive. In turn, when the remaining persisters are transferred to an environment without antibiotic a large fraction is able to revert to the growing state, allowing the population as a whole to survive. Remarkably, in this subsequent phase of growth roughly the same small fraction of persisters is retained as before the treatment. Deterministic models based on short-term optimal growth cannot explain how part of a population adopts a slow-growing state: they would predict that each cell should adopt the growing phenotype in absence of antibiotics. Cells could have a mechanism to detect the presence of unfavorable environmental conditions and adopt the persister phenotype as a response, but there are several experimental observations not explained by such a mechanism [518]: (1) a fraction of persisters exists prior to antibiotic treatment; and (2) not all cells, although genetically identical, adopt the persister phenotype. We will see that a more parsimonious description of persistence involves an optimization of long-term rather than short-term growth, which differs when environmental conditions fluctuate.

Bacterial persistence is an example of *bet-hedging*, which more generally refers to the benefit of spreading resources across multiple behavioral phenotypes to reduce the risk associated with investing all resources into any single phenotype. Returning to the example of bacterial persistence, a natural question one may ask is: what determines the precise fraction of persister cells (risk-avoiding, potentially low-reward phenotype) compared to growing cells (risky, potentially high-reward phenotype) within a given population? This question echoes a central question in financial investment: how should investors diversify their portfolio to maximize their capital in the context of uncertain returns? We will see that some of the same mathematical arguments of optimality under uncertainty can be used to analyze these two problems, showing how the optimal fraction of persisters is expected to depend critically on the probability to experience different environmental states. The terms of the analogy are presented in Table 13.1.

A pure bet-hedging strategy assumes the absence of any direct information on the current environmental state. Biologically, cells may sense signals or cues that encode varying degrees of information on their current environment. For instance, in some populations, a larger proportion of persisters is found in nutrient-poor environments compared to nutrient-rich, implying a direct relationship between shifts in environment and switches between phenotypes. These sensing or signaling mechanisms can come with associated costs however, imparted by the investment of cellular resources in, for example, the gene expression machinery. Thus, optimal cellular behavior in the face of uncertainty may be expected to involve a trade-off between stochastic (e.g., bet-hedging) and deterministic (e.g., signaling) mechanisms that balance benefit to cost in a manner that depends on evolutionary context. Other trade-offs may also exist regarding reward versus risk associated with a particular cellular decision making strategy. Analogously, financial investors face trade-offs when using incomplete information on the current state of the market and developing an investment strategy based on
* Pavlov and Ehrenberg. Optimal control of gene expression for fast proteome adaptation to environmental change. *Proceedings of the National Academy of Sciences USA*, 2013 [309]. Example of the use of optimal control to understand temporal adaptation of gene expression in response to a change in the environment.

* Giordano, Mairet, Gouzé, Geiselmann, and de Jong. Dynamical allocation of cellular resources as an optimal control problem: Novel insights into microbial growth strategies. *PLoS Computational Biology*, 2016 [306]. Example of dynamic optimization applied to a coarse-grained model of microbial growth, including a review of related work.

**Dynamic flux balance analysis of metabolic networks:**

* Mahadevan, Edwards, and Doyle. Dynamic flux balance analysis of diauxic growth in *Escherichia coli*, *Biophysical Journal*, 2002 [329]. Classical paper proposing two different approaches for dynamic flux balance analysis, with an example of diauxic growth in *E. coli*.

* Hjersted and Henson. Steady-state and dynamic flux balance analysis of ethanol production by *Saccharomyces cerevisiae*, *IET Systems Biology*, 2009 [325]. Comparison of steady-state and dynamic flux balance analysis of variably-sized metabolic models for screening metabolic engineering strategies affecting ethanol production by *Saccharomyces cerevisiae*.

* Köbis, Bockmayr and Steuer. Time-optimal adaptation in metabolic network models, *Frontiers in Molecular Biosciences*, 2022 [354]. Overview of various approaches to dynamic optimization focused on time optimality.

## Problems

Computer exercises for this chapter can be found on the book website.

### Problem 11.1 Interpretation of abrupt switches

The solutions of some of the optimization problems in section 11.3 exhibit abrupt switches in enzyme concentrations ("bang-bang behavior"), which in reality are not possible.

(a) What would be possible adjustments of the models to make the predictions more realistic? (Hint: consider explicit modeling of enzyme synthesis and enzyme degradation or dilution by growth.)

(b) What predicted behavior would you expect for these modified models?

(c) In these extended models, which possibilities would the cell have to speed up the adaptation of enzyme concentrations? Under what circumstances could this provide an actual advantage?

### Problem 11.2 Accounting for metabolites

The resource allocation model in Section 11.4 defines biomass as being composed of proteins only, neglecting notably the (small) contribution of metabolites. This has the disadvantage of putting no constraints on metabolite concentrations, which is not realistic from a biological point of view.

(a) What would be a possible adjustment of the model to integrate metabolites into the biomass composition, under the assumption that the total biomass density remains constant? (Hint: consider the definition of growth rate in Eq. (11.24).)

(b) How would the objective function for this model (Eq. (11.25)) need to be adapted accordingly?

(c) How would you expect the resulting constraint on metabolite concentrations to affect the predicted behavior of the microbial self-replicator?
<figure><img src="image_82.png" /><figcaption>Figure 13.1: Bacterial persistence as an example of a cellular strategy to cope with uncertainty in environmental conditions. **A)** Cells in a genetically identical population can display one of two distinct phenotypes that are associated with growth (pink cells) or dormancy (blue cells) in the absence of antibiotics. Only the dormant cells survive (persist) when exposed to antibiotics, and can transition back to the growth phenotype so that the population as a whole resumes growth in the absence of antibiotic. **B)** In this simplified model of bacterial persistence, the strategy *u* over two responses (phenotypes) $R_{growth}, R_{dormant}$ depends on environmental states $E_{low}$ and $E_{high}$, corresponding to low and high levels of the antibiotic, respectively. The occurrence of the states $E_{low}$ and $E_{high}$ is governed by probabilities $p(E_{low})$ and $p(E_{high})$, respectively. **C)** The multiplicative rates $f(R|E)$ associated with phenotypes $R_{growth}, R_{dormant}$ depend on environmental conditions, so that $f(R|E)$ can be represented in matrix form. The resulting optimal strategy $u_d$ corresponding to the fraction of dormant cells in the population in turn depends on the probabilities of the environmental state *E*. An analogy with Kelly betting is illustrated on the right-hand side, where the probabilities of a horse winning a race, the odds provided by a bookmaker and the optimal betting strategy are identified with $p(E)$, $f(R|E)$ and $u(R|E)$, respectively, as displayed in Table 13.1.</figcaption></figure>

the level of risk they are willing to incur.

13.2. Modeling cells growing in uncertain environments

We begin with a simple model of persistence before introducing a more general framework. This simple model assumes that bacterial cells experience an alternation of low and high antibiotics environments and can adopt two physiological states, growing or dormant (Fig. 13.1). The dormant cells are unable to replicate but persist in either high- or low-antibiotics environments while growing cells always divide when antibiotics are low in concentration but die when they are high. Mathematically, this is described by $f(R, E)$, the number of descendants of a cell with phenotype $R$ in environment $E$: $f(R = \text{dormant}, E = \text{low}) = f(R = \text{dormant}, E = \text{high}) = 1$, while $f(R = \text{growing}, E = \text{low}) = 2$ and $f(R = \text{growth}, E = \text{high}) = 0$. In absence of sensing mechanism, we consider that the fraction of dormant cells, $u_d \equiv u(R = \text{dormant})$, is a fixed quantity that only possibly evolves on very long time scales. The population thus grows by a global factor $A_{high} = f(R = \text{dormant}, E = \text{low})u_d$ if the environment is high antibiotics and by a factor $A_{low} = f(R = \text{dormant}, E = \text{low})u_d + 2f(R = \text{growing}, E = \text{low})(1 - u_d)$ if it is low antibiotics. Finally, the environment fluctuates randomly, with a probability $p_a$ to have high antibiotics and a probability $1 - p_a$ to have low antibiotics. Over a large number $T$ of generations, a population therefore experiences in average $p_a T$ periods of high
Chapter 12

Diversity of metabolic fluxes in a cell
population

Andrea De Martino and Marcelo Rivas-Astroza

Chapter overview

Even in clonal populations, cells appear to be strongly heterogeneous in terms of, e.g., protein levels, RNA levels, sizes at birth or division, interdivision times and elongation rates. Part of this variability is likely due to the inherent stochasticity of gene expression at the level of single cells. It is how- ever known that heterogeneous populations may possess an evolutionary advantage, for instance in variable environments or under stress. Despite appearing to be at odds with the idea of optimal- ity presented in the previous chapters, metabolic diversity can be described and modeled within the constraint-based framework introduced in the previous chapters. Specifically, a statistical represen- tation of heterogeneous populations can be obtained by defining suitable probability distributions on the flux polytope. This chapter addresses

○ the different sources of variation that affect microbial metabolism along with the mechanisms that may favor higher variability,

○ the methods devised to represent heterogeneous microbial populations within the framework of constraint-based models, and

○ how these approaches connect to the optimality scenario presented in the previous chapters.

**12.1 Introduction**

The theory of cellular metabolism developed up to this point through constrained-based models (CBMs) relies crucially on some type of optimality assumption: among all viable flux states encoded in the flux polytope by mass-balance, thermodynamic and regulatory constraints, cells strive for those that maximize a physiologically motivated objective function. For *E. coli* cells growing on carbon-limited substrates, for instance, it is reasonable to take such a function to be the growth yield. At the very least, these optimal states provide reference points to gauge cellular behavior. In this respect, having a good grasp of what makes a configuration of fluxes through the network ‘optimal’ with respect to a certain objective is rather important from a theoretical viewpoint. On the other hand, it is not easy to prove directly in an experiment that a certain function is *actually* being optimized (in any physical system, let alone in a microbe or a microbial population). An optimality assumption can usually be corroborated *a posteriori*, e.g. by comparing optimality predictions to experimentally measured fluxes or growth rates [355, 356], or indirectly, e.g. by showing that, in a given growth medium, certain metabolic enzymes are expressed at just the level ensuring maximal growth [148]. By looking at the behavior of individual cells in a population, however,
<table><caption>Table 13.1: Analogy between bet-hedging in biological populations and diversification strategies in Kelly's gambling and finance. The common problem in each case is an uncertain environment that makes it impossible to anticipate which phenotype or investment is optimal for future growth. In finance, the "population" is constituted by the capital which is distributed across different options (different horses of a race or different stocks of a stock market). The main limitation of the analogy is that information is not processed centrally in biological populations but at the level of each individual, with therefore no equivalent to a gambler or investor. The notations are introduced in the main text for the biological problem and in Box 2 for the gambling problem.</caption><thead><tr><th>Biology</th><th>Gambling</th><th>Finance</th></tr></thead><tbody><tr><td>Individual</td><td>Capital unit</td><td>Currency unit</td></tr><tr><td>Environment p(E)</td><td>Race results p(x)</td><td>Market state</td></tr><tr><td>—</td><td>Gambler</td><td>Investor</td></tr><tr><td>Phenotype decisions u(R)</td><td>Bets b(x)</td><td>Investment strategy</td></tr><tr><td>Multiplicative rate f(R, E)</td><td>Odds o(x)</td><td>Immediate return</td></tr><tr><td>Environmental cue P(S|E)</td><td>Side information P(y|x)</td><td>Side information</td></tr><tr><td>Population growth rate Λ</td><td>Long-term return W</td><td>Long-term return</td></tr><tr><td>Extinction probability</td><td>Probability of bankruptcy</td><td>Probability of bankruptcy</td></tr><tr><td>Growth rate variance σ<sup>2</sup></td><td>Growth rate variance σ<sub>W</sub><sup>2</sup></td><td>Volatility</td></tr><tr><td>Population size N<sub>t</sub></td><td>Capital C<sub>t</sub></td><td>Capital</td></tr></tbody></table>

antibiotics and $(1-p_a)T$ periods of low antibiotics. As further explained below, the population size $N_T$ after $T$ generation is hence expected to globally grow as $N_T = (A_{high})^{p_a T} (A_{low})^{(1-p_a)T} N_0$. This corresponds to an exponential growth (or decay) of the form $N_T = e^{\Lambda T} N_0$ with a long-term growth rate $\mathcal{L}$ given by $\Lambda = p_a \ln u_d + (1-p_a) \ln(u_d + 2(1-u_d))$.

Two bacterial populations which have different "strategies" $u_d$ will then have different growth rates $\mathcal{L}(u_d)$. The optimal strategy which maximizes $\mathcal{L}(u_d)$ is therefore when the probability $u_d$ to adopt the dormant state is

$$u_d = \begin{cases} 2p_a, & \text{if } 0 < p_a \le 1/2. \\ 1, & \text{if } 1/2 < p_a \le 1. \end{cases}$$

The interesting case is when $p_a < 1/2$, otherwise antibiotics is so often high that the population cannot grow. In this case, we find that a limited fraction of the population should be in the dormant state and that this optimal fraction depends on the frequency $p_a$ at which high antibiotics occurs.

This example can be extended to an arbitrary number of environmental states $E$ and phenotypic (response) states $R$ and to the presence of cues collected from the environment. In general, the states and cues may take discrete (as in the above example) or continuous values. The "strategy" of a cell may then be described by its probability $u(R)$ to adopt a particular phenotype $R$. This strategy depends on the environment if some signal $S$ is perceived, in which case the strategy takes the form of a conditional probability $u(R|S)$ satisfying

$$\sum_R u(R|S) = 1, \quad \text{with} \quad u(R|S) \ge 0$$

for each possible signal $S$. For the example of bacterial persistence, $u(R = \text{dormant}|S)$ may be the fraction of cells adopting a dormant phenotype within the population of cells with intracellular antibiotics concentration $S$. The fraction of growing cells would then be given by $u(R = \text{growing}|S) = 1 - u(R = \text{dormant}|S)$. By comparison, Figure 13.1B illustrates a model where $u(R|E)$ depends directly on the environmental state $E$. In finance, $u(R|S)$ would correspond to the fraction of the capital that an investor allocates to asset $R$ when receiving incomplete information $S$ on the current market state $E$. More generally, we may also consider that the probability to adopt a phenotype $R_t$ at time $t$ depends on the phenotype $R_{t-1}$ adopted at time $t-1$ by the cell or its parent, which would be described by $u(R_t|S, R_{t-1})$ or $u(R_t|S_t, R_{t-1})$ to indicate that the signal $S_t$ is obtained at time $t$.

The model also needs to specify the temporal dynamics of the environment and the relation between $S$ and $E$. The simplest assumption is that successive environmental states are uncorrelated, and occur with probability $p(E)$ and that
one cannot help but notice a salient feature: their diversity. Individual cells are macroscopically heterogeneous in terms of parameters like interdivision times, elongation rates, sizes at birth or division, etc. This suggests that a corresponding diversity is present at the level of intracellular processes like cell cycle, gene expression and, of course, metabolism. Quantitative experiments probing populations at single-cell resolution (see Experimental Methods Box 12.A) can nowadays characterize such a diversity in some detail. Among the remarkable outcomes of these studies is that, when analyzed through a lens that accounts for diversity, bacterial growth displays signatures of universality [357, 358, 359], suggesting the existence of general, system- and condition-independent control mechanisms (e.g. of cell division and growth) that do not change with specifics like strain, quality of medium, etc. Identifying these mechanisms yields robust insight (and predictive capacity) into the physiology of microbial systems (see also Chapter 14).

<figure><img src="image_112.png" /></figure>

Experimental methods 12.A Quantitative methods for single-cell analysis

At the very minimum, quantitative experimental characterization of cell-to-cell diversity in microbial populations requires (i) the possibility of achieving steady-state cell growth in controlled environments, and (ii) the possibility of identifying individual cells within a population. The two setups that are most important for the present chapter (and most widely used in general for the study of cell-to-cell heterogeneity in microbial systems) are the following.

O High-resolution optical microscopy of bacteria growing on agarose pads. Optical microscopy is the first and still most used technique to address cellular individuality [360]. Besides giving direct information about the macroscopic growth dynamics of individual cells [361, 359], it can be used in conjunction with gene expression reporters like fluorescent proteins to quantify diversity in gene expression levels [362] and dynamics [363]. Optical means usually allow to reliably follow the expression of a relatively small number of genes. In addition, however, they can also provide information about many other aspects of bacterial physiology, like motility, chemotaxis or the spatial self-organization of colonies.

O Microfluidic 'lab-on-a-chip' devices. In essence, these techniques allow to confine single cells or small lineages thereof in controlled environments for long-term data acquisition [364]. A well-known example is the 'mother machine' [365]. In a mother machine cells grow in narrow (ca. 1 µm) microfluidic dead-end channels such that (a) all cells in the same channel are daughters of a mother cell stuck at the closed side of the channel; (b) a main feeding channel carries away cells that grow out of the length of dead-end channels (which suffice to contain a few cells, usually 5 to 10); and (c) nutrient in-flow and waste out-flow from the feeding channel ensure a constant medium in all dead-end channels via diffusion. This setup effectively keeps the population size fixed. Growing bacteria can then be imaged and analyzed by standard means like time-lapse microscopy to obtain the statistics of quantities like the interdivision time or the size at birth at stationarity [358].

The setup of mother machines has the advantage that cells can be followed for many more generations than on agarose pads, since the latter tend to become overcrowded after a limited number of rounds of divisions. On the other hand, agarose pads offer a more natural environment for cell division. In addition to these, a host of other techniques are being increasingly refined and used to probe single-cell properties and behavior in bacterial populations, including single-cell metabolomics by mass spectrometry [366], nanoscale secondary ion mass spectrometry (nanoSIMS) [367], and single-cell transcriptomics [368].

It is not hard to guess why a bunch of identical cells sharing the same medium would, say, elongate at different rates. For one, gene expression has a stochastic component, from e.g. the random diffusion of transcription factors to targets to the thermal noise driving the on/off dynamics of transcription events. We also know that the cell cycle can be highly variable [369]. And other 'natural' sources of variance can be found in the dynamics of expression in genetic circuits, aging, asymmetric partitioning of cellular resources at division, inter-cellular interactions, and epigenetic modifications [370]. In other terms, a degree of variability across a population is to be expected. The question, however, is, how can variability be reconciled with the optimality picture? And related to this: can we explain cell-to-cell differences in terms of some other, perhaps more involved, optimality criterion? Are there cases in which variability is optimized? Can we describe quantitatively a microbial population in ways that account for inter-cellular diversity? Note that cell-to-cell variability is inherently a population-level concept. Addressing it therefore requires a framework that is capable of clearly distinguishing single-cell properties from population-level ones.

It is definitely possible to explain cell-to-cell variability within an optimality framework (see Chapter 5). For example, one could say that, in appropriate conditions, all microbes in a population are optimal, but the optima are slightly different for different cells. As a matter of fact, optimal states in CBMs need not be iso-
signals are derived from a conditional probability $p(S|E)$, as illustrated in Figure 13.1B where $p(S|E) = \delta(S|E)$ is equivalent to $S \equiv E$. This is sufficient to demonstrate bet-hedging or discuss the value of signaling and in the examples below we therefore make this simplifying assumption by default. More generally, to address issues of inheritance where $R_t$ depend on $R_{t-1}$, we may assume a discrete-time Markov process where the state of the next environment depends only of the previous one, with transition probabilities $p(E_t|E_{t-1})$ where $E_t$ denotes the state of the environment at time $t = 1, 2, \dots$. Even more generally, we may also want to account for the feedback that the population exerts onto its environment and consider that $E_t$ depends on the size and composition of the population.

Finally, we need to specify the dynamics of the population itself. Between time points $t$ and $t+1$, a cell adopting phenotype $R$ in the context of environment $E_t$ either dies or survives and may additionally produce offsprings. This is summarized by a quantity $f(R, E_t) \ge 0$ that indicates the mean number of descendants at time $t+1$ of an individual with phenotype $R$ in environment $E_t$ (possibly including the individual itself). Given that $u(R|S_t)$ denotes the fraction of cells or probability of the organism adopting phenotype $R$ based on sensed state $S_t$, a population is therefore expected to globally increase (or decrease) in size by a factor

$$A_t = \sum_R f(R, E_t)u(R|S_t) \qquad (13.1)$$

that depends both on the strategy $u$ and the current environmental state $E_t$. This factor $A_t$ is a stochastic variable as it depends on the stochastic variables $E_t$ and $S_t$. More explicitly, if $N_t$ denotes the size of the population at time $t$, this size will increase or decrease to $N_{t+1} = A_t N_t$ at time $t+1$ (in average). We can in this way account for the dynamics of population growth and then ask what is an “optimal” strategy $u(R|S)$ that leads to, for example, the largest population size over a given time interval. Compared to the deterministic setting, however, this is not yet a well-formulated problem as the population size varies with time and therefore generally depends on the particular sequence of environments $E_0, \dots, E_t$, which is in turn stochastic. Thus, we need to extend the concept of optimality to the stochastic regime. We examine this question in the next section.

## 13.3. Optimization in uncertain environments

In the previous section, we used notation $A_t$ to denote the fractional increase or decrease in population size given that strategy $u(R|S_t)$ is adopted in environment $E_t$. An alternate name for this quantity is the instantaneous growth rate. It follows from recursion that, given an initial population size of $N_0$ at time $t=0$, the population size at time $t$ is given by

$$N_t = A_t A_{t-1} \cdots A_1 N_0 \qquad (13.2)$$

where $A_t$ depends on the environmental state $E_t$ and is therefore a stochastic variable when the environment varies stochastically. Here the choice of an objective criterion is fundamentally linked to the time scale at which growth is considered.

### 13.3.1. Long-term versus short term optimization

At the shortest time scale, maximization of population growth over a single time step corresponds to adopting the distribution $u(R|S_t)$ that maximizes the arithmetic mean $\mathbb{E}[A]$, where $A$ denotes the random variable whose realization at time $t$ is $A_t$ (Box 1). This maximum is typically achieved by a population where all individuals adopt the same optimal phenotype – the phenotype $R$ maximizing $\mathbb{E}[f(R, E_t)u(R|S_t)] = \sum_{E,S} P(S|E)P(E)f(R, E)u(R|S)$. In the example of persistent cells, this strategy would correspond to having all cells in a growing state if the most likely environment is an absence of antibiotics. This strategy is extremely risky if these growing cells cannot survive an episode of antibiotics, which would therefore lead to extinction of the population. Taking into account the rare but important events of high antibiotics concentration requires taking a long-term perspective. Remarkably, in the long-term the problem becomes effectively deterministic due to the law of large numbers. The best known example of a law of large number applies to
<figure><img src="image_113.png" /><figcaption>Figure 12.1: Single optimum versus multiple optima in the flux polytope – (A) A two-dimensional flux polytope (shaded area) with non-negative fluxes and the defining constraints shown as dashed lines. (B) The linear objective function represented by the blue line has a unique maximum (red dot). (C) The linear objective function, represented by the blue line parallel to one of the constraints, has a continuous set of maxima points which coincides with the segment shown in red.</figcaption></figure>

lated points belonging to the flux polytope. There can in fact be infinitely many flux vectors that maximize an objective function (this happens, for instance, when an objective function attains its maxima on one of the edges or faces of the polytope, see Figure 12.1). This implies that identical cells subject to the same constraints and sharing the same objective may end up having different metabolic profiles despite carrying the same value for the objective function. In this scenario, diversity is induced by a very special feature of the objective function and, unless some other ingredient is brought into the game to lift the degeneracy, all optimal states would be equally likely for cells. If having an objective function of this type seems unlikely in a high-dimensional setup such as metabolism, one may imagine a scenario in which all cells optimize the same objective but with slightly different constraints (i.e. in a slightly different polytope, e.g. due to small variations in regulatory constraints, energy demands, or nutrient uptakes). In this case, each cell would solve its own optimization problem, ending up having, along with a different metabolic profile, a slightly different value of the objective function. Metabolic diversity is therefore induced by variability in the constraints. But it is also possible that, if cells are subject to fluctuating exogenous constraints (e.g. variable nutrient levels), they would prefer to maximize their, say, growth rate *averaged over conditions*, especially if fluctuations occur on faster timescales than those over which metabolic reactions equilibrate. In such a case, the average growth rate would be maximum (given the external variability), but other than that every cell could carry a different growth rate and a different metabolic profile. In this respect, one can say that diversity is now being optimally adapted to external conditions, or one may even think that different cells have different objective functions. This scenario, possibly unrealistic for growing microbial populations but not for other cell types (think for instance of the mixture of neurons with high energy demands and glia with low energy demands in the brain), would also lead to heterogeneous flux profiles and objectives. And so on.

It is clear from these examples that, in order to represent heterogeneity within CBMs, one must, first and foremost, clarify the origin of heterogeneity as much as possible. Next, it is necessary to shift from the language of individual flux vectors belonging to the flux polytope to that of *ensembles* of flux vectors or, more reasonably for large populations, of *probability densities* defined on the flux polytope. This transition is less trivial and more momentous than it sounds and, together with the causes of variability, is the core subject of the present chapter. We shall begin by giving a more precise characterization of the different types and sources of diversity that can be considered when modeling metabolic networks. Next, we shall introduce probability densities on the flux polytope and briefly discuss a few simple examples. We shall then address the general problem of using probability densities to represent heterogeneity and uncertainty, most notably that seen in empirical data. Finally, we will show how these ideas can be used to generalize the notion of optimality to heterogeneous populations.

## 12.2 Sources of variability and uncertainty in metabolism

Metabolic heterogeneity is widespread among clonal populations of prokaryotic and eukaryotic cells. Populations of *Escherichia coli* display diverse cell-to-cell conversion yields of glucose into final products, such
the sum $A_1 + \cdots + A_t$ of $t$ random variables $A_i$, which almost certainly behaves as $t\mathbb{E}[A]$ as $t \to \infty$. Here, the problem involves a product of random variables and a similar but different law of large number applies: the product $A_1 \times \cdots \times A_t$ does not typically behave as $(\mathbb{E}[A])^t$ but instead as $\exp(t\mathbb{E}[\ln A])$ where $\mathbb{E}[\ln A]$ is known as the geometric mean (Box 1). This corresponds to the intuition that population size typically grows exponentially in the long run, $N_t \sim e^{\mathcal{L}t}N_0$, with a well-defined long-term growth rate

$$ \mathcal{L} = \mathbb{E}[\ln A] = \sum_{E,S} p(S|E)p(E) \ln \left( \sum_{R} f(R, E)u(R|S) \right), \quad (13.3) $$

that is predictable despite the stochasticity of the environment.

Biologically, therefore, maximizing the geometric mean is equivalent to maximizing the long-term growth rate of the population. This is the relevant measure of fitness in the long-term from an evolutionary point of view, because of two populations with growth rates $\mathcal{L}_1$ and $\mathcal{L}_2$, the one with $\mathcal{L}_1 > \mathcal{L}_2$ will almost certainly exponentially outnumber the other.

### Mathematical details 13.A: Arithmetic versus geometric mean and logarithmic utility functions

Additive random processes are governed by the law of large numbers: the sum of many random variables scales with their arithmetic mean. In finance and biology, returns are compounded and growth is a multiplicative process. This is fundamentally different: the typical outcome is no longer described by the arithmetic mean but by the geometric mean [519]. A simple example illustrates this difference. Imagine a succession of environments in which the population either doubles or is reduced by 2/3, with same probability. This corresponds formally to a population size increasing as $N_t = A_t \cdots A_1 N_0$ where $A_t = 2$ (doubling) with probability 1/2 and $A_t = 1/3$ (2/3 dying rate) with probability 1/2. The arithmetic mean is 7/6 which is > 1 and suggests that the population will grow. But as each outcome has the same probability, the typical growth over $t$ generation is actually given by $2^{t/2}(1/3)^{t/2} = e^{t\mathcal{L}}$ with $\mathcal{L} = (1/2)\ln(2/3)$ which is < 0: the population will in fact most likely go extinct. Mathematically, taking the log turns the product into a sum to which the central limit theorem applies. More intuitively, the arithmetic mean is dominated by very rare events. Historically, the importance of the geometric mean for estimating risk was first understood by Daniel Bernoulli in the context of games [520, 521]. Later, it has been the subject of many debates in finance [520], reflecting the fact that alternative *utility functions* over which to optimize may be more appropriate when considering a short temporal horizon or when accounting for different degrees of risks.

The simple example of persistence that we introduced previously illustrates well how maximizing the long-term growth rate is different from optimizing the instantaneous growth rate. The arithmetic mean $\mathbb{E}[A]$ is indeed maximized by $u_d = 0$ when $p_a < 1/2$, which leads to certain extinction unless $p_a = 0$. This remains true for general models including multiple environmental states and sensing that conveys information about the environment through conditional probability $p(S|E)$. Using the long-term growth rate $\Lambda$ as a measure of fitness, it is then possible to quantify the value of information $S$ by comparing the optimal growth rate that can be achieved in presence of $S$ to that in its absence. Remarkably, for special limits of the model, corresponding to Kelly's horse-race model (Box 2), this value is given by some of the same quantities that appear in Shannon's theory of communication (Box 3).

### Mathematical details 13.B: Kelly's model

In 1956, Kelly [522] extended the work of Shannon on communication to the field of gambling. This classic model has important implications for investment strategies in finance and beyond. In the context of biology, Kelly's paper led to a clarification of the notion of value of information which is described in Box 3.

Let us recall the basic elements of Kelly's horse race. The odds paid by the bookmaker when the horse $x$ wins is $o(x)$, and the probability for this to happen is $p(x)$. A gambler can distribute his/her bets on the different horses, and $b(x)$ is the fraction of the bet set on horse $x$. Thus, a *strategy* of the gambler is defined by a vector
<figure><img src="image_114.png" /><figcaption>Figure 12.2: Minimal metabolic network with multiple optima – (A) Toy network where the top metabolite is imported by reaction $v_1$ and processed by reactions $v_2$ and $v_3$, which convert it into the bottom metabolite that is then excreted via $v_4$. When fixing $v_1 = 10$, the maximization of $v_4$ –a proxy of biomass growth rate– under mass balance results in $v_4 = 10$. There are however infinitely many flux vectors (defined by the condition $v_2 + v_3 = 10$) that are coherent with this solution, including those indicated by red blue and green flux values. (B) The subspace of optimal solutions forms a line (dark purple) in the space of feasible flux vectors (Problem 12.1). The orange-shaded triangle represents the flux polytope for $0 \le v_1 \le 10$.</figcaption></figure>

as fatty acids and tyrosine [371]; not surprisingly, the intracellular concentration of co-factors, including ATP, also vary significantly between cells [372]. *Saccharomyces cerevisiae* metabolic states have been observed to change over time for each cell. For instance, a single budding yeast does uptake oxygen before duplicating its genetic material, but it changes to an anaerobic metabolism once DNA duplication starts in order to prevent mutations related to free radicals [373]. Animal cells within a single tissue also show heterogeneous metabolisms. Non-small cell lung cancer display a remarkable diversity of preferred carbon sources. Within the tumor, some cells consume glucose and produce lactate, whereas others divert their metabolism to consume lactate as a carbon source [374].

The root cause of this metabolic heterogeneity is manifold, including uneven distribution of nutrients in the environment, asymmetric cell partitioning at division, and noise in gene expression [375, 376]. These effects are stochastic, and prevent the determination of a cell metabolic state in advance. This type of uncertainty is rooted in the nature of metabolism itself. We refer to it as *objective uncertainty*.

There is however another type of uncertainty at play, one that comes from our models of metabolism. In any metabolic network reconstruction, there can be missing reactions [377], errors or lack of knowledge about the directionality of certain reactions under *in vivo* conditions [378], and errors in the experimental estimates of certain parameters –such as exchange fluxes, or the weights of the biomass reaction [379]. Even when using a *bona fide* metabolic network conditioned by precisely measured parameters, optimality principles can lead to a reduction of the viable polytope as opposed to the identification of a single ‘optimal’ state [380] (see Chapter 5). This is exemplified in the network of Fig. 12.2.A, where the maximization of $v_4$ only reduces the viable polytope to a subspace defined by the line shown in Fig. 12.2.B. Uncertainties that stem from modeling uncertainties can be categorized as *subjective*, as they arise solely from an observer’s imperfect knowledge.

As we will see in the following, although objective and subjective uncertainties have different sources, both can be modeled using probability theory.

## 12.3 Probability densities over the flux polytope

In what follows, we shall denote the convex flux polytope by $\mathcal{P}$ and a generic flux configuration in $\mathcal{P}$ by $v = \{v_i\}_{i=1}^N$. A probability density $p$ defined on $\mathcal{P}$ is any non-negative function such that

$$ \int_{\mathcal{P}} p(\mathbf{v})dv_1 \cdots dv_N \equiv \int_{\mathcal{P}} p(\mathbf{v})d\mathbf{v} = 1. \qquad (12.1) $$
of bets b of M components b(x). At every race, the gambler invests his/her entire capital on all horses, so that $\sum_{x=1}^{M} b(x) = 1$, always betting a non-zero amount on all horses. Since no bet is zero, there is a well-defined vector of the inverse of the odds paid by the bookmaker denoted r. When the odds are fair, the bookmaker does not keep any of the invested capital and as a result $\sum_{x=1}^{M} r(x) = 1$.

At each time t, one horse, which we call x, wins with probability p(x). As a result, the capital at time t + 1 is updated according $C_{t+1} = \frac{b_x}{r_x}C_t$. As explained previously, this multiplicative process is best studied by considering instead the log of the capital, $\log\text{cap}(t) \equiv \log C_t$, which satisfies the assumptions of the law of large numbers when races are independent. In these conditions, $\log\text{cap}(t) \equiv \log C_t$ converges on long times towards the growth rate $W(\mathbf{b}, \mathbf{p})$ where

$$W(\mathbf{b}, \mathbf{p}) = \sum_{x} p(x) \log o(x) b(x). \qquad (13.4)$$

This growth rate can be rewritten using an information theoretic measure between two probability distributions, p and q, called the Kullback-Leibler divergence and defined by

$$D_{\text{KL}}(\mathbf{p}, \mathbf{q}) = \sum_{x} p(x) \log \frac{p(x)}{q(x)}. \qquad (13.5)$$

One can show that this quantity is a non-negative measure between the two probability distributions. With this notation, the growth rate can be rewritten as

$$W(\mathbf{b}, \mathbf{p}) = D_{\text{KL}}(\mathbf{p}\|\mathbf{r}) - D_{\text{KL}}(\mathbf{p}\|\mathbf{b}), \qquad (13.6)$$

It follows from this equation that the strategy $\mathbf{b}^* = \mathbf{p}$ is optimal. This strategy, known as Kelly's strategy or proportional betting, overtakes any other strategy in the long-term as illustrated in Fig. 13.2.

This formulation shows that the growth rate is the difference between the distance of the bookie's estimate from the true distribution and the distance of the gambler's estimate from the true distribution. Hence, the gambler makes money if they have a better knowledge of the winning probabilities than the bookie. The optimal long term growth rate is the positive quantity :

$$W^*(\mathbf{b}, \mathbf{p}) = D_{\text{KL}}(\mathbf{b}\|\mathbf{r}). \qquad (13.7)$$

Kelly's horse race model is formally a particular case of the model introduced in the main text when considering that one, and only one phenotype $R = R(E)$ can grow in any given environment $E$, such that $f(R, E) = f(E)$ if $R = R(E)$ and 0 otherwise. Horses $x$ may then be interpreted as both the environments $E$ and their associated phenotypes $R(E)$ so that $u(R) = b(x)$ and $f(E) = o(x)$. In biology, but also in finance where $R$ is interpreted as an asset, there is generally no one-to-one correspondence between environments $E$ and phenotypes $R$ and multiple phenotypes (assets) may grow (have non-zero return) in any given environment. The optimal strategy is then no longer necessarily proportional betting as illustrated in the example of persistence presented in the main text and as also shown in Ref. [523].

### 13.3.2. Trade-offs at intermediate time scales

So far we considered two extreme limits of immediate and infinite time scales under one important assumption: the population is always large enough to escape extinction. Eq. (13.2) is indeed valid only for large $N_t$ and does not apply anymore when $N_t \sim 1$, in which case the population size is subject to stochastic effects, called demographic noise in population biology. In our analogy with finance, the eventuality of $N_t = 0$ with no possible recovery corresponds to a risk of bankruptcy.

When considering long time scales, a population with $L > 0$ will either become extinct or grow exponentially. In this
Notice that the integral over $\mathcal{P}$ implicitly encodes two types of constraints: mass-balance equations (i.e. $\mathbf{S}\mathbf{v} = \mathbf{0}$) and ranges of variability of the form $v_{i,\min} \le v_i \le v_{i,\max}$ (see Chapter 5). The quantity $\int_{\mathcal{P}} dv$ represents therefore the a priori volume of $\mathcal{P}$ (which, understandably, is far from simple to calculate for high-dimensional polytopes like those corresponding to genome-scale metabolic network reconstructions [381]). As usual, $p(\mathbf{v})$ can be interpreted as the relative likelihood of flux configuration $\mathbf{v}$: if we imagine that a cell is assigned a flux configuration by "randomly sampling it from $\mathcal{P}$" using the rule described by $p$, then $p(\mathbf{v})d\mathbf{v}$ represents the probability that the cell's flux configuration will lie in a small volume $d\mathbf{v}$ around $\mathbf{v}$. It is clear then that probability densities on $\mathcal{P}$ provide a mathematically convenient way of describing the metabolic state of large populations (or ensembles) of cells at a given time, provided one can assume that cells have the same metabolic network and are subject to the same constraints, so that $\mathcal{P}$ is the same for all of them. For the population of cells described by $p$, the probability density clearly contains all the statistics of metabolic fluxes, from mean values to variances to correlations. For instance, by integrating $p$ over all fluxes except the $i$-th, one obtains the marginal probability density of flux $v_i$, i.e.

$$ \int_{\mathcal{P}} p(\mathbf{v})d\mathbf{v}_{\setminus i} = p_i(v_i), \qquad (12.2) $$

where the subscript $\setminus i$ corresponds to 'except for the flux of index $i$' (so $d\mathbf{v}_{\setminus i} = dv_1 \cdots dv_{i-1} dv_{i+1} \cdots dv_N$). And from $p_i$ we can immediately retrieve the statistical features of flux $v_i$ (e.g. mean value, variance, etc).

Let us make a few simple examples.

* If we assume that all cells in the population maximize the same objective function, and that there is no degeneracy in the optimal state, then

$$ p(\mathbf{v}) = \delta(\mathbf{v} - \mathbf{v}^*), \qquad (12.3) $$

where $\mathbf{v}^*$ denotes the (unique) objective-maximizing flux vector and $\delta(x)$ denotes Dirac's $\delta$-distribution.

Σ Math box 12.B Dirac's $\delta$-distribution

For our purposes, the defining property of the $\delta$-distribution in one dimension is the following: if a variable $x$ is $\delta$-distributed around the finite value $x^*$, then, for any continuous function $f$,

$$ \int_{-\infty}^{+\infty} f(x)\delta(x - x^*)dx = f(x^*). \qquad (12.4) $$

This means that, intuitively, $\delta(x - x^*) = 0$ everywhere on the real axis except at $x^*$, where its value is $+\infty$. Such a function only makes sense within an integral. In this respect, (12.3) should be seen as an abuse of notation, albeit a convenient one. There are however several ways to represent the $\delta$-distribution that comply with the above requirement. For example, one can define

$$ \begin{aligned} \int_{-\infty}^{+\infty} f(x)\delta(x - x^*)dx &:= \lim_{\sigma \to 0} \int_{-\infty}^{+\infty} f(x) \frac{1}{\sqrt{2\pi}\sigma^2} e^{-\frac{(x-x^*)^2}{2\sigma^2}} dx \\ &= \lim_{\sigma \to 0} \int_{-\infty}^{+\infty} f(x^* + \sigma y) \frac{1}{\sqrt{2\pi}} e^{-\frac{y^2}{2}} dy = f(x^*). \end{aligned} \qquad (12.5) $$

The generalization to $n > 1$ dimensions is obtained by straightforwardly assuming $\delta(\mathbf{x} - \mathbf{x}^*) = \prod_{i=1}^n \delta(x_i - x_i^*)$, so that

$$ \int_{\mathbb{R}^n} f(\mathbf{x})\delta(\mathbf{x} - \mathbf{x}^*)d\mathbf{x} = f(\mathbf{x}^*). \qquad (12.6) $$

Because the $\delta$-distribution effectively has non-zero probability mass only at a single point, it is reasonable to expect (12.6) to hold also if the integral is carried out over a compact domain $D$, provided $\mathbf{x}^*$ belongs to $D$. This is indeed the case, although the proof requires some work. For a quick guide to the many other interesting and useful properties of the $\delta$-distribution that are beyond our current scopes, see [382].

* If we can make no assumption on the cells' metabolic activity other than it has to be compatible with
<figure><img src="image_83.png" /><figcaption>Figure 13.2: Evolution of the log-capital (left plot) or of the capital itself (right plot) as a function of the number of races for Kelly's optimal strategies (red curve) and for a non-optimal strategy (yellow curve). On left plot, the straight lines have the slope of the corresponding growth rate for each strategy. Note that the fluctuations in Kelly's strategy can in fact be quite large, when plotted in normal scale instead of log-scale.</figcaption></figure>

later case, demographic noise is eventually negligible and our approach valid. At intermediate time scales, however, population sizes $N_t$ may deviate substantially from $N_0e^{Lt}$ predicted by exponential growth, and may become extinct ($N_t = 0$) as a result. To quantify these deviations, note that for the model defined in the main text where there are no correlations of the instantaneous growth rate $A_t$, the central limit theorem imposes that the quantity

$$ \Delta_t = \frac{1}{\sigma\sqrt{t}} \left( \log \frac{N_t}{N_0} - t\Lambda \right), \qquad (13.8) $$

converges on long times towards a Gaussian distribution of unit variance, where $\sigma$ is the standard deviation of the instantaneous growth rate. It follows from this property that

$$ \sigma^2 = \frac{1}{t} \mathrm{Var} \left( \log \frac{N_t}{N_0} \right), \qquad (13.9) $$

measures the deviation from exponential growth. This quantity is therefore a natural measure of risk, known in finance under the name of volatility. To understand at which time scale this risk is important, we consider Eq. 13.8, assuming $\Delta_t$ is of the order one. Risk will be important, when the term associated with fluctuations, which is of the order of $\sigma\sqrt{t}$ will be larger than the term associated with exponential growth, which is $t\Lambda$. This will happen when $t \ll (\sigma/\Lambda)^2$: the risk is relevant at intermediate time scales, long-enough for the central limit theorem to apply but not too long for deviations from exponential growth to become negligible.

This measure of risk has well known drawbacks in finance: it is symmetrical with respect to losses and gains, which does not conform to the intuitive notion of risk, and furthermore typical fluctuations are often non-Gaussian. Nevertheless, the volatility is still an important notion in the study of optimization of portfolios [527]. In this context, Markowitz introduced plots of the volatility $\sigma$ as a function of the mean growth rate, which define the so-called “efficient frontier”. This representation illustrates graphically a fundamental trade-off that exists between the maximization of the mean return and the minimization of the variance (or risk). The point of zero volatility is a risk-free strategy, which corresponds to dormant states in biology.

This trade-off is naturally present in Kelly's model introduced in Box 2. Indeed, Kelly's strategy is based on the maximization of the long-term growth rate, but at intermediate times the capital can deviate significantly from the expected exponential growth as shown in Figure 13.2. Prominent economists, such as Samuelson, strongly opposed the use of Kelly's criterion in finance precisely for that reason [528]. In practice, however, investors can mitigate this risk by using Kelly's criterion for only a fraction of the bets [529]. The resulting strategy has reduced fluctuations, and at the same time, a reduced growth rate. Another consequence of the trade-off is that the risk near the optimal
the constraints encoded by $\mathcal{P}$, then any flux vector $\mathbf{v} \in \mathcal{P}$ is equally likely to occur in a population. This means that $p$ is constant on $\mathcal{P}$. Specifically, its value must be equal to the inverse of the volume of $\mathcal{P}$:

$$p(\mathbf{v}) = \left( \int_{\mathcal{P}} d\mathbf{v}' \right)^{-1} \quad (\mathbf{v} \in \mathcal{P}) . \qquad (12.7)$$

For any given flux polytope, this distribution can be sampled at least in principle using the methods described in Chapter 4.

* Imagine having a dataset derived from a ¹³C labeling experiment (mass spectrometry) that gives the mean value $\bar{v}_i$ of every flux in the network (the average being over the population of cells used in the experiment), together with an experimental error $\sigma_i$ (which likely conflates different sources of uncertainty of which we may know very little, if anything at all), such that the experimental population-level estimate of $v_i$ is $\bar{v}_i \pm \sigma_i$. Let us assume that we know enough about the experiment to be able to define a flux polytope for the cell type ($\mathcal{P}$), and that all empirically measured averages and errors are in $\mathcal{P}$. Then, if we want to describe the population by a probability density in $\mathcal{P}$ that is uniform over the domain defined by experimental estimates, we can set

$$p(\mathbf{v}) = \prod_{i=1}^{N} \frac{\theta(\bar{v}_i + \sigma_i - v_i)\theta(v_i - \bar{v}_i + \sigma_i)}{2\sigma_i} \quad (\mathbf{v} \in \mathcal{P}), \qquad (12.8)$$

where $\theta(x)$ denotes the Heaviside (step) function defined as (Problem 12.2)

$$\theta(x) = \begin{cases} 1 & \text{for } x > 0 \\ 0 & \text{for } x < 0 \end{cases}. \qquad (12.9)$$

* (Boltzmann distribution) Let $f(\mathbf{v})$ denote a generic function of the flux vector, such as $f(\mathbf{v}) = \sum_{i=1}^{N} c_i v_i$, with $c_i$ prescribed constants. The Boltzmann distribution is defined as

$$p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta f(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}), \qquad (12.10)$$

where $\beta$ is a constant and $Z$ is a factor ensuring normalization (i.e. (12.1)), namely $Z(\beta) = \int_{\mathcal{P}} e^{\beta f(\mathbf{v})} d\mathbf{v}$. The behavior of $p$ is simple to grasp in three limits.

1. For $\beta \to 0$, (12.10) reduces to (12.7): in other words, $p$ becomes uniform over $\mathcal{P}$ (and therefore insensitive to $f$).

2. For $\beta \to +\infty$, $p$ effectively concentrates on the flux vector $\mathbf{v}^*$ that maximizes $f$ (which for simplicity we assume to be unique). To see this at a heuristic level, it suffices to notice that, for any $\mathbf{v} \neq \mathbf{v}^*$, the ratio

$$\frac{p(\mathbf{v}^*)}{p(\mathbf{v})} = e^{\beta[f(\mathbf{v}^*) - f(\mathbf{v})]} \qquad (12.11)$$

increases exponentially as $\beta$ increases. Because densities are normalized, when this ratio becomes large, $p(\mathbf{v})$ must become very small. Hence, when integrated over $\mathcal{P}$, the larger is $\beta$, the closer to $\mathbf{v}^*$ must flux vectors be in order to give a significant contribution to the integral. For $\beta \to +\infty$, the only relevant contribution comes from $\mathbf{v}^*$, so that, effectively, $p(\mathbf{v}) \approx \delta(\mathbf{v} - \mathbf{v}^*)$. This conclusion can be reached more precisely using Laplace's method (a.k.a. saddle-point approximation) to evaluate integrals of the form $\int_{\mathbb{R}^n} e^{\beta g(\mathbf{x})} d\mathbf{x}$ in the limit $\beta \to \infty$ for fixed $n$ (see e.g. [383], Ch. 27).

3. By a similar reasoning, for $\beta \to -\infty$ the only relevant contribution to integrals involving $p$ comes from the (unique, by assumption) flux vector $\mathbf{v}_*$ that minimizes $f$, so that, effectively, $p(\mathbf{v}) \approx \delta(\mathbf{v} - \mathbf{v}_*)$.

When $\beta$ varies, things depend strongly on the form of $f$ and can become rather complicated when $f$ is non-linear, especially when terms that involve the product of two or more fluxes ('high-order interactions') are present. However, in the simple case in which $f$ is linear (as outlined above), then the probability density gradually morphs from a uniform distribution over $\mathcal{P}$ to a $\delta$-distribution around the maximum of $f$ as $\beta$ increases from 0 to $+\infty$ as shown in Fig. 12.3 (and likewise when $\beta$ decreases from
Mathematical details 13.C: The value and cost of information for growing populations

To see how uncertainty may be quantified by Shannon entropy, first consider a model where $f(R, E)$ is non-zero only for one phenotype $R$ best adapted to each particular environment $E(R)$. As seen in Box 2, the optimal strategy in the long-term is proportional betting, $u(R(E)) = p(E)$. To quantify the cost of uncertainty, it must be compared to a situation where full information is available, in which case all the cells can systematically adopt the optimal phenotype, leading to an ideal growth rate $L^{**} = \sum_{E} p(E) \ln f(R(E))$. The cost of uncertainty is $L^{**} - L^{*} = - \sum_{E} p(E) \ln p(E)$, which is nothing but the Shannon entropy of the environment $H(E)$. This has a simple interpretation: the more unpredictable the environment, the larger its entropy and the lower the maximal growth rate of the population.

The reasoning can be extended to the presence of partial information, modeled by $p(S|E)$. The case of perfect information is indeed the limit case where $S = E$. The optimal strategy with partial information is a generalization of proportional betting that takes into account $S$ and the difference of growth rate is now given by the mutual information $I(S, E)$ (Problem 13.1). The mutual information is minimal when the signal $S$ is uncorrelated to $E$, in which case $I(S, E) = 0$, and maximal in presence of perfect information, in which case $I(S, E) = H(E)$ [524].

These results were first derived by Kelly [522]. They have been generalized to more general forms of $f(R, E)$ as well as to more general environmental processes in the context of financial investment in which case the cost of uncertainty and value of information are no longer equal but bounded by information theoretic quantities [524, 525].

This is illustrated in Problem 13.2 with an extension of the model of persistence presented in the main text.

Information is generally costly as it implies producing and operating an accurate sensor, which may come at the expense of growth rate. Taking into account this cost introduces a trade-off between the cost and value of information that may justify an imperfect sensor, or even explain an absence of sensor (Problem 13.1). This trade-off has for instance been invoked to explains that bacteria subject to infrequent periods of antibiotics evolved to stochastically switch their phenotype rather than to sense the presence of antibiotics [526].

While the problems of information processing in biology and in finance share many analogies, it is also important to recognize an important difference: in biology, information processing is distributed at the level of each cell, which may perceive different signals, while in finance, information is processed by an investor who centralizes the information. The value of information is bounded by information theoretic quantities only in the second case, or more generally when the same common information is available to all the cells [525]. If information processing is stochastic at the single cell level, the value of information is effectively higher (Problem 13.3).

strategy (Kelly's strategy) can be reduced significantly provided one is ready to sacrifice a small amount of growth rate, an important lesson for gamblers and investors. In order to build systematically improved gambling strategies with a reasonable amount of risk in Kelly's model, one can introduce an objective function that is a linear combination of the growth rate with the volatility of Kelly' model, $\sigma_W$, weighted by a risk aversion parameter $\alpha$ [530]. The method is illustrated in Problem 13.4 for the two-horse version of Kelly's model. By optimizing this objective function, one builds the Pareto diagram shown in Fig. 13.3 when varying the parameter $\alpha$.

A general inequality characterizes this trade-off mathematically for an arbitrary number of horses. For Kelly's gambling model with fair odds defined in the box 13.B, it has the form

$$ \sigma_W \geq \frac{W}{\sigma_q}, \qquad (13.10) $$

where $\sigma_W$ is the volatility of Kelly's model, $W$ the average growth rate (the equivalent of $\Lambda$) and $\sigma_q$ is the standard deviation of a distribution, $q(x)$ defined by $q(x) = r(x)/p(x)$. This distribution compares the probability of races outcomes described by $p(x)$ with the risk-free strategy described by $b(x) = r(x)$, for which $\sigma_W = W = 0$ [530]. Recently, a similar bound has been derived for other well-known financial models such as the Black-Scholes and the Heston models [531].
<figure><img src="image_115.png" /><figcaption>Figure 12.3: Boltzmann distribution on the flux polytope – The Boltzmann distribution, Eqn (12.10), morphs from a uniform probability density to a $\delta$-distribution concentrated on the flux vector that maximizes the function $f$ as $\beta$ varies from 0 to $+\infty$.</figcaption></figure>

0 to $-\infty$). In this respect, the parameter $\beta$ can be seen simply as a 'degree of optimization': the closer a population is to optimizing $f$, the higher the value of $\beta$. For reasons that will become clear in the next section, the Boltzmann distribution plays an especially important role in this chapter (Problem 12.3).

○ In Constrained Allocation FBA [356] (see Chapter 10), one considers an ensemble of growth-rate maximization problems constructed by sampling (from a prescribed probability density) a family of random variables representing the proteome fraction to be invested in each metabolic enzyme per unit flux of the corresponding reaction. The idea in CAFBA is that different sets of parameters effectively correspond to different cells, reflecting the cell-to-cell variability in e.g. transcription levels and protein abundances. The population-level behavior is then obtained by averaging over different choices of these parameters (i.e. over a population of heterogeneous cells). An alternative interpretation is however possible, namely that different parameters reflect the different environmental conditions that a species can encounter over its life process history. By averaging over parameters one obtains a growth strategy that levels out this environmental variability. Such a strategy may be the one that cells prefer to implement e.g. when environmental fluctuations are fast (faster than regulatory timescales). In either case, in CAFBA, randomness in a family of parameters related to the optimization problem induces randomness in the solutions, and therefore a probability density over the feasible space. This probability is unfortunately hard to write down explicitly in the case of CAFBA due to the complexity of the optimization problem. Its marginal distributions are however easy to calculate numerically. Two of them, specifically for the single-cell growth rate and acetate excretion fluxes, are shown in Fig. 2 in [356].

We could provide more examples but the key message of this section should already be visible: probability densities on the flux polytope are useful (a) when one wants to explicitly represent how uncertainties, experimental knowledge (with errors), or variability in parameters impact our knowledge of what part of the flux space $P$ is occupied by the metabolic states that occur in a true microbial population; and (b) when one is interested in representing an *optimal* (in some sense) population in a way that explicitly accounts for heterogeneity. If one has data (with errors), a probability density can provide a representation of the data, as in (12.8). It can likewise describe the solution to a population-level optimization problem, and therefore a purely theoretical prediction, as in (12.3). Or the solution to an optimization problem with uncertainty, i.e., partial knowledge or variability in some of the parameters, in which case it represents an ‘informed’ theoretical prediction (as in the CAFBA example, where the ‘information’ injected into the problem comes from the probability density from which parameters are sampled). Or it can simply be a tool to interpolate between extreme cases when we are unsure about how well a certain function is being optimized (as in (12.10)). Notice how, in our examples, different motivations activate different theoretical routes, all of which lead to working with probability densities that have *a priori* different origins and meanings even though they can be formally the same.

The two broad motivations for working with probability densities on $P$ outlined above [i.e. (a) representing uncertainty and (b) representing optimal heterogeneous populations], pose fundamentally different
<figure><img src="image_84.png" /><figcaption>Figure 13.3: Pareto diagram showing the growth rate W versus the fluctuations measured by the standard deviation of the growth rate σ<sub>W</sub> (which is the volatility for this model) in the simple case that only two horses are present. The curve can be obtained by varying a risk aversion parameter α, which enters in the definition of an objective function (see Problem 13.4 for details). The point of maximum growth rate (red square) corresponds to Kelly's strategy and divides a trade-off branch (blue solid line) from a non-trade-off branch (red solid line) (adapted from [530]).</figcaption></figure>

Let us now illustrate the implications of this trade-off for a biological population using a simple bet-hedging model with only two phenotypes. Individuals in the population can switch from phenotype A to phenotype B with a transition probability π₁, and with probability π₂ from B to A, assuming no sensing. The population grows in an environment that fluctuates between two values 1 and 2. We denote the population vector, which describes the number of individuals in each phenotype at a given time t by N(t) = (N<sub>A</sub>(t), N<sub>B</sub>(t))<sup>T</sup>, where T denotes the transpose. The subpopulation of individuals with phenotype A grows when placed in the environment i with the growth rate k<sub>Ai</sub>, while the other subpopulation with phenotype B grows with rate k<sub>Bi</sub>. The population is assumed to be large, there is no population noise, the dynamics of the system is deterministic in each separate environment. The population dynamics of the model can be described by the vector equation :

$$ \frac{d}{dt} \mathbf{N}(t) = M_{S(t)} \mathbf{N}(t), \qquad (13.11) $$

with matrices

$$ M_{S_1} = \begin{pmatrix} k_{A1} - \pi_1 & \pi_2 \\ \pi_1 & k_{B1} - \pi_2 \end{pmatrix} \quad \text{and} \quad M_{S_2} = \begin{pmatrix} -\pi_1 + k_{A2} & \pi_2 \\ \pi_1 & k_{B2} - \pi_2 \end{pmatrix}. \qquad (13.12) $$

The finite time averaged population growth rate is defined as

$$ \Lambda_t = \frac{1}{t} \ln \frac{N(t)}{N(0)}, \qquad (13.13) $$

in terms of the total population N(t) = N<sub>A</sub>(t) + N<sub>B</sub>(t), and the long term population growth rate is

$$ \Lambda = \lim_{t \to \infty} \Lambda_t. \qquad (13.14) $$

This optimal long term growth rate Λ can be obtained analytically in this model [532], but approximations are needed to evaluate the fluctuations of the growth rate, which is the equivalent of the volatility σ² of Eq. 13.9 [533]. One can then study the trade-off that exists between the average growth of the population (either measured instantaneously or over a long time) and the fluctuations of the growth rate, using the same Pareto plot used for Kelly's model in Figure 13.3. This "efficient frontier" is shown in Fig. 13.4, and as in the case of Kelly's model, in the region of fast growth rate, it is advantageous for a population to trade growth for less risky fluctuations. In this model, σ² correlates with the probability that the population N(t) goes below a certain threshold, where the population is considered as extinct. The
modeling challenges. In the first case, the key question is one of model selection: given some empirical knowledge, what is the probability density on $\mathcal{P}$ that best represents our residual uncertainty? For instance: how good of a choice for $p$ is (12.8) given the data we had? Are there criteria that can guide our choice of a probability density? We will briefly consider these issues in the upcoming Sec. 12.4. When attempting to model optimal heterogeneous populations at the theoretical level, instead, one basically has to generalize the problem tackled by CBMs like FBA to the case in which an optimal probability density is searched for instead of an optimal flux configuration. We will see how this can be done in Sec. 12.5.

## 12.4 Representing heterogeneity and uncertainty

### 12.4.1 Maximum Likelihood, Maximum a Posteriori and Bayesian inference

We have seen that probability densities on $\mathcal{P}$ can represent, under certain assumptions, populations of microbes whose metabolism can be described by the same flux polytope, and that different probability densities can be surmised to model the distribution of $v \in \mathcal{P}$ when some external information (e.g. experimental data) is available. Here, we will address the following question: how can one choose the $p$ that best represents our knowledge about the metabolic state of a population in presence of these external data?

To summarize the huge and highly involved set of problems behind the above (very general) question [383] in a way that is useful for the purposes of this chapter, we can start by assuming we have a priori chosen a form of $p$ that depends on certain free parameters and ask how to tailor parameters so that $p$ 'optimally' matches the empirical evidence. To be concrete, let us denote by $\psi$ the vector of parameters of $p$, and by $\mathbf{W} = \{\mathbf{w}^1, \mathbf{w}^2, \dots, \mathbf{w}^R\}$ a set of $R$ experimental samples of $v$. Each measurement, $\mathbf{w}$, is a vector of metabolic fluxes that ideally should include all the reactions of a metabolic network. In practice, a vector $\mathbf{w}$ typically spans only a subset of all the reactions of the metabolic network, e.g. those that are amenable to $^{13}$C labeling (TCA, glycolysis, and pentose phosphate pathways) or that correspond to exchange fluxes that can be reliably measured (glucose and oxygen consumption, or lactate and ethanol, to name a few). According to Bayes' rule (we assume all variables to be continuous), the quantities

* $p(\psi|\mathbf{W})$: the conditional probability density of the parameters given the observations (a.k.a. the posterior);

* $p(\mathbf{W}|\psi)$: the conditional probability density of the observations given the parameters (a.k.a. the likelihood);

* $p(\psi)$: the prior probability density of parameters (a.k.a. the prior);

* $p(\mathbf{W})$: the (marginal) probability density of observations (a.k.a. the evidence)

are related by the formula

$$p(\psi|\mathbf{W}) = \frac{p(\mathbf{W}|\psi)p(\psi)}{p(\mathbf{W})} \quad (12.12)$$

Ideally, what one would like to know in order to 'optimally' set the parameters of $p$ is how likely a parameter set is given the data, i.e. the full posterior $p(\psi|\mathbf{W})$, as it allows to quantify our uncertainty on the model itself. One may however also consider different (less ambitious) ways to choose parameters. The three best known methods are the following:

* Maximum Likelihood (ML) inference aims at finding the parameter vector that maximizes the likelihood:

$$\psi_{\text{ML}} = \arg \max_{\psi} p(\mathbf{W}|\psi) \quad (12.13)$$

In standard cases, this produces a single 'optimal' vector $\psi$ (hence it is called a 'point estimator'), resulting in a $p$ that models -in a context-specific manner- the metabolic heterogeneity within the cellular population.

* Maximum a Posteriori (MAP) inference aims instead at finding the parameter vector that maximizes the posterior:

$$\psi_{\text{MAP}} = \arg \max_{\psi} p(\psi|\mathbf{W}) = \arg \max_{\psi} p(\mathbf{W}|\psi)p(\psi) \quad (12.14)$$
<figure><img src="image_85.png" /><figcaption>Figure 13.4: Pareto diagram showing the population growth rate versus the fluctuations of that growth rate in a simple model of a biological population evolving in a stochastic environment with no sensing according to Eq. 13.11 [533]. In this figure, the time scale of environment fluctuations is comparable to that of phenotypic fluctuations. The inset shows the probability $P_{ext}$ that the population goes below a certain extinction level versus the risk aversion parameter $\alpha$ which measures the distance along the Pareto plot. Colored bullets represent different points on the Pareto front (adapted from [533]).</figcaption></figure>

probability of extinction is not monotonic along the Pareto front, which explains why in the region of low growth rate, it is more advantageous to prioritize instead the increase the growth rate to avoid extinction.

In the context of ecology, besides the probability of extinction, a quantity of interest is the chance for a population to grow from rarity in the presence of other species. In agreement with the above trade-off, it was found that this chance can not be predicted only from the mean growth rate, and that the mean growth rate and its variance should be both used for such a prediction [534]. In summary, the similarity of the Pareto plots (called efficient frontier in finance) obtained in Kelly's model and in models of biological populations in fluctuating environments [533], and evidences from various works in ecology, suggest that the trade-off discussed here is broadly applicable in various fields ranging from biology and ecology to economics.

## 13.4. Strategies in correlated environments

So far we considered two time scales: the time scale at which phenotypic changes occur and at which instantaneous growth is defined ($t = 1$ in our discrete-time model, which may be taken to correspond to one generation), and the longer time scale $t \sim (\sigma/\mathcal{L})^2$ beyond which population growth is effectively exponential, with growth rate $\mathcal{L}$. We saw that the choice of an optimization criterion depends fundamentally on the time horizon relative to these time scales.

Additional time scales are relevant when environmental states are correlated in time, for instance through a Markov chain $P(E_t|E_{t-1})$. This is for instance the case if conditions of high nutrient or high stress extend over several generations. As a consequence, strategies $u(R_t|S_t, R_{t-1})$ that depend on past internal states $R_{t-1}$ in addition or instead of externally driven signals $S_t$ may become valuable, since the fact that phenotype $R_{t-1}$ survived in environment $E_{t-1}$ indirectly carries information on the current environment $E_t$. We may then recognize that $R_t$ plays two distinct roles: on one hand, it determines survival and growth via $f(R_t, E_t)$ and, on the other, it provides information to determine the next state $R_{t+1}$ via $u(R_{t+1}|S_{t+1}, R_t)$. This corresponds to the fundamental distinction between phenotype and genotype in biology: the genotype $\gamma$ is what is transmitted from one generation to the next while the phenotype $\phi$ is what determines instantaneous growth. Formally, $R_t = (\phi_t, \gamma_t)$ with $f(R_t, E_t) = f(\phi_t, E_t)$ and $u(R_t|S_t, R_{t-1}) = u(R_t|S_t, \gamma_{t-1})$, by definition of $\phi_t$ and $\gamma_t$. The “central dogma” of molecular biology states that information flows from the genotype to the phenotype but not reciprocally, which corresponds here to assuming that $u(\phi_t, \gamma_t|\gamma_{t-1})$ factorizes as $d(\phi_t|S_t, \gamma_{t-1})h(\gamma_t|\gamma_{t-1})$, where
Σ Math box 12.C Gaussian model

Maximum Likelihood (ML) is the most commonly used point estimation method. As said above, the estimated parameters, $\hat{\psi}$, are computed as the argument that maximizes the likelihood of the observed data, i.e.

$$ \hat{\psi} = \arg \max_{\psi} p(\mathbf{W}|\psi) = \arg \max_{\psi} \prod_{i=1}^{R} p(\mathbf{w}^{(i)}|\psi) = \arg \max_{\psi} \sum_{i=1}^{R} \log[p(\mathbf{w}^{(i)}|\psi)], \quad (12.15) $$

where in the last step we used the fact that, as far as the solution is concerned, maximizing $p(\mathbf{W}|\psi)$ is equivalent to maximizing its logarithm. ML takes a familiar form if one follows, for instance, Theorell et al. [384] in modeling data according to a multivariate normal distribution:

$$ \mathbf{w} \sim N(\mathbf{w}|\psi) = \frac{1}{\sqrt{(2\pi)^N |\Sigma|}} \left(-\frac{1}{2}(\bar{\mathbf{v}} - \mathbf{w})^T \Sigma^{-1} (\bar{\mathbf{v}} - \mathbf{w})\right) \quad (12.16) $$

The parameters encompass the mean values, $\bar{\mathbf{v}}$, and the covariance matrix, $\Sigma$. That is, $\psi = [\bar{\mathbf{v}}, \Sigma]$. Accordingly,

$$ p(\mathbf{w}^{(i)}|\psi) = N(\mathbf{w}^{(i)}|\psi) = \frac{1}{\sqrt{(2\pi)^N |\Sigma|}} e^{(-\frac{1}{2}(\bar{\mathbf{v}}-\mathbf{w}^{(i)})^T \Sigma^{-1} (\bar{\mathbf{v}}-\mathbf{w}^{(i)}))}, \quad (12.17) $$

and

$$ \hat{\psi} = \arg \max_{\bar{\mathbf{v}}, \Sigma} \sum_{i=1}^{R} \left[ -\frac{1}{2}(\bar{\mathbf{v}} - \mathbf{w}^{(i)})^T \Sigma^{-1} (\bar{\mathbf{v}} - \mathbf{w}^{(i)}) - \log \left( \sqrt{(2\pi)^N |\Sigma|} \right) \right], \quad (12.18) $$

which leads to the well-known weighted least squares estimators of mean values ($\hat{\bar{\mathbf{v}}}$) and variances ($\hat{\Sigma}$). With $\hat{\bar{\mathbf{v}}}$ and $\hat{\Sigma}$, the frequency of any vector $\mathbf{w}$ can be computed from $N(\mathbf{w}|\hat{\bar{\mathbf{v}}}, \hat{\Sigma})$. Standard techniques, such as confidence intervals, can be applied to assess the precision of $\hat{\psi}$. Generally speaking, the larger the number of samples, $R$, the smaller the uncertainty in $\hat{\psi}$.

where the last equality follows from the fact that $p(\mathbf{W})$ does not depend on $\psi$. As for ML, the MAP estimator is a point estimator.

* Bayesian inference aims finally at computing the full posterior distribution $p(\psi|\mathbf{W})$. It is therefore a 'distribution estimator' rather than a point estimator.

Problem 12.4 should clarify the way in which point estimators differ from (and are less informative than) distribution estimators in practice.

### 12.4.2 MaxEnt inference

According to the principle of Maximum Entropy (MaxEnt) [399], among all probability densities that are consistent with given prior knowledge or data, the one having the largest value of the entropy

$$ H[p] = - \int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v}) d\mathbf{v} \quad (12.19) $$

is the one that best represents our knowledge about the system. A classical intuitive justification of the MaxEnt principle is most easily given for discrete variables [400].

Consider $N$ cells, each of which can be found in any of $K$ states (what precisely defines a state is immaterial for this reasoning). Let an assignment $\mathbf{n} = \{n(i)\}$ be given, such that $n(i)$ denotes the number of cells in state $i$ (with $\sum_{i=1}^K n(i) = N$). Because we can always exchange the states of two cells without changing $\mathbf{n}$, there are multiple 'microscopic' ways to realize an assignment $\mathbf{n}$. Combinatorics tells us that the number of different microscopic realizations of an assignment $\mathbf{n}$ is given by

$$ \mathcal{N}(\mathbf{n}) = \frac{N!}{\prod_{i=1}^{K} n(i)!} \quad (12.20) $$
<figure><img src="image_86.png" /><figcaption>Figure 13.5: Example of optimal strategies in correlated environments – **A**. We consider here a model where the environment $E_t$ is a continuous variable following a Gaussian process, $P(E_t|E_{t-1}) = \exp(-(E_t - aE_{t-1})^2/(2\sigma_X^2))/(2\pi\sigma_X^2)^{1/2}$ with two parameters $a$ and $\sigma_X^2$ that control the overall amplitude of the fluctuations $\sigma_E^2 = \sigma_X^2/(1-a^2)$ and their time scales $\tau_E = -1/\ln a$, as illustrated by the different time series. **B**. An individual inherits a genotype $\gamma_{t-1}$ which determines its phenotype $\phi_t$ with probability $d(\phi_t|\gamma_{t-1}) = \exp(-(\phi_t - \gamma_{t-1})^2/(2\sigma_D^2))/(2\pi\sigma_D^2)^{1/2}$ where $\sigma_D^2$ thus represents phenotypic noise. $\gamma_{t-1}$ also determines the genotype $\gamma_t$ of the progeny with probability $h(\gamma_t|\gamma_{t-1}) = \exp(-(\gamma_t - \gamma_{t-1})^2/(2\sigma_M^2))/(2\pi\sigma_M^2)^{1/2}$ where $\sigma_M^2$ thus represents mutational noise. The number $\xi$ of offsprings is a random variable whose mean $f(\phi_t, E_t) = k \exp(-(\phi_t - E_t)^2/2)$ depends on the phenotype $\phi_t$ as well as the current state $E_t$ of the environment. A population of such individuals grows with a long-term growth rate $\bar{L}$ that can be computed analytically [535]. **C**. Values of $\sigma_D^2$ and $\sigma_M^2$ that optimize $\bar{L}$ define four phases as a function of the environmental parameters $\tau_E$ and $\sigma_E^2$. For nearly constant environments, the optimal strategy is to maintain constant phenotypes ($\sigma_D^2 = 0$) and genotypes ($\sigma_M^2 = 0$) ("no variation"). For strongly varying but poorly correlated environments, the optimal strategy is to introduce phenotypic variations ($\sigma_D^2 > 0$) but no genotypic mutations ($\sigma_M^2 > 0$) ("phenotypic switching"). For highly correlated environments, the optimal strategy is instead to introduce genotypic mutations ($\sigma_M^2 > 0$) while canalizing the phenopype ($\sigma_D^2 = 0$) ("inherited variations"). A phase also exists where both types of variations are beneficial ("mixed"). This model thus identifies environmental variations for which bet-hedging (phenotypic switching) is expected to evolve, namely variations of environmental of sufficient amplitude but with limited temporal correlations across generations.</figcaption></figure>

$d(\phi_t|S_t, \gamma_{t-1})$ may be interpreted as a developmental kernel and $h(\gamma_t|\gamma_{t-1})$ as an inheritance kernel, with no dependence on $S_t$ (no Lamarckism). The mathematical framework that we introduced can be used to study to which extent this particular decomposition is indeed a good "strategy" [535]. The answer generally depends on the nature and amplitude of the environmental fluctuations.

Similarly, the model can be analyzed to understand the conditions under which it is advantageous to introduce phenotypic variations that are not transmitted – as in bet-hedging – versus genotypic variations that are transmitted – as with genetic mutations. Stochasticity may indeed be introduced either in the mapping from $\gamma_{t-1}$ to $\phi_t$ or the mapping from $\gamma_{t-1}$ to $\gamma_t$, or in both of them – a problem with no equivalent in finance. This is illustrated in Fig. 13.5 with a simple solvable model showing how the optimal strategy depends on the nature of the fluctuations of the environment. In particular, bet-hedging strategies where stochasticity is purely phenotypic are found to be optimal for environmental fluctuations of sufficient large amplitude but low temporal correlations from one generation to the next.

Historically, the notions of genotype and phenotype were introduced much before the molecular mechanisms that underlie them were uncovered. In general, the genotype, defined as inherited information, should not be confused with the notion of genetic information: along with DNA, a range of epigenetic states, including metabolic states, are also transmitted from cell to cell which represent genotypic information. In other words, the physiological state of a cell, which we analyzed in most of this book from the standpoint of a phenotype determining current growth, may also represent valuable genotypic information for future generations.
Σ Math box 12.D Inference in metabolic network modeling

In metabolic network modeling v is usually a vector of fluxes. Unfortunately, the number of samples is usually very small [385], which may lead to $\hat{\psi}$ over-fitted to the sample set. One way to overcome limited sample sizes is to regularize the estimation procedure by incorporating prior information on $\psi$ via the MAP estimation method (12.14). The evidence $p(\psi)$ in MAP can be used to encode the distribution of v values observed in previous experiments or formulated as a plausible non-informative probability distribution. For example, Heinonen et al. [386] formulated $p(\psi)$ as a multivariate normal distribution with mean values equal to zero, and variances for each flux adjusted to prevent fluxes extending beyond their lower and upper bounds defined in $P$. MAP estimation can be considered as an ML estimation whose objective function has been augmented by the prior distribution of $p(\psi)$. In this sense, MAP estimation is a 'regularized' ML estimation, which helps prevent overfitting.

MAP estimation however does not exploit the capacity of Bayes' theorem to explore the full set of values that the parameters can achieve. By producing a distribution estimation of the parameters, Bayesian inference allows quantifying the parameters' variability. Compared to point estimation methods, though, Bayesian inference is computationally expensive as it requires to assess how different values of $p(\psi)$ affect $p(W|\psi)$. Fortunately, some families of $p$ are susceptible to methods such as Gibbs sampling or Markov Chain Monte Carlo that offer an efficient way to compute the posterior numerically [387]. This is the case, for instance, for the truncated multivariate normal distributions that Heinonen et al. [386] used for the likelihood and prior functions appearing in (12.12).

The posterior can then be used to derive statistical features of quantities that depend on $\psi$, e.g. metabolic fluxes. In practice, most parameters underlying the mechanisms that govern cellular metabolism -e.g., enzymes' allosteric regulation or the local conditions within cells' organelles- remain unknown. Various hypotheses can be advanced to close this knowledge gap. Alas, it is not uncommon to have conflicting scenarios. For instance, to explain overflow metabolism in S. cerevisiae and E. coli [388, 389, 390], numerous plausible explanations have been pushed forward, including ATP savings for the production of non-oxidative enzymes (which by being smaller, compared to their oxidative counterparts, require less ATP in their synthesis) [391, 392], limited uptake rates capacity [393], and an upper limit on the dissipation of Gibbs energy [394]. (See [395] for an excellent review of optimization-based explanations.) Because each mechanism can be encoded through a different prior, it is clear that the choice of the prior is a delicate matter in Bayesian inference. Generally speaking, the choice of the prior becomes less and less problematic the more data we have, i.e. the better sampling we have of the state space of the system. However, if data is scant, the prior will leave an important imprint on the resulting posterior. In these cases, a careful selection of the prior is paramount. Among the methods most commonly employed are (a) the construction of empirical priors (namely priors that encode previous knowledge about parameters), (b) the use of so-called "non-informative priors" (i.e. priors that reflect 'vague knowledge' about parameters, like the fact that a certain parameter is non-negative) [396], and (c) the selection of priors based on the Maximum Entropy principle (see below) [397, 398].

If all $n(i)$'s are large enough, we can use Stirling's approximation ($n! \approx (n/e)^n$) to see that

$$ \mathcal{N}(\mathbf{n}) \approx e^{NH(\mathbf{n})} \quad , \quad H(\mathbf{n}) = -\sum_{i=1}^{K} \frac{n(i)}{N} \ln \frac{n(i)}{N} \equiv -\sum_{i=1}^{K} p(i) \ln p(i) \equiv H(\mathbf{p}) \quad , \quad (12.21) $$

where $p(i)$ denotes the fraction of cells in state $i$ (or, equivalently for us, the probability to find a cell in state $i$). $H$ is the entropy of the assignment $\mathbf{n}$, and is in essence a measure of the microscopic degeneracy that underlies a macroscopic arrangement. The distribution $\mathbf{p} = \{p(i)\}$ carrying the largest entropy subject to certain constraints is therefore the one having the largest underlying microscopic degeneracy given those constraints. So, if one were to randomly pick a microscopic state given those constraints, the most likely macroscopic state would be the maximum entropy distribution. In other terms, the MaxEnt distribution is the least biased distribution compatible with the constraints, as any other distribution satisfying the same constraints would correspond to a smaller underlying degeneracy, thereby neglecting some feasible (i.e. constraint-satisfying) microscopic configurations. In this respect, a MaxEnt distribution requires the least information besides prior knowledge (i.e. constraints). (A more detailed justification for using the MaxEnt principle as an inference tool is given e.g. in [400].) If for instance cells are assigned to states in a completely random way, the MaxEnt distribution is the solution of

$$ \max_{\mathbf{p}} - \sum_{i=1}^{K} p(i) \ln p(i) \quad \text{subject to} \quad \sum_{i=1}^{K} p(i) = 1 \quad , \quad (12.22) $$
13.5. Perspectives

We presented optimal strategies that biological populations may exploit for coping with uncertain environments and drew
analogies with problems of gambling and financial investments. Optimality assumes a measure of performance which,
however, is not readily defined when environments are changing stochastically. In particular, the time scale over which
the problem is considered is critical. This difficulty has led to multiple debates over the concept of fitness in biology
which partly mirror those over the concept of utility in economics.

While the analogy with finance is instructive, its limitations should also be kept in mind. Most importantly, the states
that individuals of a biological population adopt are not centrally controlled by a gambler or an investor. This raises a
question that is absent in finance but central in evolutionary biology: is a strategy that is optimal for the population
but detrimental to some of its members – as for instance the persister cells that “sacrifice” their current growth for the
sake of future growth – evolutionary stable? A strategy that is optimal for a population may indeed never be achieved
through evolution as natural selection at the individual level may favor non-cooperating individuals – an issue known
as a “conflict between levels of selection” which implies that a strategy may be optimal at the population level but
not evolutionarily stable. To address this question, we may extend our model to treat strategies as variables that are
themselves subject to evolution (Problem 13.5). For the model discussed in this chapter, the results show that strategies
that optimize the long-term growth rate are indeed evolutionarily stable (but this is no longer necessarily the case when
considering, for instance, sexually reproducing populations [536]).

The same extension of the model to evolving strategies shows that knowledge of the statistics of the environment (p_d for our example) is not required a priori but can effectively be learned through an evolutionary process. This solves a problem that appears also in gambling and finance where the statistics of the environment must be inferred from past experience. The question has been particularly studied in finance, where optimal learning strategies known as universal portfolios have been proposed [537]. In the simpler case of Kelly's model, the gambler may for instance record previous race results and use them together with Bayesian inference to predict the probability of the race outcomes [538]. With biological populations, however, learning must be performed at the individual level. One theoretical proposal that goes beyond random mutations is for instance that biological populations may use a reinforcement mechanism akin to Hebb's rule in neural learning [539].

Finally, we note that the models that we presented rely on a strongly simplifying assumption: the environmental changes
occur independently of the population. In fact, the environment is often also changing as the population grows, for
instance through the consumption of nutrients. Even more generally, the environment may comprise other individuals
from the same or other populations with which they may interact. This ecological dimension is the subject of other
chapters.

Recommended readings

**Persistence** This reference describes the phenomenon of bacterial persistence. Nathalie Q Balaban, Jack Merrin, Remy Chait, Lukasz Kowalik, and Stanislas Leibler. Bacterial persistence as a phenotypic switch. Science, 305(5690):1622–1625, 2004. doi: 10.1126/science.1099390

**Information theory for decision making under uncertainty** This book is a classic text on the use of information theory in problems from finance. Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. John Wiley & Sons, Hoboken, 2005. doi: 10.1002/047174882X

**Model for information sensing** The Kussel-Leibler model is one of the first to incorporate information theory into cellular behavior and signaling. Edo Kussell and Stanislas Leibler. Ecology: Phenotypic diversity, population growth, and information in fluctuating environments. Science, 309(5743):2075–2078, 2005. ISSN 00368075. doi: 10.1126/science.1114383.
<figure><img src="image_116.png" /></figure>

Economics analogy 12.E Maximum Entropy economic equilibrium

Most of economic theory relies on the assumption that markets are capable of allocating resources optimally, i.e. so that the utilities of each of the participating agents is maximized (an assumption can be seen as the analog of each cell in a population maximizing its growth rate). In order to achieve optimal states (called 'equilibria' in economics), agents endowed with a priori different preferences, resources and goals identify the actions that maximize their utilities (e.g. transactions, trade, or production) and carry them out. This process however can become more and more demanding as the number of agents that take part in the market gets larger and larger, because (in short) the set of viable transactions for each agent can become exceedingly large. How can one describe the equilibria that arise from these situations?

A possible approach, used at least since [401], is based on the Maximum Entropy principle. The idea, in short, is the following. Once every agent has somehow chosen their preferred actions (i.e. once a system-wide 'configuration of individual actions' has been selected), the market as a whole presents a set of transactions to be carried out that aggregate the choices of individual agents. When looked at the aggregate level, though, each set of transactions can correspond to more than one configuration of individual actions. (This can happen, for instance, because agents have a degree overlap in their characteristics which makes them indistinguishable from an economic perspective.) If one assumes that agents choose their actions at random from their set of viable transactions, then some sets of transactions are bound to be more likely than others, simply because they can be realized in more 'microscopic' ways (for instance, by interchanging agents of the same type). It is then reasonable to think that the likelihood of any particular set of transactions will be larger, the larger the number of microscopic ways in which it can be realized. Taking entropy as a measure of multiplicity, the most likely set of transactions, then, is the one that maximizes the entropy.

A model of market where the above program is worked out in detail is found in [401]. The 'statistical equilibrium' theory that follows from the use of the Maximum Entropy principle generalizes the standard competitive equilibrium discussed in microeconomics by providing a description of optimality in large markets with heterogeneous participants. This line of work has also inspired further developments that explicitly included agents' heterogeneity into the theory of competitive equilibria [402, 403, 404]. To the best of our knowledge, a similar approach has not yet been used to model heterogeneous microbial systems.

which can be found via the method of Lagrange multipliers (Problem 12.5). If other constraints are imposed, though, the MaxEnt distribution will clearly change (Problem 12.6). For our purposes, the continuous case with entropy given by (12.19) can be seen as a straightforward generalization of the discrete one.

To get some grasp of the scenario that the MaxEnt rule provides within CBMs, let us work out one especially noteworthy case, namely the MaxEnt probability density of flux configurations with a given mean value of a generic function $f$ of the fluxes. This probability density is the solution of

$$ \max_{p(\mathbf{v})} - \int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v}) d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v}) d\mathbf{v} = 1 \quad \text{and} \quad \int_{\mathcal{P}} f(\mathbf{v}) p(\mathbf{v}) d\mathbf{v} = \bar{f} . \qquad (12.23) $$

To find it, we construct the functional

$$ \mathcal{L}[p] = H[p] + \alpha \left[ \int_{\mathcal{P}} p(\mathbf{v}) d\mathbf{v} - 1 \right] + \beta \left[ \int_{\mathcal{P}} f(\mathbf{v}) p(\mathbf{v}) d\mathbf{v} - \bar{f} \right], \qquad (12.24) $$

where $\alpha$ and $\beta$ are Lagrange multipliers for the normalization and the mean-value-of-$f$ constraints, respectively. Variation of $\mathcal{L}$ with respect to $p$ yields the maximum condition

$$ -1 - \ln p(\mathbf{v}) + \alpha + \beta f(\mathbf{v}) = 0 . \qquad (12.25) $$

Solving for $p$ results in

$$ p(\mathbf{v}) = \frac{e^{\beta f(\mathbf{v})}}{e^{1-\alpha}} . \qquad (12.26) $$

The normalization condition however determines the value of $\alpha$, as one must have

$$ \int_{\mathcal{P}} e^{\beta f(\mathbf{v})} d\mathbf{v} = e^{1-\alpha} \equiv Z(\beta) . \qquad (12.27) $$
## Problems

**Problem 13.1 Kelly strategy with partial information** In analogy with Kelly's problem of betting on horse races, assume that different environments $E$ occur with independent probabilities $p(E)$ at each generation with a single phenotype $R = E$ permitting growth by a factor $f(E)$. In absence of any information, the optimal strategy $u(E)$ for long-term growth is proportional betting, $u(E) = p(E)$ (Box 2). Now assume that an information $S$ is available to each member of the population that relates to $E$ through a transition probability $q(S|E)$, i.e., $q(S|E)$ is the probability of perceiving $S$ given $E$.

(1) Show that the long-term growth rate can be written in the form

$$ \Lambda = \sum_{S} p(S) \left[ \sum_{E} p(E|S) \ln(f(E)u(E|S)) \right] \quad (13.15) $$

where $p(S)$ is the probability to perceive $S$ averaged across all environments and $p(E|S)$ is the probability that environment is $E$ given that $S$ is perceived. Write $p(E|S)$ as a function of $p(E)$ and $q(S|E)$.

(2) Justify that the optimal strategy is $u(E|S) = p(E|S)$.

(3) Compare the optimal long-term growth rate in presence of information to the optimal growth rate in absence of information and show that the difference is given by the mutual information

$$ I(E, S) = \sum_{E,S} q(S|E)p(E) \ln \frac{q(S|E)}{p(S)} \quad (13.16) $$

The mutual information $I(E, S)$ therefore quantifies the value of information $S$ in this particular context.

(4) Acquiring information is generally costly. If the presence of the information channel $q(S|E)$ reduces the long-term growth rate by $c$, what are the conditions on $p(E)$ for the presence of this channel to be beneficial?

(5) The cost $c$ may be expected to depend on the precision of the sensor. Consider for instance a channel that reveals the correct environment with probability $1 - \epsilon$ and otherwise does not reveal anything (so-called erasure channel). Given a cost $c(\epsilon)$ that increases when $\epsilon$ decreases, which value of $\epsilon$ provides an optimal trade-off between the value and the cost of information?

**Problem 13.2 Value of information beyond Kelly's model** Consider the model of bacterial persistence introduced in the main text where cells can adopt two phenotypes, one growing irrespectively of the environment and the other growing only in absence of antibiotics.

(1) Express the long-term growth rate $\mathcal{L}$ in presence of an information $S$ modeled by an information channel $q(S|E)$.

(2) What is the optimal strategy given $S$?

(3) Show by comparing to a situation with no information that the value of information can be strictly lower than $I(S, E)$.

**Problem 13.3 Stochastic sensing at the level of individual cells** In the two previous problems, the information $S$ is assumed to be common to each member of the population. Here we assume instead that each individual has its own sensor $q(S|E)$ so that $S$ may differ from one individual to the next.

(1) Justify that in this case the long-term growth rate takes the form

$$ \Lambda = \sum_{E} p(E) \ln \left( \sum_{R,S} f(R, E) u(R|S) q(S|E) \right) \quad (13.17) $$

(2) Use the concavity of the logarithm (Jensen's inequality) to justify that the same information channel $q(S|E)$ has more value at the individual level than at the population level.
<figure><img src="image_117.png" /><figcaption>Figure 12.4: MaxEnt modeling of single-cell growth rate distributions – Empirical distributions (blue points, left) are reproduced by a MaxEnt assumption where the mean growth rate is constrained (right), leading to a Boltzmann distribution over the flux polytope (Eq. (12.30)).</figcaption></figure>

One is then left with

$$p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta f(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}). \qquad (12.28)$$

The value of $\beta$ must be determined from the constraint on the mean value, namely from

$$\frac{1}{Z(\beta)} \int_{\mathcal{P}} f(\mathbf{v}) e^{\beta f(\mathbf{v})} d\mathbf{v} = \bar{f}. \qquad (12.29)$$

Notice that the result is nothing but Boltzmann's distribution (12.10). We have therefore found that (12.10) is the MaxEnt distribution for a given mean value of the function $f$. This means that if we have a dataset returning the empirical mean value of an observable $f$ over a population of cells, our knowledge is best represented by assuming that $p(\mathbf{v})$ is of the form (12.10), with $\beta$ ensuring the matching of empirical and theoretical means.

This suggests a possible way to represent single-cell growth-rate distributions [405], such as the E. coli populations growing in rich media studied e.g. in [358, 359] (see Figure 12.4).

Let us assume that all cells in the population can be described by the same flux polytope $\mathcal{P}$ and let $\lambda(\mathbf{v})$ denote the growth rate associated to flux configuration $\mathbf{v}$. We can ask the following question: what is the $p(\mathbf{v})$ on $\mathcal{P}$ that best represents our knowledge that the mean growth rate of cells is $\bar{\lambda}$ (empirical)? The answer is

$$p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta \lambda(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}), \qquad (12.30)$$

where $Z(\beta) = \int_{\mathcal{P}} e^{\beta \lambda(\mathbf{v})} d\mathbf{v}$, and where $\beta$ is set so that the empirical mean growth rate ($\bar{\lambda}$) matches the theoretical mean, i.e.

$$\frac{1}{Z(\beta)} \int_{\mathcal{P}} \lambda(\mathbf{v}) e^{\beta \lambda(\mathbf{v})} d\mathbf{v} = \bar{\lambda}. \qquad (12.31)$$

We can therefore solve the above equation (numerically) and analyze the resulting distribution. One sees from (12.30) that $\beta$ has a 'natural' unit given by $\lambda_{\max}^{-1}$, the inverse maximum growth rate achievable in $\mathcal{P}$ (which is easily computed by LP). In the populations analyzed in [405], the value of $\beta$ that ensures the matching condition ranges from $190/\lambda_{\max}$ to $300/\lambda_{\max}$, suggesting that indeed the degree of optimization of $\lambda$ is significant. The most remarkable result, however is that the marginal distribution of the growth rate computed from (12.30), namely

$$p(\lambda) = \int_{\mathcal{P}} \delta(\lambda - \lambda(\mathbf{v})) p(\mathbf{v}) d\mathbf{v}, \qquad (12.32)$$

matches the overall empirical growth-rate distributions. In other words, if one adjusts the parameter of (12.30) so that the theoretical mean growth rate and the experimental one coincide, then (12.32) reproduces the entire empirical growth-rate distribution. This observation confirms the empirical evidence that the variance of single-cell growth-rate distributions is a function of the mean, such that, if growth rates
**Problem 13.4** Pareto front for Kelly's model Let us consider Kelly's model with fair odds for two horses. Let the probability that the first horse wins be $p$, the bet and the odd on the first horse be $b$ and $1/r$.

(1) Write the expression of the mean growth rate $\langle W \rangle$, and of the volatility $\sigma_W$ for this problem. Show that there is a risk free strategy when $b = r$.

One introduces the objective function

$$J = \alpha \langle W \rangle - (1 - \alpha) \sigma_W. \qquad (13.18)$$

(2) From the optimization of $J$ show that the optimal strategy has the two branches shown in Fig. 13.3. Show that the optimal bets on these two branches are of the form $b^{\pm} = p \pm \gamma\sigma$, where $\gamma = (1 - \alpha)/\alpha$ and $\sigma = \sqrt{p(1-p)}$.

(3) Show that the slope of the Pareto border has the form

$$\frac{d\sigma_W}{d\langle W \rangle} = \frac{\sigma}{p-b}. \qquad (13.19)$$

What happens to this slope near Kelly's point and near the risk free strategy ?

**Problem 13.5** Evolution of an optimal strategy Here we consider evolving the strategy itself.

(1) Implement numerically the model of bacterial persistence introduced in the main text for a large but finite population. To this end, consider $N$ individuals (e.g., $N = 1000$), each with an attribute $R$. For each individual, draw a random number $\xi$ of descendants, with mean $f(R, E_t)$ where $E_t$ drawn from $P(E)$ is common to all individuals. Assign a $R$ to each of these descendants with probability $u(R)$. If the total number of descendants $N_t$ is non-zero, record the ratio $N_t/N$ and re-sample at random the population to bring back its size to $N$. Show that provided that $N$ is large enough and $N_t$ does not reach 0 then $(\sum_t \ln(N_t/N))/t$ provides a good approximation to the growth rate $\Lambda$ in the limit of large $t$.

(2) Extend the model to make $u_d u(R = \text{dormant})$ an attribute of each individual. Assume that $u_d$ is transmitted from one parent to one of its offspring as $u_d = \min(1, \max(0, u_d + \mu))$ where $\mu$ is normally distributed with variance $\sigma_M^2$. Show that provided that $\sigma_M^2$ is small enough, the distribution of $u_d$ evolves to be centered around the optimal $u_d$.
are re-scaled by the mean, distributions roughly collapse on 'universal curves' [358, 359]. In addition, the analysis of [406] has shown that predictions for individual fluxes obtained from (12.30) (i.e. mean values plus standard deviations) provide a better fit to experimentally measured fluxes than growth-rate maximizing fluxes obtained from FBA. This is especially important as it suggests that, despite the relatively high degree of optimization, the cell-to-cell variability underlied by (12.30) is biologically relevant.

In the following section we will use this observation as a springboard for the analysis of optimal heterogeneous populations.

## 12.5 Representing optimal populations

Let us start from a rather abstract question. Suppose that an organism is actually maximizing a certain function $F$, unknown to us, which depends on metabolic fluxes $v$ as well as on a set of other variables $w$ that are not part of metabolism: $F \equiv F(v, w)$. We shall denote by $v^*, w^*)$ the (supposedly unique) configuration of variables where $F$ attains its maximum. Let's furthermore say that we have a guess for what the organism's objective function might be, and that this guess is only a function of metabolic fluxes, which we denote by $f \equiv f(v)$. If we trust our guess, and if $f$ is maximized by the (supposedly unique) flux vector $\hat{v}$, our prediction for the fluxes would be $\hat{v}$. Question: what is the probability that $\hat{v}$ is the true optimum, i.e. that $\hat{v} = v^*$? Note that $f(v^*) \equiv f^* < \hat{f} \equiv f(\hat{v})$ (i.e. at the 'true' optimum the value of $f$ is bound to be smaller than the maximum value of $f$).

The answer goes like this: according to the MaxEnt principle, the probability density $p(v)$ for any flux configuration $v$ to be the true state of the system (i.e. the true optimum) should be undetermined other than by our knowledge that the real optimum has some value of $f$ below $\hat{f}$. What is the correct constraint to enforce (besides normalization) if we are to look for such a $p(v)$? We could impose that allowed configurations strictly have some fixed value of $f < \hat{f}$. This choice would lead to a uniform density over all states with a given value of $f$. In this way, though, we are imposing that states with a different value of $f$ are strictly inaccessible, which is not part of our knowledge. However, if we impose that only the mean value of $f$ is constrained, MaxEnt will return a probability density with the exact same mean value as the uniform density just described (by construction) but a much larger entropy, just because -intuitively- it will assign a non-zero probability to all states. Hence, as long as we have no other information, the best prediction we can make for $p(v)$ is given by the probability density that maximizes the entropy $H[p]$ subject to the constraint $\langle f \rangle \equiv \int_P p(v)f(v)d\mathbf{v} = f^*$, i.e. by the solution of

$$ \max_p - \int_P p(v) \ln p(v) d\mathbf{v} \quad \text{subject to} \quad \int_P p(v) d\mathbf{v} = 1 \quad \text{and} \quad \int_P p(v) f(v) d\mathbf{v} = f^* . \quad (12.33) $$

We now know the result to be given by (12.10), i.e.

$$ p(v) = \frac{1}{Z(\beta)} e^{\beta f(v)} \quad (v \in \mathcal{P}) , \quad (12.34) $$

where $\beta$ is the Lagrange multiplier enforcing the constraint $\langle f \rangle = f^*$. What this means in practice is this: if one is modeling a microbe's metabolism and is unsure about the objective function but has a guess ($f$), information theory suggests that the best one can do is to assume that metabolic flux configurations are selected according to (12.34). Ideally, the value of $\beta$ for which one obtains the best agreement between predictions based on sampling (12.34) and experiments is the 'degree' to which the system optimizes $f$. If $f$ is the true objective function, then the agreement between theory and experiments will get better and better as $\beta$ increases. It is important to keep in mind that (i) while we have assumed that the organism is actually maximizing something, we didn't really use the fact that $F$ is maximized at $(v^*, w^*)$ (only that the true state of the system has some value of $f$ below $\hat{f}$); (ii) this is a totally ideal situation (for instance, experimental data have errors, so whether comparisons between theory and experiments are informative doesn't only depend on the theory but also on the quality of the data).

The fact that (12.34) is 'optimal' in a rather fundamental sense (a priori different from the sense in which $f$-maximizing populations are optimal) encourages to view distributions described by (12.30) through a
<figure><img src="image_118.png" /><figcaption>Figure 12.5: Fitness-information bound (general form) – The orange line encodes the maximum mean growth rate achievable for any given value of the information content I (Eq. (12.35)) of a metabolic flux distribution (or the minimum value of I required to achieve any given mean growth rate).</figcaption></figure>

different lens. When we maximize the entropy at fixed mean growth rate, in practice, we are looking for the 'broadest' probability density (i.e. the most variable population) on $\mathcal{P}$ that is compatible with the given mean. In other terms, we are saying that, given a mean growth rate, the optimal population is the one that has the largest possible variability. To quantify variability in a more readily understandable way, it is convenient to transform it into a measure of the amount of information encoded in $p$. One can reason as follows: if no prior information is available about the population, uncertainty is maximal and all flux vectors in $\mathcal{P}$ must be considered to be equally likely. This means that, for such a population, the probability density over $\mathcal{P}$ is uniform (see (12.7)). We shall denote the entropy of the uniform distribution over $\mathcal{P}$ by $H(0)$. When we inject information into the problem (e.g. the fact that the population has a certain mean growth rate), then the probability density is no longer uniform but given, say, by (12.30). The uncertainty is therefore reduced by $H(0) - H(\beta)$, where $H(\beta)$ is the entropy of (12.30). (Clearly, $H(0)$ is just the entropy of (12.30) for $\beta = 0$.) The quantity

$$I = \frac{H(0) - H(\beta)}{\ln 2} \qquad (12.35)$$

denotes the amount of information (in bits, hence the factor $\ln 2$) injected by a non-zero value of $\beta$. Reformulating our population-level optimization, we can say that, for any fixed mean growth rate $\langle \lambda \rangle$, the optimal population is the one carrying the smallest value of $I$. A short calculation (Problem 12.7) shows that $\langle \lambda \rangle$ and $I$ are related by

$$\beta \langle \lambda \rangle = I \ln 2 + \int_{0}^{\beta} \langle \lambda \rangle d\beta', \qquad (12.36)$$

where it should be noted that $\langle \lambda \rangle$ is an increasing function of $\beta$ (as $\beta$ increases, the density gets more and more concentrated around the growth-rate maximizing flux vector, thereby leading to an increase of $\langle \lambda \rangle$).

The curve $\langle \lambda \rangle$ versus $I$ described by (12.36) can therefore be computed numerically for any metabolic network reconstruction (as the only ingredients required are encoded in the flux polytope $\mathcal{P}$) [405]. The resulting line (see Figure 12.5) separates the $(\langle \lambda \rangle, I)$ plane in a viable (achievable) region and a forbidden region where the mean growth rates are too large for the amount of information encoded in the population. This 'phase diagram' yields, first and foremost, a general prediction linking the mean growth rate (fitness) of a microbial population to its metabolic heterogeneity: all populations must have fitness-heterogeneity values in the viable region. Recent work relying on an advanced statistical inference framework has shown that actual microbial populations indeed lie in the viable part of the plane [407]. In addition, it provides a quantitative definition of an optimal population that accounts for variability: optimal populations have fitness-heterogeneity pairs that lie on the boundary between the viable and the forbidden region. In this respect, results from [405, 406, 407] can be summarized by saying that heterogeneous, faster-growing E. coli populations (mean growth rate larger than roughly 1/h, richer growth media) are very close to optimality, while slower-growing ones tend to have mean growth rates and information contents that get more and more sub-optimal the less rich is the growth medium. (Of course, this notion of optimality refers to the growth rate and information content as the key parameters to evaluate a population's performance. It may well be, and this is an issue definitely worth exploring, that slower-growing population are optimal
Chapter 14

Economy of organ shapes and function

Christophe Goupil, Éric Herbert, Cyril Karamaoun, Benjamin Mauroy and Frédérique Noël

Chapter highlights

* The central concept of economy presented in the previous chapters are extended to the level of the organ and the living organisms.

* The development of organs in pluricellular living beings is constrained by several requirements. Nutrients, energy and form are few examples of developmental constraints that are analyzed, in the point of view of the economy of the organ function.

* In mammals, the respiratory system is submitted to a high level of constraints, mainly energetic and morphometric, that did shape the lung through evolutionary processes.

* The lung is the central organ of respiration, connecting the outer atmosphere to the cellular compartment through the ventilation process.

* We show that the constraints on this major organ imply a high level of complexity of the organ's shape and a precise control of the ventilation.

* The control of growth and shape of the lung is highlighted by previous works on the scaling laws that exist and govern the development and function of this organ throughout the entire mammal class.

* Through several examples, we demonstrate how these scaling i.e., *allometric* laws control the ventilation, and the respiratory processes in general.

14.1. Optimization of organs and systems

In the previous chapters, the central model of the cell has been deeply explored. On another scale, the integration of cells into larger structures such as tissues, organs, and entire systems in multicellular organisms requires an extension of the main concepts presented in this book. Nevertheless, the completion of the multiple functions of an organ follows the same general principles as for individual cells, including the economic aspects.

14.1.1. Organs and constraints

The way that a cell population aggregates itself into a high-level structure, as part of a pluricellular organism, has been determined through evolutionary processes, following a more general path of specialization of structures and functions. Each organ evolved to fulfill its functions in the most optimized manner. This observation leads us to interrogate the concept of optimization for such a large structure – from a cellular point of view.
with respect to some other parameter(s).) At any rate, the above definition of optimality coincides with the standard one (growth-rate maximization) for $\beta \to +\infty$, in which case variability goes strictly speaking to zero as all cells collapse on the same flux configuration. And we now understand how it generalizes it: by stressing the way in which heterogeneous populations can be optimal despite growing at sub-maximal rates.

For later convenience, note that, because the entropy is a convex functional, the solution to the MaxEnt problem is the same as the solution to

$$ \max_p \int_{\mathcal{P}} p(\mathbf{v})\lambda(\mathbf{v})d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v})d\mathbf{v} = 1 \quad \text{and} \quad -\int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v})d\mathbf{v} = H^*. \quad (12.37) $$

The above problem has perhaps a more direct interpretation: the optimal population is the one that has the largest mean growth rate at fixed variability (entropy) or, equivalently, at fixed information content.

Before moving on, we notice that, in the above setting, optimality of heterogeneous populations has a rather simple mechanistic interpretation in terms of how populations 'occupy' the flux polytope. If one considers the uniform distribution on $\mathcal{P}$, Eq. (12.7), and calculates the marginal distribution for the growth rate (i.e. (12.32)), one finds that the growth-rate landscape in which populations grow is extremely skewed towards slow growth rates: the overwhelming majority of metabolic flux configurations corresponds to slow-growing cells, i.e. with growth rates roughly two orders of magnitude below $\lambda_{max}$. This implies that, whatever flux vector we are in, a small random change to it is overwhelmingly more likely to decrease our growth rate than increase it. In this respect, slow states have an 'entropic' advantage over fast states. On the other hand, by definition, fast-growing flux configurations replicate faster than slow-growing ones, and therefore have a replicative advantage. It is therefore tempting to interpret the probability density (12.30) as resulting from the balance between these two tendencies. One can for instance imagine that a microbial population grows and evolves in time in $\mathcal{P}$ due to (i) replication events, and (ii) small random changes of the flux vector (due e.g. to gene expression noise). Ref. [405] has indeed shown that such a population evolves toward a distribution very close to (12.30), where the role of $\beta$ is played by the inverse rate of diffusion of the population in $\mathcal{P}$, that is, by the inverse of the rate at which small random changes occur: fast rate implies small $\beta$, and vice versa. (As the mathematical analysis of this scenario requires the toolbox of non-linear Fokker-Planck equations, it is beyond the scopes of this chapter.)

The above theory can be extended in various directions. We shall limit ourselves to one example here, namely that of optimal populations in fluctuating environments [408]. The basic assumption we make is that the growth rate $\lambda$ is a function of both the flux vector $\mathbf{v}$ and of a single (for simplicity) exogenous variable $s \ge 0$ representing the stress level to which the population is subject: $\lambda \equiv \lambda(\mathbf{v}, s)$. We furthermore assume that $s$ is a random variable with probability density $P(s)$. For any value of $s$, $\lambda$ will be maximized by a certain flux vector $\mathbf{v}^* \equiv \mathbf{v}^*(s)$. If fluctuations of $s$ are sufficiently slow, then cells may be able to perfectly adapt their metabolic response to every value of $s$ they encounter. But this is unlikely to be possible in rapidly fluctuating environments. In the latter case, it is instead reasonable to assume that cells will try to maximize their average growth rate, where the average is taken over the distribution of $s$. The relevant quantity is now the probability density to observe a certain flux configuration $\mathbf{v}$ given that the state of the environment is $s$ ($p(\mathbf{v}|s)$), while the objective function (to be maximized over $p(\mathbf{v}|s)$) is just

$$ \langle \lambda \rangle = \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s) \lambda(\mathbf{v}, s) d\mathbf{v}. \quad (12.38) $$

We should now specify the constraints. One is simple and concerns normalization: $\int_{\mathcal{P}} p(\mathbf{v}|s) d\mathbf{v}$ should be equal to one for all $s$. To introduce the second one, we note that, because one expects $\mathbf{v}$ to encode information about the environment, it is convenient to constrain the mutual information of $\mathbf{v}$ and $s$, i.e.

$$ I(\mathbf{v}; s) = \int ds P(s) \int_{\mathcal{P}} p(|s|) \log_2 \frac{p(\mathbf{v}, s)}{p(\mathbf{v})P(s)} d\mathbf{v}, \quad (12.39) $$

where $p(\mathbf{v}, s) = P(s)p(\mathbf{v}|s)$ is the joint distribution of $\mathbf{v}$ and $s$, whereas $p(\mathbf{v}) = \int ds P(s)p(\mathbf{v}|s)$. Clearly, $I = 0$ if $p(\mathbf{v}, s)$ factorizes over $\mathbf{v}$ and $s$ and it gets larger and larger as $\mathbf{v}$ and $s$ become more and more
In the context of organ function, optimization can be defined through the processes by which the functions are fulfilled as best as possible while minimizing the associated cost variables.

Among those variables, energy plays a central role. Thus, one possibility for constraining the organ would be to maintain its function at an optimal level while minimizing its cost in energy.

This effect can be expressed mathematically. Let us define the cost in energy $\mathcal{E}$ which depends on one or more variables $x \in \mathbb{R}^n$ ($n \ge 1$). Furthermore, let us define one or more equality constraints to our problem, $c(x) = 0$, where $c: \mathbb{R}^n \to \mathbb{R}^m$. The optimization problem under constraints comes down to finding an optimal value for $x$ that minimizes the function $\mathcal{E}(x)$ while $x$ satisfies $c(x) = 0$. This results in

$$ \min_{x \in \mathbb{R}^n} \mathcal{E}(x), \text{ such that } c(x) = 0. \qquad (14.1) $$

The optimization under constraint problem can be solved using the Lagrangian function,

$$ \mathcal{L}(x, \lambda) = \mathcal{E}(x) - \sum_{k=1}^{m} \lambda_k c_k(x), $$

where the $\lambda_k$ are Lagrange multipliers. Indeed, if we assume that $x_*$ is the optimal solution to the problem (14.1), then thanks to the Lagrange multiplier theorem, there exists a unique Lagrange multiplier $\lambda_*$ such that,

$$ \nabla \mathcal{E}(x_*) = \lambda_*^T J_c(x_*), $$

where $J_c$ is the Jacobian of the function $c$. It implies that the optimal solution $x_*$ is a stationary point of $\mathcal{L}$, satisfying the condition of minimal energy expenditure.

A study of the constraints on the cardiac system offers an excellent example of energy optimization, due to high consumption of the heart. The cardiac pump delivers deoxygenated blood to the lung through pulmonary circulation and brings oxygenated blood to the whole body through systemic circulation [540].

Interestingly, blood pressure developed in both ventricles are not of the same order of magnitude, with a left ventricular pressure approximately ten times larger than the right ventricular one [540]. This makes sense from an energetic point of view; the heart requires a non-negligible amount of energy to fulfill its role of blood pumping. Furthermore, as with any mechanical system, only a fraction of the energy consumed (mainly in the form of ATP) is converted in mechanical work – around 25% [541], the rest being dissipated as heat. Thus, the pumping work tends to be optimized from an energy consumption point of view.

The pressure needed to irrigate the pulmonary circulation is low; the lung presents a small value of resistance to perfusion, and its apex is located only centimeters above the heart position. To the contrary, the pressure developed in the systemic circulation must allow the oxygenated blood to irrigate all the organs, including high-energy consumers – muscles, brain – that located further above heart position, developing a hydrostatic pressure that the blood flow has to overcome [540]. It is to be noted that this energetic optimization is also connected to the metabolism requirements, with a pumping work closely related to the body's $O_2$ consumption [542], which ensures an optimized adaptation of the cardiac output to the body energy requirements.

Although often considered as a major aspect, this energetic constraint is far from being the only condition for a proper functioning of the organ. Other variables such as nutrient consumption, metabolic integration or physical constraints participate in shaping the organ function.

Brain development in primates, and humans especially, is a prime example of effect that combination of several constraints has on energetic and nutrients availability. The underlying mechanisms that determine the evolution towards a large and complex brain structure in humans are still debated [543, 544]. However, it is evident that the development and normal function of this organ is dependent on adequate and specific energetic and nutrients inputs.

From the energetic point of view, brain metabolism largely depends on glucose consumption. However, in case of high
correlated. Putting these pieces together, we can write the cell's optimization problem as

$$ \begin{gathered} \max_{p(\mathbf{v}|s)} \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s) \lambda(\mathbf{v}, s) d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v}|s) d\mathbf{v} = 1 \quad (\forall s) \\ \text{and} \quad \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s) \log_2 \frac{p(\mathbf{v}|s)}{p(\mathbf{v})} d\mathbf{v} = I^* . \quad (12.40) \end{gathered} $$

A comparison with (12.37) should clarify how the above generalizes the previously discussed optimization framework. Again using the method of Lagrange multipliers one finds that the optimal probability density is now given by (Problem 12.8)

$$ p(\mathbf{v}|s) = \frac{p(\mathbf{v})}{Z(s, \beta)} e^{\beta \lambda(\mathbf{v}, s)}, \quad (12.41) $$

where

$$ Z(s, \beta) = \int_{\mathcal{P}} d\mathbf{v} p(\mathbf{v}) e^{\beta \lambda(\mathbf{v}, s)}, \quad (12.42) $$

while $\beta$ is a Lagrange multiplier.

The meaning of (12.41) is straightforward: when $\beta \to 0$, the metabolic flux configuration $\mathbf{v}$ becomes independent of $s$, implying $I = 0$. As $\beta$ increases, $\mathbf{v}$ and $s$ get more and more correlated, while $p^*(\mathbf{v}|s)$ tends to get more and more sharply peaked around $\mathbf{v}^*(s)$. In the limit $\beta \to +\infty$ cells respond to each value of $s$ by selecting the exact flux configuration that maximizes $\lambda$. To achieve this, maximal $I$ is required. A detailed study of the optimal probability density emerging in this case within a highly coarse-grained model of metabolism has been carried out in [408], showing how complex metabolic strategies (including the coexistence of slow-growing, persistent states with fast-growing ones) arise as optimal responses to a fluctuating environment.

## 12.6 Concluding remarks

Metabolic variability in cell populations has, as we have discussed, multiple origins, both rooted in unavoidable stochastic effects and (possibly) in the fact that, in certain cases, being heterogeneous can be optimal for a microbial population. Models can account for variability by representing (sufficiently large) populations via probability densities defined on the flux polytope. Two main (different) goals can be achieved. First, one can look for the probability density that yields the best (in a precise sense) description of a set of empirical data. Methods like Maximum Likelihood and Maximum Entropy provide different, albeit related, approaches to this task. Second, one can formulate optimization problems for populations, whose general solution is a probability density rather than a single flux configuration. Solutions to these problems can highlight how fitness and variability are related in optimal populations, providing useful theoretical benchmarks for real microbial systems. While possibly more demanding from a mathematical viewpoint (and certainly more demanding from a computational viewpoint), these approaches expand the scope of CBMs, including in terms of predictive power. In addition, they can refine the notion of optimality and provide insights into the fundamental principles that govern the organization of metabolism across populations. The question of whether variability confers an advantage to microbial populations is however very general, and goes beyond the metabolic level of CBMs on which we focused here. A broader discussion of these aspects is presented in Chapter 13.

## Recommended readings

* Gardiner, C. (2009) *Stochastic Methods: A Handbook for the Natural and Social Sciences*. Springer Series in Synergetics (Springer, Berlin)

* Jaynes, E. T. (1957) Information theory and statistical mechanics. Physical Review, 106:620

* MacKay, D. J. (2003) *Information theory, inference and learning algorithms* (Cambridge University Press)

* De Martino, A., & De Martino, D. (2018) An introduction to the maximum entropy approach and its
consumption and/or deprivation, ketones metabolism takes place in order to furnish a fast and rich source of energy for the organ. Ketones are catabolized mainly in the liver [545], and have the important property of being able to cross the blood-brain barrier, to the contrary of long chains of saturated fatty acids [546].

In parallel, proper brain development and function require a large input of specific nutrients that are not common in every food source [547]. Among those, iodine [548] and iron [549] appear to be essential for the brain, and exert a strong constraint on its adequate growth and functioning. The notable presence of iodine-enriched food sources close to the sea shores, compared to traditional terrestrial food sources, is thought to have favored the recent development of the so-called *shore-based paradigm* of human brain evolution i.e., that the access to seafood produce supported and enhanced brain development in early hominid populations, leading to increased brain mass and cognitive functions in those populations [547].

Among these considerations, let us remind that any organ has to develop and function in specific localization and body environment. Thus, the constraints applied to an organ and its development and function can also be of morphometric nature.

14.1.2. The example of the lung

As an example of an organ submitted to geometric limitations, the lung has to face, from its early development to its mature state, multiple constraints on its morphometry and proper functioning. The principal role of this organ is, as well known, to establish the connection between the respiratory gases in the atmosphere and these in circulation in the body i.e., O₂ as a reactive agent, and CO₂ as a by-product that has to be eliminated from the organism.

To fulfill its role, the lung has evolved in a manner that maximizes the gas exchange surface – as diffusion is a surface phenomenon – in a reduced thoracic volume. This surface-to-volume requirement has forged the lung structure has it is known; an intricate dichotomic bronchial tree that conducts the air inwards and outwards, to and from the alveolar sacs, respectively. This semi-fractal, space-filling structure, presents the advantages of an extremely wide exchange surface enclosed in a relatively small volume [540].

The mechanisms of development of the lung branching structure in a closed environment is still a debated question [550, 551]. Indeed, the tree structure presents a series of specific characteristics necessary for a proper functioning of the organ. Among these, the space-filling aspect of the bronchial tree is remarkable, as it solves the problem of the surface-to-volume constraint of the organ. In addition, the whole bronchial tree is a self-avoiding structure, as no bronchus enters in contact with other ones in its local environment, which ensures a proper circulation of the air in the structure. It is striking that these properties, which can be found in fractal geometries, are observed in any well-functioning lung structure, leading to important developmental questionings.

The pattern of branching of the bronchi, although strongly stereotyped in the first generations starting from the trachea, appears to follow a space-filling procedural development rather than a deterministic branching pattern [550, 552, 553]. Accordingly, some authors have developed a set of hypotheses that tend to explain these mechanisms. A group of restricted genes would encode the steps of branching and growth of the bronchi during the organ development, ensuring a proper structure of the adult lung [551], following procedural steps somehow encoded in genes or groups of genes coding for periodicity, bifurcating and rotating routines [554]. However, to the best of our knowledge, these genes have not been determined nor a proper molecular mechanism of stereotypical branching.

Among the questions raised by the *programmed morphogenesis* approach, the link between the molecular dimensions and the organ world are still elusive. Another path for branching procedure, which could reconcile the deterministic point of view with the problem of the transfer of information along different orders of magnitude is the *self-organized morphogenesis* approach. Several authors [550, 555, 556] suggested that the branching routine of the bronchial tree is less stereotyped than thought, especially in the central and distal generations. This hypothesis is supported by the observation that modelling approaches using stochastic space-filling routines, constructed based on a stereotyped proximal tree, are capable of generating tri-dimensional branched structures that satisfy the constraints of a morphometric adult
application to inference problems in biology. Heliyon, 4:e00596

O Scharfenaker, E., & Yang, J. (2020). Maximum entropy economics. The European Physical Journal Special Topics, 229:1577-1590

## Problems

Computer exercises for this chapter can be found on the book website.

### Problem 12.1 Multiple optima in the model from Fig. 12.2

Use the model in Fig. 12.2 to study different objective functions, specifically combinations of fluxes. Can you find other cases in which the optimum is not unique?

### Problem 12.2 Relationship between Dirac's $\delta$ and Heaviside's $\theta$

Show that, for a real variable $x$, a continuous function $f$ and upon integration over $\mathbb{R}$, $\frac{d}{dx}\theta(x) = \delta(x)$. Hint: Use the fact that $\frac{d}{dx}[\theta(x)f(x)] = \theta'(x)f(x) + \theta(x)f'(x)$.

### Problem 12.3 Well-defined versus ill-defined flux spaces

Using the sampling methods introduced in Chapter 4. and a linear objective function of your choice, write a program that will sample a toy two-dimensional flux polytope according to (12.10), and check the outcome for a few values of $\beta$. Then try changing the shape of the polytope in different ways by changing the constraints. What features of the polytope can make sampling harder and/or less accurate (i.e. require a larger number of samples)? Can you work out a modification of the sampling algorithms that alleviates these problems?

### Problem 12.4 MAP inference versus Bayesian inference

Consider a Bernoulli random variable with parameter $\psi$, i.e. such that the probability of having $k$ successes in $n$ trials given $\psi$ is

$$p(k|\psi) = \binom{n}{k} \psi^k (1-\psi)^{n-k}, \qquad (12.43)$$

and assume that the prior for $\psi$ is a $\beta$-distribution with parameters $a$ and $b$, i.e.

$$p(\psi) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \psi^{a-1} (1-\psi)^{b-1}. \qquad (12.44)$$

Calculate the full posterior $p(\psi|k)$ and the MAP estimator for $\psi$ as a function of $k$, $n$, $a$ and $b$. Then assume $a = b = 2$ and compare the following situations: (i) a Bernoulli process that returned 2 successes in 3 trials; (ii) a Bernoulli process that returned 20 successes in 33 trials. Show that the MAP estimator for $\psi$ is 60% for both (i) and (ii) (so the two processes are indistinguishable to MAP), while the posterior is different. Knowing the posterior, which process would you pick if you were asked to point to the one that is more likely to have $\psi = 0.6$?

### Problem 12.5 Uniform distribution as the solution of (12.22)

Show that the solution of the maximization problem (12.22) is the uniform distribution $p(i) = 1/K$ for all $i$.

### Problem 12.6 MaxEnt distribution in different cases

Assume that a certain real variable $x$ takes values $x(i)$ in the $K$ states (one can for instance think of $x(i)$ as the growth rate of cells in state $i$). Show that the MaxEnt distributions for constraints imposed on (i) normalization of the distribution, (ii) normalization and mean value of $x$, (iii) normalization, mean value of $x$ and second moment of $x$, and (iv) normalization and mean of the logarithm of $x$, are, respectively, uniform, exponential, Gaussian, and power-law.

### Problem 12.7 Fitness-information relationship

Retrieve formula (12.36).
<figure><img src="image_87.png" /><figcaption>Figure 14.1: Proposed mechanism for the morphogenesis of biological branched structures. In this approach, the gradient of concentration of a key molecular activator guides the growth of specific tissue layer through the activation of the associated receptor. In this example, the local gradient of concentration of FGF10 (black arrows) activates the budding of the epithelium layer (branching mechanism – left). As the tissue curvature flattens, the local concentration vanishes and the growth stops, preventing the tissue overlap (avoiding mechanism – right) [550].</figcaption></figure>

lung [557, 558].

The core concept of the self-organized morphogenesis approach relies on the observation that key molecular components are necessary and sufficient for proper growth and branching of the bronchial epithelium. Among these, the fibroblast growth factor 10 encoded by the *fgf10* gene has been demonstrated to play a central role in epithelial proliferation [559], whose activity is highly regulated [550, 560].

In 2012, Clément et al. [550] proposed a scenario for the spontaneous emergence of a tree structure. This scenario is based on the sole diffusion of a protein promoting cell proliferation, such as FGF10, in an environment with two layers that mimic the bronchial epithelium and the lung mesothelium. In addition, the layers present a resistance to folding and are growing as a function of the received flow of proteins. This scenario forms a model for the lung development and has been studied using mathematical and numerical tools.

To mimic the diffusion process from the outer layer of the organ (mesothelium) to the inner layer (bronchial epithelium), Clément et al. solved the Laplace equation applied on the protein concentration *c*:

$$ \Delta c = 0 $$

Then, they considered that each layer was growing according to a function of the local protein gradient:

$$ \frac{dx}{dt} = f_m (||\nabla c(\mathbf{x})||) \quad \text{for } \mathbf{x} \text{ in the mesothelium} $$

$$ \frac{dx}{dt} = f_e (||\nabla c(\mathbf{x})||) \quad \text{for } \mathbf{x} \text{ in the bronchial epithelium} $$

The functions $f_m$ and $f_e$ are increasing functions, typically with a sigmoid shape. To avoid the epithelium to catch up with the mesothelium, $f_e$ is kept smaller than $f_m$. A smoothing of the layers based on a fixed characteristic size is then performed in order to mimic the layers resistance to folding.

With this model, Clément et al. observed the spontaneous formation of branching patterns similar to those observed during bronchial development, as depicted in Figure 14.1 and, based on an extended model, in Figure 14.2.

Hence, this compact modelling approach is sufficient for observing *de novo* branching and growth patterns in a simulated tissular environment. Since then, this self-organized morphogenesis approach has been used as a framework for other organs [562] and other branched systems [563]. To date, the question of the mechanisms of development of branching
**Problem 12.8 Optimal flux distribution in fluctuating stress conditions**

Retrieve formula (12.41) (hard).
<figure><img src="image_88.png" /><figcaption>Figure 14.2: Tridimensional spontaneous emergence of a tree (*budding*) based on the model of Clément et al. (2014) [561]. An eighth of the budded sphere has been sliced out to show the branching core (blue). Notice the self-avoiding and space-filling branching, which are commonly found in biological tree structures.</figcaption></figure>

organs is not clearly elucidated. However, the link between the molecular and cellular components requires further investigation, in order to unveil the determinants at the scale of the tissues and organs.

However, these examples of constrained organ development raises several issues that need to be discussed in details. Among these, one can notice that the shape of the system appears as central in the developmental considerations, especially under constraints.

**Philosophical remarks 14.A: The origin of shape?**

How growth and organ specialization define the shape and structure of the mature organ is a long-debated scientific question, which has not unveiled all its secrets and mysteries, as we will discuss in Section 14.3.

**Box 14.B: Interlude – Energy**

When energy is transferred to a system, its response manifests itself at the microscopic level by the excitation of its individual degrees of freedom, and at the global level when collective excitations are possible. In generic terms, a thermodynamic machine is defined as a system where an incident energy flow *dispersed* is converted into an energy flow *aggregated* and a loss flow. This conversion is performed by a thermodynamic *working fluid* which, carrying entropy, leads to a coupling between the respective potentials through the equations of state.

In the case of thermal machines, the dispersed form of energy is called heat and its associated potential is temperature, while the aggregated form is called work and its associated potential is, for example, pressure. Temperature and pressure are linked by one or more equations of state. The system response results from the collective response of the microscopic degrees of freedom of the working fluid. Thus, part of the energy received by the working fluid can be made available to a load on a global, and possibly macroscopic, scale for a given purpose as useful work, the remainder being redistributed (dispersed) at the microscopic level and dissipated due to internal friction and any other dispersion processes imposed by the boundary conditions [564]. Conversion efficiency is therefore closely related to the proportion of energy allocated to the system's collective modes.

Living organisms are open, out-of-equilibrium and dissipative systems, as they continuously exchange energy and
# Chapter 13

## Cells in the face of uncertainty

David Lacoste, Olivier Rivoire, and David Tourigny

### Chapter overview

* Organisms that grow and survive in uncertain environments may need to change their physiological state as the environment changes.

* When the environment is uncertain, one strategy known as bet-hedging is to make these changes randomly and independently of the environment, to ensure that at least part of the population is well adapted.

* Organisms that collect information from their environment may also use this information to modulate their changes of physiological states.

* We review these different strategies and point out parallels with the theory of optimal financial investments.

### 13.1 Introduction

To a large extent, the content of this textbook prior to the current chapter has dealt with models of microorganisms under the implicit assumption that the dynamics of both environmental factors and intracellular components are deterministic, and that behavior is optimized uniformly across cells in a population. On longer time scales however, natural selection also acts on populations and these populations may encounter environments that fluctuate across both time and space. Under these conditions, natural selection may not favor a homogeneous deterministic cellular response across the population, but rather select for a certain level of population diversity and heterogeneity, including behaviors arising from mechanisms that are fundamentally stochastic. Stochasticity is inherent to intracellular processes such as gene expression and signal transduction due to the small number of molecules that they involve. It is often referred to as “noise”, but this terminology can be misleading because may also fulfill an essential role in cellular function and survival, for example during growth in uncertain environmental conditions. The purpose of this chapter is to highlight this role, introduce the mathematical models necessary for understanding it, and draw a new economic analogy with problems of investment in finance.

Before expanding upon the role that uncertainty plays in shaping cellular behavior, we briefly point out some general limitations of deterministic models based on optimal regulation of behavior in time as described in Chapter 11. In that chapter, it was assumed that microorganisms have evolved, under selective pressures exerted by the environment, to optimize a specific objective criterion or combination of objective criteria that were shared by all cells of a population. This assumption was then incorporated into an optimal control framework to explain how cellular behavior (e.g., enzyme expression) is optimally regulated in time depending on deterministic interactions between a microbial population and its environment. In particular, we consider optimal control strategies across a prescribed time window. Defining in such case
matter with their environment [565, 566]. Unlike classical thermodynamic engines, for which equilibrium models can be constructed using extremal principles, such a possibility does not exist in the case of living organisms due to the absence of truly identifiable equilibrium states and the absence of principle for non-equilibrium systems. Nevertheless, assuming a global system close to equilibrium, the development of a tractable thermodynamic model of metabolism can be based on notions from classical equilibrium thermodynamics. In this approach, the working fluid acts as a conversion medium, characterized by its thermoelastic properties, or chemicoelastic coefficient for chemical systems.

Despite the complexity of biological systems, it is possible to apply the Onsager's phenomenological approach of locally linearized non-equilibrium thermodynamics Onsager [567]. Through this approach, it is possible to identify the non-equilibrium processes that link the degradation of the chemical potential of food by its digestion into a macroscopic form of energy made available for muscular work. By applying Onsager's approach and integrating it with macroscopic systems, we can describe the behavior of certain thermodynamic conversion machines under mixed boundary conditions [564, 568, 569]. In the case of Dirichlet boundary conditions the system is driven by the potential differences, meanwhile in Neumann boundary conditions the system is driven by the fluxes. Mixed conditions are located between these two extreme configurations. These lead to feedback effects and the emergence of complex dynamic behaviors [570].

a)

<figure><img src="image_89.png" /></figure>

b)

<figure><img src="image_90.png" /><figcaption>Figure 14.3: a) Illustration of muscle as an energy converter. The incoming energy flow $\Phi_+$ is converted into mechanical power $P = F_M I_M$ and a waste fraction $\Phi_-$, $I_M$ is the so-called *metabolic intensity*. b) Plot of the system's response under varying metabolic intensities $I_M$. The response extends from the basal resting point to the point of exhaustion, via the point of maximum work production.</figcaption></figure>

If we apply this description to the case of living organisms that have been reduced to chemical conversion machines, we obtain a thermodynamic formalism (see Figure 14.3) that regains the phenomenological description of the muscular response proposed by Hill [571, 572]. In Hill's phenomenology, the metabolic force $F_M$ and the contracting velocity $v$ are linked by three constants represented by the equation $F_M = \frac{c}{v+b} - a$. The thermodynamic formulation gives $F_M = \frac{(F_{iso}+R_{fb})I_T}{I_M+I_T} - (R_{fb}I_T + R_M I_M)$ where $I_M \propto v$. The thermodynamic approach gives us access to the physical meaning of the parameters i.e., $F_{iso}$ is the isometric force of the muscle, $I_T$ defines the threshold of acceptable metabolic intensity, $R_M$ is the viscous resistance to displacement and $R_{fb}$ the feedback resistance induced by the mixed conditions previously mentioned.

A proxy for the flow released by the muscle is the quantity of oxygen breathed in during ventilation [573]. To achieve an effort of a given intensity, the level of $O_2$ adjusts accordingly. Naturally, this quantity cannot grow indefinitely, and is limited by the absolute size of the organ that enables this exchange and by the relative size of this organ.
assumes the organism has perfect information on how the environment will change (including in response
to actions taken) over time. In an uncertain environment, this information is simply not available. An alter-
native is instantaneous optimization of growth rate at each time point but this is a shortsighted strategy
that excludes any partial information on future environmental states that the organism may have acquired
over the course of evolution. Such deterministic models may be suitable for deterministically changing en-
vironments, but cannot account for stochastic behaviors that may be advantageous to population growth
in uncertain environments.

In this chapter, it will be shown how principles of optimality can be formulated to study the behavior of organisms growing under uncertainty. Unlike the deterministic setting however, optimality will instead need to be defined in terms of probabilities and expected returns. Analogous to the general unification of deterministic models for cellular behavior using an optimal control theory framework, models including uncertainty are unified by the subject of *stochastic optimal control*. Beyond biology, this subject has wide-reaching applications to engineering but the most relevant analogy is with finance where stochastic strategies of portfolios diversification mirror stochastic strategies of cellular diversification. This will add a new economic analogy to the economic analogies of previous chapters.

13.2 Strategies to cope with uncertainty: a financial analogy

We will use the topic of as a recurring example throughout this chapter (Figure 13.1). When a clonal pop-
ulation of bacteria is exposed to an antibiotic, not all cells within the population are killed – a small sub-
population, although genetically identical to the rest, may nevertheless be in a distinct phenotypic state
that is growth-dormant and resistant to treatment (Figure 13.1A). While the peers of this dormant sub-
population previously grew well in the absence of antibiotic, upon exposure to treatment these growing
cells are killed, and only the dormant cells (the persisters) remain alive. In turn, when the remaining persis-
ters are transferred to an environment without antibiotic a large fraction is able to revert to the growing
state, allowing the population as a whole to survive. Remarkably, in this subsequent phase of growth
roughly the same small fraction of persisters is retained as before the treatment. Deterministic mod-
els based on short-term optimal growth cannot explain how part of a population adopts a slow-growing
state: they would predict that each cell should adopt the growing phenotype in absence of antibiotics.
Cells could have a mechanism to detect the presence of unfavorable environmental conditions and adopt
the persister phenotype as a response, but there are several experimental observations not explained by
such a mechanism [409]: (1) a fraction of persisters exists prior to antibiotic treatment; and (2) not all cells,
although genetically identical, adopt the persister phenotype. We will see that a more parsimonious de-
scription of persistence involves an optimization of long-term rather than short-term growth, which differs
when environmental conditions fluctuate.

Bacterial persistence is an example of , which more generally refers to the benefit of spreading resources
across multiple behavioral phenotypes to reduce the associated with investing all resources into any single
phenotype. Returning to the example of bacterial persistence, a natural question one may ask is: what
determines the precise fraction of persister cells (risk-avoiding, potentially low-reward phenotype) com-
pared to growing cells (risky, potentially high-reward phenotype) within a given population? This question
echoes a central question in financial investment: how should investors diversify their portfolio to max-
imize their capital in the context of uncertain returns? We will see that some of the same mathematical
arguments of optimality under uncertainty can be used to analyze these two problems, showing how the
optimal fraction of persisters is expected to depend critically on the probability to experience different
environmental states. The terms of the analogy are presented in Table 13.1.

A pure bet-hedging strategy assumes the absence of any direct information on the current environmental
state. Biologically, cells may sense signals or cues that encode varying degrees of information on their
current environment. For instance, in some populations, a larger proportion of persisters is found in
nutrient-poor environments compared to nutrient-rich, implying a direct relationship between shifts in
environment and switches between phenotypes. These sensing or signaling mechanisms can come with
associated costs however, imparted by the investment of cellular resources in, for example, the gene ex-
pression machinery. Thus, optimal cellular behavior in the face of uncertainty may be expected to involve
compared to the size of the individual.

For an individual, this is an intrinsic limitation on the ability to produce effort. So, depending on the size of the individual, which constrains its volume, the respiratory system must be optimized to maximize the flow of O₂. By comparing inter-species data and using a generic description, it is then possible to find an allometric law, as we shall see in this chapter.

In the next section, the respiratory system, and the lung as its central organ, will be studied in details in light of the concepts of organ optimization. Indeed, the lung, its structure, its functioning, its efficiency, are all the result of a series of optimization under constraints that shaped the organ through evolution.

## 14.2. The lung as a model organ for optimization under constraints

At the core of the respiration process, the lung is the organ that connects the ambient air to the blood, allowing to transport oxygen from the ambient air to blood and carbon dioxide from blood to the ambient air. The needs of the body in oxygen and carbon dioxide, the respiratory gases, determines the lung function, which is based on a complex geometrical structure and on several physical and chemical processes.

### 14.2.1. Lung morphology, a complex structure

A basic description of the lung structure would consist in dividing it in two parts: the bronchial tree and the exchange surface with blood.

The function of the bronchial tree is limited to the transport of the respiratory gases and no exchange occurs in this part of the lung. The bronchial tree forms a cascade of bifurcating airways with cylindrical shapes. There is an average of seventeen successive bifurcations in the human lung. The trunk of the tree is called the trachea; it is connected to the ambient air through the tracheo-pharyngeal pathway. The leaves of the tree are called the terminal bronchioles; they are connected to the exchange surface with blood. At each bifurcation the size of the airways is decreasing, with a tracheal diameter of about 2 cm and a diameter of the terminal bronchioles of about 0.3 to 0.5 mm.

The exchange surface with blood consists in a foam-like structure that is an assembly of exchange units called the acini. Each acinus is also shaped as a bifurcating airway tree, but the size of the airways is conserved at the bifurcations. There is an average of six successive bifurcations in a typical acinus. The acinar airways are called the alveolar ducts and their walls are garnished with bubble-like structures, the alveoli. The alveoli walls are mainly blood capillaries, called pulmonary capillaries, and they are the location of the respiratory gas exchanges. Each terminal airway of the bronchial tree feeds an average of two acini. The auto-similar, multi-scaled structure of the bronchial tree and of the acini allows the lung to contain a very large exchange surface that is folded in the thorax. In a typical human, the exchange surface is about 70-100 m² [540].

Since the morphology of the lung is complex, it becomes necessary to make assumptions in order to have a simple model while conserving the principal geometrical properties. Our model is then based on the assembly of self-similar trees with cylindrical branches and symmetric bifurcations that mimic the two functional zones (see Figure 14.4). To account for the core geometrical properties of the lung, we assume that the dimensions of the branches in the conductive tree decreases from one generation to the next with a ratio $h = (\frac{1}{2})^{\frac{1}{3}}$ [574, 575, 576], while in the acinus we assume that the size of the bronchi remains constant [574]. Note that the airways spatial distribution such as the branching angles or the orientations of the branching planes is not taken into account in our model since it is not really relevant for the computation of oxygen transport in the lung.
<figure><img src="image_119.png" /><figcaption>Figure 13.1: Bacterial persistence as an example of a cellular strategy to cope with uncertainty in environmental conditions – A) Cells in a genetically identical population can display one of two distinct phenotypes that are associated with growth (pink cells) or dormancy (blue cells) in the absence of antibiotics. Only the dormant cells survive (persist) when exposed to antibiotics, and can transition back to the growth phenotype so that the population as a whole resumes growth in the absence of antibiotic. B) In this simplified model of bacterial persistence, the strategy *u* over two responses (phenotypes) $R_{growth}, R_{dormant}$ depends on environmental states $E_{low}$ and $E_{high}$, corresponding to low and high levels of the antibiotic, respectively. The occurrence of the states $E_{low}$ and $E_{high}$ is governed by probabilities $p(E_{low})$ and $p(E_{high})$, respectively. C) The multiplicative rates $f(R|E)$ associated with phenotypes $R_{growth}, R_{dormant}$ depend on environmental conditions, so that $f(R|E)$ can be represented in matrix form. The resulting optimal strategy $u_d$ corresponding to the fraction of dormant cells in the population in turn depends on the probabilities of the environmental state $E$. An analogy with Kelly betting is illustrated on the right-hand side, where the probabilities of a horse winning a race, the odds provided by a bookmaker and the optimal betting strategy are identified with $p(E)$, $f(R|E)$ and $u(R|E)$, respectively, as displayed in Table 13.1.</figcaption></figure>

a trade-off between stochastic (e.g., bet-hedging) and deterministic (e.g., signaling) mechanisms that balance benefit to cost in a manner that depends on evolutionary context. Other trade-offs may also exist regarding reward versus risk associated with a particular cellular strategy. Analogously, financial investors face trade-offs when using incomplete information on the current state of the market and developing an investment strategy based on the level of risk they are willing to incur.

## 13.3 Modeling cells growing in uncertain environments

We begin with a simple model of persistence before introducing a more general framework. This simple model assumes that bacterial cells experience an alternation of low and high antibiotics environments and can adopt two physiological states, growing or dormant (Fig. 13.1). The dormant cells are unable to replicate but persist in either high- or low-antibiotics environments while growing cells always divide when antibiotics are low in concentration but die when they are high. Mathematically, this is described by $f(R, E)$, the number of descendants of a cell with phenotype $R$ in environment $E$: $f(R = \text{dormant}, E = \text{low}) = f(R = \text{dormant}, E = \text{high}) = 1$, while $f(R = \text{growing}, E = \text{low}) = 2$ and $f(R = \text{growth}, E = \text{high}) = 0$. In absence of sensing mechanism, we consider that the fraction of dormant cells, $u_d \equiv u(R = \text{dormant})$, is a fixed quantity that only possibly evolves on very long time scales. The population thus grows by a global factor $A_{high} = f(R = \text{dormant}, E = \text{low})u_d$ if the environment is high antibiotics and by a factor $A_{low} = f(R = \text{dormant}, E = \text{low})u_d + 2f(R = \text{growing}, E = \text{low})(1 - u_d)$ if it is low
<figure><img src="image_91.png" /><figcaption>Figure 14.4: Illustration of the lung model used in this work. The tree in beige mimics the bronchial tree, where oxygen and carbon dioxide are only transported along the branches. The tree in blue mimics the acini, where the respiratory gases are transported along the branches. They are also captured by the alveoli that cover the walls of the branches.</figcaption></figure>

### 14.2.2. Lung dynamics: where physics enters the play

The transport of the respiratory gases to and from blood involves a combination of physical processes which ensure that the needs of the body in respiratory gases are fulfilled.

**Diffusion: no energy costs, but too weak.** As blood entering the pulmonary capillaries has an oxygen partial pressure lower than the oxygen partial pressure in the alveolar air, oxygen flows to the blood by the process of diffusion that tends to balance the partial pressures between blood and the alveolar air. For the lung's point of view, the blood acts as an oxygen sink. The transport of carbon dioxide in the lung relies on the same processes that that for oxygen, except that blood flowing in the alveoli membranes acts as a source of carbon dioxide. The diffusion process is passive in the lung i.e., no energy is spent by the organ to perform the transport. Notice that this is not true from the pulmonary blood circulation point of view, as blood has to be incessantly renewed to maintain the respiratory gas partial pressure difference between the alveolar air and the blood. However, at the metabolic time scale, the diffusion process has a limited range in the airway tree.

Were the transport of the respiratory gas only based on diffusion, the lung could not maintain the respiratory gas flow at a level compatible with the mammals metabolisms. The reason behind this limitation stands in the size of the airway tree. The pathways from the ambient air to the respiratory zone are too long and narrow for the diffusion to provide gas flows compatible with the metabolism of mammals. In the case of the human lung, the typical length of these pathway is of about $L_p = 30$ cm [575]. The characteristic time $t_p$ for an oxygen molecule to travel by diffusion through all such a pathway can be estimated using a dimensional analysis. Using $L_p$ and the diffusion coefficient of oxygen in air $D = 0.2$ cm²·s⁻¹ [577], the order of magnitude of $t_p$ can be estimated with:

$$t_p = \frac{L_p^2}{D} \approx 4500 \text{ s} = 1 \text{ hour and } 15 \text{ minutes!}$$

Hence, a pure diffusive transport of the respiratory gas cannot fit the mammals needs. Actually, in human, the order of magnitude of the length $L_D$ traveled by diffusion during the typical time of inspiration i.e., $t_i = 2$ seconds, is $L_D = \sqrt{D \times t_i} \approx 6.3$ mm. Thus, in the resting human, diffusion can transport oxygen from the terminal bronchioles to the nearby exchange surface. However, at a time scale compatible with the metabolism, diffusion cannot reach the upper part of the bronchial tree. It cannot either reach the deeper parts of the respiratory zone, which is non active at
<table><caption>Table 13.1: Analogy between bet-hedging in biological populations and diversification strategies in Kelly's gambling and finance. The common problem in each case is an uncertain environment that makes it impossible to anticipate which phenotype or investment is optimal for future growth. In finance, the "population" is constituted by the capital which is distributed across different options (different horses of a race or different stocks of a stock market). The main limitation of the analogy is that information is not processed centrally in biological populations but at the level of each individual, with therefore no equivalent to a gambler or investor. The notations are introduced in the main text for the biological problem and in Box 13.B for the gambling problem.</caption><thead><tr><th>Biology</th><th>Gambling</th><th>Finance</th></tr></thead><tbody><tr><td>Individual</td><td>Capital unit</td><td>Currency unit</td></tr><tr><td>Environment p(E)</td><td>Race results p(x)</td><td>Market state</td></tr><tr><td>-</td><td>Gambler</td><td>Investor</td></tr><tr><td>Phenotype decisions u(R)</td><td>Bets b(x)</td><td>Investment strategy</td></tr><tr><td>Multiplicative rate f(R, E)</td><td>Odds o(x)</td><td>Immediate return</td></tr><tr><td>Environmental cue P(S|E)</td><td>Side information P(y|x)</td><td>Side information</td></tr><tr><td>Population growth rate Λ</td><td>Long-term return W</td><td>Long-term return</td></tr><tr><td>Extinction probability</td><td>Probability of bankruptcy</td><td>Probability of bankruptcy</td></tr><tr><td>Growth rate variance σ<sup>2</sup></td><td>Growth rate variance σ<sub>W</sub><sup>2</sup></td><td>Volatility</td></tr><tr><td>Population size N<sub>t</sub></td><td>Capital C<sub>t</sub></td><td>Capital</td></tr></tbody></table>

antibiotics. Finally, the environment fluctuates randomly, with a probability $p_a$ to have high antibiotics and a probability $1 - p_a$ to have low antibiotics. Over a large number $T$ of generations, a population therefore experiences in average $p_a T$ periods of high antibiotics and $(1 - p_a) T$ periods of low antibiotics. As further explained below, the population size $N_T$ after $T$ generation is hence expected to globally grow as $N_t = (A_{high})^{p_a T} (A_{low})^{(1-p_a)T} N_0$. This corresponds to an exponential growth (or decay) of the form $N_T = e^{\Lambda T} N_0$ with a long-term growth rate $\Lambda$ given by $\Lambda = p_a \ln u_d + (1 - p_a) \ln(u_d + 2(1 - u_d))$.

Two bacterial populations which have different "strategies" $u_d$ will then have different growth rates $\Lambda(u_d)$. The optimal strategy which maximizes $\Lambda(u_d)$ is therefore when the probability $u_d$ to adopt the dormant state is

$$u_d = \begin{cases} 2p_a, & \text{if } 0 < p_a \le 1/2. \\ 1, & \text{if } 1/2 < p_a \le 1. \end{cases}$$

The interesting case is when $p_a < 1/2$, otherwise antibiotics is so often high that the population cannot grow. In this case, we find that a limited fraction of the population should be in the dormant state and that this optimal fraction depends on the frequency $p_a$ at which high antibiotics occurs.

This example can be extended to an arbitrary number of environmental states $E$ and phenotypic (response) states $R$ and to the presence of cues collected from the environment. In general, the states and cues may take discrete (as in the above example) or continuous values. The "strategy" of a cell may then be described by its probability $u(R)$ to adopt a particular phenotype $R$. This strategy depends on the environment if some signal $S$ is perceived, in which case the strategy takes the form of a conditional probability $u(R|S)$ satisfying

$$\sum_R u(R|S) = 1, \quad \text{with} \quad u(R|S) \ge 0$$

for each possible signal $S$. For the example of bacterial persistence, $u(R = \text{dormant}|S)$ may be the fraction of cells adopting a dormant phenotype within the population of cells with intracellular antibiotics concentration $S$. The fraction of growing cells would then be given by $u(R = \text{growing}|S) = 1 - u(R = \text{dormant}|S)$. By comparison, Figure 13.1B illustrates a model where $u(R|E)$ depends directly on the environmental state $E$. In finance, $u(R|S)$ would correspond to the fraction of the capital that an investor allocates to asset $R$ when receiving incomplete information $S$ on the current market state $E$. More generally, we may also consider that the probability to adopt a phenotype $R_t$ at time $t$ depends on the phenotype $R_{t-1}$ adopted at time $t-1$ by the cell or its parent, which would be described by $u(R_t|S, R_{t-1})$ or $u(R_t|S_t, R_{t-1})$ to indicate that the signal $S_t$ is obtained at time $t$.
rest. Actually, this last phenomenon, called the screening effect [577], plays a crucial role in the lung. It is described in details later in this chapter.

More generally, the limited spatial range of diffusion has many consequences on the living systems. An emblematic example is its role on the size limitation of insects [578], where diffusion in the tracheal tubes is the only mean of respiratory gas transport. It participates to the explanation of why the increased atmospheric oxygen concentration during the Palaeozoic era allowed insects to be larger than today as, following Fick's law, the diffusive flow is proportional to the gradient of partial pressure between the ambient air and the inner body.

**Convection : the rescuer.** We have seen that the diffusion process is too weak to transport the respiratory gas through the whole airway tree. In the absence of other transport mean, the oxygen partial pressure in the lung would decrease and the flow of oxygen to blood would drop. Similarly, the carbon dioxide partial pressure would increase and prevent the exchanges with blood to occur. Consequently, the air in the lung has to be renewed in order to expel the excess of carbon dioxide and to refresh the inhaled air volume with oxygen. This phenomenon is called the ventilation. The ventilation is a dynamic i.e., time-dependent, process based on the succession of inhalation and exhalation of a volume of air, the tidal volume, at a given rate, the breathing frequency.

Ventilation is performed thanks to a set of muscles that surround the lung and modify its volume. At rest regime, the main acting muscle is the diaphragm, located at the base of the lung. By first pulling onto the lung, this muscle deforms the lung tissues, creating a negative pressure drop and the transport a volume of ambient air inside the lung; this is the inspiration phase. At rest, the elastic energy stored in the tissues during the inspiration phase allows for a passive recoil of the lung and a volume of air equal to the volume inhaled is expelled; this is the expiration. Then the cycle repeats following the same procedure, at least at rest. Since the duration of a breath cycle for a resting human is about four to five seconds, a human performs, on average, about six to seven hundred millions breaths during her/his lifetime.

**modelling the oxygen transport.** The transport of oxygen in the lung is then driven by three phenomena: diffusion, convection by the airflow and exchange with blood through the alveoli walls in the alveolar ducts. The partial pressure of oxygen averaged over the lumen area is transported along the longitudinal axis *x* of the airway. Hence, in each airway of our idealized lung, the mean partial pressure of oxygen *P* over the airway section follows,

$$ \frac{\partial P}{\partial t} - D \frac{\partial^2 P}{\partial x^2} + u \frac{\partial P}{\partial x} = \beta (P_{\text{blood}} - P), \qquad (14.2) $$

where *D* is the oxygen diffusion coefficient, *u* is the velocity of the airflow, *β* is a reactive term and *P*<sub>blood</sub> is the partial pressure of oxygen in the capillary blood. The reactive term *β* mimics the exchanges with blood through the alveolar membrane. This coefficient depends on the diffusion coefficient of oxygen in water, on the solubility coefficient of oxygen in water, on the thickness of the alveolar-capillary membrane, and on the radius of the alveolar duct. It is equal to zero in the bronchial tree since no exchange with the blood happens in this part of the lung and is positive in the acini. The oxygen partial pressure in blood is determined by assuming that the flow of oxygen leaving an alveolar duct through the alveolar-capillary membrane is equal to the flow of oxygen captured by blood, accounting for the oxygen captured by hemoglobin and for the oxygen dissolved in plasma [540]. Finally, all generations are linked through bifurcations by assuming continuity between generations and conservation of the quantity of oxygen at each bifurcations.

Experimental methods 14.C: Conditions for the numerical simulations

Our model takes as input the ventilation parameters: the tidal volume *V*<sub>T</sub> (in mL) and the breathing frequency *f*<sub>b</sub> (in min<sup>-1</sup>) and outputs the mean amount of oxygen exchanged with blood over a respiratory cycle. To validate our model, we performed computations at rest by assuming that a human breathes around 12 times per minute and inhales around 500 mL of air for each breathing cycle. With these parameters, our transport model gives an oxygen flow exchanged with blood of 230 mL·min<sup>-1</sup>, which is close to the average physiological value of 250 mL·min<sup>-1</sup> [540].