Bacterial metabolic heterogeneity: from stochastic to deterministic models

1.   Introduction

Described for the first time by Monod ^[1], 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 (Figure 1). 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 ^[2,3,4]. 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 ^[5]. This fact was called "metabolic heterogeneity". Takhaveev and Heinemann ^[6] 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 that consumes 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. ^[7]. 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.

2.   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). Letting $K$ increase to infinity allows to make the link between the stochastic model and the deterministic limit model for large populations

Let us consider two different substrates, Sugar 1 which is preferential and Sugar 2, and a stochastically described population of bacteria split into two compartments constituted respectively of Sugar 1 consumers and of Sugar 2 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 per unit of volume we introduce the rescaling

The mass concentration of each sugar is described by a continuous process

which corresponds in the same order to Sugar 1 and Sugar 2. 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 Sugar 1 and Sugar 2, and acetate as the metabolite: Sugars are consumed sequentially, and when the preferential sugar glucose is abundant the xylose consumers switch to consume it; the reverse transition is more complex and requires that an activated xylR protein binds to the DNA of each switching individual. The activation of this protein is caused by the presence of xylose but inhibited by glucose due to the catabolic repression.

We describe the complete culture medium by the Markov process

evolving as follows.

● Demography. An individual growing on Sugar 1 divides at rate $b_1(R^K_1(t), A^K(t))$ due to Sugar 1 and metabolite co-consumption. Likewise, an individual growing on Sugar 2 divides at rate $b_2(R^K_2(t), A^K(t))$ due to Sugar 2 and metabolite co-consumption. This results in the jumps transitions

● State transitions. An individual growing on Sugar 1 switches its state in order to consume Sugar 2 at rate $\eta_1(R^K(t))$, which depends on both resources since this is inhibited by the catabolic repression due to Sugar 1, the preferential sugar. Likewise, an individual growing on Sugar 2 switches its metabolic state to consume Sugar 1 at rate $\eta_2(R^K_1(t))$ which depends only on the abundance of Sugar 1. This results in the jumps transitions

● Resource dynamics. These are linked to the biomass and metabolite synthesis by the biochemical reactions that happen continuously inside the batch. An individual growing on Sugar 1 consumes a small amount $\mu_1(R^K_1(t), A^K(t))/q_1K$ of Sugar 1 and produces a small amount $\theta_1\mu_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 Sugar 2 consumes a small amount $\mu_2(R^K_2(t), A^K(t))/q_2K$ of Sugar 2 and produces a small amount $\theta_2\mu_2(R^K_2(t), A^K(t))/K$ of metabolite per unit of time, before dividing as a result of this consumption. Note that it could happen that the metabolite inhibits the consumption of sugars as is the case for acetate in glucose-xylose consumption, see Enjalbert et al. ^[7]. Finally, each individual consumes a small amount $\mu_3(A^K(t))/q_3K$ of metabolite per unit of time by a separate independent pathway 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

Note that the terms of inhibition (due to the metabolite on each growth rate, and due to Sugar 1 on the switching rate) and their incidence are antagonistic. Indeed, decreasing any of the coefficients $\lambda$ and $k_i$ increase the repression on the corresponding mechanisms. In addition, since the sugars and the metabolite are co-consumed by independent pathways, the birth rates defined in the transitions (2.2) chosen as

ensure a conservation law on average:

To illustrate this model, we use the parameters described in Table 1 taken from recent batch experiments ^[8] for the sugars glucose and xylose.

The subpopulation densities per unit of volume (or biomass concentrations) and the resources mass concentrations are expressed in $g/L$. For the typical case that will be illustrated in simulations, the volume is one liter and the default initial conditions are given by

