Deterministic and stochastic populationlevel simulations of an artificial lac operon genetic network
 Michail Stamatakis^{1}Email author and
 Kyriacos Zygourakis^{1}
DOI: 10.1186/1471210512301
© Stamatakis and Zygourakis; licensee BioMed Central Ltd. 2011
Received: 10 January 2011
Accepted: 26 July 2011
Published: 26 July 2011
Abstract
Background
The lac operon genetic switch is considered as a paradigm of genetic regulation. This system has a positive feedback loop due to the LacY permease boosting its own production by the facilitated transport of inducer into the cell and the subsequent derepression of the lac operon genes. Previously, we have investigated the effect of stochasticity in an artificial lac operon network at the single cell level by comparing corresponding deterministic and stochastic kinetic models.
Results
This work focuses on the dynamics of cell populations by incorporating the above kinetic scheme into two Monte Carlo (MC) simulation frameworks. The first MC framework assumes stochastic reaction occurrence, accounts for stochastic DNA duplication, division and partitioning and tracks all daughter cells to obtain the statistics of the entire cell population. In order to better understand how stochastic effects shape cell population distributions, we develop a second framework that assumes deterministic reaction dynamics. By comparing the predictions of the two frameworks, we conclude that stochasticity can create or destroy bimodality, and may enhance phenotypic heterogeneity.
Conclusions
Our results show how various sources of stochasticity act in synergy with the positive feedback architecture, thereby shaping the behavior at the cell population level. Further, the insights obtained from the present study allow us to construct simpler and less computationally intensive models that can closely approximate the dynamics of heterogeneous cell populations.
Background
Since the introduction of the operon concept by Jacob et al. [1], the lac operon genetic switch has been considered as a paradigm for genetic regulation. Several experimental studies of this system over the past several decades have elucidated the underlying biomolecular interactions and a plethora of mathematical models have integrated the complex interplays of the key biochemical species in order to predict the behavior of the system [see, for example, [2]].
Most of these models, however, pertain to the single cell behavior [see, for example the models reviewed in [3]] with a limited number of studies focusing on cell populations, or taking comparative approaches. For instance, in a minireview article, Vilar et al. [4] compared different models pertaining to single cells and cell populations, in order to illustrate the performance and limitations of different methodologies. At the single cell level Vilar et al. [4] simulated four phenomenological ordinary differential equations (ODEs) for the concentrations of permease, inducer, and βgalactosidase. A stochastic single cell model was then developed by transforming the phenomenological deterministic rates into propensities. However, the latter transformation is not unique, contrary to the case of a mechanistic model where the reaction rate expressions are derived on the basis of statistical mechanical assumptions [see for instance chapter X in ref. [5], and refs [6–8]]. Thus, the stochastic model simulated by Vilar et al. [4] is based on ad hoc assumptions. In addition, the Vilar et al. models do not appear to account for the dilution of the concentrations due to cellular growth [9], for cell division and for stochastic partitioning or DNA duplication effects, which are important sources of extrinsic noise [10–12].
Furthermore, van Hoek and Hogeweg [13, 14] studied the effect of intrinsic and extrinsic sources of stochasticity, as well as genetic mutations and spatial heterogeneity, at the singlecell and the population level from an evolutionary perspective. Their deterministic model consists of an intracellular and an extracellular part, with the latter capturing the influxes and effluxes of glucose and lactose into and from the cells, as well as the diffusion of these species over a grid, the points of which are assumed to be vacant or occupied by single cells. The intracellular part of the model employs ten differential equations for the concentrations of mRNA, βgalactosidase, permease, lactose, allolactose, glucose, glucose6phosphate, cAMP, ATP, and the size of the cell. The rate expressions used in this model contain Hill functions to model saturation and cooperativity effects, and account for dilution due to cell growth. Hence, stricto sensu the model is not a mechanistic one, and the incorporation of stochastic effects invokes ad hoc assumptions rather than a formal approach based on statistical mechanics. In particular, van Hoek and Hogeweg [13] use the deterministic equation for the mRNA species, but interpret the mRNA concentration as the probability of a single mRNA molecule being present in the cell. They subsequently use this probability to infer the frequency of translational bursts.
Thus, there remain several open questions regarding the emergence of population behavior from the complex interplay between reaction dynamics and stochasticity at the single cell level. Investigating this connection and ultimately understanding cell population dynamics is significant for two reasons. First, typical biology experiments involve cell populations rather than single cells and, thus, phenotypic distributions obtained for instance from flow cytometry pertain to the cell population rather than the lifetime of a single cell. Second, there have been mathematical modeling studies for simple genetic networks, suggesting that the behavior of the single cell is very different from the behavior of the cell population [12, 15, 16]. It would therefore be interesting to investigate whether this is the case in the more complex lac operon system, or whether one can formulate average single cell models that can adequately capture the behavior of the cell population.
In a previously published article, Stamatakis and Mantzaris [17] have presented a kinetic scheme that captures the salient features of an artificial lac operon system, which can be constructed in the lab by introducing mutations to the wildtype system. In the present paper we will incorporate that kinetic scheme into two Monte Carlo (MC) frameworks that simulate cell population dynamics.
The first framework assumes stochastic reaction occurrence and takes into account stochastic DNA duplication, division and partitioning. McAdams and Arkin [18–20] were the first to employ the Gillespie MC algorithm [7, 8] to simulate gene induction and protein synthesis at the single cell level. This approach was later extended to account for cell growth and division [21–23]. However, all these and many other studies simulated the dynamics of biochemical pathways at the single cell level. In a recent study we generalized and expanded the chemical master equation (CME) to the cell population level [10]. The resulting cell population master equation (CPME) is simulated with a Monte Carlo algorithm, and captures stochasticity in the intracellular reaction, stochastic DNA duplication, and division events, as well as in the partitioning of the content of a mother cell to the two daughters. This novel framework is applicable to cell populations since every single daughter cell is tracked, thereby making possible the calculation of any statistical property of the population.
In order to better understand how stochastic effects shape cell population distributions, we develop here a second framework that assumes deterministic reaction dynamics and stochastic DNA duplication, cell division and partitioning. Throughout this work, the single cell models are derived from the reaction network developed in ref. [17] by invoking standard assumptions based on statistical mechanics and reaction rate theory [7, 8, 24] rather than ad hoc techniques [4]. We also explicitly take into account the dilution effect due to cell growth in the deterministic case [9]. Furthermore, stochasticity in the times of division and DNA duplication, as well stochastic partitioning effects are effects that were neglected in previous studies [4], but are explicitly taken into account in our cell population frameworks.
Comparisons of the predictions of the two frameworks reveal that stochasticity can create or destroy bimodality, and may enhance phenotypic heterogeneity by creating heavy tailed distributions, a phenomenon that has not being shown before and can be investigated experimentally. The insights provided from the new study also allow us to construct simpler and less computationally intensive models that can closely approximate the dynamics of heterogeneous cell populations. Specifically, for the case of deterministic reaction dynamics, we use the continuum model formulation which assumes that all cells in the population occupy a continuous and expanding biotic phase [9]. Under this assumption, a set of mass balances for a representative "average" cell in the population provides a lumped description of cell population dynamics. We demonstrate an excellent agreement between the continuum and the cell population models, which is encountered for the first time and contradicts previous studies. Extending these results in the stochastic case, we show that the calculation of distributions of intensive quantities (such as species concentrations) can be performed on the basis of single cell simulations, instead of computationally demanding populationlevel ones.
Methods
We are now ready to present the frameworks that will be used for the simulating cell populations in this study. Both frameworks treat the occurrence of division and DNA duplication events as stochastic processes. Their difference lies in the treatment of reactions and the partitioning of species between the daughter cells: the first framework treats these as stochastic, whereas the second one as deterministic processes. We will refer to the former as the "population model with stochastic reaction dynamics", and the latter as the "population model with deterministic reaction dynamics".
Population Model with Stochastic Reaction Dynamics
In an earlier study [10], we presented the formulation of a cell population master equation (CPME) that describes cell population dynamics and takes into account the major sources of heterogeneity: stochasticity of intracellular reactions, DNAduplication, cell division, and random partitioning of species contents into the two daughter cells. That formulation also takes into account cell growth and respects the discrete nature of the molecular contents and cell numbers. Our approach assumes that each cell of a population can be completely described by a state vector z = (X, V) where X is a vector with n entries for species copy numbers and V is the volume of the cell. Additional morphometric characteristics like cell membrane area or length can be easily incorporated into our framework. We also classify the chemical species into nonchromosomal DNA species and chromosomal DNA species that may exist in various states.
where the terms F_{ R } , F_{ S } , F_{ G } and F_{ D } describe, respectively, the stochastic dynamics of intracellular reactions, DNA duplication, cell growth, and cell division. The form of these terms is presented in Eq. 6 of Additional file 1. Our earlier publication [10] provides the details about the development of the CPME and the Monte Carlo algorithm that simulates the stochastic processes it describes (this information is also included in Section S1 of Additional file 1 for convenience). Here, we just note that the partitioning of cellular volumes upon division is governed by a symmetric beta distribution. Further, the partitioning of nonchromosomal DNA species follows the binomial distribution whereas for chromosomal DNA species partitioning is governed by a multivariate hypergeometric distribution. These choices result in division events where the daughters may have different sizes and species contents; however, on average any daughter inherits half the contents and has half the size of the mother cell. In a more general setting, it is possible to use sums of nonsymmetric beta distributions to describe asymmetric (biased) division (see Section S2 in Additional file 1). In this case, one daughter will consistently inherit more of the mother's molecular content.
Population Model with Deterministic Reaction Dynamics
where the terms , F_{ S } , F_{ G } and describe, respectively, the deterministic dynamics of intracellular reactions and cell growth, and stochastic DNA duplication and cell division. Section S1 of Additional file 1 provides the details for the derivation of Eq. 2, which is essentially the deterministic limit of Eq. 1 for all species being present in large numbers inside the cell. In this limit, the stochastic fluctuations due to the individual reaction events are suppressed. As a consequence, reaction dynamics become deterministic (see Eqs. 11, 14 and pertinent discussion in ref. [6]) and the partitioning of species is now governed by Dirac delta distributions (see Section S3 of Additional file 1).
Structured Continuum Model: a Lumped Description of Cell Population Dynamics
The Artificial lac Operon Reaction Network
To simulate the artificial lac operon genetic network, we will use the reaction network of Stamatakis and Mantzaris [17]. This is a minimal model that neglects effects such as the σ70 dependence of the lac promoter and assumes that only one lacO operator is functional; thus, DNA looping is not accounted for. The existence of three operator sites has been previously shown to result in stronger repression and higher sensitivity in induction [27], as well as lower sensitivity to changes in repressor molecules and lower transcriptional noise [28]. Computational studies of mutations and deletions in these operators have been able to successfully reproduce experimental results [29, 30].
Here, our intention has been to model an artificial lac operon that can be constructed in the lab rather than the natural lac operon system. This artificial system incorporates the positive feedback from the LacY permease resulting in bistable behavior, and thus, allows us to isolate the contributions of stochasticity and bistability in shaping the behavior of the cell population. A more complicated pathway incorporating the aforementioned interactions would be intractable within the cell population frameworks that we have developed, and it is not within the scope of this work to investigate these effects.
Reactions and Propensity Functions for the Stochastic lac Operon Model
Reaction  Propensity Function^{1, 2, 3}  

