Bacterial metabolic heterogeneity: from stochastic to deterministic models.

We revisit the modeling of the diauxic growth of a pure microorganism on two distinct sugars which was first described by Monod. Most available models are deterministic and make the assumption that all cells of the microbial ecosystem behave homogeneously with respect to both sugars, all consuming the first one and then switching to the second when the first is exhausted. We propose here a stochastic model which describes what is called "metabolic heterogeneity". It allows to consider small populations as in microfluidics as well as large populations where billions of individuals coexist in the medium in a batch or chemostat. We highlight the link between the stochastic model and the deterministic behavior in real large cultures using a large population approximation. Then the influence of model parameter values on model dynamics is studied, notably with respect to the lag-phase observed in real systems depending on the sugars on which the microorganism grows. It is shown that both metabolic parameters as well as initial conditions play a crucial role on system dynamics.


Introduction
Described for the first time by Monod [8], the diauxic growth consists in a biphasic growth in a bacterial population consuming two different sugars in a closed medium. The corresponding curve of biomass density at the macroscopic scale shows two distinct exponential phases separated by a "plateau" called lag-phase. The explanation proposed by Monod is that the preferred sugar (which is in some sense "easier" to metabolize) is consumed first while the metabolic pathway allowing the consumption of the second one is suppressed. When the concentration of the first sugar becomes low enough, this repression is lifted. Then, the microorganism may produce the enzymes necessary to metabolize the second sugar: this is the lag-phase. The second exponential growth is observed until the second sugar is eventually consumed. Until recently, it was admitted that the explanation given above was homogeneous within the cell population in the sense that each individual adopted exactly the same behavior at the same time: each cell first consumed the sugar that was "easiest to metabolize" first, then the other one after a duration corresponding to the lag-phase. Such an assertion implies that the latency time would simply be a constant depending only on the sugars involved. In order to better understand this phenomenon and test hypotheses, many models of diauxic growth have been proposed in the literature ( [7,11,10]). All such models have in common to make the hypothesis that each cell of the microorganism under consideration exhibits the same behavior with respect to both substrates at a given time. In addition, most approaches make use of deterministic models that are not suited for low biomass densities.
However, recent investigations suggest that lag phases are controlled by the inoculum history and organized with heterogeneity among individual cells (Bertrand [4]). This fact was called "metabolic heterogeneity". Takhaveev and Heinemann [9] suggested that this heterogeneity could be induced by mechanisms linked to ecological factors, gene expression, and other inherent dynamics, or by interaction between individuals, which all also depend on environment changes.
In this paper, following the idea that there is an intrinsic heterogeneity of cells within the ecosystem -which yields a metabolic heterogeneity -we develop a stochastic model of diauxic growth. More precisely, we propose a model of a batch culture of a pure strain growing on two different sugars. In this model, the metabolic heterogeneity is modeled via the possible emergence of a subpopulation able to consume the second sugar while the first one is not yet totally consumed. In other words, all cells do not exhibit the same behavior with respect to each substrate at a given time. To be as close as possible to the observations, the model accounts for the fact that in such situations the acetate produced -which is a growth-inhibiting metabolite -is co-consumed by each cell, as shown by Enjalbert et al. [5]. One goal of this work is to link both stochastic and deterministic approaches in order to explain the observations available at different scales, and to study the main parameters that control the length of the lag-phase.
The paper is organized as follows. First, the stochastic model is presented. Secondly, its behavior for large populations is approximated, allowing us to write a model consisting in a set of deterministic differential equations. Then, the model is used to investigate the role of a number of model parameters and of initial conditions on the substrate consumption dynamics and on the length of the lag-phases. Eventually the main conclusions and perspectives are drawn. An appendix provides some additional information on proofs and simulations.