Figure 2 shows that this model is able to predict the diauxic growth observed by Monod (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.

3.   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/\sqrt{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.

Theorem 3.1. Let us assume that

and that $\sup_{K}\mathbb{E}\Bigl(\bigl\lVert \bigl(n^K(0), R^K(0), A^K(0)\bigr) \bigr\rVert \Bigr) < +\infty$. Let $(n(t), r(t), a(t))_{t\geq 0}$ be the unique solution with initial condition $(n^0, r^0, a^0)$ of the differential system\\

Then the stochastic process $(n^K(t), R^K(t), A^K(t))_{t\geq 0}$ is approximated for large $K$ by $(n(t), r(t), a(t))_{t\geq 0}$ 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/\sqrt{K}$. Under adequate assumptions on the initial conditions, this can be made rigorous through a functional central limit theorem: The process

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\ge0}$ 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 ^[9] is a classic book on the subject, and Anderson and Kurtz ^[10] and Bansaye and Méléard ^[11] 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$.

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

Furthermore, we are interested in how the metabolic parameters and initial conditions can influence the length of the lag-phase in the deterministic model. We perform a sensitivity analysis by considering an approximation of the lag duration defined as the time elapsed between the instant when the intake of preferential sugar and metabolite reaches a minimum threshold $\epsilon_1 > 0$, and the instant when a new phase of population growth is detected beyond another thresold $\epsilon_2 > 0$. Note that as in Eq (2.6) for the birth rate (2.5) in the stochastic model, we have the following conservation law for the deterministic limit

Then, since Sugar 2 is almost constant during the lag-phase, the intake of the preferred sugar and metabolite is proportional to the quantity

We then set

where

For all simulations of the deterministic model, we use the classical Runge-Kutta numerical scheme of order 4. Figure 5 illustrates this construction of the interval $[t_{\epsilon_1}, t_{\epsilon_1, \epsilon_2}]$ approximating the lag phase.

We present several simulations in which we vary the values of the maximal switching rate $\bar{\eta}_1$, the inhibition coefficient $k_i$ and the initial conditions in different intervals and plot the corresponding lag duration approximation. Figure 6 reveals that as the switching rate $\bar{\eta}_1$ increases, the lag phase become considerably shorter. This is particularly noticeable for small values of this parameter. This is also the case for the inhibition coefficient $k_i$ that affects the lag phase considerably for small values. Hence, the lag phase sensitivity is significant for strong inhibition. Finally, we consider varied initial conditions for the subpopulations which all have the same total population size, and observe that the lag phase duration seems to be well correlated to the initial proportion of sugar1 consumers. This situation is interesting in order to understand how the transplantation from a first culture medium to another could affect the lag.

5.   Conclusions

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 Escherichia 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 $\bar{\eta}_1$ from Sugar 1 to Sugar 2 consumption and the inhibition coefficient $k_i$ of the Sugar 1 to Sugar 2 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.

Acknowledgments

We warmly thank our biologist colleagues Manon Barthe, Muriel Cocaign and Brice Enjalbert who have brought this problem to our attention and shared fruitful discussions with us. This work has been supported by the Chair "Modélisation Mathématique et Biodiversité" of Veolia Environnement-École Polytechnique-Muséum national d'Histoire naturelle-Fondation X and by the ANR project JANUS (ANR-19-CE43-0004-01).

Conflict of interest

The authors declare no conflict of interest.

A. Appendix

A.1. Proof of main results

Let us comment the proof of Theorem 3.1 which can easily adapted from the results in ^[9,10,11]. Let us firstly note that we can express the stochastic process $(n^K(t), R^K(t))_{t\geq 0}$ as

where the processes $M^K_{1}$ and $M^K_{2}$ are square integrable martingales such that

The proof firstly consists in showing that the sequence of laws of the stochastic processes $(n^K(t), R^K(t), A^K(t), t\geq 0)_{K}$ is relatively compact. It is based on 2-moments estimates, uniform on finite time intervals and on $K$ and on a well known criterion of uniform tightness (cf. for example ^[11]). Then there exists at least one limiting probability measure (on the path space). Using the fact that the jump amplitudes are going to $0$ when $K$ tends to infinity, uniformly on finite time intervals, we deduce that these probability measures only charge continuous trajectories. Moreover, the moment estimates and Eq (A.1) allow to prove that the martingale part converges in probability to $0$ when $K$ tends to infinity. Therefore, it is easy to deduce that the limiting probability measures only charge the solutions of the dynamical system (3.1). The last step consists in proving the uniqueness of such a solution, which is due to a Cauchy-Lipschitz Theorem.

A.2. Numerical simulations

In order to simulate the Markov process $(X^K(t))_{t\ge0} = (n^K(t), R^K(t), A^K(t))_{t\geq 0}$ defined in Eq (2.1) for various sets of parameters, we propose an algorithm simulating numerically the differential system satisfied by the resources in between the jump instants, while the jump instants and the jump amplitudes are simulated directly in terms of the past. The ideas are based on first principles according to the Markov property.

The jump structure of $(X^K(t))_{t\ge0}$ can be described locally at each state $x$ by the value $\alpha(x)\ge0$ of a jump rate function $\alpha$ and if $\alpha (x)>0$ by a probability measure $\pi(x, dh)$ for drawing the amplitudes of the jumps. More precisely, there are overall $p\ge1$ possible non-null jump amplitudes $h_1, \cdots, h_p$, taken at each state $x = (n, r, a)$ at respective rates $\alpha_1(x)\ge0, \cdots, \alpha_p(x)\ge0$, 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\ge0}$ has not jumped yet at time $t>0$ is given in terms of the rate function $\alpha$ by

