Effect of Intrinsic Noise on the Phenotype of Cell Populations Featuring Solution Multiplicity: An Artificial lac Operon Network Paradigm

Heterogeneity in cell populations originates from two fundamentally different sources: the uneven distribution of intracellular content during cell division, and the stochastic fluctuations of regulatory molecules existing in small amounts. Discrete stochastic models can incorporate both sources of cell heterogeneity with sufficient accuracy in the description of an isogenic cell population; however, they lack efficiency when a systems level analysis is required, due to substantial computational requirements. In this work, we study the effect of cell heterogeneity in the behaviour of isogenic cell populations carrying the genetic network of lac operon, which exhibits solution multiplicity over a wide range of extracellular conditions. For such systems, the strategy of performing solely direct temporal solutions is a prohibitive task, since a large ensemble of initial states needs to be tested in order to drive the system—through long time simulations—to possible co-existing steady state solutions. We implement a multiscale computational framework, the so-called “equation-free” methodology, which enables the performance of numerical tasks, such as the computation of coarse steady state solutions and coarse bifurcation analysis. Dynamically stable and unstable solutions are computed and the effect of intrinsic noise on the range of bistability is efficiently investigated. The results are compared with the homogeneous model, which neglects all sources of heterogeneity, with the deterministic cell population balance model, as well as with a stochastic model neglecting the heterogeneity originating from intrinsic noise effects. We show that when the effect of intrinsic source of heterogeneity is intensified, the bistability range shifts towards higher extracellular inducer concentration values.