(11) 


(12) 


(13) 


(14) 


(15) 


(16) 


(17) 


(18) 


(19) 


(110) 


(111) 


(112) 


(113) 


(114) 


(115) 


(116) 


(117) 


(118) 


(119) 


(120) 


(121) 


(122) 


(123) 


(124) 


(125) 


Rate Equations for the Deterministic lac Operon Model
(21) 
 

(22) 
 
(23) 
 
(24) 
 
(25) 
 
(26) 
 
(27) 
 
(28) 
 
(29) 
 
(210) 

Parameters of the lac operon models
Symbol  Value  Units  Description 

R _{ 0,E. coli }  0.4  μm  E. coli radius 
L _{ E. coli }  2.3  μm  Representative E. coli length 
A _{ E. coli }  5.8  μm^{2}  E. coli membrane area (for the above length) 
V _{ E. coli }  1.0  fL  E. coli volume (for the above length) 
O _{ T }  1  (copy number)  operator molecular content 
k _{ sMR }  0.23  nM·min^{1}  lacI transcription rate 
k _{ sR }  15  min^{1}  LacI monomer translation rate constant 
k _{ 2R }  50  nM^{1}·min^{1}  LacI dimerization rate constant 
k _{ 2R }  10^{3}  min^{1}  LacI dimer dissociation rate constant 
k _{ r }  960  nM^{1}·min^{1}  association rate constant for repression 
k _{ r }  2.4  min^{1}  dissociation rate constant for repression 
k _{ dr1 }  3·10^{7}  nM^{2}·min^{1}  association rate constant for 1^{st} derepression mechanism 
k _{ dr1 }  12  min^{1}  dissociation rate constant for 1^{st} derepression mechanism 
k _{ dr2 }  3·10^{7}  nM^{2}·min^{1}  association rate constant for 2^{nd} derepression mechanism 
k _{ dr2 }  4.8·10^{3}  nM^{1}·min^{1}  dissociation rate constant for 2^{nd} derepression mechanism 
k _{ s1MY }  0.5  min^{1}  lacY transcription rate constant 
k _{ s0MY }  0.01  min^{1}  leak lacY transcription rate constant 
k _{ sY }  30  min^{1}  lacY translation rate constant 
k _{ p }  0.12  nM^{1}·min^{1}  LacYinducer association rate constant 
k _{ p }  0.1  min^{1}  LacYinducer dissociation rate constant 
k _{ ft }  6·10^{4}  min^{1}  TMG facilitated transport constant 
h _{ t }  1.55·10^{6}  dm·min^{1}  TMG passive diffusion permeability constant 
λ_{ MR }  0.462  min^{1}  lacI mRNA degradation constant 
λ_{ MY }  0.462  min^{1}  lacY mRNA degradation constant 
λ_{ R }  0.2  min^{1}  LacI monomer degradation constant 
λ_{ R2 }  0.2  min^{1}  LacI dimer degradation constant 
λ_{ Y }  0.2  min^{1}  LacY degradation constant 
λ_{ YIex }  0.2  min^{1}  LacYinducer degradation constant 
λ_{ I2R2 }  0.2  min^{1}  repressorinducer degradation constant 
g  0.0231  (min^{1})  cell growth rate parameter 
n _{ d }  25  (dim/less)  division propensity sharpness exponent 
V _{ d,crit }  15  (fL)  critical volume for division 
q  80  (dim/less)  beta distribution sharpness exponent 
n _{ s }  25  (dim/less)  DNA duplication propensity sharpness exponent 
V _{ s,crit }  10  (fL)  critical volume for DNA duplication 
There is considerable uncertainty for some parameters. For instance, the equilibrium constant for repression, k_{r}/k_{r} is in the range of 10^{13} to 10^{11} M [31–33]. The halflife for dissociation of operator DNA fragments from the repressor has been reported as 30 ~ 49 s [31, 34]; thus, we have taken the dissociation rate k_{r} = 2.4 min^{1} and k_{r} = 960 nM^{1}·min^{1} which results in k_{r}/k_{r} = 2.5·10^{12} M. This value for k_{r} turns out to be much higher than the experimentally measured one of 0.6 nM^{1}·min^{1}, which is deduced by the 59 s time needed for a repressor to find an operator [35], assuming V_{ E.coli } = 10^{15} L. Yet, simulations of the single cell deterministic model with this value of k_{r} produce quantitatively identical bifurcation structures with the nominal parameter set, provided that the value of k_{r} has been adjusted to keep the thermodynamic constant for repression the same. Since the timescale for repression will be different in this case, the single cell stochastic model is expected to exhibit bistability for different induction levels. Previous work [17] showed that slower repressoroperator association and dissociation results in wider bistable regions. In general, however, we expect all qualitative features reported here, to be valid for slower repressoroperator dynamics as well.
Results and discussion
Comparison of Deterministic and Stochastic Reaction Dynamics
Deterministic Reaction Dynamics: Phenotypic Distributions and Statistics
The lac operon system is well known for its ability to exhibit bistable behavior when artificial inducers are used [37]. In our model bistability arises from the positive feedback that LacY exerts on its own expression. Other studies suggest that DNA loop structures also play a key role in the bistable behavior of the lac operon system [refer to [38, 39]]. In this work, however, we focus solely on the aforementioned autocatalytic architecture.
The oscillations shown in Figure 2 are due to the fluctuations of the promoter concentration during growth and DNA duplication. In particular, at the beginning of the life cycle of a newborn cell, a single promoter exists and the cell has a small volume. Thus, the promoter concentration is high and drops as the cell grows. Right after the DNA duplication event, the concentration of the promoter doubles and subsequently drops again as the cell continues to grow up to the point it divides. Upon cell division, the promoter concentration could also change if division is not symmetric, and thus the daughters inherit half the promoter content but they have volumes that are not equal to V_{ mother } /2. Of course, the concentrations of all other species (that are modeled deterministically) remain the same during DNA duplication and cell division (see section S3 in Additional file 1).
For the simulations of Figure 3c and Figure 3d we have biased the process in order to observe such transitions. Still, we must emphasize the fact that we did not impose conditions that are totally unrealistic. Both extremely asymmetric divisions and long delays in division are possible, albeit highly improbable. Thus, they are not observed within the reasonable simulated time intervals of 6000 min. These observations are in line with results obtained with the models of Vilar et al. [4] and van Hoek and Hogeweg [13], which also predicted no switching between the two equilibrium states, namely, induced and uninduced.
Effect of Stochasticity on Cell Population Behavior
Having analyzed the case where reaction dynamics follow deterministic laws, we will now investigate the case of stochastic reaction occurrence since the small copy numbers of molecules encountered in this system are expected to result in significant intrinsic stochasticity.
For low extracellular TMG concentrations, [I_{ ex } ], stochasticity creates a heavy tail in the NDF, whereas the simulation with deterministic reactions exhibits a narrow peak in the total LacY concentration (panel a). This heavy tail can be attributed to the autocatalytic mechanism present in the lac operon system, an assertion that could be used as a starting point for formulating experimentally testable hypotheses. Heavy tails have also been observed in the absence of positive feedback [40]. In that case heavy tails result from the additive noise due to partitioning of species, as well as due to the exponential evolution of protein number during a division cycle, when protein numbers are large, or the fluctuations in protein production rate, for small protein numbers. These effects exist in our model as well; note, however, that in this study we observe heavy tails in the intensive quantities (species concentrations) as opposed to the extensive quantities (number of protein molecules) of ref. [40]. Since the increase of the protein numbers inside the cell is accompanied by an increase in volume, focusing on the concentration would offset the effect of cell growth. Further, heavy tails have been observed for the single cell lac operon model in the absence of growth [17], further supporting our assertion that this behavior is due to the autocatalytic dynamics of the system.
For intermediate inducer concentrations [I_{ ex } ], the deterministic reaction NDF is clearly bimodal, in contrast to that for stochastic reaction, which exhibits a much less prominent higher mode (panel b). Thus, intrinsic noise in this case suppresses bistability; yet, a heavy tail indicative of the positive feedback dynamics is still observed. Note that intrinsic noise does not always result in the suppression of bistability: in a previous study pertaining to the single cell level [17] it was shown that intrinsic stochasticity alone can also result in extending the bistable region to parameter values for which a deterministic model would predict monostable behaviour. Finally, for high [I_{ ex } ] the NDF for deterministic reaction is unimodal (panel c). However, the NDF for stochastic reaction appears bimodal, with a sharp peak at total LacY concentration equal to zero and a wide peak at high [Y]_{T}. In this case, stochasticity in reaction occurrence appears to extend the region where the NDF is bimodal and widen the upper mode of the distribution.
In all cases (Figure 5), the heterogeneity exhibited by the population model with stochastic reaction occurrence is significantly larger than that exhibited by the model with deterministic reaction occurrence. This observation agrees with previous results by van Hoek and Hogeweg [13], whose model also predicted that, in the absence of spatial and genetic heterogeneity, the population variability exhibited by the stochastic simulation is much larger than that of the deterministic one.
Structured Continuum Model for Deterministic Reaction Dynamics
Simulating in detail the intra and intercellular processes at the population level is computationally demanding. Therefore, it is natural to ask whether one can adequately predict the dynamics of the population average with the use of a lumped model that neglects heterogeneity. For this purpose, we will use the structured continuum model formulation [9] that treats all cells as a lumped and expanding (due to cell growth) biotic phase. Section S4 in Additional file 1 presents all the transient mass balances that constitute the structured continuum model defined by Eqs. 3 and 4. In our case, μ is equal to g (Eq. 4) because we have assumed exponential growth, which implies a constant average specific growth rate. Note again that these balances are written for the average intracellular concentrations (intensive quantities) of the chemical species of interest over the entire cell population. The structured continuum model is valid for all times. The assumption here is that the concentrations in any cell at any time remain close to the concentrations predicted by this continuum model. Note that we did not include the dilution effect for the operator species O, since the DNA duplication process continuously regenerates O.
Simulations with the structured continuum model require estimates for the average membrane area over the volume 〈A/V〉 and the average total operator concentration 〈[O]_{ T }〉. These average quantities can either be obtained directly from the population simulation or estimated if we know approximately the sizes at which the cells duplicate their DNA or divide.
Using the parameters of Table 3, this estimated average equals 5.69 μm^{1}. Cell population simulations give an average equal to 5.76 μm^{1}. Thus, using these heuristic arguments, we estimated the average ratio 〈A/V〉 with a remarkably low error (1.2%).
Using the parameters of Table 3, we obtain an estimate of the average operator concentration equal to 2.40 nM. Cell population simulations give an average equal to 2.48 nM (or 3.2% error).
The qualitative agreement between the steady states of the structured continuum model and the population averages is excellent. The agreement between the two models is remarkably good when the populationaverage values for A/V and [O] _{ T } are used in the structured continuum model. When we use the estimated values from Eqs. 8 and 9, the total LacY concentration at maximal induction is slightly underestimated due to the underestimation of the average total operator concentration. Thus, the level of agreement between the steady states of the structured continuum model and the averages of the population model depends on how good the estimates for the average cell characteristics are.
The above results pertain to the time invariant behaviors of the cell population. It is interesting, however, to also compare the dynamical behavior of the structured continuum and the cell population model.
The dashed and dashdotted curves on panels (a) and (b) of Figure 7 present predictions of the structured continuum model. For the dashdotted curves, the average values for the cellular characteristics 〈A/V 〉 and 〈[O]_{ T }〉 were taken from the cell population simulations, while for the dashed curves these values were estimated using Eqs. 8 and 9. For panel (a), the structured continuum models are simulated for [I_{ ex } ] = 60 μM, using as initial conditions the steady state concentrations calculated for [I_{ ex } ] = 0 μM. The dynamics of the opposite switch are presented in panel (b).
The agreement between the transient behaviors predicted by the structured continuum model and the population model is excellent even in the case where the intracellular dynamics are significantly slower than the proliferation rates of the cell as shown in panels (c) and (d) of Figure 7. Panel (c) shows that for slow lacY mRNA dynamics the switching from the low to the high state takes approximately 300 min, which is much longer than the 30 min average doubling time (panel d). Even so, the structured continuum model yields excellent predictions of the population average. Similar observations were previously reported by Vilar et al. [4], who showed that an approach taking into account single cell stochasticity and population level heterogeneity may not be needed, depending on the system studied and the conditions of interest. In particular, Vilar et al. [4] noted that, when the lac operon is under the influence of high inducer concentrations, averages taken over independent cells versus the overall population give very similar results.
This agreement between the steady state and the transient behavior of the structured continuum and the cell population models can be explained as follows. In this system, the only coupling between the cells is due to stochastic partitioning that can generate variability but no bias on the total LacY concentration, in the sense that one does not observe consistently higher LacY_{T} concentrations in one daughter versus the other. In fact, the two newborn daughters may have different LacY_{T} contents, but they always have equal LacY_{T} concentrations, which are identical to their mother's LacY_{T} concentration just before division. Similarly, unsynchronized DNA duplication and division events also contribute to the observed variability in LacY_{T}, but they cannot consistently bias the LacY concentration.
In general, the key to understanding this effect lies in the fact that in the cell population, division results in the removal of an "old" mother cell and the addition of two "young" daughter cells. If the properties of the "old" cells are different than those of the "young" daughters, then this effect results in the properties of the "young" ones being overrepresented in the population. On the other hand, such an effect is absent in a cell chain, where each division results in the removal of a mother cell and the addition of a single daughter, and thus, the properties of both subpopulations are equally represented in the cell chain probability distribution.
Let us now suppose that the property we chose to investigate is a species concentration (intensive variable). The concentrations of mother and daughter cells are on average the same, due to the symmetric properties of the binomial and beta distributions that model division. Consequently, the cell chain and cell population distributions will be practically indistinguishable. On the other hand, if we chose an intensive property to investigate, the two distributions would be different in general, as we have shown in our previous work [Figure 6 in ref. [10]]. This point is rigorously proven for deterministic reaction dynamics in an earlier publication [41], where we establish that a population balance equation transforms into a continuum model for the species concentrations if certain conditions are satisfied. The conditions, in summary, call for equal species concentrations in mother and daughter cells, sizeindependent intensive reaction rates, sizedependent but concentrationindependent growth and division rates and partitioning probability density function. All conditions hold true to a good approximation in our present simulations and they are also expected to be biologically plausible for a range of systems of interest.
These observations contradict previously published results by Mantzaris for different biological systems [12, 16, 42, 43]. In these studies, comparisons of the steady states of structured continuum models with the stationary averages of cell population balances (CPBs) showed that the regions in which the two models exhibit specific types of behavior, such as bistability, were vastly different. These differences were interpreted as effects of cell population heterogeneity. However, we believe that the differences shown by Mantzaris [12, 16, 42, 43] are not genuine effects of heterogeneity, but rather stem from the modeling assumptions employed.
More specifically, these studies incorporated into the CPB single cell models that are written for species concentrations. However, the CPB accepts single cell models written for cellular contents (amounts), which are extensive quantities, and not concentrations which are intensive [9]. One immediate consequence of this fact is by partitioning the concentrations upon cell division one would violate mass conservation. Moreover, the cell volume is not taken into account in any of these studies. Cells in the exponential phase, however, are continuously growing and dividing, thereby changing their volume roughly twofold during one division cycle. In turn, the change in volume shifts the thermodynamic equilibrium point and also affects the dynamics if multimolecular reactions are present in the reaction scheme, which is the case in all of the systems considered in these studies [12, 16, 42, 43]. Therefore, neglecting the change of the cell volume as time progresses changes creates a priori an inconsistency between the single cell expressions incorporated into the population balance and the structured continuum model.
Furthermore, some of these studies used unequal partitioning to artificially generate complex behavior, such as oscillations in a reaction network with 0^{th} and 1^{st} order reactions. However, the asymmetry in E. coli division is negligible: cells may stochastically divide in two unequal daughters, but consistent generation of one large and one small daughter has not been experimentally observed. In fact, it has been shown that the distribution of daughter cell sizes has a mean corresponding to equal partitioning and a small coefficient of variation [44, 45].
Simulation of Cell Chains for Stochastic Reaction Dynamics
Given the high computational expense of simulating the population model with stochastic reaction dynamics, we pose the question of whether there is a simpler method for obtaining good approximations for the distributions of phenotypic characteristics. In contrast to the case of deterministic reaction dynamics, we cannot use a continuum model for comparison purposes here. The main reason is that in chemical systems far from the thermodynamic limit, the size influences noise strength. In the case of stochastic reaction occurrence, therefore, growth results in dilution of molar contents and also suppression of stochastic fluctuations. Consequently, we cannot create a lumped model that accounts for growth with just the incorporation of a dilution term as was done in the deterministic case.
However, we can take a different approach. We have already observed that cell population dynamics emerge from single cell behavior and that, on average, the two daughters share the same concentrations as their mother cell. Therefore, instead of tracking the NDF in a cell population, we could compute the PDF in a cell chain. Simulation of a cell chain tracks only one daughter after each division event. Thus, instead of focusing on the expected number of cells of the population that exist in state z, we turn our attention to the probability of finding a single cell of the cell chain at state z. Note, however, that we will be comparing the PDF and NDF of intensive quantities, in particular the total LacY concentration.
Conclusions
This study generalized the deterministic and stochastic single cell lac operon models of an earlier publication [17] using two cell population frameworks that account for stochastic and deterministic dynamics of the biochemical reactions. We subsequently performed simulations to investigate the behavior of the cell population and demonstrated the effect of stochasticity in reaction dynamics. We concluded that this source of stochasticity can amplify population heterogeneity by introducing heavy tails to the phenotypic distributions and can create or destroy bimodality.
We also carried out a systematic comparison of predictions obtained by a structured continuum model and a detailed cell population model with deterministic reaction dynamics. These comparisons showed that the structured continuum model gives satisfactory results for the average LacY_{T} concentration of the population, even in the case where the reaction dynamics are much slower than the proliferation dynamics. This agreement between the two models was attributed to the similar intensive properties (such as species concentrations) between mother and daughter cells.
Finally and in the case of stochastic reaction dynamics, we demonstrated that by simulating the dynamics of a cell chain we can obtain very accurate approximations of the cell population dynamics. The PDFs obtained by cell chain simulations for the LacY_{T} concentration were in excellent agreement to the cell population NDF computed with the Monte Carlo algorithm implementing the CPME of Eq. 1.
Our study shows that for cell populations in which the cells interact weakly, through division only, it is possible to accurately model and explain the population behavior in terms of the single cell dynamics. For such systems, the key parameters for describing the behavior of the population are the kinetic constants of the underlying pathway of interest, and the physiological functions that express the single cell growth rate, DNAduplication and division propensity, as well as the partitioning mechanism. Intrinsic noise can be inherently accounted for, once the cell size has been specified and thus no additional parameters are needed for this purpose. Such a description is expected to be of great importance in bioinformatics studies focusing on population variability, since, for cells interacting though division only this variability can be explained in terms of a limited number of parameters.
Finally, deviations between the experimentally observed population dynamics and the behavior predicted by our framework may indicate the presence of more complex effects that are not accounted for in this framework. For instance, the cells could be nonisogenic, or coupled with strong interaction mechanisms. Another source of complexity is the existence of multiple compartments in the cell, which would invalidate the assumption of a single wellmixed intracellular space. Such effects would have to be incorporated to the framework in order to obtain a more accurate description of the system of interest.
Author's Contributions
MS developed the mathematical models and performed the simulations, analyzed the results and drafted the manuscript. KZ conceived of the study, participated in its design and coordination, and refined the manuscript. Both authors read and approved the final manuscript.
Declarations
Acknowledgements
The financial support of NIHNIGMS through grant R01 GM071888 is gratefully acknowledged. This work was also supported in part by the Shared University Grid at Rice funded by NSF under Grant EIA0216467, and a partnership between Rice University, Sun Microsystems, and Sigma Solutions, Inc.
Authors’ Affiliations
References
 Jacob F, Perrin D, Sanchez C, Monod J: L'operon: groupe de gènes à expression coordonnée par un opérateur. Comptes Rendus Hebdomadaires des Séances de L'Académie des Sciences Serie D: Sciences naturelles 1960, 250: 1727–1729.
 Santillán M, Mackey MC: Quantitative approaches to the study of bistability in the lac operon of Escherichia coli. J R Soc Interface 2008, 5(Suppl 1):S29S39.PubMed CentralView ArticlePubMed
 Stamatakis M: Stochasticity and Cell Population Heterogeneity in an Artificial lac Operon Genetic Network. PhD thesis. Houston: Rice University 2009.
 Vilar JMG, Guet CC, Leibler S: Modeling network dynamics: the lac operon, a case study. J Cell Biol 2003, 161(3):471–476. 10.1083/jcb.200301125PubMed CentralView ArticlePubMed
 van Kampen NG: Stochastic processes in physics and chemistry. New York, Amsterdam: NorthHollandPersonalLibrary; 1992.
 Gillespie DT: The chemical Langevin equation. J Chem Phys 2000, 113(1):297–306. 10.1063/1.481811View Article
 Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys 1976, 22(4):403–434. 10.1016/00219991(76)900413View Article
 Gillespie DT: Exact stochastic simulation of coupled chemical reactions. J Phys Chem 1977, 81(25):2340–2361. 10.1021/j100540a008View Article
 Fredrickson AG: Formulation of Structured Growth Models. Reprinted from Biotechnology and Bioengineering, Vol. XVIII, No. 10, Pages 1481–1486 (1976). Biotechnol Bioeng 2000, 67(6):720–725. 10.1002/(SICI)10970290(20000320)67:6<720::AIDBIT10>;3.0.CO;2IView ArticlePubMed
 Stamatakis M, Zygourakis K: A Mathematical and Computational Approach for Integrating the Major Sources of Cell Population Heterogeneity. J Theor Biol 2010, 266(1):41–61. 10.1016/j.jtbi.2010.06.002PubMed CentralView ArticlePubMed
 Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic Gene Expression in a Single Cell. Science 2002, 297(5584):1183–1186. 10.1126/science.1070919View ArticlePubMed
 Mantzaris NV: From SingleCell Genetic Architecture to Cell Population Dynamics: Quantitatively Decomposing the Effects of Different Population Heterogeneity Sources for a Genetic Network with Positive Feedback Architecture. Biophys J 2007, 92: 4271–4288. 10.1529/biophysj.106.100271PubMed CentralView ArticlePubMed
 van Hoek M, Hogeweg P: The effect of stochasticity on the lac operon: An evolutionary perspective. PLoS Comput Biol 2007, 3(6):e111. 10.1371/journal.pcbi.0030111PubMed CentralView ArticlePubMed
 van Hoek MJA, Hogeweg P: In silico evolved lac operons exhibit bistability for artificial inducers, but not for lactose. Biophys J 2006, 91(8):2833–2843. 10.1529/biophysj.105.077420PubMed CentralView ArticlePubMed
 Mantzaris NV: Singlecell geneswitching networks and heterogeneous cell population phenotypes. Comput Chem Eng 2005, 29(3):631–643. 10.1016/j.compchemeng.2004.08.009View Article
 Mantzaris NV: Stochastic and deterministic simulations of heterogeneous cell population dynamics. J Theor Biol 2006, 241(3):690–706. 10.1016/j.jtbi.2006.01.005View ArticlePubMed
 Stamatakis M, Mantzaris NV: Comparison of Deterministic and Stochastic Models of the lac Operon Genetic Network. Biophys J 2009, 96(3):887–906. 10.1016/j.bpj.2008.10.028PubMed CentralView ArticlePubMed
 McAdams HH, Arkin A: Stochastic mechanisms in gene expression. Proc Nat Acad Sci USA 1997, 94(3):814–819. 10.1073/pnas.94.3.814PubMed CentralView ArticlePubMed
 McAdams HH, Arkin A: Simulation of prokaryotic genetic circuits. Annu Rev Biophys Biomol Struct 1998, 27: 199–224. 10.1146/annurev.biophys.27.1.199View ArticlePubMed
 McAdams HH, Arkin A: It's a noisy business! Genetic regulation at the nanomolar scale. Trends Genet 1999, 15(2):65–69. 10.1016/S01689525(98)01659XView ArticlePubMed
 Gibson MA, Bruck J: Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. J Phys Chem A 2000, 104(9):1876–1889. 10.1021/jp993732qView Article
 Lu T, Volfson D, Tsimring L, Hasty J: Cellular growth and division in the Gillespie algorithm. Syst Biol 2004, 1(1):121–128. 10.1049/sb:20045016View Article
 Swain PS, Elowitz MB, Siggia ED: Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Nat Acad Sci USA 2002, 99(20):12795–12800. 10.1073/pnas.162041399PubMed CentralView ArticlePubMed
 McQuarrie DA: Stochastic Approach to Chemical Kinetics. J Appl Probab 1967, 4(3):413–478. 10.2307/3212214View Article
 Ramkrishna D: Population Balances: Theory and Applications to Particulate Systems in Engineering. San Diego, CA: Academic Press; 2000.
 Crudu A, Debussche A, Radulescu O: Hybrid stochastic simplifications for multiscale gene networks. Bmc Syst Biol 2009, 3: 89. 10.1186/17520509389PubMed CentralView ArticlePubMed
 Saiz L, Vilar JMG: Proteinprotein/DNA interaction networks: versatile macromolecular structures for the control of gene expression. IET Syst Biol 2008, 2(5):247–255. 10.1049/ietsyb:20080091View ArticlePubMed
 Vilar JM, Leibler S: DNA looping and physical constraints on transcription regulation. J Mol Biol 2003, 331(5):981–989. 10.1016/S00222836(03)007642View ArticlePubMed
 Saiz L, Vilar JMG: Ab initio thermodynamic modeling of distal multisite transcription regulation. Nucleic Acids Res 2008, 36(3):726–731.PubMed CentralView ArticlePubMed
 Vilar JMG: Accurate Prediction of Gene Expression by Integration of DNA Sequence Statistics with Detailed Modeling of Transcription Regulation. Biophys J 2010, 99(8):2408–2413. 10.1016/j.bpj.2010.08.006PubMed CentralView ArticlePubMed
 Goeddel DV, Yansura DG, Caruthers MH: Binding of Synthetic Lactose Operator DNAs to Lactose Repressors. Proc Nat Acad Sci USA 1977, 74(8):3292–3296. 10.1073/pnas.74.8.3292PubMed CentralView ArticlePubMed
 Riggs AD, Newby RF, Bourgeois S: lac RepressorOperator interaction II. Effect of galactosides and other ligands. J Mol Biol 1970, 51(2):303–314. 10.1016/00222836(70)901440View ArticlePubMed
 Riggs AD, Newby RF, Bourgeois S: lac RepressorOperator interaction I. Equilibrium studies. J Mol Biol 1970, 48(1):67–83. 10.1016/00222836(70)902196View ArticlePubMed
 Lin SY, Itakura K, Rosenberg JM, Wilcox G, Bahl C, Wu T, Narang SA, Dickerson RE, Riggs AD: Molecular Mechanisms in the Control of Gene Expression. Edited by: Nielich DP, Rutter WJ, Fox CF. New York: Academic Press; 1976:143–158.View Article
 Elf J, Li GW, Xie XS: Probing transcription factor dynamics at the singlemolecule level in a living cell. Science 2007, 316(5828):1191–1194. 10.1126/science.1141967PubMed CentralView ArticlePubMed
 Trueba FJ, Woldringh CL: Changes in cell diameter during the division cycle of Escherichia coli. J Bacteriol 1980, 142(3):869–878.PubMed CentralPubMed
 Dreisigmeyer DW, Stajic J, Nemenman I, Hlavacek WS, Wall ME: Determinants of bistability in induction of the Escherichia coli lac operon. IET Syst Biol 2008, 2(5):293–303. 10.1049/ietsyb:20080095View ArticlePubMed
 Oehler S, Alberti S, MüllerHill B: Induction of the lac promoter in the absence of DNA loops and the stoichiometry of induction. Nucleic Acids Res 2006, 34(2):606–612. 10.1093/nar/gkj453PubMed CentralView ArticlePubMed
 Narang A: Effect of DNA looping on the induction kinetics of the lac operon. J Theor Biol 2007, 247(4):695–712. 10.1016/j.jtbi.2007.03.030View ArticlePubMed
 Krishna S, Banerjee B, Ramakrishnan TV, Shivashankar GV: Stochastic simulations of the origins and implications of longtailed distributions in gene expression. Proc Nat Acad Sci USA 2005, 102(13):4771–4776. 10.1073/pnas.0406415102PubMed CentralView ArticlePubMed
 Stamatakis M: Cell population balance, ensemble and continuum modeling frameworks: Conditional equivalence and hybrid approaches. Chem Eng Sci 2010, 65(2):1008–1015. 10.1016/j.ces.2009.09.054View Article
 Mantzaris NV: Effects of Population Heterogeneity on the Dynamics of Cell Populations. In 7th IFAC Symposium on Dynamics and Control of Process Systems. Cambridge, MA; 2004.
 Mantzaris NV: A cell population balance model describing positive feedback loop expression dynamics. Comput Chem Eng 2005, 29(4):897–909. 10.1016/j.compchemeng.2004.09.012View Article
 Harvey RJ, Marr AG, Painter PR: Kinetics of Growth of Individual Cells of Escherichia coli and Azotobacter agilis. J Bacteriol 1967, 93(2):605–617.PubMed CentralPubMed
 Marr AG, Harvey RJ, Trentini WC: Growth and division of Escherichia coli. J Bacteriol 1966, 91(6):2388–2389.PubMed CentralPubMed
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.