This allows to construct the first jump instant as follows. If the non-decreasing continuous process $(\Lambda(t))_{t\ge0}$ and its left-continuous inverse $(\Lambda^{-1}(t))_{t\ge0}$ 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 = \Lambda^{-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 $\pi(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\ge0}$ up to the first jump instant. In the present situation this consists in simulating the components $(R^K(t), A^K(t))_{t\geq 0}$ of the Markov process (2.1) by solving the differential system (2.4) in which the other components of (2.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\geq 0}$ we must be able to compute the integral $\Lambda(t)$ and its inverse $\Lambda^{-1}(t)$ defined in Eq (A.2). 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 $\tilde{\alpha}$ such that $\alpha \le \tilde{\alpha}$ and that the corresponding $\tilde{\Lambda}$ and that $\tilde{\Lambda}^{-1}$ defined similarly to Eq (A.2) are simpler to compute than $\Lambda$ and $\Lambda^{-1}$. We simulate the process $(X^K(t))_{t\ge0}$ by an acceptance-rejection method which proposes a jump from state $x$ at rate $\tilde{\alpha}(x)$ and accepts it with probability $\alpha(x)/\tilde{\alpha}(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 $\tilde{\alpha}(x)- \alpha(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 $\tilde{\alpha}$ is a constant. Then the true jump instants of $(X^K(t))_{t\ge0}$ constitute a thinning of a Poisson process of constant intensity $\tilde{\alpha}$, 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 $\alpha(x)/\tilde{\alpha}(x)$.

Let us come back to our model and denote by $\mathcal{L}(X^K(0))$ the distribution of the initial random vector $X^K(0)$, by $\mathcal{E}\left(\lambda\right)$ the exponential law with parameter $\lambda>0$ and by $\mathcal{U}([0, 1])$ the uniform law on~$[0, 1]$. If we moreover denote by $(\phi(x, t-t_0))_{t\geq t_0}$ the flow of the process $(X^K(t))_{t\geq 0}$ from an initial condition $X^K(t_0) = x$ until the next jump time, the above description can be summarized in the following algorithm.

Algorithm:

Simulate $x_0\thicksim \mathcal{L}(X^K(0))$

$T_0 \longleftarrow 0$;

$k \longleftarrow 0$;

Repeat

  Simulate $\epsilon_{k+1}\thicksim\mathcal{E}\left(\tilde{\alpha}\left(x_k \right)\right)$;

  $T_{k+1} \longleftarrow T_{k} + \epsilon_{k+1}$;

  Follow the flow $\left(\phi(x_k, t-T_k)\right)_{t\geq T_k}$ for resources, until the moment $T_{k+1}\wedge T$;

  $x_{k+1} \longleftarrow \phi(x_k, T_{k+1}\wedge T-T_k)$;

  If $T_{k+1} < T$, then

   Simulate $U_{2k}\thicksim\mathcal{U}([0, 1])$;

   If $U_{2k}\tilde{\alpha}(x_{k+1})\leq \alpha(x_{k+1})$, then

    $i \longleftarrow 1$;

    Simulate $U_{2k+1}\thicksim\mathcal{U}([0, 1])$;

    $s \longleftarrow \alpha_1(x_{k+1})$;

    While $i < p$ and $U_{2k+1}\alpha(x_{k+1}) > s$, do

      $i \longleftarrow i + 1$;

      $s \longleftarrow s + \alpha_i(x_{k+1})$;

    End_While.

    $x_{k+1} \longleftarrow x_{k+1} + h_i$;

   End_If.

  End_If.

  $k \longleftarrow k+1$;

Until $T_k\geq T$.