Introduction
The phenotype of a cellular population is not exclusively the result of single-cell level complex chemical networks; cells interact with each other leading to phenotypic variations amongst the individual members of isogenic populations, a phenomenon commonly known as cellular heterogeneity. The literature reporting cellular heterogeneity is large and here we cite some representative examples, e.g., the variations of phage burst size [1], the transcriptional states heterogeneity in sporulating cultures of Bacillus subtilis [2], and the lysogenic states of phageinfected bacteria [3,4]. The effect of heterogeneity has been studied in transcriptomics [5,6], metabolomics [7], pathogens [8][9][10][11][12], as well as in mitochondrial activity [13][14][15]. It is also noteworthy to report that the design of contemporary biomedical therapies and of synthetic circuits with robust performance incorporates the effects of heterogeneity [16][17][18].
For an isogenic cell population residing in a uniform extracellular environment, there exist two fundamentally different sources of heterogeneity [19]: The first one originates from unequal partitioning of the mother intracellular content to its offsprings during division [20,21]. The unevenly distributed regulatory molecules lead to different phenotypes, and the phenomenon is repeated due to the operation of the cell cycle. This type of heterogeneity is called extrinsic [19]. The regulatory molecules, which control the network of intracellular reactions and determine the cells phenotype exist in small amounts [22][23][24], and even small fluctuations can lead to an uncontrolled-uncertain outcome (phenotype). Thus, cells with approximately the same amount of regulatory molecules can feature utterly different phenotypic behaviour; this type of heterogeneity is called intrinsic [19].
Several models simulating heterogeneous populations have been developed in order to elucidate the effect of the different sources of heterogeneity. Shah et al. [25] were the first to model the stochastic behaviour of cell populations by developing a Monte Carlo algorithm for the dynamics of the cell mass distribution. Hatzis et al. [26] extended this algorithm to describe the dynamics of a growing population of phagotrophic protozoa. However, these models are computationally expensive due to the exponentially growing number of simulated cells of the population. To overcome the extensive requirements in CPU time, Constant-Number Monte Carlo (CNMC) algorithms are used [27,28] simulating a constant number of cells that are assumed to be a representative sample of the studied population. More recent studies include the work of Shu et al. [29] in which the population balance models incorporate extrinsic heterogeneity and intracellular stochastic processes through Itô stochastic differential equations; a chemical master equation for the population level, which models uncertainty of intracellular reactions, DNA duplication and content partitioning has been presented in [30][31][32]. Zechner et al. [33] used low-order moments through the moment closure approach to approximate intrinsic and extrinsic distributions; Toni and Tidor [17] employed van Kampen's O-expansion for the approximation of intrinsic stochastic dynamics and incorporated extrinsic heterogeneity through variability of kinetic parameters and initial conditions. Finally, we report agent-based modelling approaches of cell, which have been presented in [34][35][36][37].
In this work, we apply a CNMC algorithm developed by Mantzaris [27] modelling the dynamics of an isogenic population. The algorithm takes into account the random nature of cell division, and unequal partitioning of intracellular content at cell division modelling extrinsic heterogeneity. In this model, interactions between individual cells are not taken into consideration. In addition, a Langevin approximation [19] of the reaction dynamics at the single-cell level is used to incorporate the effect of intrinsic heterogeneity. In our case study, all cells carry the lac operon genetic network [38,39]; it is an artificial genetic network with a positive feedback architecture, featuring solution multiplicity within a range of extracellular inducer (IPTG, TMG or lactose) concentration values at the single-cell level. Bistability is also present at the population level, however the range of solution multiplicity is significantly altered. This has been demonstrated in [40,41] by solving deterministic cell population balance models, which incorporate the effect of extrinsic heterogeneity. In order to quantify the effect of intrinsic heterogeneity, we need to apply stochastic modelling (here the CNMC algorithm), since deterministic models disregard the discrete nature of cell populations and cannot incorporate the intrinsic source of heterogeneity.
The disadvantage of stochastic modelling is its inefficiency to perform a systems level analysis of a population's long time behaviour for a wide range of parameter values, by executing exclusively temporal simulations. In a recent publication [21], we demonstrated the application of a multiscale framework, which enables the application of well-established numerical algorithms utilizing short bursts of dynamic stochastic simulations. This framework is known as "equation-free" method [42][43][44]; it wraps around fine scale models interchanging information between microscopic (single-cell) and macroscopic (population) level, in order to perform numerical tasks including the computation of steady-state solutions, stability, and bifurcation analysis. The interchange of information between the different levels of description is feasible by constructing a discrete time-mapping, the coarse time-stepper, which reports the evolution of macroscopic quantities of interest at discrete time instances.
Here, we apply the equation-free method and wrap it around the CNMC model [27], which describes the dynamics of an isogenic cell population. When each individual of the population carries the lac operon genetic network, the existence of a range of extracellular inducer concentration values (IPTG) is expected, within which the population can feature utterly different phenotypic behaviour, i.e., high or low expression levels of the lacY gene can be both observed for the same parameter value. We apply bifurcation analysis, by means of the pseudo arc-length parameter continuation technique [45], in order to accurately determine the limits of the solution multiplicity region, which is then compared with the results obtained from the deterministic cell population balance models. This comparison, which quantifies the effect of intrinsic heterogeneity on the phenotype of a cell population, reveals some interesting findings. In particular, when the intrinsic heterogeneity effect is strengthened, the bistability interval is located at higher extracellular inducer IPTG concentration values. This bistability interval can be even expunged when we consider high asymmetry in the partitioning mechanism and sharper rates for the cells division. The accurate determination of the bistability interval is important for the understanding of possible phenotypic switching when the cell population operates at the proximity of the limits of this interval.
The paper is organized as follows: In the next section we present a simplified single-cell reaction rate expression, which describes the dynamics of the lac operon genetic network. In order to incorporate intrinsic noise effects, the reaction rate is augmented with a gaussian noise term (Langevin approach). We then briefly describe the homogeneous model, which neglects all sources of heterogeneity, the deterministic cell population balance model, which incorporates extrinsic heterogeneity, and finally the stochastic CNMC algorithm, which takes into account both extrinsic and intrinsic source of heterogeneity. Upon description of the stochastic CNMC model, we present the multiscale equation-free methodology, which enables the performance of coarse steady-state, bifurcation and stability analysis. In the Results section, we present steady state solutions of cell populations as a function of the IPTG concentration for different levels of intrinsic and extrinsic heterogeneity. Finally, we outline the main findings of this work and propose potent future research directions.

