Stochastic models coupling gene expression and partitioning in cell division in Escherichia coli

BioSystems 193-194 (2020) 104154 Available online 28 April 2020 0303-2647/© 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Stochastic models coupling gene expression and partitioning in cell division in Escherichia coli Ines S.C. Baptista, Andre S. Ribeiro * Laboratory of Biosystem Dynamics, Faculty of Medicine and Health Technology, Tampere University, 33014, Tampere, Finland ARTICLE INFO Keywords: Stochastic models Single gene expression Genetic circuits Partitioning in cell division Escherichia coli ABSTRACT Regulation of future RNA and protein numbers is a key process by which cells continuously best fit the environment. In bacteria, RNA and proteins exist in small numbers and their regulatory processes are stochastic. Consequently, there is cell-to-cell variability in these numbers, even between sister cells. Traditionally, the two most studied sources of this variability are gene expression and RNA and protein degradation, with evidence suggesting that the latter is subject to little regulation, when compared to the former. However, time-lapse microscopy and single molecule fluorescent tagging have produced evidence that cell division can also be a significant source of variability due to asymmetries in the partitioning of RNA and proteins. Relevantly, the impact of this noise differs from noise in production and degradation since, unlike these, it is not continuous. Rather, it occurs at specific time points, at which moment it can introduce major fluctuations. Several models have now been proposed that integrate noise from cell division, in addition to noise in gene expression, to mimic the dynamics of RNA and protein numbers of cell lineages. This is expected to be particularly relevant in genetic circuits, where significant fluctuations in one component protein, at specific time moments, are expected to perturb near-equilibrium states of the circuits, which can have long-lasting consequences. Here we review stochastic models coupling these processes in Escherichia coli, from single genes to small circuits. 1. Introduction Statistical fluctuations in biochemical processes have long been considered to explain rare events or mechanisms (Delbrück, 1940). For example, stochastic models were shown to best explain the single-cell distributions of numbers of virions in infected bacteria, where cell-to-cell variability was observable (Delbrück, 1945). However, noise in RNA and protein numbers was not usually considered until recently, since observations of mean gene expression levels of large cell populations (usually highly expressing genes, such as the Lac gene (Lutz et al., 2001)) using technologies such as qPCR, plate reading, Western blot, abortive initiation reactions (McClure, 1980), etc., were well explained by continuous models (Ackers et al., 1982). Subsequently, with the introduction of fluorescent proteins (Shimomura et al., 1962), single-cell fluorescence microscopy (Sanderson et al., 2014), synthetic genetic constructs (Lee et al., 2016), and flow cytometry (Adan et al., 2017) allowed observing significant cell-to-cell variability in protein numbers, even in isogenic cell populations in homogenous environments. In bacteria, one important cause for this variability are the small numbers of RNA in optimal growth conditions, usually less than 10 per gene (Taniguchi et al., 2010). Noise in prokaryotic gene expression is considered to be a key source of bacterial phenotypic heterogeneity. This heterogeneity is expected to promote survival of cell populations and cell lineages in fluctuating environments (Acar et al., 2008; Charlebois and Bal� azsi, 2016; Healey et al., 2016) by enhancing adaptability via specialization of biological functionalities (Ackermann, 2015), by allowing phenotype selection (Süel et al., 2006) and by providing negative frequency-dependent selection, i.e. the emergence of rare, lesser fit phenotypes in, e.g. environments with diverse resources (Charlebois and Bal� azsi, 2016; Healey et al., 2016). One of the earliest observations of a stochastic decision-making mechanism of gene expression occurred in a study of the immunity phase-shift in cells lysogenic for λCI857, conducted by Neubauer and Calef (1970). Two genes, CI and Cro, form a toggle switch by repressing one another. This circuit is thus expected to be in one of two states: either CI is ‘ON’ and Cro is ‘OFF’, or vice versa. Once reaching one of these states, the circuit is expected to remain stable thereafter. However, * Corresponding author. Arvo Ylp€ on katu 34, P.O.Box 100, 33014, Tampere University, Finland. E-mail address: [email protected] (A.S. Ribeiro). Contents lists available at ScienceDirect BioSystems journal homepage: http://www.elsevier.com/locate/biosystems https://doi.org/10.1016/j.biosystems.2020.104154 Received 10 September 2019; Received in revised form 3 April 2020; Accepted 16 April 2020 BioSystems 193-194 (2020) 104154 2 at rare moments, the circuit was observed to switch to the other state, which is consistent with stochastic dynamics. To show this, Arkin et al. (1998) design a kinetic model of the circuit using the stochastic formulation of chemical kinetics (Gillespie, 1977, 1992) and found that, unlike deterministic models, the stochastic model predicted the statistics of the circuit. Since then, several stochastic models of gene expression have been proposed for various organisms, particularly for the model organism E. coli. In E. coli, and similar organisms in general, little regulation appears to be exerted on RNA degradation, in that this process seems to be largely independent from the RNA sequence (Bernstein et al., 2002; Deutscher, 2006; Chen et al., 2015). In detail, while there is some evidence for regulation of bulk RNA concentrations as a function of cell growth rates (Esquerr� e et al., 2014), so far we lack evidence for RNA sequence dependent regulation, unlike in transcription initiation, where several regulatory mechanisms are encoded at the promoter region of each gene (e.g., binding sites for transcription factors). Thus, most stochastic models of gene expression set the key regulatory mechanisms at the stage of transcription initiation (for reviews see, e.g. (Gibson and Mjolsness, 2001; de Jong, 2002; Ribeiro, 2010)). As time-lapse microscopy and single-molecule tracking became common practices, cell-to-cell diversity was observed to also emerge in cell division due to, e.g., morphological asymmetries in division (for a review see (Kysela et al., 2013)), or asymmetries in the partitioning of components between sister cells, e.g., due to non-random, protein spatial distributions (Llopis et al., 2010). Interestingly, these asymmetries in cell division can be affected by the partitioning scheme of the components being partitioned (Huh and Paulsson, 2011), and by the presence of macromolecules (including the nucleoid), which can enhance/reduce heterogeneities in the spatial distribution of the components (Gupta et al., 2014a). Further, asymmetries in the location of the plane of division (Gupta et al., 2014b) can also enhance asymmetries in cell division. Meanwhile, at a higher scale, if the components asymmetrically partitioned are integrated into a genetic circuit, the asymmetry in one component may propagate to other proteins of the circuit (Munsky et al., 2012). Here we review recently proposed stochastic models of the combined effects of noise in gene expression and stochastic partitioning of RNA and proteins in cell division, in the context of bacterial cell populations. Given the myriad of sources of cell-to-cell diversity in gene expression components, it is not simple to identify the effect of each source when analyzing experimental data. Consequently, their modelling can also be difficult. Therefore, models usually account for the most common sources of noise (e.g. in gene expression) and then add a specific source studied (for reviews, see (Hellweger et al., 2016; Gupta and Mendes, 2018)). As such, we start by reviewing stochastic models of gene expression and of partitioning schemes in cell division, separately. Next, we review models that include both of these sources of noise. Finally, we describe freely available simulators for these models that use the common language of representation of chemical reactions in order to suit users with various scientific backgrounds. For reviews of models of gene expression, models of genetic circuits, and a comparison of delayed versus non-delayed stochastic models of gene expression, we refer to, respectively, (de Jong, 2002; Karlebach and Shamir, 2008; Hasty et al., 2001), and (Ribeiro, 2010). 2. Stochastic models of gene expression Several models of stochastic gene expression have been proposed, usually based on observations of viral and bacterial gene expression systems. One of the first stochastic models was proposed by Ko et al. (Ko, 1991, 1992), from empirical evidence that the expression dynamics of individual genes did not match their mean behavior (Ko et al., 1990). Several models followed (Arkin et al., 1998; Mcadams and Arkin, 1997; Gibson and Bruck, 2000; Sasai and Wolynes, 2003; Ozbudak et al., 2002). For reviews, see (Karlebach and Shamir, 2008; Gibson and Mjolsness, 2001). Later on, to consider the complex, multi-step nature and/or significant length in time of gene expression, delayed stochastic models were introduced (Bratsun et al., 2005; Barrio et al., 2006; Roussel and Zhu, 2006; Ribeiro et al., 2006; Gaffney and Monk, 2006). Usually, the delays were introduced to account for the time-length of protein folding and activation, and/or of transcription initiation. We start by describing a stochastic gene expression model with neither time delays nor regulatory processes. For this, we assume a constitutive gene, i.e. always active (in a ‘ON’ state, represented as ‘Gene ON ’). We model the production of RNA and proteins separately, both as 1-step processes (reactions (2.1) and (2.2), respectively). This separation of gene expression into transcription and translation allows for the separate regulation of the kinetics of the two events. Also included are separate, constant exponential decays of RNA and proteins (reactions (2.3) and (2.4), respectively) (Peccoud and Ycart, 1995; Munsky et al., 2012): GeneON ��! k1RNA (2.1) RNA��! k2P(2.2) RNA��! k3∅(2.3) P��! k4∅(2.4) Here, k 1 and k 2 are the RNA and protein (P) production rates, respectively, while k 3 and k 4 are the RNA and protein degradation rates by 1step decay processes, respectively. Next, one can introduce a regulatory process (Jacob and Monod, 1961). The most common is repression by transcription blocking, i.e. a molecule binds near or at the promoter region, blocking access to the RNA polymerase (RNAP). In this model, a gene can either be ‘ON’, with transcription occurring at a constant rate, or ‘OFF’ (Gene OFF ), if repressed. Transitions between the ON and OFF states (at exponentially distributed intervals (Gardiner, 2004)) occur at the rates k ON and k OFF (Kepler and Elston, 2001), as follows: GeneOFF ⇄ kON kOFF GeneON (2.5) Usually, the two events depicted in reactions (2.5) require the intervention of one or more regulatory molecules (e.g. transcription factors.). As examples, in reaction (2.6), we model a gene whose activation requires the binding of molecule ‘A’, while in reaction (2.7) we model a gene whose repression requires the binding of molecule ‘B’ (similarly to above, the explicit representations of A and B allows the control of the numbers of these molecules in time and/or in individual cells): GeneOFF þA⇄ kON kOFF GeneON:A(2.6) GeneOFF:B⇄ kON kOFF GeneON þB(2.7) By designing more complex models than those above, it is possible to enhance their realism. For example, in accordance to empirical data, and unlike what is represented in reaction (2.1), transcription initiation is a multi-step process (McClure, 1985). This is modeled in reactions (2.8) where, first, the RNAP has to find a transcription start site (TSS), usually located in the promoter region of the gene. The binding between the RNAP and the TSS that forms a ‘closed complex’ is reversible (McClure, 1980). In detail, in the first step of the set of reactions (2.8), k cc is the rate at which the RNAP finds the promoter region and performs a 1D diffusion process along the DNA until successfully binding to the TSS (Bai et al., 2006; Wang and Greene, 2011), while k -cc is the opposite event. Due to being reversible, in general, the RNAP forms several closed I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 3 complexes until one of them successfully progresses into a stable, open complex, following the melting of the DNA and formation of a transcription ‘bubble’, at the rate k oc (Bai et al., 2006). The kinetics of this process (e.g. the number of closed complex formations per open complex formation) can be tuned by changing k cc and k oc . The remaining steps (here regulated by k escape ) represent the overall rate of promoter escape and clearance (after completion of the open complex), accounting for the rate of abortive initiation (Duchi et al., 2016; Liang et al., 1999). In normal conditions, ‘promoter escape’ is much faster than the two previous rate-limiting steps (k escape >> k cc and k oc ) (McClure, 1985; Duchi et al., 2016). Subsequently, as the RNAP clears the promoter region, it forms an elongation complex that will synthesize the RNA. Overall, the multi-step nature of transcription initiation can be well modeled by reactions (2.8), with Gene ON being an active, unrepressed gene, R being the RNA polymerase, Pro CC the closed complex, and Pro OC the open complex: GeneON þR⇄ kcc k cc ProCC ��! koc ProOC ��! kescape GeneON þRþRNA (2.8) In this model, the rate of transcription elongation is not accounted for in the kinetics of RNA production. This is because it does not affect the mean time length between consecutive RNA production events, only the variability of this time length. The modeling of the multi-step nature of transcription initiation (Walter et al., 1967; Kierzek et al., 2001; Saecker et al., 2011; Chen et al., 2019, 2020) can be of importance, particularly the first two rate-limiting steps, as these can be significantly long-lasting. Also, by tuning k cc and k oc individually, one can alter significantly the shape of the distribution of time intervals between RNA production events, even when not changing the mean rate of RNA production (Startceva et al., 2019), which is not possible if using a one-step model of RNA production without time delays (Lloyd-Price et al., 2016). Second, regulatory mechanisms in E. coli may act on only one of these steps, or in both by different degrees (Lutz et al., 2001; M€ akel€ a et al., 2017), thus requiring their direct representation. As an example of how differences in the rates of the two main rate limiting steps of transcription initiation can affect the shape of the distribution of time intervals between consecutive RNA production events and, thus, of noise in RNA production, we performed example simulations of the model (2.8), setting different values for k cc and k oc , while maintaining (1/(R �k cc )þ1/k oc ) constant. We also consider RNA degradation (reaction (2.3)), to ensure realistic RNA numbers (as well as the contribution from noise of the degradation process). Models, parameter values of the rate constants, and initial values of each component are shown in Fig. 1. We consider the following conditions. In (A), the formation of the closed complex is faster than the formation of the open complex (k cc > k oc ). In (B), the rates of formation of the closed and open complex are identical (k cc ¼k oc ). Finally, in (C) the formation of the closed complex is slower (k cc <k oc ). The expected time between transcription events is the same in all conditions. To compare the kinetics of RNA production in the three conditions, we performed simulations and extracted the number of RNAs in a period of time when the RNA numbers are expected to be near equilibrium (for Fig. 1. Example simulations of a model of a 2-rate limiting step process of RNA production by an active gene and RNA degradation. Shown are the reactions of RNA production and RNA degradation (top left table, reactions 1 and 2, respectively). Also shown are the initial quantities of each component (top right table), and parameter values different and identical to each condition (mid tables, left and right, respectively). Three conditions are considered: (A) the rate of closed complex formation is higher than the rate of open complex formation (k cc >k oc ); (B) k cc ¼k oc ; and (C), k cc <k oc . In all conditions, the mean expected RNA production rate (1/(R �k cc )þ1/ k oc ) is identical. Also identical are k escape and k d . Finally, the bottom figures show the probability distributions of time intervals (Δt) between RNA production events (right figures), and the number of RNAs in individual cells at any given moment (left figures). Each figure also shows the mean ( μ ) and standard deviation ( σ ) of the distribution. Note that the axes ranges differ. We performed 2200 simulations of each model, each 10 6 s long, and extracted the number of RNAs once per second, between moments 2 �10 5 s to 1 �10 6 s, during which time the RNA numbers are expected to be near equilibrium. I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 4 details, see legend of Fig. 1). We also obtained the time intervals (Δt) between consecutive RNA production events. For this, we modeled probes (not shown in reaction 2.8) that appear at the same time as an RNA is produced, but do not degrade. As such, the intervals can be obtained from when a new probe appears (avoiding interference from RNA degradation). The distributions of RNA numbers and Δt values of each condition are shown in Fig. 1, along with their respective mean ( μ ) and standard deviation ( σ ). As expected, since (1/(R �k cc )þ1/k oc ) does not differ between conditions, both μ (Δt) and μ (no. RNAs) do not differ. However, one finds that the standard deviation ( σ ) is minimal when k cc and k oc are identical (condition B), while it does not differ significantly between conditions A and C since the absolute difference between k cc and k oc is identical in the two conditions. To validate these conclusions, we performed two-sample t-tests between the pairs of distributions of extracted values of σ of the distributions of Δt values, with the null hypothesis that the data comes from normal distributions with equal means and unequal variances. We found that the p-values were approximately 0 when comparing A and B and, B and C (for p-values smaller than 0.01, we reject the null hypothesis). Meanwhile, comparing A and C, we obtained a p-value of 0.13 and, thus we cannot reject the null hypothesis. Similar conclusions were obtained when comparing the distributions of RNA numbers. Overall, using the two-step model (2.8), it is possible to regulate the mean and variability in RNA numbers independently from one another (to a certain degree), similarly to models with an ON-OFF mechanism whose kinetics is independent from the transcription process (e.g. a model combining reactions 2.1 and 2.5). Rather than model the events in transcription initiation explicitly, one can instead consider how long they take to be completed once initiated. This may result in more accurate dynamics of RNA numbers. For example, simulations of reaction (2.8) usually assume that the time taken by each step follows an exponential distribution (particularly when simulated using the Stochastic Simulation Algorithm (Gillespie, 1977)). However, this is not necessarily the case (e.g. the opening of the DNA for reading is not expected to have this kinetics). To correct for this, one can introduce time delays ( τ ) between the start and completion of an event as in (Ribeiro et al., 2006). The delays can be constant or, to be more accurate, be randomly extracted from a desired distribution each time a reaction occurs (the simulator proposed in (Lloyd-Price et al., 2012) can be used for this purpose). For this, transcription can be modeled as: RþPro��! k1Proð τ 1Þ þ RNAð τ 2Þ þ Rð τ 2Þ(2.9) In (2.9), the RNA polymerase (R) binds to the promoter (Pro) as previously but, instead of releasing the products once the reaction occurred, it does so only after the time delays (after τ 1 seconds in the case of the promoter region, and after τ 2 , with τ 2 > τ 1 , in the case of the RNA and RNAP, to account for transcription elongation and termination). One advantage of this model is that the delays can be directly extracted from empirical distributions each time the reaction occurs, to obtain accurate distributions of intervals between RNA production events. To exemplify the ability of this modeling strategy to tune single-cell RNA numbers, consider the model (2.9) and three conditions, differing in the distribution from which τ 1 values are randomly drawn from, each time a transcription event occurs (see Fig. 2 for details). For simplicity, we consider three conditions (A, B, and C), with τ 1 values being drawn from Gaussian distributions differing in mean and standard deviation (represented as Gaussian ( μ , σ 2 )). The conditions also differ in the rates of RNA degradation (k d ), which are tuned so that the conditions differ little in mean expected RNA numbers. Example results in Fig. 2 show that the distribution of time intervals between the production of consecutive RNAs and the distribution of Fig. 2. Example models of a 2-rate limiting step process of RNA production from an active gene with time delays in the release of the promoter (Pro), the RNA, and the RNAP (R). Shown are the reactions of RNA production along with RNA degradation (top left table), the initial quantities of each component (top right table), as well as the general parameter values and the parameter values specific to each condition (mid tables). Time delays are randomly drawn from Gaussian distributions each time an event occurs. Three conditions are considered: (A), (B), and (C), which differ in the mean and standard deviation of the Gaussian distribution from which time delays are drawn from. Finally, the bottom figures show the probability distributions of time intervals (Δt) between RNA production events and of the number of RNAs in individual cells at any given moment (once reaching near equilibrium). These figures also show the mean ( μ ) and standard deviation ( σ ) of the distributions. Note that the axes ranges differ. We performed 1200 simulations of each model, each 10 5 s long, and extracted the number of RNAs once per second, between moments 2 �10 4 and 1 �10 5 , during which time the RNA numbers are expected to be near equilibrium. I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 5 RNA numbers can be tuned by the time delays in transcription initiation, even when tuning, e.g. RNA degradation rates, so that mean RNA numbers are not affected. To validate these conclusions (by establishing that the dynamics differs significantly between the conditions), we performed two-sample t-tests between the pairs of distributions of extracted values of σ of the distributions of single-cell RNA numbers values, with the null hypothesis that the data comes from normal distributions with equal means and unequal variances. We found that the pvalues were approximately 0 for all possible pairs of distributions (for pvalues smaller than 0.01, we reject the null hypothesis). Similarly, using this methodology, translation can be modeled as: Ribosome þRBS��! k2RBSð τ 3Þ þ Ribosomeð τ 4Þ þ Pð τ 5Þ(2.10) where RBS stands for the ribosome binding site region of the RNA since, in E. coli, once this region is free, ribosomes can bind to initiate translation of the RNA (which allows for several proteins to be produced simultaneously from the same RNA). Meanwhile, τ 5 stands for the time that it takes for protein folding and activation. Finally, proteins can be degraded as in reaction (2.4) (its noted that this model assumes a process of protein degradation as simple as the one for RNA, which may not suffice in realism). These models can be expanded to include several other phenomena involved in gene expression. As an example, we describe how one can introduce the influence of positive supercoiling buildup, known to affect highly expressed genes (see e.g. (Chong et al., 2014)). A set of reactions similar to reactions (2.5) can be used to model transcription locking due to positive supercoiling buildup (PSB), which causes transcriptional bursting (Golding et al., 2005; Chong et al., 2014) (Fig. 3). Namely, transcription locking resulting from the accumulation of PSB from the activity of neighboring genes can be modeled by reaction (2.11) (Chong et al., 2014). Escape from this state requires the intervention of Gyrase (Reece and Maxwell, 1991; Drlica, 1992). This is modeled by reaction (2.12), whose kinetics depends on the rates of Gyrase association and dissociation from the DNA (Reece and Maxwell, 1991) and the number of Gyrases in the cell. In this regard, one could also add a time-delay, to account for the time taken by Gyrase to resolve sufficient coils to release the promoter from a locked state (usually this is much faster than the time scale of intervals between transcription events and, thus, we opted for not including it). ProON ��! kPSB ProLocked (2.11) ProLocked þG��! kGProON þG(2.12) Promoter locking/unlocking due to PSB, represented in (2.11) and (2.12), can be coupled with the models of transcription. E.g., it can be coupled with the model in (2.8), by setting that the promoter (Pro), when in the ON state, being able to either be bound by RNAP (via reaction (2.8)) or to become locked due to supercoiling (via reaction (2.11)). As such, one could set a model where there is competition for Pro ON , each time it is available. A simpler version of this model (transcription modeled as a 1-step process) was used in (Chong et al., 2014; Golding et al., 2005). Finally, an additional reaction (e.g. reaction (2.3)) would allow for RNA to degrade (which is not necessary, if the model only intends to predict intervals between consecutive RNA production events, rather than RNA numbers). As an example, Fig. 4 shows the time intervals between RNA production events of such a model, absent of RNA degradation and of regulation by transcription factors, for simplicity (reactions (2.6) and (2.7) could be added for implementing activation/repression mechanisms, respectively). Also shown are the reactions simulated, their parameter values, and initial numbers of each component. From the results, the variability of the time intervals between consecutive RNA production events increases as Gyrase numbers decrease, even when the mean of the interval between consecutive RNA production events is kept near constant (here, this was achieved by compensating increasing OFF periods due to positive supercoiling buildup with increased transcription initiation rates). In summary, stochastic modelling of gene expression can be made increasingly complex, depending on which processes are modeled and to which level of detail. The level of detail desired usually should follow the precision and nature of the measurements from which the empirical data is obtained: e.g., from mean RNA and protein production rates in cell populations, to single-cell level, to single-RNA (Golding et al., 2005) and single-protein numbers (Taniguchi et al., 2010), and finally, to the single nucleotide and codon level (Herbert et al., 2008). The models should also consider, e.g., what perturbations are they expected to mimic the effects of. For example, when the measurements explore the effects of overexpressing Gyrase, it will be relevant to add reactions modeling positive supercoiling buildup and represent Gyrases and their actions explicitly. 3. Stochastic partitioning in cell division Some of the first stochastic models of single-cell distributions of protein levels in growing cell populations were developed for E. coli (Berg, 1978; Rigney, 1979), whose division is morphologically near-symmetric (Marr et al., 1966; Trueba, 1982). These models usually assumed that proteins were homogenously distributed in the mother cell, at the moment of division, or were evenly distributed by the two geometric sides of the cell. The latter assumption is usually accurate for many proteins existing in large numbers, unless they are clustered Fig. 3. Schematic representation of transcriptional bursting. The gene can be in the OFF state or in the ON state, during which time there is a significant increase in the numbers of mRNA transcribed by the gene. I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 6 (Kuwada et al., 2015). In this regard, a genome wide study of protein localization (Kuwada et al., 2015) using cells from the ASKA þlibrary (Kitagawa et al., 2006) reported that only 40% of the (869) proteins observed were partitioned asymmetrically. Further, ~50% of these (~250 proteins) formed aggregates, likely due to inclusion bodies, which ‘artificially’ promote asymmetry, due to their preferential localization at the cell poles (Lindner et al., 2008), caused by a phenomenon of nucleoid exclusion from midcell (Woldringh et al., 1994; Winkler et al., 2010). E.g. assume a cell with two such clusters and that, once the cell divides, each daughter cell inherits only one cluster (Gupta et al., 2014a). Since these clusters will likely differ in size, this already generates asymmetries in division. Further, in the subsequent cell division events, the partitioning of the clusters will be necessarily uneven, as only one of daughter cells will inherit the only protein cluster of the mother cell. Recent observations of these phenomena were reported in, e.g. (Neeli-Venkata et al., 2016), where both natural and synthetic protein aggregates, as well as the nucleoids, were tracked, which allowed the development of more complex models of partitioning in cell division. Given the above (and since the localization of the proteins was widely diverse, suggesting that fusion with GFP did not play a major role (Kuwada et al., 2015)), the actual fraction of proteins that are asymmetrically partitioned in cell division in natural conditions may be less than 25%. This implies that, for some proteins the asymmetric partitioning is selectively advantageous while for others it is not (Erjavec et al., 2008; Kysela et al., 2013). It further suggests that E. coli may have evolved ingenious asymmetric protein partitioning processes (Ventura and Sourjik, 2011). One case where it could be advantageous is the asymmetric partitioning of unwanted protein aggregates, which could be a means to ‘renew’ most cells of a lineage, at the cost of a few individuals that would inherit most aggregates (Stewart et al., 2005; Lindner et al., 2008). In general, there are at least two means by which partitioning in cell division can be asymmetric. Either the cell components are unevenly split even though division is morphological symmetric (Fig. 5, bottom left) or, division is morphologically asymmetric, causing the larger daughter cell to inherit more components, provided that they are homogenously distributed (Fig. 5, top right). Finally, both asymmetries (positional and morphological) can occur, which could further enhance functional asymmetries between sister cells. To introduce these asymmetries in stochastic models of gene expression of growing cell populations, Huh and Paulsson developed several models of partitioning schemes (Huh and Paulsson, 2011), shown in Fig. 5. They range from ‘perfectly homogeneous’ (top left scheme) to ‘entirely heterogenous’ (one-gets-all, bottom right scheme). The intermediate schemes considered here are based on known mechanisms by which cellular components (e.g. non-functional proteins) can form larger components or, instead, on cell components that preferentially locate in certain positions (e.g. the Tsr protein complexes preferentially locate at the cell pole due to nucleoid exclusion (Neeli-Venkata et al., 2016) as well as due to the action of Tol-Pal complexes (Santos et al., 2014)). Finally, the authors considered morphological asymmetries (see the bottom right model in Fig. 5). These strategies have since been applied to various models of cellular processes (Jahn et al., 2015). Fig. 4. Example models of a 2-rate limiting step process of RNA production from an active gene (Pro ON ) subject to locking (Pro Locked ) due to positive supercoiling buildup and unlocking by the action of Gyrases (G). Shown are the reactions of RNA production along with promoter locking and unlocking (top left table), the initial quantities of each component (top right table), and the general parameters and the parameters specific to each condition (mid and bottom tables, respectively). Three conditions are considered: (A), (B), and (C), differing in the numbers of Gyrases (G) and in the kinetics of the closed complex formation (k 1 and k -1 ). Finally, the bottom right figures show example probability distributions of time intervals (Δt) between consecutive RNA production events, along with the mean ( μ ) and standard deviation ( σ ) of the distributions of each condition. Note that the axes ranges differ. We performed 1200 simulations of each model, each 10 6 s long, and extracted the time intervals between consecutive RNA production events. I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 7 Assuming the use of, e.g., the simulator SGNS2 (Lloyd-Price et al., 2012), one can introduce such partitioning schemes with ease. For example, assume a growing population with cells expressing a given protein. To introduce partitioning of this protein by pair formation, one can introduce reactions by which these proteins bind and unbind from one another. The rates of these events will define the fraction of paired proteins. One can also set reactions for pairs to bind to a third protein, etc., creating larger clusters. To introduce uneven partitioning at the single protein level, one can define uneven probabilities of the two daughter cells inheriting a protein. Similarly, the effects of morphological asymmetries in division can also be modeled by setting uneven probabilities of inheritance and then defining different rate constants in the two daughter cells (to account for differences in their volumes). It is also possible to model spatial heterogeneities, such as due to the presence of the nucleoid at midcell. E.g., in (Gupta et al., 2014a) the cell was modeled as a 1-dimensional structure, divided into N equal sections. Then, to simulate the effects of the presence of a nucleoid at midcell on the spatial dynamics of large molecules, the model tuned different probabilities of cell components moving between next nearest neighbor sections as a function of their position. Initially, the model assumed equal probabilities of a molecule in section n to move to sections nþ1 or n-1 (except at the extreme sections). Next, to introduce nucleoid occlusion, which heavily reduces the chances of large cellular components to occupy the midcell region, it was decreased the chance of large molecules of moving towards the geometric center of the cell in accordance with the shape of an inverted Gaussian function, centered at midcell (i.e. the closest from midcell, the less likely it is to further approach the cell center). This model has significant flexibility. E.g. it allows introducing single-cell variability in the nucleoid center or size by introducing distinct Gaussian functions in individual cells. A similar strategy can be used to model nucleoid replication. E.g. one can implement two, lesser wide Gaussian functions, that first emerge at midcell and then move over time towards each becoming centered in between midcell and the cell extremities. 4. Models combining stochastic genetic circuits and partitioning in cell division The first models considering cell-to-cell variability in protein numbers in growing cell populations (Berg, 1978; Rigney, 1979) made evident that cell division alters the shape of the single-cell distribution of protein numbers. It was shown that, even when assuming binomial partitioning in division, single-cell distributions of protein numbers could not be well described by a Poisson distribution, if proteins numbers are small. One of the reasons for this is that, unlike noise from gene expression, noise from cell division is introduced at specific points in time, causing a temporary surge of cell-to-cell variability in protein numbers, which then “dissipates” over time, until the next division event. Consequently, the effects of stochastic partitioning should differ with the rate of cell division (Bertaux et al., 2018; Scott et al., 2010; Schwabe and Bruggeman, 2014; Klumpp et al., 2009). Also, the models suggested that the variability in RNA and protein numbers introduced by cell division can be regulated by the timings of cell division (e.g. the degree of synchronization, the nature of the partitioning process such as its bias, and the kinetics of gene expression events (Wang et al., 2015, 2018; Soltani and Singh, 2016)). In this regard, it has been considered that, in real cells, when ‘implementing’ different growth rates, cells exhibit different mRNA lifetimes (Esquerr� e et al., 2015) and RNAP and ribosomes abundance (Klumpp et al., 2009), among other. This suggests that the cells execute different global transcriptional programs, which likely differ in noise levels (Bar-Even et al., 2006; Esquerr� e et al., 2014). Meanwhile, Bertaux and colleagues (Bertaux et al., 2018), using a 2-step gene expression model (transcription followed by translation), showed that a single-cell distribution of protein numbers can be kept independent from the growth rate, by tuning the translation rates. As examples of how partitioning in cell division affects distributions of gene expression products in an isogenic cell population, we simulated the dynamics of single-cell protein distributions in cell lineages subject to different partitioning schemes. Models and results are shown in Fig. 6. We considered: (A) perfect proteins’ partitioning; (B) random size partitioning and (C) all-or-nothing partitioning. Each cell was set to grow at a constant rate and divide when reaching a specific length. In division, the length of the daughter cells relative to the mother cell was set according to a binomial distribution with p ¼0.5 (except for case B, where it is based on a beta distribution with α ¼β ¼2). The gene expression model in each cell is depicted in Fig. 6 (reactions (2.9), (2.3), (2.10) and (2.4) for RNA and proteins production and degradation, respectively). The rate constants and parameters of the partitioning schemes are also shown in Fig. 6. One important information required to implement cell division is whether a component is replicated (e.g. a gene) or partitioned (e.g. proteins) in division. In the case of SGNS2, the partitioning of the components of the mother cell is a process by which a component either “remains” in one daughter cell, or it “moves” to “the other” daughter cell. This is implemented as follows (see an example complete implementation in Supplementary Material). Let “split” be a cell division event represented as a chemical reaction, with substrates, a rate constant, and products, as in (4.1): split:Protein@Cell –[1]–>@Cell þPromoter@Cell þ:Protein@Cell; (4.1) Here, the cell component “Protein” of the mother cell is consumed during a split (as it is set to be a substrate of the reaction) and then appears as a product (of the splitting process). Consequently, it is not replicated but will appear in one of the daughter cells. Meanwhile, the component “Promoter” is not consumed, since it is not a substrate, but is produced. As such, a “Promoter” will be replicated (in the sense that it Fig. 5. Schematic representation of partitioning schemes of cell components due to asymmetric positioning and/or morphologically asymmetric divisions. (Top left) Perfect partitioning, (Top center) Pair formation, (Top Right) Random size partitioning, (Bottom left) Preferential partitioning, (Bottom center) Cluster formation, and (Bottom right) All or nothing partitioning. I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 8 (caption on next page) I.S.C. Baptista and A.S. Ribeiro BioSystems 193-194 (2020) 104154 9 will “remain” in the original cell and will also appear in the “other daughter cell” that was produced in this process). For purposes of implementation alone, note that the product “Protein” is preceded by a “:”, while the product “Promoter” is not. For details on the various possible schemes of partitioning of the components that are not replicated, see (Lloyd-Price et al., 2012). Finally, in the case of SGNS2, there is the possibility of delayed events. To handle this (i.e. for components in the waitlist to not be lost during cell division), the simulator will always implement the same rules of partitioning in cell division to both the components in and the components out of the waitlist. However, those in the waitlist will remain in the waitlist of the daughter cells until the waiting time is over, as which point they will other appear in both daughter cells or in only one of the daughter cells, depending on the rules set, as described above. From Fig. 6, first, comparing conditions A, whose proteins are distributed evenly between sister cells (see inset in Fig. 6A), and B, whose proteins partitioning follows an unbiased binomial distribution, one finds that the mean difference between protein numbers of sister cells following division is near-identical (aside from stochastic fluctuations). However, the latter shows a larger standard deviation and kurtosis, due to the randomness in protein partitioning. Meanwhile, in condition C, there is an absolute asymmetry in the partitioning, causing the mean and standard deviation of the difference in protein numbers between sister cells to be much higher than in A and B (while the kurtosis is lower). Regarding small genetic circuits (Segall et al., 1986; Arkin and Ross, 1994; Arkin et al., 1998; Samoilov et al., 2002; Wolf and Arkin, 2002, 2003; Levchenko et al., 2004; Klumpp et al., 2009; Uriu, 2016; Bertaux et al., 2018; Rosenfeld et al., 2005), the effects of cell division are expected to differ significantly between circuits. For example, in genetic switches, since their stability depends on the number of proteins of the active gene, the higher the asymmetry in the partitioning of the ‘dominant’ protein during cell division, the more likely it will be that state switching will occur in the daughter cell receiving less proteins. Meanwhile, in genetic clocks, the same perturbation is expected to cause the circuit to lose ‘track of time’ but not ‘change state’, since the system only has one stable ‘noisy attractor’. It is also expected a relationship between the period of oscillation of the Repressilator and the effects of tuning cell division rates. The effects of the ‘uncertainty’ introduced by stochastic partitioning and division times on the dynamics of small genetic circuits has recently been explored in several studies, assuming various initial conditions and models. For example, Stamatakis and Mantzaris (2010) studied the effects of intrinsic noise and division cycle on the dynamics of a model stochastic oscillator. In some conditions, coherence resonance emerged between the cell cycle and the oscillatory dynamics of the oscillator, which enhanced the robustness of the oscillation. Meanwhile, Gonze (2013) studied the perturbations introduced by the cell division cycle on the oscillatory dynamics of a Repressilator. The study focused on the effects of noise due to the limited number of proteins and cell-to-cell variability in these numbers generated by noisy partitioning in cell division and noisy division times. Interestingly, they reported that, within the range of parameter values analyzed, cell-to-cell heterogeneity was a stronger source of variability in the oscillations of the Repressilator than the division events. In (Lloyd-Price et al., 2014), noise in partitioning in division was also found to affect the long-term dynamics of toggle switches and Repressilators, for various partitioning schemes. The behavior of these circuits was analyzed for binomial partitioning, pair formation (an ordered partitioning scheme) and random accessible volume (a disordered partitioning scheme) as defined in (Huh and Paulsson, 2011). It was found that changing the partitioning scheme sufficed to alter significantly the circuits’ long-term dynamics. Another study (Tourigny, 2014) showed that the oscillations can suffer phase shifts after cell division although, if the oscillations have faster periods than cell division, the oscillations will likely be robust to the perturbations. Meanwhile, Ahmed et al. (2015) studied the dynamics of populators of genetic oscillators assuming that cell division acts as a ‘resetting’ process and showed that, in this scenario, cell division causes phase drifts. When studying genetic oscillators, Veliz-Cuba and colleagues (Veliz-Cuba et al., 2015) reported that there is a high cell-to-cell variability in the oscillation amplitude, even though sister cells were strongly correlated for many minutes after cell division. They noted that, since fluctuations due to intrinsic noise have short-term effects (Rosenfeld et al., 2005), the long-term variability observed is likely from extrinsic noise sources, namely, the partitioning of cellular components in cell division. They also noted that time delays in gene expression can enhance/decrease the impact of cell division on the circuits’ dynamics. Meanwhile, Jaruszewicz-Bło� nska and Lipniacki (2017) explored the influence of the cell cycle on the dynamics of model genetic switches (with one of the genes close to the origin of replication and the other close to the terminus of replication site). It was found that the rate of the cell cycle can affect the state of the toggle switch and that this influence depends on the locations of the genes of the switch and their state at the moment of division. Finally, results in (Bertaux et al., 2018) suggest that the regulation of cell size when changing cell division rates may allow enhancing/reducing the effects on gene expression noise. Oppositely to the studies above, Paijmans and colleagues (Paijmans et al., 2016) considered that in live cells, since most small genetic circuits conduct complex operations (such as counting time and decision making, as described above), their dynamics needs to be reliable (i.e. stable), even following DNA replication. To investigate potential mechanisms evolved to ensure their stability in cell division, they started by using models to identify which events during cell division most perturb the circuits. Next, they proposed means by which bacteria (specifically, the cyanobacterium Synechococcus elongatus) protect the circuits from those perturbations. Relevantly, they showed that one such mechanism could be the multi-copy presence of the circuits and their asynchronous replication, as a means to minimize changes in gene expression rates and, thus, protein numbers. Subsequently, in (Paijmans et al., 2017), the authors investigated sources of robustness of the repressilator (Elowitz et al., 2002) and the dual-feedback oscillator (Stricker et al., 2008) in growing cells populations, and suggested that the perturbations due to cell division can couple the period of oscillation to the period of the cell cycle. Finally, recent models have become more detailed regarding cell volume growth and cell growth stage dynamics (Song et al., 2015; Marguet et al., 2019; Mura et al., 2019), and focused on how these variables affect the relationship between noise from division and noise Fig. 6. Probability distributions of differences in protein numbers between sister cells following cell division, when splitting proteins according to ‘perfect partitioning’ (A, blue), ‘random size partitioning’ (B, green) and ‘all or nothing partitioning’ (C, in red) (cell lengths and/or RNA and protein numbers). Shown are the distributions as well as their mean ( μ ), standard deviation ( σ ), and kurtosis (k), and respective standard errors. The insets show example time series of protein numbers in a single cell line (vertical black lines represent moments of cell division). Simulations using SGNS2, assuming the model depicted in the top left table. Parameter values common to all conditions are depicted in the Table ‘General parameters’, while those differing between conditions are shown in the middle table describing conditions A, B, and C. All cells have a constant growth rate and divide when reaching the same specific length. However, after division, sister cells can differ in size, implying that they will take different amounts of time to perform their division. Here, ‘L’ stands for cell length while ‘N’ stands for number of RNA or proteins, prior or after cell division. We performed 100 simulations per condition (6 division events from one original cell, per simulation). Each cell lifetime was 1800 s. Proteins number were extracted once per second. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) I.S.C. Baptista and A.S. Ribeiro