The stochastic model
First and foremost, let us introduce the parameter K > 0 that scales the initial number of individuals. Indeed, the population size varies widely between different kinds of bacterial cultures, and may range from a few individuals in microfluidics (then K is very small) to billions or more in fermenters (then K is very large).
Let us consider two different substrates, Sugar1 which is preferential and Sugar2, and a stochastic population of bacteria split into two compartments constituted respectively of Sugar1 consumers and of Sugar2 consumers. Let N K 1 (t) and N K 2 (t) denote respectively the numbers of individuals in each compartment and Here we are introducing the specific scaling in which K can be seen as proportional to the carrying capacity of the medium and 1/K as proportional to the individual biomass, and in order to capture the two subpopulation densities we introduce the rescaling The mass concentration of each sugar is described by a continuous process , t ≥ 0 , which corresponds in the same order to Sugar1 and Sugar2. We also take into account the mass concentration A K (t) of a metabolite produced during the consumption of sugars by each individual and co-consumed with them. As an illustration, we may consider a mixed medium with glucose and xylose as Sugar1 and Sugar2, and acetate as the metabolite.
We describe the complete culture medium by the Markov process evolving as follows.
• Demography. An individual growing on Sugar1 divides at rate b 1 (R K 1 (t), A K (t)) due to Sugar1 and metabolite co-consumption. Likewise, an individual growing on Sugar2 divides at rate b 2 (R K 2 (t), A K (t)) due to Sugar2 and metabolite co-consumption. This results in the jumps • State transitions. An individual growing on Sugar1 switches its state in order to consume Sugar2 at rate η 1 (R K (t)), which depends on both resources since this is inhibited by the catabolic repression due to Sugar1, the preferential sugar. Likewise, an individual growing on Sugar2 switches its metabolic state to consume Sugar1 at rate η 2 (R K 1 (t)) which depends only on the abundance of Sugar1. This results in the jumps • Resource dynamics. These are linked to the biomass and metabolite synthesis by the biochemical reactions that happen inside the batch. An individual growing on Sugar1 consumes continuously a small amount µ 1 (R K 1 (t), A K (t))/q 1 K of Sugar1 and produces continuously a small amount θ 1 µ 1 (R K 1 (t), A K (t))/K of metabolite per unit of time, before dividing as a result of this consumption. Similarly, an individual growing on Sugar2 consumes continuously a small amount µ 2 (R K 2 (t), A K (t))/q 2 K of Sugar2 and produces continuously a small amount θ 2 µ 2 (R K 2 (t), A K (t))/K of metabolite per unit of time, before dividing as a result of this consumption. Finally, each individual consumes a small amount µ 3 (A K (t))/q 3 K of metabolite per unit of time before dividing as a result of this consumption. This leads us to describe the resource dynamics by the dynamical system The typical situation we will consider is the following. The initial conditions satisfy in which (n 0 , r 0 ) = (n 0 1 , n 0 2 , r 0 1 , r 0 2 ) is fixed. The above rate functions involve Monod-type and classic inhibition functions and are of the forms In addition, the choice b j (r j , a) = µ j (r j , a) + µ 3 (a) , j = 1, 2 , ensures a conservation law on average: Parameter Biological signification Default value for simulations µ 1 Maximal growth rate on Sugar1 6.50e-01 κ 1 Monod constant on Sugar1 3.26e-01 λ Inhibition coefficient due to the metabolite 4.70e-01 µ 2 Maximal growth rate on Sugar2 5.70e-01 κ 2 Monod constant on Sugar2 4.68e-01 µ 3 Maximal growth rate on the metabolite 1.47e-01 κ 3 Monod constant on the metabolite 6.45e-01 η 1 Maximal switching rate from Sugar1 to Sugar2 2.04e-03 k 1 Regulation coefficient of the Sugar1 to Sugar2 transition 1.20e-02 k i Inhibition coefficient of the Sugar1 to Sugar2 transition 1.03e-03 η 2 Maximal switching rate from Sugar2 to Sugar1 6.60e-01 k 2 Regulation coefficient of the Sugar2 to Sugar1 transition 4.50e-02 q 1 Individual yield on Sugar1 5.50e-01 θ 1 Metabolite yield on Sugar1 6.00e-01 q 2 Individual yield on Sugar2 4.50e-01 θ 2 Metabolite yield on Sugar2 5.60e-01 q 3 Individual yield on the metabolite 2.50e-01 V Bioreactor volume 1.0 n 0 Initial subpopulation densities (2.80e-01, 0.0) r 0 Initial sugar concentration (8.15, 9.05) To illustrate this model, we use the parameters described in Table 1 taken from recent batch experiments (Barthe et al. [3]) for the sugars glucose and xylose. Figure 2 shows that this model is able to predict the diauxic growth observed by Monod (see Figure 1). This can be observed even for a small number of individuals. We additionally observe that the trajectories oscillate randomly for small K and become smoother as K becomes larger. This observation will be developed in the next section, in which the stochastic model will be shown to be approximated by a deterministic model when K increases to infinity.