Single-Cell Lac Operon Dynamics
In this work, we study the dynamics of an isogenic population carrying the lac operon genetic network. It consists of a promoter (lacP), an operator (lacO) and three genes (lacZ, lacY, lacA) which encode the necessary proteins for the metabolism of lactose. The three lac operon genes are inhibited by a lacI repressor. The lacI repressor binds to the operator site, lacO, and prevents binding of the RNA polymerase inhibiting the transcription of the three genes' DNA to the corresponding mRNA.
In the presence of lactose or its analogues, TMG or IPTG, the inducer is transported into the cell, binds to the repressor lacI through a bimolecular reaction and the operator lacO becomes free of lacI, hence initiating the transcription. Upon expression of lacY, the transport of the inducer is facilitated resulting in further expression of the three lac operon genes. Thus, the expression of lacY gene promotes its further expression, and the network functions as an autocatalytic system or a positive feedback loop [21,40,41].
A simplified description of the reaction steps has been described in [2,40,46] as follows: O 0 and O 1 denote the free and occupied operator sites. The repressor R binds to and unbinds from the occupied and the free operator sites at a rate proportional to k r and k −r , respectively. The extracellular IPTG inducer, I ex , is transported through the cell membrane at a rate proportional to k t . Binding of the intracellular IPTG, I, to the repressor molecules occurs at rate proportional to k r2 , and the rate constant of unbinding is k −r2 . We also consider the rate of expression of the lac operon genes to produce lac permease, Y, to be proportional to the amount of free operator sites O 0 ; Lac permease facilitates the transport of IPTG, playing the role of the enzyme and I ex playing the role of the substrate in a scheme following the Michaelis-Menten kinetics. Finally, we assume that Y and I degrade following first-order kinetics at a rate constant λ Y and λ I , respectively. Values of the kinetic constants are given in Table 1.
By assuming that the total number of operator sites and repressor molecules remain constant, i.e., O T = O 0 + O 1 , and R T = R + I 2 R, one can deduce a deterministic model describing the single-cell lac operon dynamics: The mathematical description of the lac operon dynamics can be further simplified by assuming that [2,40] 2. the repressor inactivation reaction 2I + R Ð I 2 R is in equilibrium with equilibrium constant: Based on these assumptions, the lac operon dynamics can be described through the following set of ordinary differential equations for Y, and the intracellular IPTG: By defining the following dimensionless quantities: and setting the dimensionless time: t ¼ t=t, the dimensionless lac permease: x ¼ ½Y=ŷ, and the dimensionless intracellular IPTG: v ¼ ½I=Î, Eqs (3) become: with κ a dimensionless degradation parameter and: By substituting with the values reported in Table 1: l Y l Y þk t % 3 Â 10 À3 << 1, which suggests that the dynamics of the dimensionless IPTG amount, v, are much faster compared to the dynamics of the dimensionless lacY amount, x. Furthermore, σ % κ, which yields: x % v, and the lac operon dynamics can be described by the following equation for the dimensionless amount of lac permease: In order to numerically verify the validity of this reduction, we present in Fig 1 a comparison between the full deterministic model (Eqs (2)), and the reduced model (Eq (7)) for two different external IPTG concentration values leading to an uninduced (low lac permease amount) and an induced (high lac permease amount) state. In particular, we present the dynamics of the lac Y concentration for [I ex ] = 10 and 40 μM, showing good agreement between the full and reduced deterministic models. The lac permease concentration of roughly 400 nM at the fully induced state is in agreement with experimental data reported in [49].
Langevin approximation of the single-cell reaction rate The model described by Eq (7) is deterministic and it does not account for intrinsic noise effects, related with the small amount of regulatory molecules and slow operator fluctuations [19]. In order to take into account intrinsic noise effects, the reaction rate is re-formulated by adopting a Fokker-Planck approximation of the chemical master equation proposed in [19]. The approximation is based on the assumption that the operator fluctuations occur on a faster time scale compared to the production and degradation rate of the monomer. The monomer concentration which is denoted by ρ i ((x), t) satisfies the equation: where K is the transition matrix containing the reaction rates for transitions between the operator's chemical states and L is the diagonal matrix of the form: g is a matrix the j th column of which contains the net production rates of the q monomer species when the operator is in the j th chemical state. Likewise, h is the matrix the j th column of which contains the diffusion coefficients of the monomer species in the j th chemical state. Kepler and Elston [46] derived the following Fokker-Planck equation:  (2) with the reduced model (dashed line), which is described by Eq (7)  where p(x, τ) is the probability density function of the random variable, x and: yields the following expressions for the drift and diffusion terms, A(x) and B(x), respectively: BðxÞ and the Langevin stochastic differential equation (SDE) for the reaction rate is given by: where ξ(τ) is a Gaussian white noise process. The parameter K is a measure of the rate of the operator fluctuations and y Ã is the reference number of molecules, quantifying the two main sources of intrinsic heterogeneity, i.e., slow operator fluctuations and small numbers of molecules. The Langevin approximation is computationally advantageous over Monte Carlo simulations, by avoiding simulations of the full process, and sampling paths of the process (generated from Eq (15)) in shorter time [46]. We note that for very fast operator fluctuations (K ! 1) and/or large number of molecules (y Ã ! 1), the Langevin Eq (15) reduces to the deterministic expression of the reaction rate, Eq (7). Special treatment should be provided in cases, where the number of reference molecules is small enough (y Ã < 20), and the Langevin approach becomes unreliable to provide the probability of large deviation from the typical evolution (rare events) [51,52]. However, in this work the results presented correspond to sufficiently large amount of regulatory molecules and sufficiently slow operator fluctuations, which can be accurately modelled by the efficient Langevin approximation. In Fig 2, we illustrate the accuracy of the Langevin approximation by comparing the dynamics of Eq (15) for y Ã = 20 and K = 200, with Monte Carlo simulations of the chemical Master-equation for the reaction system described by Eqs (1), using the direct Gillespie's algorithm [53,54]. In particular, we present the average evolution of 100 simulations of the single-cell lac operon dynamics for [I ex ] = 15μM, and [I ex ] = 45μM leading to an uninduced ([Y] 1 % 15.5nM) and an induced ([Y] 1 % 400nM) state, respectively.
The stochastic single-cell model (Eq (15)) is more realistic compared to the corresponding deterministic one (Eq (7)), since it can capture the intrinsic noise effects, however we are interested in studying the behaviour at the cell population level incorporating also the effect of heterogeneity.

