Coarse-grained analysis of stochastically simulated cell populations with a positive feedback genetic network architecture
aa r X i v : . [ q - b i o . M N ] D ec Journal of Mathematical Biology manuscript No. (will be inserted by the editor)
Coarse-grained analysis of stochastically simulated cellpopulations with a positive feedback genetic networkarchitecture.
I. G. Aviziotis · M. E. Kavousanakis · I.A. Bitsanis · A. G. Boudouvis the date of receipt and acceptance should be inserted later
Abstract
Among the different computational approaches modelling the dy-namics of isogenic cell populations, discrete stochastic models can describewith sufficient accuracy the evolution of small size populations. However, for asystematic and efficient study of their long-time behaviour over a wide rangeof parameter values, the performance of solely direct temporal simulations re-quires significantly high computational time. In addition, when the dynamicsof the cell populations exhibit non-trivial bistable behaviour, such an analysisbecomes a prohibitive task, since a large ensemble of initial states need tobe tested for the quest of possibly co-existing steady state solutions. In thiswork, we study cell populations which carry the lac operon network exhibit-ing solution multiplicity over a wide range of extracellular conditions (inducerconcentration). By adopting ideas from the so-called “equation-free” method-ology, we perform systems-level analysis, which includes numerical tasks suchas the computation of coarse steady state solutions, coarse bifurcation analysis,as well as coarse stability analysis. Dynamically stable and unstable macro-scopic (population level) steady state solutions are computed by means ofbifurcation analysis utilising short bursts of fine-scale simulations, and the
I.G. AviziotisNational Technical University of Athens, School of Chemical Engineering, Athens, 15780,Greece E-mail: [email protected]. KavousanakisNational Technical University of Athens, School of Chemical Engineering, Athens, 15780,GreeceTel.: +30-772-3290E-mail: [email protected]. BitsanisInstitute of Electronic Structure and Laser, Foundation for Research and Technology Hellas,GR-711 10 Heraklion, Greece E-mail: [email protected]. BoudouvisNational Technical University of Athens, School of Chemical Engineering, Athens, 15780,Greece E-mail: [email protected] I. G. Aviziotis et al. range of bistability is determined for different sizes of cell populations. Theresults are compared with the deterministic cell population balance (CPB)model, which is valid for large populations, and we demonstrate the increasedeffect of stochasticity in small size populations with asymmetric partitioningmechanisms.
Keywords lac operon · Monte Carlo simulations · equation-free method · bistability · cell heterogeneity Mathematics Subject Classification (2000) · · Molecular biology, genomics, transcriptomics and proteomics have providedus the appropriate powerful tools for the investigation and understanding ofthe immensely complex processes, which occur at the single-cell level. How-ever, the biological behavior is not a matter of solely intracellular processes;it also depends on the interactions, which take place between the cells ofan isogenic population, leading to a phenotypic variability amongst them.This phenomenon is known as cell population heterogeneity and has been ob-served in an abundance of biological systems, e.g., phage burst size variations(Delbr¨uck 1945), transcriptional states heterogeneity in sporulating cultures of
Bacillus subtilis (Chung and Stephanopoulos 1995), heterogeneity in endothe-lial cell surface markers (Oh et al. 2004) and in various isogenic
Escherichiacoli systems (Elowitz et al. 2002).Until the previous decade, the biological paradigms and modeling frame-works for the design and control of biochemical processes did not consider cellpopulation heterogeneity (Chung and Stephanopoulos 1995; Fedorrof and Fontana 2002).Their key assumption was that all cells of an isogenic population behavelike the average cell, and simple ordinary differential equations were usedto describe the population behavior (Avery 2006; Davidson and Srette 2008).However, such an approach can result in incorrect predictions as shown in(McAdams and Arkin 1998; Mantzaris 2005; Kavousanakis et al. 2009) and onemust account for the heterogeneous nature of cell populations. Furthermore,much of our knowledge is based on ensemble measurements, and despite thefact that cell-to-cell differences are always present to some degree, the collectivebehaviour of a population may not represent the behaviour of the individuals(Altschuler and Wu 2010).The accurate prediction of the collective behaviour of cell populations isof great importance in the majority of biotechnological applications, wherethe objective is to maximise their productivity, rather than increasing the effi-ciency of each individual cell. Furthermore, phenotypic heterogeneity may belinked to the viability of a cell population and its ability to adapt to abruptchanges of their environment (McAdams and Arkin 1999; Sumner and Avery 2002;Veening et al. 2008). Thus, the resistance of certain infectious bacteria to an-tibiotics could be explained on the basis of the existence of small sub-populations oarse-grained analysis of stochastically simulated cell populations 3 that survive the medical treatment and resume growing after the antibiotic hasbeen removed (Booth 2002).From the discussion above, it naturally emerges that a mathematical de-scription that incorporates heterogeneity in cell population dynamics is ofgreat importance for a variety of biological systems. Fredrickson and coworkersintroduced a special class of models in the 1960s for the prediction of the behav-ior of heterogeneous cell populations (Eakman et al. 1966; Tsuchiya et al. 1966;Fredrickson et al. 1967), known as Cell Population Balance (CPB) models.These models are partial integro-differential equations, which describe thephysiological state of cells, i.e., a vector whose components can be the in-tracellular content of chemical species, as well as morphometric characteristicsof the cells, e.g., size. CPB models incorporate the type of heterogeneity orig-inating from the unequal partitioning of most intracellular components (withthe exception of DNA) between the daughter cells (Block et al. 1990). Due tothe operation of the cell cycle, this phenomenon repeats itself, thus leading tofurther variability. However, the validity of CPB models is limited only whenthe continuum assumption is justified, i.e., for large size populations. Further-more, the deterministic approach does not account for stochastic effects duringdivision, which can play a key role in small size cell populations. For these rea-sons, and in order to provide an -as much as possible- accurate description ofthe system behaviour, we resort to fine-scale, stochastic modelling approaches.The first modelling approach for capturing the stochastic behavior of cellpopulations can be found in (Shah et al. 1976), who developed a Monte Carloalgorithm to describe the dynamics of the cell mass distribution. In Hatzis etal. (1995), the Monte Carlo algorithm was extended to describe a system ofincreased complexity describing the growth of phagotrophic protozoa. Thesealgorithms are computationally intensive due to the exponential growth intime of the number of cells in the population. Constant-number Monte Carlo(CNMC) algorithms (Smith and Matsoukas 1998; Mantzaris 2006) overcomethe extensive requirements in CPU time, by simulating a constant number ofcells that are assumed to be a representative sample of the overall population.Despite the fact that stochastic modelling provides a more realistic descrip-tion of the physical problem, their computational efficiency is limited when acoarse-level analysis is required. For example, if we are interested in studyingthe asymptotic behaviour of a cell population with respect to certain parame-ters (e.g., extracellular inducer concentration), then a large number of dynamic-long time- stochastic simulations is required. Furthermore, for cases where thesystem is expected to feature solution multiplicity within a range of parametervalues, the proper initial states need to be chosen in order to converge to allpossible steady state solutions, with no guarantee of success since the limitsof bistability are not known a priori .For a systematic and efficient analysis of the cell population coarse be-haviour, we adopt ideas from the so-called “equation-free” methodology (Gear et al. 2002;Kevrekidis et al. 2003; Kevrekidis et al. 2004), a computer-assisted multiscaleframework that enables models at a fine (microscopic) level of description toperform numerical tasks at a coarse (macroscopic) level. The backbone of the
I. G. Aviziotis et al. methodology is the development of a computational super-structure, whichwraps around the microscopic simulator and reports the evolution of coarsequantities of interest at discrete time instances. This discrete time-mapping,which is called the coarse time-stepper, allows the communication betweenthe single-cell and population level. Using as a tool the coarse time-stepper,we implement well established numerical techniques for the performance ofsystems-level analysis tasks, such as the computation of coarse steady-statesolutions, coarse bifurcation analysis and coarse stability analysis.Here, we demonstrate the efficiency of this multiscale methodology, whenthe fine-scale, microscopic simulator is a CNMC model (Mantzaris 2006), whichsimulates the evolution of an isogenic cell population carrying the lac operongenetic network. Within a range of extracellular inducer concentration val-ues (IPTG), the lac operon exhibits bistability, i.e., for the same parametervalue asymptotic cell phenotypes of high or low expression levels of the lac
Ygene can be both observed. At the single-cell level, the solution multiplicity isattributed to the nonlinear expression of the rate of reaction in which the intra-cellular content ( lac
Y) participates (Mantzaris 2007). The same behaviour isalso inherited to the population level, when solving a deterministic CPB model,which incorporates cell heterogeneity(Kavousanakis et al. 2009). In fact, cellheterogeneity has a significant impact on the average population phenotypeshifting the bistability region towards higher IPTG concentration values. Theaccurate determination of the bistability region limits is realised by applyingbifurcation analysis, and in particular by means of the pseudo arc-length pa-rameter continuation technique (Keller 1977). In this work, we perform thesame parametric analysis using a stochastic model, in order to study the effectof random events during division on the limits of the bistability range.The performance of such systems-level analysis is feasible by means ofthe equation-free methodology. In particular, we perform coarse steady statecomputations and coarse bifurcation analysis, utilising CNMC simulations(Mantzaris 2006). This study enables the determination of the bistability rangefor different population sizes, which is compared with results obtained fromthe numerical solution of deterministic CPB model. Major differences are ob-served for small size cell populations, with highly asymmetric partitioningmechanism; the computed bistability region for small size populations is lo-cated within the range of higher IPTG concentration values compared to thedeterministic analysis findings.The accurate determination of the bistability region becomes a valuabletool for the understanding of possible phenotypic switching induced by randomevents (Libby and Rainey 2011). When the cell population operates at theproximity of the bistability interval limits, small fluctuations can drive thesystem to large phenotypic changes. Here, we illustrate switching from high tolow (and reverse) levels of lac
Y gene expressions, when the parameter value(here the IPTG concentration) is in the vicinity of critical turning points,which are computed from the coarse bifurcation analysis.This paper is organised as follows: In Section 2, we present the CNMCmodel, which simulates the evolution of an isogenic cell population. We briefly oarse-grained analysis of stochastically simulated cell populations 5 describe the lac operon network and present results of direct stochastic sim-ulations for different IPTG concentrations. The equation-free methodology isdescribed in Section 4; in Section 5, we present results of the proposed coarse-grained analysis, studying the effect of partitioning asymmetry, and of thedivision rate sharpness for different sizes of cell populations. Finally, we sum-marise the main findings of this study in Section 6 and propose future researchdirections. lac operon), and we denote with S τ the randomstate of the cell population at a given dimensionless time τ : S τ ≡ { X i ( τ ) = x i , i = 1 , , ..., N } , (2.1)where X i ( τ ) is the intracellular content of cell i , and N is the constant numberof cells of the simulated population. Each cell undergoes division at a rate Γ ( x )and a partition probability density function P ( x, x ′ ) determines the intracellu-lar content of newborn cells. Between division events, the content of each cellchanges according to an ordinary differential equation, the exact formulationof which depends on the genetic network the cells carry:d x i d τ = R ( x i ) , i = 1 , , ..., N. (2.2) Computation of time between division events
The time between division events, T , is a random variable which depends on S τ ; its cumulative distribution function is given from (Mantzaris 2006): F τ ( z | τ ) = 1 − exp " − Z z N X i =1 Γ ( x i ( τ + z ′ ))d z ′ . (2.3)In order to determine the time between two successive division events, a ran-dom number p is then generated from a uniform distribution, with p ∈ [0 , P Ni =1 R T Γ ( x i ( τ + z ′ )) d z ′ ln [1 − p ] + 1 = 0 . (2.4)Equation (2.4) is solved iteratively with the Newton-Raphson method; thetrapezoid rule is applied for the numerical computation of the integral. The I. G. Aviziotis et al. evolution of cell’s i intracellular content, x i , during this intervening time periodis computed by integrating in time Eq. (2.2) with the explicit forward Eulerscheme. Alternative time integration schemes can be also applied, howeverwith no significant effect on the results of simulations. Determination of the cell undergoing division
After the computation of time T between division evens, the state vector S τ is updated by integrating Eq. (2.2) (with the forward Euler scheme) and thecell undergoing division is selected. The dividing cell (denoted with index k )is randomly selected from the conditional distribution function:Pr { k = j | S τ + T } = Γ ( x j ( τ + T )) P Ni =1 Γ ( x i ( τ + T )) . (2.5) Determination of the content of the two daughter cells
Upon division of cell k , the content of mother cell is partitioned to the new-born cells according to a partition probability density function P ( x, x ′ ), whichdetermines the mechanism through which a mother cell of content x ′ producesa daughter cell with content x and a second cell with content x ′ − x . In the re-sults presented below, we choose to work with the simple discrete partitioningmechanism, formulated as: P ( x, x ′ ) = 12 f δ ( f x ′ − x ) + 12(1 − f ) δ ((1 − f ) x ′ − x ) , (2.6)where δ is the Dirac function and f is an asymmetry parameter. In practice, thedivision of the mother cell with intracellular content x ′ results in two daughtercells with content f x ′ and (1 − f ) x ′ , respectively. From this definition, it isclear that f ∈ [0 , . f corresponding to more asymmetricintracellular content partitioning. When f = 0 . Restoration of the sample
Each division event leads to the generation of two newborn cells from themother cell. Upon determination of the content of each of the daughter cells,the mother cell is replaced by the first daughter cell. In order to maintain thesize of population constant at a pre-specified value, N , the second daughtercell replaces a randomly chosen cell from the population, with all cells havingequal probability of being selected and replaced. Finally, the dimensionlesstime is updated, τ → τ + T , and the algorithm is repeated until a pre-specifiedstoppage time is reached. oarse-grained analysis of stochastically simulated cell populations 7 ∂n ( x, t ) ∂t + ∂∂x [ R ( x ) n ( x, t )] + Γ ( x ) n ( x, t ) == 2 Z x max x Γ ( x ′ ) P ( x, x ′ ) n ( x ′ , t )d x ′ − n ( x, t ) Z x max Γ ( x ) n ( x, t )( d ) x, (2.7)where n ( x, t ) is the cell density function and denotes the number of cells withcontent x at time t , divided by the total number of cells at the same time.The maximum intracellular content value is denoted with x max . The boundaryconditions imposed to (2.7) require that the cells of the population do not growoutside the domain [0 , x max ]: n (0 , t ) = n ( x max , t ) = 0 . (2.8)If we take the first-order moment of (2.7) and apply the mass conservationof the intracellular component at cell division then the average intracellularcontent h x i is computed from:d h x i d t = Z x max R ( x ) n ( x, t )d x − h x i Z x max Γ ( x ) n ( x, t )d x. (2.9)For the numerical solution of the deterministic cell population balancemodel (partial integro-differential equation) a moving boundary algorithm wasdeveloped and applied in Kavousanakis et al. (2009). The dynamics of (2.7)-as well as of the stochastic model- are fully determined through the definitionthe division rate Γ ( x ), the partitioning mechanism P ( x, x ′ ) and the single-cellreaction rate R ( x ), which contains the intracellular species described in themodel. These three functions are generally known as intrinsic physiologicalstate functions (IPSF).For the division rate, we consider the following normalised power law(Dien 1994): Γ ( x ) = (cid:18) x h x i (cid:19) m , (2.10)where m controls the sharpness of Γ . Large m values correspond to sharperdivision rates, i.e. a larger separation of rates under which cells below andabove a certain intracellular content can undergo division. The partitioning ofintracellular content upon division is discrete following the expression (2.6).The expression of the single-cell reaction rate R ( x ) depends on the geneticnetwork carried by the cells of the population. In this work, we study a classof regulatory genetic networks with positive feedback architecture, and in par-ticular the lac operon, which constitutes one of the most well-studied generegulatory networks (Beckwith and Zipser 1970; Miller and Reznikoff 1978). I. G. Aviziotis et al.
The lac operon consists of the promoter lac
P, the operator lac
O and threegenes, which encode the proteins required for the metabolism of lactose (seeFig. 1). The gene lac
Z encodes the enzyme β - galactosidase and the lac A en-codes transacetylase. The gene lac
Y is responsible for the encoding of theprotein lac permease, which is the potent transport facilitator.
Lac permeasecontributes to the transport of lactose or an analogue, such as TMG or IPTGin the interior of the cell. In the absence of lactose or of an extracellular in-ducer (TMG, IPTG), the expression of lac operon genes is turned off by theinhibitive action of lac I. Fig. 1
Sketch of the positive feedback loop network, lac operon with IPTG induction
The inhibitor lac
I binds to the operator site, lac
O, prevents binding of theRNA polymerase thus inhibiting the transcription of the genes’ DNA into thecorresponding mRNA. On the other hand, in the presence of lactose, TMGor IPTG, the inducer is transported into the cell, binds to the repressor lac
Ithrough a bimolecular reaction and the operator lac
O becomes free of lac
I,hence initiating the transcription. Upon expression of lac
Y, further transportof the inducer occurs at a higher rate resulting in further expression of thethree lac operon genes. Thus, the expression of lac
Y gene promotes its furtherexpression, and the network functions as an autocatalytic system or a positivefeedback loop.Genetic networks with genes carrying the positive feedback architecture ex-hibit bistable behavior at the single-cell level (Mantzaris 2007; Kavousanakis et al. 2009).A simple mathematical model, which captures the basic features of the posi-tive feedback loop architecture is described in Mantzaris (2005).. A simplifieddescription of the reaction steps has been proposed by Kepler and Elston(2001): oarse-grained analysis of stochastically simulated cell populations 9 O k → Y (3.1) O k → Y (3.2) Y λ → ∅ (3.3) O + Z φ ⇄ αφ O (3.4) Y + Y χ ⇄ βχ Z. (3.5)The O , O symbols denote the fraction of free and occupied operatorsites; Y is the monomer produced by the expression of lac Y gene in eitherthe occupied or the unoccupied state of the operator lac
O. The rate of lac
Yexpression in the unoccupied state ( k ) is significantly lower than that in theoccupied state ( k ). Z is the dimer of the gene product, which binds to the freeoperator site leading to an occupied state. Assuming that the production ratesof the monomer product are proportional to the fractions of unoccupied andoccupied operator sites and that the degradation of Y is a first order reaction,the single-cell monomer dynamics are described by:d Y d τ = k O + k O − λY, (3.6)where λ is the degradation rate constant. The number of operation sites isconserved: O + O = 1 . (3.7)It is further assumed that the unoccupied and occupied states are in equi-librium and the same holds for the dimerization reaction: O Z = αO , (3.8) Y = βZ, (3.9)where α and β are the equilibrium constants of the operator transition anddimerization reactions, respectively. Substituting (3.7) and (3.8) into (3.6)yields the following expression for the reaction rate:d Y d τ = k αβ + k Y αβ + Y − λY. (3.10)In order to obtain the dimensionless form of the reaction rate R ( x ), wenon-dimensionalize the intracellular content Y and time τ as follows: x = YY ∗ τ = tt ∗ . (3.11) Setting: k t ∗ /y ∗ = 1, π = k o /k , ρ = αβ/y ∗ , δ = λt ∗ (Mantzaris 2005),and substituting to (3.10) the following non-dimensional form for the rateof change of the dimensionless lac Y amount (the intracellular content x ) isobtained: d x d t ≡ R ( x ) = πρ + x ρ + x − δx, (3.12)where π is the relative rate of expression when the operator lac O is eitherfree or occupied ( π <<
1, due to significantly lower rate of gene expressionwhile being in the unoccupied state); ρ is a parameter inversely proportionalto the extracellular inducer concentration (IPTG) and δ is the dimensionlessdegradation rate.The main feature of lac operon is that its dynamics exhibit a bistablebehaviour, where two stable steady stats co-exist with an unstable steady statewithin a significant region of the parameter space ( π, ρ, δ ). If we examine thesimplest case of a cell population, where all individuals behave like the averagecell (homogeneous population), then the number cell density is expressed as n ( x, t ) = δ ( x − h x i ), where δ is the Dirac function. By introducing the numbercell density expression for the homogeneous cell population to (2.9), then thedynamics of h x i are described from the following ordinary differential equation:d h x i d t = R ( h x i ) − h x i = πρ + h x i ρ + h x i − δ h x i − h x i . (3.13)In Fig. 2, we present the dependence of the average intracellular contentof a homogeneous cell population, h x i , as a function of the inverse propor-tional of the IPTG concentration, ρ . The upper branch (I) corresponds tostable steady state solutions with high expression levels of lac Y, and the lowerbranch (III) to solutions with low expression level of lac
Y. The intermediatebranch of solutions corresponds to dynamically unstable steady states, i.e.,an infinitesimally small perturbation drives the system towards one of the co-existing stable steady state solutions (branch (I) or (III)). There exists a wideparametric region, ρ ∈ [0 . , . lac operonnetwork and present direct temporal simulation results in the next paragraph. oarse-grained analysis of stochastically simulated cell populations 11 ρ h x i stableunstableturning point (I) (II) (III)Homogeneous Cell Population Fig. 2
Dependence of average intracellular lac
Y content, h x i , of a homogeneous cell pop-ulation on the inverse IPTG inducer concentration, ρ . Parameter set values: π = 0 .
03, and δ = 0 . lac operon genetic network. Theinitial condition is a shifted Gaussian distribution with mean value µ = 0 . σ = 0 .
05 ( N (0 . , . )). The density function n ( x, τ )denotes the number of cells of dimensionless content x ( lac Y) at dimension-less time τ , divided by the total number of cells. The results are producedby incorporating: (a) Eq. (3.12), which describes the lac Y rate of change, (b)Eq. (2.10) which describes the division rate of each cell, and (c) the discretepartitioning mechanism (2.6) with f = 0 .
5, in the CNMC model (describedin Sec. 2). During the initial stage of the simulation the Gaussian distribu-tion splits into a three-humped distribution ( τ = 0 . τ = 2 . ρ . In order to computethe possibly co-existing steady state solutions through direct temporal sim-ulations, one needs to choose appropriate initial conditions. Indicatively, wepresent in Fig. 4 two asymptotic (steady) states (normalised with respect totheir average expression h x i ) for the same set of parameter values: f = 0 . m = 2, π = 0 . δ = 0 .
05, and ρ = 0 .
10, and different initial states. The solidline distribution has an average intracellular content of h x i = 0 .
61, when theinitial condition is a shifted Gaussian distribution N ( µ = 0 . , σ = 0 . ).The asymptotic distribution depicted by the dashed line, has an average intra- x n ( x ) t=0 x n ( x ) t=0.4 x n ( x ) t=1.0 x n ( x ) t=2.5 x n ( x ) t=4.0 x n ( x ) t=10.0 Fig. 3
Snapshots of the time evolution of the density function of a population of 10,000cells simulated with the CNMC model. Parameter set values: f = 0 . m = 2, π = 0 . δ = 0 .
05, and ρ = 0 . cellular content of h x i = 0 .
04, when the initial condition is a shifted Gaussiandistribution with a smaller mean value of µ = 0 .
05 and the same standard de-viation, σ = 0 .
05. Thus, solution multiplicity is also predicted by performingstochastic simulations of a cell population. However, the approach of explor-ing the solution space of steady states over a wide range of parameter valuesthrough direct temporal simulations is not computationally efficient. The com-putational load is significantly high, since it is required to perform long timeinterval executions of direct temporal simulations for different sets of param-eter values. Furthermore, the accurate determination of the bistability regionis not guaranteed (if not impossible), since stochastic noise at the vicinityof the bistability region (turning points) can cause unpredictable phenotypicswitches.In order to perform such an analysis for a system the description of whichis available at a fine-microscopic level, it is required to resort to alternativemethods of multiscale analysis. In this work, we adopt ideas from the equation-free methodology framework (Kevrekidis et al. 2003), which enables the per-formance of coarse-level numerical tasks, utilising information which originatesfrom short bursts of microscopic simulations. A description of this methodol-ogy and its application to the studied system of the isogenic cell population isprovided in the next section. oarse-grained analysis of stochastically simulated cell populations 13 x/ h x i n ( x / h x i ) µ = 0 . , σ = 0 . µ = 0 . , σ = 0 . Fig. 4
Asymptotic states (average of 100 copies for noise reduction purposes) of a popula-tion of 10,000 cells simulated with the CNMC model, and different initial conditions. Bothinitial conditions are shifted Gaussian distributions with mean µ and standard deviation σ = 0 .
05. The initial condition with µ = 0 .
25 converges to the steady state solution de-picted with the dashed line ( h x i = 0 .
61) and the initial state with µ = 0 .
05 to the steadystate distribution depicted with the solid line ( h x i = 0 . f = 0 . m = 2, π = 0 . δ = 0 .
05, and ρ = 0 . We are interested in studying the macroscopic behaviour of the heterogeneouscell population, utilising simulations which are performed at a microscopic /cell-level. In particular, our target is to compute the coarse asymptotic be-haviour of cell populations phenotype carrying the lac operon genetic net-work , and explore its dependence on the extracellular inducer concentration(IPTG). When closed descriptions of the macroscopic variables under studyare available, the performance of systems-level tasks, such as the computa-tion of steady state solutions, and bifurcation analysis is feasible through anarsenal of analytical and numerical techniques. In our case study, such de-terministic descriptions are available (CPB models), however their validity islimited only when large size populations are simulated. In order to study thesteady state behaviour of smaller size populations, one would need to performan extensive number of direct stochastic simulations, for adequately large timeintervals, so as to ensure the convergence of the system to a “coarse” steadystate solution. However, the computational requirements for such an approachare prohibitively large; in addition, and since the studied system of the lac operon network exhibits bistability, noise effects can lead to abrupt pheno-typic changes when the parameter value is at the proximity of critical turningpoints.
Here, we adopt a computational framework, which enables the perfor-mance of systems-level tasks, by utilising short bursts of appropriately ini-tialised stochastic simulations (instead of performing long time executions).The backbone of this framework, commonly known as the “equation-free”approach (Gear et al. 2002; Siettos et al. 2003; Kevrekidis et al. 2004), is the“coarse time-stepper”.4.1 Coarse time-stepperSuppose we are interested in studying the dynamics and/or asymptotic be-haviour of a coarse variable f . In our case, f is taken to be the cumulativedistribution function of cells, which are simulated by the CNMC model. Anupscaling of the discrete model should yield explicit equations governing thedynamics of f , e.g., the deterministic CPB models, which are valid understrong assumptions (continuum assumption). An alternative approach, whichcan bypass the derivation of such closures is to apply the so-called coarsetime-stepper , a schematic of which is shown in Fig. 5. The coarse time-stepperemploys two steps that link the fine-scale simulator with coarse-level compu-tations: restriction and lifting Fig. 5
Schematic of the coarse time-stepper for the model of an isogenic cell populationcarrying the lac operon network.
Restriction
Suppose that the stochastic simulation of N cells yields a set of intracellular lac Y contents: x = { x i } , i = 1 , .., N . Then the cumulative distribution func- oarse-grained analysis of stochastically simulated cell populations 15 tion (CDF) f is trivially computed by sorting into an ascending order the N -dimensional x vector ( x i ≤ x i +1 ) and plot the graph, which consists of thepoints: ( x , p ) = ( x i , p i = ( i − . /N ) , i = 1 , ..., N . This provides our restric-tion of microscopic data, U , to the macroscopic description, f , through anoperator M , i.e., f = M U . In case of many repetitions of the same stochasticsimulation, we compute the average restriction. Lifting
We choose to work with the CDF rather than the -natural to think- probabilitydensity function (PDF), n , because it is the derivative of a numerically speci-fied function (the CDF) and numerical differentiation is poorly conditioned. Inpractice, it is easier to work with the inverse CDF (ICDF), x ( p ), which givesthe intracellular content x i = x ( p i ) of a given cell, i . Assuming, that the ICDFis smooth enough, we can use a low-dimensional description of it based on thefirst few of an appropriate sequence of orthogonal polynomials (Gear 2001).Orthogonal polynomials are simple approximations, yet they are not useful forthe CDF which has a possible support from −∞ to + ∞ . However, the inverseCDF, x ( p ) has a finite support ( p ∈ [0 , j -degree orthogonal polynomials φ j ( p ), which must be monotone,non-decreasing and must lie between 0 and 1: x ( p ) ≈ q X j =0 α j φ j ( p ) . (4.1)Computationally, the φ j , j = 0 , ..., q polynomials are represented by theirvalues on the point set p and it is seldom to find it necessary to use morethan the first six basis functions ( q = 5). In Fig. 6(a)-(b), we present the upto 5 th order degree basis functions, which are used for the coarse descriptionof the heterogeneous cell population. They can be pre-computed (they dependonly on N ) and form a ( q + 1) × N matrix, Φ . Then the coefficients α = α i , i = 0 , ..., q in Eq. (4.1) are computed from α = Φ x , and the approximationto x can be computed from x ≈ Φ T α . Thus, we can evaluate the intracellularcontent of each cell using a small set of α values. This constitutes the liftingstep, which maps the macroscopic data to consistent microscopic realisationsthrough an operator, µ , i.e., U = µ f . This choice of lifting and restrictionoperators ensures that lifting from macroscopic to the microscopic and thenrestricting down again has no effect (except round-off), i.e., µ M ≈ I .Given the vector of α values we generate through lifting a consistent mi-croscopic realisation (intracellular content values x for N cells). We performa number of CNMC simulations for the N -cell population and after a shorttime interval τ we obtain its restriction (average CDF, and updated α values).Through lifting and restriction steps, we sidestep the derivation of a closure forthe evolution of the macroscopic α variables; instead we construct a discretetime-mapping: α ′ = G τ ( α ) , (4.2) p p o l y n o m i a l b a s i s f un c t i o n s , φ φ φ φ p p o l y n o m i a l b a s i s f un c t i o n s , φ φ φ φ (a) (b) Fig. 6
Graphical representation of the orthogonal polynomial basis functions (up to q = 5order) for the coarse description of cell populations. which can be utilised for the performance of black-box based numerical tasks.One of these tasks is the computation of the coarse steady state α ∗ , by solvingthe following set of non-linear equations: R ≡ α ∗ − G τ ( α ∗ ) = , (4.3)with the Newton-Raphson method, which requires during each iteration thesolution of the linearised system: ∂ R ∂α δα = − R → (cid:20) I − ∂ G τ ∂α (cid:21) δα = − R . (4.4)We re-iterate that an explicit expression of the operator G τ is not available;however, we can still approximate it numerically by appropriate initialisationof the coarse time-stepper. For example, the ( i, j ) element of matrix ∂ G τ ( α ) ∂α isapproximated by the forward finite difference scheme: (cid:18) ∂ G τ ∂α (cid:19) i,j ≈ ( G τ ( α j + ǫ )) i − ( G τ ( α j )) i ǫ , (4.5)where ǫ is a small number (perturbation number).Furthermore, an extension of this application is the performance of bifurca-tion analysis by means of pseudo-arc-length parameter continuation techniques(Keller 1977). Such techniques enable the computation of both stable and un-stable coarse steady state solutions, and the accurate computation of bistabil-ity limits. In the computations presented below, the continuation parameter isthe (dimensionless) inverse concentration of the extracellular inducer IPTG, ρ . In addition, the stability of each computed solution can be quantified by oarse-grained analysis of stochastically simulated cell populations 17 determining the eigenvalues of the matrix, ∂ G τ ( α ∗ ) ∂α . When at least one eigen-value crosses the unit circle in the complex plane, the solution is characterisedas dynamically unstable (Strogatz 1994). In this section, we apply the equation-free methodology in order to perform aparametric analysis of the asymptotic behaviour of heterogeneous cell popu-lation as a function of the extrcellular inducer IPTG. In the results presentedbelow, we choose as continuation parameter the inverse of the extracellularIPTG inducer ( ρ ), which can be adjusted experimentally. First, we validate ourcomputational methodology by comparing the findings of coarse bifurcationanalysis when applied on a large size cell population against results obtainedfrom the numerical solution of the deterministic CPB model. In particular,we examine the case of a 10,000 cell population with symmetric partitioningmechanism ( f = 0 . lac Y content as a function of thedimensionless ρ parameter. The comparison between the results obtained fromthe stochastic CNMC model and the deterministic CPB model shows verygood agreement, even when predicting unstable steady state solutions (dashedlines in Fig. 7). The computation of both stable and unstable steady statesolution is feasible through the application of pseudo arc-length parametercontinuation techniques, performed within the equation-free framework. Thisapplication enables the tracing of the entire solution space, contrary to theperformance of solely long time interval stochastic simulations, which requireappropriately chosen initial conditions in order to converge to the upper orlower branch of stable steady state solutions (see also discussion in Sec. 3.1).Furthermore, we can quantify the stability of the obtained steady statesolutions be determining the eigenvalues of the ∂ G τ ( α ∗ ) ∂α matrix. In Fig. 8we present the results of coarse stability analysis as obtained for (a) an up-per branch (b) intermediate branch and (c) lower branch steady solution. InFig. 8(a) and (c), all eigenvalues lie within the limits of the complex planeunit circle, and the corresponding steady state solutions are characterised asstable. The existence of one eigenvalue crossing the complex plane unit circlein Fig. 8(b), renders the corresponding steady state solution as dynamicallyunstable.Upon validation of our methodology, we can now investigate the effect ofnoise effects during division on the range of bistability, which is increasing forsmall size populations. In particular, we perform coarse bifurcation analysis forcell population sizes of N = 500 and N = 1000 cells and compare the resultswith the deterministic CPB model (2.7), which is numerically solved with thefinite element method (Kavousanakis et al. 2009). The bifurcation diagram ofthe steady state expression level of the average intracellular content h x i as ρ h x i deterministic CPB stochastic CNMC Fig. 7
Dependence of steady state average expression level, h x i on the inverse IPTG con-centration ρ . The black lines (solid and dashed) correspond to results obtained from thedeterministic CPB model, whereas the lines with open squares correspond to the stochasticCNMC model executed with 10,000 cells (average of 100 copies for noise reduction). Thedashed lines in both the deterministic and stochastic models correspond to unstable steadystate solutions. Parameter set values: f = 0 . m = 2, π = 0 .
03, and δ = 0 . a function of the dimensionless parameter ρ is presented in Fig. 9(a), whenequal partitioning is considered ( f = 0 . ρ values (or higher IPTG concentrationvalues) as the population size decreases. In particular, the upper ρ limit of thebistability region as computed from the deterministic CPB model is locatedat ρ = 0 . N = 1000 and 500 cells are ρ = 0 .
140 and 0.138, respectively. The lower ρ limit of the bistability regionalso shifts towards lower values for smaller size populations, ρ = 0 .
085 forthe deterministic CPB model and ρ = 0 . N = 1000 and 500cells, respectively. For noise reduction purposes, we performed 1000 copies ofCNMC simulations for the case of N = 1000 cells and 5000 copies for thecase of N = 500 cells. The short time interval τ which reports the discreteevolution of coarse variables is τ = 0 . h x i ) as obtainedfrom CNMC simulations with N = 1000 cells for the same parameter value, ρ = 0 . h x i = 0 .
62 for the upper solution branch and h x i = 0 .
038 for the lower solu-tion branch. The dashed line depicts the unstable steady state solution, whichcan be tracked only computationally through the application of the arc-lengthparameter continuation algorithm. The unstable solution has an average of h x i = 0 . oarse-grained analysis of stochastically simulated cell populations 19 −1 −0.5 0 0.5 1−1−0.500.51 Real Part I m ag i n a r y P a r t upper stable branch - ρ =0.1405 −1 −0.5 0 0.5 1−1−0.500.51 Real Part I m ag i n a r y P a r t intermediate unstable branch - ρ =0.1003 (a) (b) −1 −0.5 0 0.5 1−1−0.500.51 Real Part I m ag i n a r y P a r t lower stable branch - ρ =0.1169 (c) Fig. 8
Results of coarse stability analysis. Eigenvalues of the ∂ G τ ( α ∗ ) ∂α matrix correspondingto (a) an upper branch stable steady state ( ρ ≈ . ρ ≈ . ρ ≈ . f = 0 . m = 2, π = 0 .
03, and δ = 0 . N = 10 , N = 500 cells. One can clearly observe the shrinking of the bistabilityinterval for lower f values, suggesting that an increase of asymmetry in parti-tioning leads to the extinction of the bistability region. In particular, while thebistability region for symmetric partitioning ( f = 0 .
5) lies within the interval ρ ∈ [0 . , . ρ ∈ [0 . , . f = 0 .
2. Asimilar behavior is also observed for larger population sizes, as well as for thedeterministic CPB model.Interestingly enough, the effect of population size has increasing signif-icance for more asymmetric partitioning mechanisms. Indicatively, we referthe reader to Fig. 11, in which the asymmetry partitioning value is f = 0 . ρ ∈ [0 . , . N = 1000cells is simulated extends over the interval ρ ∈ [0 . , . f =0.5 h x i ρ deterministic1000 cells500 cells (a) x/ h x i n ( x / h x i ) upperintermediatelower (b) Fig. 9 (a) Effect of the population size on the average expression of the lac
Y gene steadystate ( h x i ), as a function of the inverse IPTG concentration ρ . The black line corresponds tothe deterministic CPB model, the line with full squares to corresponds to CNMC simulationwith N = 1000 cells, and the line with open triangles corresponds to the CNMC simulationwith N = 500 cells. Solid and dashed lines correspond to stable and unstable steady states inall cases. (b) Three coexisting (normalised with respect to h x i ) coarse steady state solutionsat ρ = 0 . N = 1000 cells (average of 1000copies). The solid line with open circles corresponds to the stable upper branch solution, thedashed line with open rectangles corresponds to the unstable intermediate branch solution,and the solid line with full triangles depicts the lower stable branch solution. Parameter setvalues: f = 0 . m = 2, π = 0 .
03, and δ = 0 . Large discrepancies between the stochastic and deterministic approachesare also observed in larger size cell populations when the partitioning mecha-nism is asymmetric. In Fig. 12(a), we present the dependence of the average lac
Y content on the parameter ρ as obtained from stochastic CNMC simula-tions with 5,000 cells and the direct CPB model for f = 0 .
2. The stochasticsimulations predict a significant shift of the bistability region towards smaller oarse-grained analysis of stochastically simulated cell populations 21 ρ h x i f =0.5 f =0.4 f =0.3 f =0.2 Fig. 10
Effect of the partitioning asymmetry factor, f , in the average expression of the lac Y gene steady state ( h x i ), as a function of the inverse IPTG concentration ρ . The linewith full circles corresponds to f = 0 .
5; the line with open squares corresponds to f = 0 . f = 0 . f = 0 . m = 2, π = 0 .
03, and δ = 0 . ρ h x i f = 0.2 deterministic1000 cells500 cells Fig. 11
Effect of the population size on the average expression of the lac
Y gene steady state( h x i ), as a function of the inverse IPTG concentration ρ when the asymmetry partitioningvalue is f = 0 .
2. The black line corresponds to the deterministic cell population balancemodel, the line with full squares corresponds to CNMC simulations with N = 1000 cells,and the line with open triangles corresponds to CNMC simulations with N = 500 cells.Solid and dashed lines correspond to stable and unstable steady state solution branches.Parameter set values: m = 2, π = 0 .
03, and δ = 0 . ρ values. However, when the partitioning mechanism is considered to be sym-metric ( f = 0 .
5) these discrepancies practically vanish (see Fig. 12(b)). ρ h x i Asymmetric Partitioning - f = 0.2 deterministicCNMC - 5,000 cells Symmetric Partitioning - f =0.5 ρ h x i deterministicCNMC - 5,000 cells (a) (b) Fig. 12
Dependence of the average lac
Y content on the inverse concentration of IPTG( ρ ) for (a) asymmetric partitioning, f = 0 .
2, and (b) symmetric partitioning, f = 0 . π = 0 . δ = 0 . m is pre-sented in Fig. 13, when simulations of N = 500 are executed with the CNMCalgorithm. Larger m values correspond to higher cell division rates. The bista-bility region clearly shrinks and shifts towards higher IPTG concentrationvalues for increasing m . The same behaviour has been also observed whenconsidering the deterministic CPB model (Kavousanakis et al. 2009). Here,the effect of cell division sharpness when we compare the deterministic modelwith the stochastic simulations is not as prominent as in the case of the effectof partitioning asymmetry parameter, f . The larger discrepancies in the phe-notypic population behaviour as predicted from the deterministic and CNMCmodelling is observed for sharp division rates. Indicatively, we present the bi-furcation diagrams obtained from the deterministic CPB and the stochasticCNMC model with N = 500 cells for m = 4 in Fig. 14. The bistability regionfor the deterministic model is located within the interval ρ ∈ [0 . , . ρ ∈ [0 . , . ρ = 0 . oarse-grained analysis of stochastically simulated cell populations 23 ρ h x i m =1 m =2 m =3 m =4 Fig. 13
Effect of the division rate sharpness (quantified by m ) in the average expression ofthe lac Y gene steady state ( h x i ), as a function of the inverse IPTG concentration ρ . Higher m values correspond to sharper division rates. The line with full circles corresponds to m = 1;the line with open squares corresponds to m = 2; the line with full triangles corresponds to m = 3 and the line with open circles to m = 4. Solid and dashed lines correspond to stableand unstable steady state solutions, respectively. The following set of parameter values isused: f = 0 . π = 0 .
03, and δ = 0 . ρ h x i m =4 deterministic1000 cells500 cells Fig. 14
Effect of the population size on the average expression of the lac
Y gene steady state( h x i ), as a function of the inverse IPTG concentration ρ when cell division rate is sharp, m = 4. The black line corresponds to the deterministic cell population balance model,the line with full squares corresponds to CNMC simulations with N = 1000 cells, and theline with open triangles corresponds to CNMC with N = 500 cells. Solid and dashed linescorrespond to stable and unstable steady state solutions, respectively. Parameter set values: f = 0 . π = 0 .
03, and δ = 0 . h x i ≈ . τ ≈ lac Y ( h x i = 0 . h x i values is presented in Fig. 16. The extracellular inducer ρ = 0 . large period, during which the cell population exhibits low level lac Y pheno-type, stochastic noise drives the system to states of high average expressionlevels. The time at which the phenotypic switch occurs is random. time , τ h x i ρ h x i upper branchlower branch ρ = 0 . N = 500 cells Fig. 15
One stochastic simulation starting from a coarse steady state solution at ρ = 0 . h x i = 0 . h x i = 0 .
03) after a long dimensionless time elapse ( τ ≈ ρ ,and with open circles the corresponding stable steady states at ρ = 0 . f = 0 . π = 0 . δ = 0 . N = 500 cells, ρ = 0 . time, τ h x i ρ h x i lower branchupper branch ρ = 0 . N = 500 cells Fig. 16
One stochastic simulation starting from a coarse steady state solution at ρ = 0 . h x i = 0 .
037 (lower branch). Due to stochastic noise the cell population jumps towardsthe upper stable branch ( h x i = 0 .
62) after a long dimensionless time elapse ( τ ≈ ρ , and with open circles the corresponding stable steady states at ρ = 0 . f = 0 . π = 0 . δ = 0 . N = 500 cells, ρ = 0 . We presented a multiscale methodology for the computation of the macro-scopic steady state behaviour of cell populations, which carry the lac operongenetic network, and are simulated with stochastic fine-scale models. In par-ticular, we employ the CNMC model presented in Mantzaris (2006), whichincorporates the effect of random events during division and describes withsufficient accuracy the behaviour of small size cell populations. Despite its ac-curacy, the stochastic modelling approach lacks efficiency when a coarse-levelanalysis is required, due to significant computational requirements. In addi-tion, for systems which exhibit solution multiplicity within a range of parame-ter values, the strategy of performing solely direct temporal simulations is notthe appropriate one, since a very large number of different initial states needto be tested in order to drive the system to different basins of attraction. Amore systematic analysis is applicable when alternative multiscale techniquesare applied.In this work, we adopt techniques from the “equation-free” framework(Kevrekidis et al. 2003; Kevrekidis et al. 2004), which enables the performanceof numerical tasks at a systems-level, such as the computation of steady statesolutions and the performance of bifurcation analysis at a coarse level of de-scription. We compute both stable and unstable macroscopic steady state so-lutions, and determine with accuracy the critical extracellular inducer con-centration values (limits of bistability), which signal transitions between thedifferent lac
Y expression levels. Bifurcation diagrams are computed for dif-ferent sets of parameter values, in order to examine the effect of partitioningasymmetry of the intracellular content during division, as well as the divisionrate sharpness on the range of bistability. When the partitioning mechanismis highly asymmetric, a shift towards higher IPTG values, as well as shrink-ing of the bistability range is observed. This effect was also demonstrated in(Kavousanakis et al. 2009) for the deterministic CPB model case.A comparison between the stochastic model and the deterministic CPBmodel shows that stochasticity tends to shift the bistability range towardshigher IPTG concentration values. The effect of random events during divi-sion becomes increasing when CNMC simulations are performed with low sizecell populations, and we demonstrated that in the case of highly asymmetricpartitioning the discrepancies between stochastic and deterministic models aresignificant. Furthermore, stochastic simulations can predict phenotypic switch-ing between different lac
Y expression levels, and we illustrated such transitionsat the neighbourhood of the bistability upper and lower limits. Such non-trivialdynamical behaviour cannot be predicted by deterministic modelsIn this work, we focused on the effect of heterogeneity, originating from theuneven distribution of intracellular components of mother cells to their off-springs and random events during division. A second source of heterogeneitystems from the fact that intracellular chemical reactions are regulated by smallconcentrations of chemical species (regulatory molecules) (Alberts et al. 1994).Experimental studies (Blake et al. 2003; Elowitz et al. 2002) have demonstrated that the process of gene expression on small synthetic genetic networks isstochastic. A modelling approach that incorporates the effect of single-cell in-trinsic noise effects is presented in Mantzaris (2007), where the deterministicreaction rate R ( x ) (see Eq. (2.2),) which describes the evolution of the intracel-lular species, is replaced by a Langevin equation (Miller and Reznikoff 1978).The region of bistability has been determined for the single-cell stochasticmodel, however a systematic analysis of the intrinsic noise effect on the population-level phenotype has not yet been performed; such an analysis can be performedby applying the equation-free methodology.Finally, we report that the presented coarse-grained analysis can be triv-ially adjusted for cell populations, which carry different genetic network ar-chitecture, such as the genetic toggle switch (Gardner et al. 2000), which alsoexhibits bistable behaviour, and switching between the co-existing stable statesis possible through thermal or chemical induction. References
Alberts et al. 1994. Alberts B, Bray D, Lewis J, Raff M, Roberts K, Wat-son JD (1994) Molecular biology of the cell, Garland Publishing, 3rd Edition,New YorkAltschuler and Wu 2010. Altschuler SJ, Wu LF (2010) Cellular Hetero-geneity: Do Differences Make a Difference? Cell, 141:559-563. doi:10.1016/j.cell.2010.04.033Avery 2006. Avery SV (2006) Microbial cell individuality and the underlyingsources of heterogeneity. Nat Rev Microbiol, 4:577-587. doi: 10.1038/nrmi-cro1460Beckwith and Zipser 1970. Beckwith JR, Zipser D (1970) The LactoseOperon, Cold Spring Harbor Laboratory, Cold Spring Harbor, New YorkBlake et al. 2003. Blake WJ, Kaern M, Cantor CR, Collins JJ (2003) Noise ineukaryotic gene expression. Nature, 442:633-637. doi:10.1038/nature01546Block et al. 1990. Block DE, Eitzman PD, Wangensteen JD, Srienc F (1990)Slit scanning of Saccharomyces cerevisiae cells: quantification of asymmetriccell division and cell cycle progression in asynchronous culture. BiotechnolProgr, 6:504-512. doi:10.1021/bp00006a015Booth 2002. Booth JR (2002) Stress and the single cell: intrapopulation di-versity is a mechanism to ensure survival upon exposure to stress, Int J FoodMicrobiol, 78:19-30. doi:10.1016/s0168-1605(02)00239-8Chung and Stephanopoulos 1995. Chung JD, Stephanopoulos G (1995) Stud-ies of transcriptional state heterogeneity in sporulating cultures of
Bacillussubtilis . Biotechnol Bioeng, 47:234-242. doi:10.1002/bit.260470215Davidson and Srette 2008. Davidson CJ, Srette MG (2008) Individual-ity in Bacteria. Annu Rev Genet, 42:253-268. doi: 10.1146/an-nurev.genet.42.110807.091601Delbr¨uck 1945. Delbr¨uck M (1945) The burst size distribution in the growthof bacterial viruses (Bacteriophages). J Bacteriol, 50:131-135. oarse-grained analysis of stochastically simulated cell populations 27
Dien 1994. Dien BS (1994) Aspects of cell division cycle related behavior of
Saccharomyces cerevisiae . Growing in batch and continuous culture: A single-cell growth analysis, Dissertation, University of MinnesotaEakman et al. 1966. Eakman JM, Fredrickson AG, Tsuchiya HM (1966)Statistics and dynamics of microbial cell populations. Chem Eng Prog, 62:37-49.Elowitz et al. 2002. Elowitz MB, Levine AJ, Siggia ED, Swain PS (2002)Stochastic gene expression in a single cell. Science, 297:1183-1186.doi:10.1126/science.1070919Fedorrof and Fontana 2002. Fedorrof N, Fontana W (2002) Small numbers ofbig molecules. Science, 297:1129-1131. doi: 10.1126/science.1075988Fredrickson et al. 1967. Fredrickson AG, Ramkrishna D, Tsuchiya HM (1967)Statistics and dynamics of prokaryotic cell populations. Math Biosci, 1:327-374. doi:10.1016/0025-5564(67)90008-9Gardner et al. 2000. Gardner TS, Cantor CR, Collins JJ (2000) Construc-tion of a genetic toggle switch in Escherichia coli. Nature, 403:339-342.doi:10.1038/35002131Gear 2001. Gear CW (2001) Projective Integration Methods for Distributions,NEC Technical Report, 2001-130, Princeton, NJGear et al. 2002. Gear CW, Kevrekidis IG, Theodoropoulos C (2002)“Coarse” integration/bifurcation analysis via microscopic simulators: micro-Galerkin methods, Comput Chem Eng, 26:941-963. doi:10.1016/S0098-1354(02)00020-0Hatzis et al. 1995. Hatzis C, Srienc F, Fredrickson AG (1995), Multistagedcorpuscular models of microbial growth: Monte Carlo simulations, BioSys-tems, 36:19-35. doi:10.1016/0303-2647(95)01524-OKavousanakis et al. 2009. Kavousanakis ME, Mantzaris NV, Boudouvis AG(2009) A novel free boundary algorithm for the solution of cell populationbalance model, Chem Eng Sci, 64:4247-4261. doi:10.1016/j.ces.2009.06.054Keller 1977. Keller HB (1977) Numerical solution of bifurcation and nonlin-ear eigenvalue problems. In: Rabinowitz P (ed) Applications of BifurcationTheory, Academic Press, New York, pp: 359-384.Kelley 1995. Kelley CT (1995) Iterative methods for linear and nonlinearequations, SIAM, Philadelphia.Kepler and Elston 2001. Kepler TB, Elston TC (2001) Stochasticity in tran-scriptional regulation: origins, consequences and mathematical representa-tions, Biophys J, 81:3116-3136. doi:10.1016/S0006-3495(01)75949-8Kevrekidis et al. 2003. Kevrekidis IG, Gear CW, Hyman JM, Kevrekidis PG,Runborg O, Theodoropoulos C (2003) Equation-free, coarse-grained multi-scale computation: enabling microscopic simulators to perform system-levelanalysis. Comm Math Sci, 1:715-762.Kevrekidis et al. 2004. Kevrekidis IG, Gear CW, Hummer G (2004) Equation-free: the computer-aided analysis of complex multiscale systems, AIChE J,50:1346-1355. doi:10.1002/aic.10106Lee and Matsoukas 2000. Lee K, Matsoukas T (2000) Simultaneous coagula-tion and break-up using constant-N Monte Carlo, Powder Technol, 110:82-89. doi:10.1016/S0032-5910(99)00270-3Libby and Rainey 2011. Libby E, Rainey PB (2011) Exclusion rules, bottle-necks and the evolution of stochastic phenotype switching, Proc R Soc B-BiolSci, 7:3574-3583. doi: 10.1098/rspb.2011.0146Mantzaris et al. 2001. Mantzaris NV, Daoutidis P, Srienc F (2001) Numer-ical solution of multivariable cell population balance models. I:Finite dif-ference methods, Comput Chem Eng, 25:1411-1440. doi:10.1016/S0098-1354(01)00709-8Mantzaris 2005. Mantzaris NV (2005) A cell population balance model de-scribing positive feedback loop expression dynamics, Comput Chem Eng,29:897-909. doi:10.1016/j.compchemeng.2004.09.012Mantzaris 2006. Mantzaris NV (2006) Stochastic and deterministic simula-tions of heterogeneous cell population dynamics, J Theor Biol, 241:690-706.doi:10.1016/j.jtbi.2006.01.005Mantzaris 2007. Mantzaris NV (2007) From single-cell genetic architecture tocell population dynamics: Quantitatively decomposing the effects of differentpopulation heterogeneity sources for a genetic network with positive feedbackarchitecture, Biophys J, 92:4271-4288. doi:10.1529/biophysj.106.100271McAdams and Arkin 1998. McAdams HH, Arkin A (1998) Simulation ofprokaryotic genetic circuits. Annu Rev Biophys Biomol Struct, 27:199-224.doi: 10.1146/annurev.biophys.27.1.199McAdams and Arkin 1999. McAdams HH, Arkin A (1999) It’s a noisy busi-ness! Genetic regulation at the nanomolar scale. Trends Genet, 15:65-69. doi:10.1016/S0168-9525(98)01659-XMiller and Reznikoff 1978. Miller JH, Reznikoff WS (1978), The Operon, ColdSpring Harbor Laboratory, Cold Spring Harbor, New York.Oh et al. 2004. Oh P, Li Y, Yu J, Durr E, Krasinska KM, Carver LA, Testa JE,Schnitzer JE (2004) Subtractive proteomic mapping of the endothelial sur-face in lung and solid tumours for tissue-specific therapy, Nature, 429:629-635. doi:10.1038/nature02580Powell 1956. Powell EO (1956) Growth rate and generation time of bacteriawith special reference to continuous culture, J Gen Microbiol, 15:492-511.doi:10.1099/00221287-15-3-492Ptashne 1987. Ptashne M (1987) A genetic switch: gene control and phagelamda, Blackwell, Cambridge (1987)Russo-Marie et al. 1993. Russo-Marie F, Roederer M, Sager B, Herzen-berg LA, Kaiser D (1993) β -galactosidase activity in single dif-ferentiating bacterial cells, Proc Natl Acad Sci, 90:8194-8198.doi:10.1073/pnas.90.17.8194Shah et al. 1976. Shah BH, Borwanker JD, Ramkrishna D (1976) MonteCarlo simulation of microbial population growth, Math Biosci, 31:1-23.doi:10.1016/0025-5564(76)90037-7Siettos et al. 2003. Siettos CI, Armaou A, Makeev AG, Kevrekidis IG (2003)Microscopic/stochastic timesteppers and coarse control: a kinetic MonteCarlo example, AIChE J, 49:1922-1926. oarse-grained analysis of stochastically simulated cell populations 29 Smith and Matsoukas 1998. Smith M, Matsoukas T (1998), Constant-numberMonte Carlo simulation of population balances. Chem Eng Sci, 53:1777-1786.Spudich and Koshland 1976. Spudich JL, Koshland Jr DE (1976) Non-genetic individuality: chance in the single cell, Nature, 262:467-471.doi:10.1038/262467a0Stocker 1949. Stocker BAD (1949) Measurements of rate of mutation of flag-ellar antigenic phase in Salmonella typhimurium. J Hyg Camb, 47:398-413.Strogatz 1994. Strogatz SH (1994) Nonlinear Dynamics and Chaos: with Ap-plications to Physics, Biology, Chemistry, and Engineering, Westview Press,Colorado.Sumner and Avery 2002. Sumner ER, Avery SV (2002) Phenotypic hetero-geneity: differential stress resistance among individual cells of the yeast