Large population approximation
The amplitude of any jump occurring in the population is bounded by a factor of the weight 1/K attributed to a single individual and hence has a variance of order 1/K 2 . Moreover, the mean number of jumps per unit of time is of order K, the order of magnitude of the number of individuals. Heuristically, for a large population the process should approach a limit deterministic continuous process dictated by the mean values, with random oscillations around this limit corresponding to variances of order 1/K and hence to standard deviations of order 1/ √ K. We can build on these heuristics and prove that the stochastic model is indeed approximated by a deterministic model in the limit of large K. This yields the following deterministic limit.

(5)
Then the stochastic process (n K (t), R K (t), A K (t)) t≥0 is approximated by (n(t), r(t), a(t)) t≥0 for large K in the sense that This theorem allows to explain on a rigorous basis the observations we have made on the simulations in the previous section. Before discussing the proof methods, let us address the important question of the range of validity of the approximation.
For small K and most notably for populations consisting of a few individuals, the deterministic system is not a good approximation of the stochastic model and does not provide a pertinent model for the population. On the contrary, when K is large enough for the approximation to be accurate, the deterministic system provides a pertinent model on which theoretical studies and numerical computations can be performed for qualitative and quantitative investigations on the population. Therefore, it is fundamental to obtain a precise evaluation of the size of K required for the approximation to be tight and to assess the error made in terms of K. The heuristics given before the theorem indicate that that the error terms should be of order 1/ √ K. Under adequate assumptions on the initial conditions, this can be made rigorous through a functional central limit theorem: the process √ K n K (t), R K (t), A K (t) − (n(t), r(t), a(t)) , t ≥ 0 , converges as K goes to infinity to a Gaussian process of Ornstein-Uhlenbeck type, with mean and covariance structure expressed solely in terms of the limit process (n(t), r(t), a(t)) t≥0 and of the variance of the jumps in a sufficiently explicit fashion to be well evaluated. This allows to evaluate the minimal size of K required for a tight approximation and to provide confidence intervals on this, as well as the possibility for intermediate sizes of K to simulate the deterministic limit process and add to it fluctuations simulated according to this Gaussian process in order to obtain a tighter approximation. The proofs of Theorem 3.1 and of the functional central limit theorem build on the heuristic explanation given before the theorem using probabilistic compactness-uniqueness methods. Ethier and Kurtz [6] is a classic book on the subject, and Anderson and Kurtz [1] and Bansaye and Méléard [2] provide pedagogical expositions well suited to the present field of application.
We illustrate these convergence results in Figure 3, by the simulations of a hundred trajectories of the total biomass for the stochastic and the limiting model, for three increasing values of the scale parameter K.  Figure 3: Ten independent stochastic trajectories, the empirical mean over a hundred independent trajectories, and the deterministic limit simulated for each of the total populations in the cases K = 10, K = 100, and K = 1000.

Heterogeneity and lag-phase sensitivity
As shown in Figure 4, the model is able to capture the heterogeneity of the population observed by biologists as well as the diauxic growth at the level of the total population size highlighted by Monod. We present several simulations, each with either one varying metabolic parameter or a varying initial condition. These varying factors have been carefully chosen after a numerical exploration showing that they have a strong influence on the lag-phase duration, and in particular that they allow to reproduce dynamics similar to the ones obtained by Monod, see

Influence of the metabolic parameters on the diauxie length
We here consider two major parameters: • The maximal switching rateη 1 from Sugar1 consumption to Sugar2 consumption, which describes individual switches on high Sugar2 medium after Sugar1 exhaustion. As this parameter increases, the transitions become more frequent and the lag phase shorter.
• The inhibition coefficient k i of the Sugar1 to Sugar2 transition, which describes the catabolic repression due to the abundance of the preferential sugar. Due to its role in the inhibition functions, as this parameter decreases the transitions become more frequent and the lag phase shorter.
Our simulations are represented in Figure 5.