Deterministic Cell Population Balance Model
The extrinsic source of heterogeneity can be captured using a special class of models, which are commonly known as cell population balance models [19,27,41,55]. In particular, we study the dynamics of the number density function n(x, t) corresponding to the number of cells per biovolume unit, which at time t have intracellular content between x and x + dx [55]: where R(x) is the single-cell reaction rate (Eq (7)), which quantifies the rate of production or consumption of the intracellular content, x; Γ(x) is the division rate of each cell with content x.
For the division rate we adopt the following normalized power law [56]: where m regulates the sharpness of the division rate, and hxi is the average content of the population. P(x, x 0 ) is a partition probability density function, which describes the mechanism of the distribution of a mother cell intracellular content among the two daughter cells. A simple formulation of the partition probability density function is: where f is the asymmetry parameter, and δ is the Dirac function. According to Eq (18), a mother cell with intracellular content x will provide a fraction fx to one of its offsprings, and the remaining (1 − f)x to the second one. The asymmetry parameter ranges within the interval [0,0.5], with lower values corresponding to asymmetric partitioning and f = 0.5 describing a symmetric partitioning mechanism. R(x), Γ(x), and P(x, x 0 ) describe processes occurring at the single-cell level and are generally known as intrinsic physiological state functions (IPSF).
A simplification of the cell population balance model is the homogoneous model, which assumes that all cells behave like the average cell, with content hxi h . The homogeneous model is derived by Eq (16), when the number density function is expressed as n(x, t) = δ(x − hxi h ) [2,19,41]. Then the dynamics of the homogenous population are given by the following differential equation: The Stochastic CNMC Model For a more realistic simulation of the dynamics of isogenic populations capturing all sources of heterogeneity, we implement the CNMC stochastic algorithm developed by Mantzaris [19], which is described below. We consider isogenic populations with cells carrying the same gene regulatory network, whose random state, S τ , at a dimensionless time instance, τ, is: is the intracellular content of cell i and N is the constant number of cells considered for simulation.

Time intervening two successive division events
The time interval, T, intervening two successive division events is a random variable depending on S τ with a cumulative distribution function given by the expression [27]: where z denotes the set of time instances during which the next division event may occur for a given state, S τ . A random variable, p 1 , is generated from a uniform distribution [0, 1] and T is computed by solving the nonlinear equation: During the time interval, T, the intracellular content of each cell evolves according to the expression of R(x, t). Intrinsic noise is modelled using the Langevin SDE, Eq (15). The Gaussian white noise ξ(τ) is treated as a standard Brownian motion or a standard Wiener process and Eq (15) is solved using the Itô interpretation of the stochastic integral [57,58]. If the parameters of intrinsic noise are neglected, then the reaction rate of the intracellular content becomes deterministic (Eq (7)) and the stochastic algorithm accounts for effects of only extrinsic source of heterogeneity. At the end of time interval, T, the cell undergoing division is selected from the conditional distribution function [27]: where k is the index of the selected cell for division. Upon division of a cell, k, its content is distributed to its offsprings following a partition probability density function, P(x, x 0 ). In this work, we select the simple discrete partitioning mechanism formulated by Eq (18), which generates a daughter cell with content fx k , and a second daughter cell with content (1 − f)x k . Finally, and in order to maintain the number of simulated cells constant, we randomly pick a cell which is replaced by the second offspring generated during division. The dimensionless time is updated, i.e., τ ! τ + T, and the algorithm is repeated until a pre-specified time, t stop , is reached.
As stated above, the long time behaviour of an isogenic population carrying the lac operon genetic network is expected to feature bistable behaviour within a range of IPTG concentrations (* 1/ρ). Despite the fact that stochastic modelling provides a more realistic description of the population dynamics, it lacks efficiency when the systematic study of their long time behaviour is required, and especially in cases where multiple solutions can co-exist for the same IPTG concentration. In particular, exploring the solution space of (coarsely) time-invariant solutions over a wide range of parameter values, by executing solely temporal stochastic simulations is computationally demanding and renders the systematic study as a practically infeasible task. Furthermore, and since bistability regions are sought, when temporal simulations are performed at the vicinity of their intervals (critical turning points), stochastic noise can drive the population to unpredictable phenotypic switches.