Influence of the initial conditions on the diauxie length
The initial condition n(0) = (n 1 (0), n 2 (0)) of the population also has a strong influence on the duration of the lag-phase. Indeed, if the population density of Sugar2 consumers is significant with respect to the density of Sugar1 consumers in the initial condition, then the lag-phase will be short. Our simulations are represented in Figure 6.

Conclusion
In this paper, we proposed a stochastic model of diauxic growth of a microorganism on two different sugars. The model assumes that the individuals preferentially consume one of the sugars while the metabolic pathway allowing the consumption of the second one is repressed until the first sugar is exhausted. To account for the fact that all individual do not behave homogeneously with respect to the consumption of sugars -which is called metabolic heterogeneity -it is supposed that some individuals can switch their metabolism in such a way they can consume the second sugar while the first one is not totally exhausted. Thus the model involves two different subpopulations: the first one which grows on the first sugar, and the second one, which emerges from the first subpopulation and consumes the second sugar. In addition, three resource variables with continuous dynamics are added: the two sugars, and the intermediate metabolite which is produced when the sugars are consumed and then re-consumed by both subpopulations. Then, the deterministic model that approximates the stochastic model dynamics is derived using a large population approximation. Using parameter values that are supposed to be close to those we can find in real experiments, for instance when E. coli grows on both glucose (the preferential sugar) and xylose, we performed a number of simulations in order to investigate the influence of the most important parameters on the model dynamics. Further, we show the importance of the weighting factor K, which allows us to understand what is the population size starting from which the deterministic model can be used to approximate the stochastic model dynamics. Finally, it is shown that several parameters, such as the maximal switching rateη 1 from Sugar1 to Sugar2 consumption and the inhibition coefficient k i of the Sugar1 to Sugar2 transition, as well as the initial conditions of the system significantly influence the lag-phase, allowing us to pave the way and to suggest strategies to minimize the lag-phase in practical experiments.
0, · · · , α p (x) ≥ 0, and The strong Markov property yields interesting consequences for the construction of the process. The future of the process after each jump is independent from its past given the new state. Thus, in order to construct the process it is sufficient to be able to do so from time 0 until the first jump instant, and then iterate the procedure by considering each jump instant as a new time origin. Moreover, starting at time 0 the probability that the process (X K (t)) t≥0 has not jumped yet at time t > 0 is given in terms of the rate function α by This allows to construct the first jump instant as follows. If the non-decreasing continuous process (Λ(t)) t≥0 and its left-continuous inverse (Λ −1 (t)) t≥0 are defined by and D is an exponential random variable of parameter 1, then Hence, we can simulate the first jump instant T 1 of the process by taking T 1 = Λ −1 (D) while simultaneously constructing the process on [0, T 1 ). If X K (T 1 −) = x then X K (T 1 ) = x + h for a jump amplitude h drawn according to π(x, dh). Using this construction directly for an actual simulation raises several issues. The first problem is that we must be able to simulate the process (X K (t)) t≥0 up to the first jump instant. In the present situation this consists in simulating the components (R K (t), A K (t)) t≥0 of the Markov process (1) by solving the differential system (2) in which the other components of (1) remain constant between jumps. This cannot be done exactly but can be approximated numerically quickly and with precision.
The second problem is that simultaneously to (R K (t), A K (t)) t≥0 we must be able to compute the integral Λ(t) and its inverse Λ −1 (t) defined in (7). This can be done numerically but is often costly in computer time and inefficient. This has a practical solution which we proceed to describe. We introduce a functionα such that α ≤α and that the correspondingΛ and that Λ −1 defined similarly to (7) are simpler to compute than Λ and Λ −1 . We simulate the process (X K (t)) t≥0 by an acceptance-rejection method which proposes a jump from state x at rateα(x) and accepts it with probability α(x)/α(x) and else rejects it. There are various ways to justify that this construction is correct. One of these is to consider the rejection as the introduction of a jump of amplitude 0 taken at the excessive rateα(x) − α(x) (the process does not actually jump, and this is called a "fictitious jump") and reason as above. The simplest situation is when the dominating functionα is a constant. Then the true jump instants of (X K (t)) t≥0 constitute a thinning of a Poisson process of constant intensityα, which can be easily simulated, in which a jump instant of this Poisson process taken when X K (T 1 −) = x is accepted with probability α(x)/α(x).