The equation-free methodology
When deterministic descriptions are available, the bistability limits can be tracked by means of bifurcation analysis applying well-established pseudo arc-length parameter continuation techniques as performed in [19,41]. However, deterministic models are only simplified approximations failing into incorporating important information, e.g., originating from intrinsic noise effects.
More realistic descriptions are provided through stochastic simulations (e.g., the CNMC model), which however suffer from severe computational limitations; since we are interested in studying the long-time behaviour of cell populations over a wide range of parameter values, it is required to perform an extensive number of stochastic simulations, which renders the systematic study of the problem as a practically infeasible process. Alternatively, one can resort to multiscale computational techniques, such as the equation-free methodology, which utilises information originating from fine-scale (microscopic) level simulations, projects it to a coarse-macroscopic level, and enables the performance of numerical tasks, such as steady-state computations, stability and bifurcation analysis. The interchange of information between micro-and macro-scopic level is enabled through a computational structure, the coarse time-stepper [43,44,59,60], which reports the evolution of macroscopic variables of interest at distinct time instances and is schematically illustrated in Fig 3; a macroscopic variable is a statistical measure of the microscopically simulated cell population, e.g., the Cumulative Distribution Function (CDF) of the intracellular content, which is denoted with f(x, t). Through the lifting step, we generate a number, N copies , of random states which are consistent with the studied coarse variable (CDF distribution). Each of these microscopic states are simulated with the CNMC model for a short time period, T, and the average CDF of the updated microscopic states is computed through the restriction step. In effect, the steps described above construct a discrete time-mapping: where the operator G T is unknown. Below, we describe in detail the Restriction and Lifting procedures used in this study.
Restriction. If we denote with x i the content of the i th cell, then the CDF is computed by sorting into ascending order the x = {x i }, i = 1, . . ., N vector, and plotting it against the vector p = {p i = (i − 0.5)/N)}, i = 1, . . ., N. This constitutes a restriction of microscopic data, U, to the macroscopic variable, f through an operator, M, i.e., f = MU. For noise reduction purposes, many repetitions of the stochastic simulation are required and then the average restriction is computed.
Lifting. For the lifting step, we choose to work with the inverse CDF (ICDF), x(p), which produces the intracellular content x i = x(p i ) of a given cell, i. Assuming, that the ICDF is smooth enough, and given that it has a finite support, p 2 [0, 1], it can be approximated using a low-dimensional description, e.g., the first few of an appropriate sequence of orthogonal polynomials [60]. These polynomials are approximated by a series of j-degree orthogonal polynomials, ϕ j (p), which are monotone, non-decreasing functions and lie within the interval [0, 1]: The ϕ j polynomials are represented by their values on the point set p and the first four basis functions (q = 3) are sufficient for the description of the corresponding fine scale state. In cases where the ICDF is not smooth enough, a larger number of basis functions are required. The wherep ¼ 0:5 À arcsinð1 À 2pÞ=p. Each of the polynomial basis functions is computed on the point set p constructing a (q + 1) × N matrix, F. The coefficients α = α i are computed from: and the microscopic state x can be approximated from: Thus, the intracellular content of each cell of the population can be constructed by a small set of α values. The lifting procedure described above constructs a microscopic description U from a macroscopic variable f through an operator, μ, i.e., U = μf. It should be noted here that the choice of lifting and restriction operators needs to guarantee that lifting from coarse level to fine scale and then restricting to coarse level again has no effect, i.e., μM % I.
Healing Time. The naturally arisen question is whether the level of description using 4 basis functions is sufficient enough, and higher order approximation is required. To test this, we perform the following computational experiment: A CNMC simulation is interrupted at time τ = τ inter and the microscopic state x original (τ inter ) is obtained. Then, we perform the lifting step using 4 basis functions and construct a lifted microscopic state x lifted (τ inter ). The next step is to perform two different CNMC simulations for a time horizon T hor using as initial microscopic states: (a) x original (τ inter ) and (b) the lifted x lift (τ inter ). Finally, we compute the evolution of the α values computed using 5 = 4+1 basis functions for the two different simulations at distinct reporting time instances, and compare their relative error, (α lift − α original )/α original .
In Fig 4(a) we show an example of this test computation, where a CNMC simulation is interrupted at dimensionless time τ = 10 to obtain the x original (10) and the lifted distribution x lifted (10) using four α coefficients (q = 3). Then the original and lifted distributions are simulated over a time interval τ 2 [10,11], and at distinct time instances with an interval of Δτ = 0.1, we compute the five first α coefficients using a 4 th order approximation. As expected, the 5 th α coefficient of the lifted distribution is significantly different from the respective α coefficient of the original distribution; however, it takes only a small time interval of Δτ = 0.1 to see that the 5 th α coefficient of both the original and lifted distributions have no difference (their relative difference drops below the stochastic noise computed from 2s a j =a original j , where σ α j is the standard deviation of the noise of coefficients calculated using a number (here 50) of direct CNMC copies [61]). The time required for higher order coefficients of lifted distributions to converge to the ones obtained from original distributions is referred to as healing time [42,43], and it is a preparation time interval for the coarse-time stepper during which the dynamics of fast variables (higher order coefficients or moments of the distribution) equilibrate quickly and get slaved by the evolution of the lower order coefficients. The first four coefficients show no difference within the entire time interval τ 2 [10,11].
We also illustrate that using a low order approximation does not affect the shape of the distribution of cells as a function of their intracellular content. In particular, in Fig 4(b) we show the original and the lifted with up to 3 rd order polynomial approximation ICDF distributions.
Clearly, lifting using a 3 rd order approximation (or 4 orthogonal basis functions) is accurate enough for a coarse-level description of the cell population.
Coarse steady state computations and parametric analysis. As reported above, the lifting and restriction steps constitute the coarse time-stepper, which maps the dynamics of the coarse variable f at discrete time instances according to the general expression of Eq (24), or if we study the coarse variable α according to the general expression: with the operator G τ being unknown. This black-box simulator can be utilised for the performance of numerical tasks, such as steady state computations with the Newton-Raphson method. A steady-state solution, α Ã satisfies: If we are interested in applying the Newton-Raphson method for the computation of the coarse steady-state solution α Ã , it is required to solve the following set of non-linear equations: In order to solve the non-linear system of equations with α Ã being the unknowns we apply the Newton-Raphson method, which requires during each iteration the solution of the linearised system: Since the explicit expression of G T is not available, its approximation is performed numerically using a simple forward finite difference scheme: where is a small perturbation number.
In Fig 5, we compare the number density function, n([Y] as a function of the intracellular content, [Y] as computed from (a) a direct temporal simulation of the full CNMC model performed until time, t % 60 hrs, and (b) the equation-free based, Newton-Raphson computation described above. In particular, we compare the steady-state distributions for two cases of cell populations featuring low and high level lac permease concentration levels, showing very good agreement in both cases. It is evidently clear that the equation-free approach, even though uses a low-level description, it produces number cell density functions which compare remarkable well with the time consuming, full simulations of the CNMC model.
Using the same idea, one can perform bifurcation analysis by means of pseudo arc-length parameter continuation enabling the computation of the entire solution space, including stable and unstable steady state solutions. As a by-product of this process, the bistability range is accurately computed. Below, we present the results of this analysis quantifying the effect of intrinsic noise effect on heterogenous cell populations, by choosing as continuation parameter the external IPTG concentration, [I ex ].

Effect of intrinsic noise
In order to decompose the effect of the different sources of heterogeneity, we perform the following steps. First, we solve the steady state solution of the homogeneous model (Eq (19)), which neglects all sources of heterogeneity and considers populations, where each individual cell carries the same intracellular amount. The extrinsic source of heterogeneity is then incorporated by solving the DCPB model (Eq (16)); for details of the numerical solution of the DCPB model, we refer the reader to Kavousanakis et al [41]. Finally, we employ the equationfree method to perform coarse bifurcation analysis of the stochastic model (CNMC) for different values of parameters K and y Ã , which quantify the effect of intrinsic noise.
In Fig 6(a), we present a bifurcation diagram showing the dependence of the average lacY expression of a population carrying the lac operon genetic network on the extracellular inducer concentration, [I ex ] for four different cases: (a) the homogeneous model, (b) the DCPB model solved with the finite element method [41], (c) the CNMC model neglecting intrinsic noise effects (K, y Ã ! 1), and (d) the CNMC model incorporating the intrinsic source of heterogeneity. The single-cell reaction rate is given from Eq (7), with π = 0.03 and κ = 0.05. The partitioning mechanism is discrete (see Eq (18)) and symmetric (f = 0.5); the division rate is given from Eq (17) with m = 2. Finally, CNMC models simulate populations consisting of N = 10,000 cells; we use 50 copies of CNMC simulations for each parameter value for stochastic noise reduction purposes.
In all cases, one can observe an S-shaped bifurcation diagram with two stable steady state solution branches (solid lines) and a branch of unstable solutions (dashed lines). In Fig 6(b), we depict the average number density distribution of 50 copies of CNMC simulations n([Y]), When the population resides in an environment, which is rich in the inducer IPTG, then the average expression level of lacY is high. By decreasing the external IPTG concentration, there exists a critical value (left turning point) which signals the abrupt transition towards low expression level of lacY. In a reverse experiment, where we start with low IPTG concentration values, the average expression of lacY is low, and by increasing IPTG there exists a critical value (right turning point) beyond which the average expression will jump towards higher values. The two turning points are distinct suggesting that transitions between populations featuring high and low lacY expression levels, through modification of the IPTG concentration are hysteretic.
The homogeneous model results show that the bistability region spans over a wide range of [I ex ] values: [I ex ] 2 [17.1,26.9]μM. Extrinsic heterogeneity shifts the bistability range towards higher [I ex ] values (DCPB and CNMC with K, y Ã ! 1). Here, the results obtained from DCPB and CNMC show very good agreement, since the size of the cell population simulated by the CNMC model is large (10,000). The range of solution multiplicity is shifted further towards higher [I ex ] values, by taking into account the effect of intrinsic noise (CNMC model with finite K and y Ã values). In particular, we report that the lower limit of the bistability region as computed from the CNMC model for K = 500 and y Ã = 50 is located at [I ex ] = 24.1μM compared to the value of [I ex ] = 22.5μM, which is computed when neglecting the intrinsic source of heterogeneity. In addition, the upper [I ex ] limit of the bistability region is shifted towards higher values for the CNMC model with intrinsic noise ([I ex ] = 35.7μM) compared to [I ex ] = 29.6μM, for the CNMC model with K, y Ã ! 1. Thus, the combined effect of extrinsic and intrinsic heterogeneity shifts the bistability region towards higher IPTG concentration values. It should be noted here that similar findings have been presented in the experimental work of Maeda and Sano [62] reporting hysteric transitions within IPTG concentration values of 2.5 − 50μM (for certain wild type promoters), as well as in Matsumoto et al. [63] reporting bistability in the range of 5-15μM.
In order to characterise a steady state solution as stable or unstable, we perform coarse stability analysis by determining the eigenvalues of the WG t ða Ã Þ Wa matrix (see Eq (32)). We present an indicative case for steady state solutions obtained from the model incorporating intrinsic heterogeneity (lines with open triangles in Fig 6(a)), and in particular for parameter value [I ex ] = 28.8μM. The spectra of eigenvalues presented in Fig 7(a)-7(c) correspond to the upper, intermediate, and lower branch, respectively. In Fig 7(a) and 7(c), the eigenvalues of the matrix, WG t ða Ã Þ Wa , lie within the limits of the complex plane unit circle, and the two solutions are dynamically stable. On the contrary, in Fig 7(b), one eigenvalue crosses the complex plane unit circle, and the corresponding steady state solution is characterised as dynamically unstable.

Effect of operator fluctuations and reference number of molecules
We now investigate the effect of the intrinsic heterogeneity "intensity", which is quantified by the K and y Ã values. In particular, the effect of intrinsic noise is strengthened by lowering the operator fluctuations (K) and the reference number of regulatory molecules (y Ã ). We demonstrate this effect in Fig 8, where the dependence of the steady state expression level of the average intracellular content, h[Y]i, on the extracellular IPTG concentration, [I ex ], is presented for different K and y Ã values. In Fig 8(a), we illustrate this effect by lowering the parameter K and keeping the parameter y Ã constant. By lowering the K value, the intrinsic noise effect is intensified and the bistability region is shifted towards higher, [I ex ], values.
In particular, the bistability range for K = 1000 lies within the inerval [I ex ] 2 [23.3,33.1]μM, whereas the corresponding interval for K = 250 starts from [I ex ] = 24.8μM and ranges up to unrealistically high (tending to infinity) values (the dimensionless ρ tends to 0). A similar trend is observed by lowering the y Ã parameter value with K fixed (see Fig 8(b)), however the degree of change is not as significant as in the variable, K, case. In this case, the bistability range for y Ã = 500 spans over the interval [  In the cases presented above, one of the two intrinsic noise parameters is kept constant and the intrinsic heterogeneity effect is strengthened by lowering the value of the other parameter. The two parameters of intrinsic source of heterogeneity act in a collaborative manner and strengthen further the effect of intrinsic noise by decreasing both K and y Ã values. In Fig 9, we compare the results obtained from the DCPB model with the CNMC model for: (a) K = 1000, y Ã = 500, (b) K = 500, y Ã = 50 and (c) K = 250, y Ã = 25. The CNMC model with the largest K and y Ã values resembles best the behaviour of the DCPB model, whereas low K and y Ã values shift the bistability limits towards higher [I ex ] values. Interestingly enough, when the effect of intrinsic noise is sufficiently intensified, the right turning point reaches infinity (ρ ! 0) and does not correspond to IPTG values with physical meaning suggesting that transitions between the high and low level expression states are non reversible: i.e., a high to low lacY expression level transition is feasible by decreasing the IPTG concentration; however the reverse transition becomes infeasible considering K = 250 and y Ã = 25, since the upper end of the bistability region is located at [I ex ] ! 1 values.

Effect of cell division sharpness and asymmetric partitioning
In addition to intrinsic noise parameters, we also examine the effect of cell division sharpness, quantified by parameter, m, and asymmetric partitioning, which is quantified by parameter, f. The effect of asymmetric partitioning is illustrated in Fig 10(a)   The impact of asymmetric partitioning becomes more significant when combined with intrinsic heterogeneity as shown in Fig 10(b) for f = 0.3. In particular, the CNMC model neglecting intrinsic noise effects (K, y Ã ! 1) agrees well with the corresponding DCPB model, whereas the CNMC model with K = 500 and y Ã = 50 predicts a shifted towards higher [I ex ] values bistability interval.
Furthermore, the effect of division rate sharpness for a population of N = 10,000 is presented in Fig 11(a). Higher single-cell division rates (larger m values) shift its upper end towards higher [I ex ] values. In particular, when m = 3 the upper end tends to infinity suggesting that transitions between high and low lacY expression levels are irreversible.   In Fig 11(b), we compare the CNMC model for m = 3 and K = 500, y Ã = 50 with the corresponding DCPB model and the stochastic CNMC model, which neglects intrinsic noise effects (K, y Ã ! 1). By incorporating the intrinsic source of heterogeneity, the effect of asymmetric partitioning becomes more intense leading to large discrepancies, compared to the ones between the DCPB model and the CNMC model, which neglects the intrinsic source of heterogeneity. In particular, the bistability regions for the CNMC with K, y Ã ! 1, and the DCPB model are [ As reported above, noise can induce rapid changes of the average phenotype of large cell populations, when the extracellular conditions are at the vicinity of critical turning point values. In Fig 12, we present a single copy simulation of a population at [I ex ] = 30.8μM, which In all cases presented above, the bistability region is shifted towards higher [I ex ] values by intensifying the effect of intrinsic noise. If we consider slower dynamics for the single cell division rate, then we observe a reverse effect. In Fig 13,   (m = 1) a low to high lacY expression level transition is expected to occur at a lower IPTG concentration value, when intrinsic noise effects are taken into consideration, whereas for higher cell-division rates intrinsic noise has always an opposite effect (low to high transitions occur at larger IPTG concentration values).

Discussion
We present a multiscale computational methodology, which enables the systems-level study of cell populations simulated by means of stochastic models. In particular, a CNMC model is employed to simulate the dynamics of heterogeneous cell populations which carry the lac operon genetic network, featuring solution multiplicity over a range of extracellular inducer concentration values. In this work, we decompose the effect of different sources of heterogeneity (extrinsic and intrinsic) and study their effect on the bistability range of IPTG values. The extrinsic source of heterogeneity has been shown to have a significant effect on the phenotype of cell populations when performing deterministic cell population balance model computations [41]. However, the effect of intrinsic heterogeneity cannot be described by deterministic modelling since it involves the computation of Langevin stochastic differential equations for the intracellular reaction network. On the other hand, stochastic simulations which incorporate intrinsic noise effects cannot be used for a systems-level study of the problem; the number of stochastic simulations, which are required to compute accurately the range of bistability can be massive, rendering this approach as a computationally infeasible one.
To bypass this impediment, we employ the equation-free methodology, which utilises finescale information and projects it to a coarse-macroscopic level, for which well established numerical algorithms can be applied. Here, we perform bifurcation analysis, using pseudo arclength continuation techniques in order to explore the dependence of the solution space on the extracellular inducer (IPTG) concentration, and accurately determine the range of bistability. Our stochastic CNMC-based computations are validated against deterministic descriptions, using large populations of cells (N = 10,000) and by neglecting the intrinsic noise effects (K, y Ã ! 1). Then, we explore the effect of intrinsic source of heterogeneity and demonstrate that when strengthened the range of bistability shrinks and shifts towards higher IPTG concentration values for sufficiently high single-cell division rates. When the effect of intrinsic noise is sufficiently strengthened, the turning point signalling the transition from low to high lacY expression levels tends to infinity, suggesting that the transition between high and low expression levels is irreversible, through modification of IPTG concentration values. A similar trend is also observed for populations dividing in a more asymmetric fashion. In addition, we show that the existence of intrinsic heterogeneity can lead to non-trivial dynamical behaviours such as phenotypic switching between co-existing stable steady states. We present such a transition at the neighbourhood of the bistability upper limit, when a sharp division rate is applied. Finally, we demonstrate the disappearance of the bistability interval for populations with low single-cell division rates and asymmetric partitioning mechanism, when intrinsic heterogeneity is also incorporated. In this case, the transition from low to high lacY expression levels occurs at a lower IPTG value compared to the case, where the intrinsic source of heterogeneity is neglected.
We also report that the equation-free framework is quite flexible to adopt for the performance of coarse-grained analysis of cell populations with genetic networks of increased complexity and different architectures. Indicatively, we report the cases of the lac operon genetic network with a promoter containing three repressor binding sites with cooperative interaction among them [64] and the genetic toggle switch [65]. The study of these networks requires new formulations for the single-cell reaction rate, the division rate and the partitioning mechanism; however, the application of the multiscale equation-free methodology does not require any modification and can be readily applied for the performance of their efficient systems level study.