Noise-driven cell differentiation and the emergence of spatiotemporal patterns

The emergence of phenotypic diversity in a population of cells and their arrangement in space and time is one of the most fascinating features of living systems. In fact, understanding multicellularity is unthinkable without explaining the proximate and the ultimate causes of cell differentiation in time and space. Simpler forms of cell differentiation can be found in unicellular organisms, such as bacterial biofilm, where reversible cell differentiation results in phenotypically diverse populations. In this manuscript, we attempt to start with the simple case of reversible nongenetic phenotypic to construct a model of differentiation and pattern formation. Our model, which we refer to as noise-driven differentiation (NDD) model, is an attempt to consider the prevalence of noise in biological systems, alongside what is known about genetic switches and signaling, to create a simple model which generates spatiotemporal patterns from bottom-up. Our simulations indicate that the presence of noise in cells can lead to reversible differentiation and the addition of signaling can create spatiotemporal pattern.


Introduction
The traditional idea of a living cell where every organelle, every reaction, and every interaction is part of a clock-like order has long been shattered by the understanding that biological systems usually struggle to function in noisy environments. One might consider life to be an uphill battle against pandemonium, where disarray is the norm and spheres of order-i.e., biological systems-are rarities that are unlikely to appear in the first place. In this view, noise is a nuisance that natural selection always attempts to eliminate. It is for the same reason that selection cannot increase the fidelity of replication beyond a certain threshold; the biological cost of increasing fidelity simply becomes too high at that point [1].
A different view has recently gained some grounds [2][3][4][5]. In this view, biological systems that regulate and utilize the noise can have higher fitness under certain circumstances. Had biological systems been utterly deterministic, adaptation-i.e., the emergence of a new phenotype or a change in the gene expression pattern to utilize a new food source-would have been impossible without the emergence of new mutations. In reality, noise in the cell can result in beneficial non-genetic diversity in otherwise genetically homogenous populations-e.g., cyanobacteria [6] and yeast [7]. But what mechanism can account for the presence of phenotypic diversity amongst daughter cells that are genetic clones of each other? Is it possible for a stochastic mechanism to explain the non-genetic diversity? Even if such stochastic explanation were offered, how could this explanation possibly account for the ordered spatiotemporal patterns in spatially-extended cell population?
The model of cell differentiation proposed in this work, henceforth referred to as the noisedriven differentiation (NDD) model, accounts for the peculiarities of this biological phenomenon by weaving noise into an explanation of cellular behaviors at the time of differentiation. While on the surface, this approach might seem lofty and even radical, the model discussed in this paper is parsimonious when it comes to the mechanisms requisite for its operation. The NDD model rests on 8 components (Table 1). Some can be regarded as facts, based on reliable empirical evidence from biological systems (components #1 and #2), while others are more accurately described as assumptions (components #3-8).
There is a plethora of phenomena within a cell that can contribute to its intrinsic noisee.g., transcription regulation, transcription factor binding to the DNA, RNA processing in eukaryotes, translation, post-translational modifications, protein complex formation, protein and RNA degradation, etc. Single-cell level measurements of gene expression further cements the notion that cells are intrinsically noisy when it comes translating its genotype into Table 1. The components of the NDD model.

Components Justifications
1 Noise, resulting from a plentitude of sources, is an inseparable part of a living cell.
2 Stochastic partitioning of cytoplasm during cell division and the random distribution of molecules in the cytoplasm determine the cytoplasmic contents of the daughter cells.
3 The fate of a cell is determined when it is born. Based on the assumption that cell-fate-determining factors are in small numbers in a cell and the stochastic distribution of these factors during cell division determines the fate of the newly-born daughter cells.
4 Cell fate is determined by a switch. Genetic switches have been observed in a variety of taxa (reviewed in [2]), and has been proposed as a model to account for cell differentiation (e.g., see [20]). 5 The interaction between the building blocks of the switch determines its bias.
Our assumption based on our knowledge of wellknown genetic switches, such as λ phage (see [21]).
6 All the information needed to construct the switch is genetic.
We assume that, while stochasticity is what drives the decision made by the switch, the information necessary to construct the switch is encoded in the genetic content of a cell. 7 The robustness of the switch is the result of a complex network of interactions.
8 Cell fate is influenced by its location and its environment. � We assume the the switch determining cell fate should, in addition to being swayed by the intrinsic factors, be influenced by its neighbors.
phenotype [23]. The displacement of the division plane relative to the middle of the cell can result in an unequal distribution of cell content between the daughter cells, even if molecules are homogeneously distributed within the cell. In fact, the central role of asymmetric cell division in the diversification of cells, from Drosophila to mammals has been known for many years [18,19]. The components #1-2 reflect the role of stochasticity in living systems based on these observation. Thus far, two types of solutions to the problem of cell differentiation have been proposed: the first category consists of models that rely on cell-cell communication (reviewed in [24]) and the second category relies on asymmetric cell division (reviewed in [25]). The research project within the confines of the former category is mainly a quest to find the building blocks of the apparatus that makes the specific kind of cell-cell communication needed for cell differentiation. The latter category, on the other hand, presumes the asymmetric cell division to result in differentiation. Hitherto unknown and often complicated mechanisms have been proposed to explain the asymmetric distribution of fate-determining factors during cell division [26,27]. Both categories rely on physical interactions at the cellular level. While we agree with the importance of the asymmetric cell division, it seems to us that a stochastic model of differentiation, like the NDD model, negates the need for new mechanisms. In this model, we adopt the view that stochastic processes result in differentiated cells due to the distribution of key proteins, instead of cells differentiating by receiving signals after they are born (component #3).
The component #4 is based on the idea that characteristics of a cell can be changed by a switch (Not a very recent idea, e.g., [28]). The notion that cell fate is determined by a switch is best illustrated by the now famous case of the λ phage. The process by which the phage decides to integrate into the host's genome-i.e., lysogenic-or to replicate copies of itself in the cell until it bursts open-i.e., lytic-can be explained by a stochastic switch which makes that portentous decision in a probabilistic fashion, while taking into account the presence of certain key factors [29]. One can assume that the bias of this switch is determined by the interactions of its building blocks (component #5). For example, upon infecting bacterial cells, λ phage proceeds to lyse the host, but as the concentration of CII protein increases, so does the likelihood of the reactions suppressing the activation of pR and pL promoters, relevant to the onset of the lytic trajectory, which in turn, tilts the scale away from lysis towards lysogeny [21]. We propose that phenotypic diversity arises from the effect of the noise on a genetic circuit that exhibits a switch-like behavior (component #6). The notion that different phenotypes are produced from the same genotype as a consequence of noise is widely observed in nature (reviewed in [30]).
How robust can be a fate-determining toggle switch in the face of new mutations? Sharifi-Zarchi [22] took advantage of the gene expression profiles of 442 mouse embryonic cells to construct a network of key transcription factors (TFs). While a regulatory circuit with two TFs could explain differentiation, They reasoned that such a simple switch is susceptible to mutations. To construct a robust switch, they built a circuit with two clusters of TFs with correlated expressions. Expectedly, the alternative switch, which involved more interactions, was much more robust. We would expect different levels of robustness for a switch, given its biological importance in evolution (component #7).
The components #1-7 are sufficient to generate a population of cells with different proportions of two phenotypes (Fig 1). While this kind of fate determination is adequate vis-à-vis primitive cells with no organization, it does not allow the emergence of multicellularity. An additional component is necessary to explain this major transition from mere phenotypic differentiation to ordered spatiotemporal patterns in the body of a multicellular organism. For self-organization to occur, we assume that the toggle switch determining cell fate should, in addition to being swayed by the intrinsic factors, be influenced by its neighbors (component #8).
To test the general veracity of the NDD model, we used a simple model of cell aggregation. In this model, a simple switch is defined that can switch between phenotypes, A and B. We look at the noise resulting from the random distribution of fate-determining proteins during cell division (component #1), the stochastic positioning of cytoplasmic cell division (component #2), and the effect of signaling from neighboring cells on other cells (component #8) on a model of cell aggregation where the fate of cell depends on a genetic switch (components #4-7). In the discussion, we juxtapose the NDD model with some of the more recent attempts at modeling cell differentiation.

Materials and methods
In the cell aggregation model, the population is made up of cells, where each cell is a circular particle defined by its state variables-e.g., spatial position, size, and phenotype. The simulation geometry is a L × L square and no flux boundaries. It is assumed that the relative amount of two key transcription factors, X and Y, controls the cell types; hence, in this model, a cell can have two phenotypes, A and B, as shown in Fig 1. The dominance of protein X leads to phenotype A and the dominance of protein Y results in phenotype B. In fact, a positive feedback loop influences the decision-making process. Two negatively coupled repressors mutually inhibit the expression of the gene that encodes the other repressor-i.e., a toggle switch (component #4). The rate of this mutual repression is represented in the form of a Hill function [31]. This positive feedback loop results in two stable steady states, hence implies non-linear dynamical equations. Non-linear differential equations govern the changes in the number of the repressor proteins, X and Y (Fig 1); Here, β is the effective rate of protein synthesis and n is the Hill coefficient, which represents the degree of competence. The number of repressors are represented in the unit of their dissociation constants, K s , and time is rescaled by degradation rate of proteins [31][32][33]. Biologically-reasonable values were chosen for the parameters used in our simulation such that Eq 1 would be bistable (following [31]). This bistable regulatory network has two attractors corresponding to its stable steady states. Based on the amount of proteins at the cell division time, the cell can be in the domain of each attractor, which determines its fate. Depending on the intensity of inhibitory effects of TFs (through the values of constants in the Hill function [31]), the two domains of attractors could be equal or not (component #5). Fig 2 shows an example of such behavior in our cell aggregation model. S1 Movie shows the changes in the distribution of TFs in cells around their attractors during the emergence of generation 12.

Population growth algorithm
Simulation starts with a single cell with phenotype A. Each iteration in the simulation can be divided into four steps: 1. Cell growth: The repressor proteins inside the cytoplasm interact with each other and their numbers, X and Y, are updated; however, because of their low copy numbers, instead of deterministic equations (Eq 1), their fluctuations are captured by the Gillespie algorithm [34] as a stochastic dynamics for discrete values. According to this algorithm, a probability of occurrence will be assigned to every biochemical reaction in the system. Every protein (X or Y) is produced with a probability according to the first term on the right hand sides of the Eq 1. As the number of protein X increases, it further represses the production of protein Y and vice versa. Every protein degrades proportional to its number. In every step of the Gillespie algorithm, one of the above reactions occurs and the time will be updated. The process continues until the number of proteins reaches a steady state. The production of X and Y, until it reaches a steady state, is simulated as a "burst", independent of cell growth. After this burst of gene expression, cells grow linearly in size.

Cell division:
Even after the number of proteins in a cell reaches the steady state, the cell continues to grow. The growth stops only after the cell reaches a critical size. At this point the cell divides into two daughter cells. The content of the mother cell is distributed among her daughters according to a uniform distribution. In reality and in the presence of active transportation, one can still expect a uniform distribution of molecules in the cytoplasm [35], making this assumption biologically reasonable. The position at which cell division occurs is randomly chosen based on a normal distribution, N ð0; s 2 Þ; which its mean value is the mid plane of the cell (also the reference point), and σ 2 as the variance around the mean value (component #2). In this study σ 2 varies up to 0.1 parentage of the radius of the cell. At the time of birth, the phenotype of each newborn cell is determined based on the cytoplasmic contents (number of key proteins, X, and Y at the time of birth) inherited from the mother cell (component #3). During the cell growth, the number of each protein has a stochastic trajectory in the domain of its attractor and finally it will reach its steady state. In this model, phenotypic change is reversible, meaning that the phenotype can change between the two possible states over generations. Since in our simulations, daughter cells have similar volumes, we consider the number of proteins distributed between them, and not their concentrations.

3.
Relaxation: After a cell divides, the cells push each other outwards to make room for the new daughter cells [36]. Simulation proceeds by repeating the steps #1-3. It is worth noting that, without considering self-organization, the process described above would result in a disordered blob of cells.

Self-organization:
To involve the self-organization phenomenon in the process of cell maturation (component #8), cells secrete some signaling molecules, with concentration C s , which affects the propensities in the Gillespie algorithm and, consequently, the production of proteins. The signaling molecules diffuse in the medium according to the following reaction-diffusion equation: Here, k sp , k sc and k sd represent, respectively, the rate of production, consumption and decay of the signaling molecules and D s is the diffusion coefficient of the signaling molecules. C A and C B respectively show the number of cells with phenotype A and B at each point of the medium. In our simulations, we used D s = 10 −11 m 2 /s, k sp = 0.01kg −1 s −1 , k sc = 0.0001kg −1 s −1 , k sd = 0.01s −1 , and the protein dissociation constant is K s = 0.01m −3 . In these simulations, the secreting cells are those with phenotype B; hence, the production of signaling molecules is proportional to the amount of B cells. Since both phenotypes consume these molecules, the consumption depends on the number of both A and B cells. When B cells emerge, they secret signaling molecules, which diffuse in their environment. The effective concentration of the signaling molecules at any location determines if a cell at that location is affected by the signal, which would decrease the production of protein X and augment the production of protein Y. Consequently, their surrounding cells would have less chance of producing protein X and their offspring is less likely to be in the domain of attraction of protein X.

Results
The overall behavior of the cell aggregation model demonstrates the principles of our framework-that is, the stochasticity results in phenotypic heterogeneity as the population grows in size (S1 Movie). To further illustrate how each source of noise affects the cell differentiation, we focused on each source separately in the simulations.

The stochastic positioning of division plane and the stochastic distribution of key proteins affect differentiation
One source of intrinsic stochasticity stems from the random positioning of the division plane. This factor would disproportionately influence the number of molecules that exist in low numbers within cytoplasm. In this work, it has been postulated that the determinants of cell fate are low in numbers and thus, greatly affected by stochasticity.
To demonstrate this phenomenon, the position of the division plane was allowed to vary with respect to the mid plane of the cell. Starting from a cell with phenotype A, in which the protein X is dominant, the population heterogeneity -i.e., emergence of phenotype B-was traced over 12 generations. The results are shown in Fig 3. When the division plane is situated in the middle of the cell, and the TFs are relatively abundant, very few cells differentiate. As the variance in the cell-division plane-i.e., σ 2 in the normal distribution used to choose the position of division plane-increases, so does the proportion of B cells. This phenomenon is dependent on the number of proteins, since such a bias is more pronounced when the number of proteins is relatively low. In fact, with large numbers of TFs in a cell, it will be more likely for its daughters to have almost the same density of TFs as their mother. Thus, they will be in the same domain as the mother in the phase space, and their fates will be identical to hers. This can be seen clearly in the lower curves in Fig 3. However, for low copy numbers of TFs, the difference between TF numbers in two daughter cells becomes more prominent and can even lead to different cell fates. Therefore, it is possible to have heterogeneity in the population in the absence of any other noise, i.e., cells with low TF numbers are heterogeneous even with no variance in division-plane displacement (Fig 3a). In Fig 3a, the distribution of TFs is deterministic across the length of the cell-e.g., if X = 10 and cell divides exactly on its equator, each daughter receives 5 TFs. We add stochasticity to the distribution of TFs by randomly drawing the number of TFs at any location along the length of the cell from a uniform distribution (U(0, L), where L is the length of the cell, provided P L l¼0 X l ¼ TF X Uð0;LÞ jjUð0;LÞjj and

PLOS ONE
Noise-driven cell differentiation P L l¼0 Y l ¼ TF Y Uð0;LÞ jjUð0;LÞjj ). This additional source of noise has a tangible, though slight, effect on the percentage of differentiated cell (Fig 3b).

Signaling can create spatial order
In the cell aggregation model, B cell can release signals in the environment. These signals diffuse at a slow rate and, consequently, have a very short radius of influence. The absorption of these signals by other cells in the population affects the number of proteins involved in the switch-that is, switching to the phenotype B during cell division becomes more likely (Fig 4). When this environmental signaling is added to the population, the cells organize in a non-random fashion, a stark contrast to the random heterogeneity observed before (S2 Movie).    change in the global clustering coefficient of two cell types A and B, as the increase in the intensity of noise. In graph theory, the clustering coefficient, C, represents a measure of cohesion of a graph, i.e., how clustered vertices are, C ¼ 3ð triangles on the graphÞ ð connected triples of verticesÞ [37,38]; here a triangle means three neighbour cells. It could be seen that in the presence of signaling, Fig 5b, the clustering coefficient of orange cells increases strongly; which means the differentiated cells created a denser population. The organization observed will increase over time and the community of orange cells will develop (S3 Movie).

Discussion
Molecular processes in the cell are noisy events that result in varying degrees of heterogeneity. Taming this inherent noise is vital for the emergence and the continuation of life. In fact, life can be characterized as a system with the capacity to control noise. The phenotype of a cell is generally stable, but during cell division, this cell can produce daughter cells with different phenotypes via symmetric or asymmetric cell division. The resulting non-genetic phenotypic diversity is a way to achieve adaptation in a fluctuating environment by producing phenotypically diverse offspring without any need for genetic change. Given the variety of sources of noise, the cell fate determination can be a stochastic process. One can imagine a few genes involved in cell fate determination, where the noise in the cell affects the proportion of daughter cells born with a certain phenotype. The ability to change the phenotypic proportion of daughter cells via a stochastic mechanism, which is also tunable, is a superb strategy to outcompete rivals bereft of such gift. Proposing a stochastic model of cell differentiation is not an entirely novel concept, e.g., see [39,40] as examples of an impressive body of work produced by Kunihiko Kaneko and his colleagues on this subject and [41,42] as similar proposals regarding the possible role of stochasticity in generating phenotypic diversity. The NDD model is an attempt to merge observations about the role of noise in living systems and the prevalence of toggle switches in biological processes to explain cell differentiation. Here, we utilized a bistable switch to demonstrate how such model would create population-level heterogeneity and how signaling can give spatial order to this heterogeneity. This approach is in keeping with the recent emphasis on the importance and the prevalence of noise in biological functions, specifically cell fate [2,4,43]. The model of cell aggregation used in this study allowed us to test all the components of the NDD model, barring components #6 and #7, which demand through investigations of their own. This model of cell aggregation provides us with a relatively realistic depiction of the process that results in phenotypic differentiation in a population. To illustrate the virtues of our approach, it is worth considering how it differs with some of the more prominent approaches to explain cell differentiation.
Olimpio et al. [44] explored the role of cell-cell communication in the emergence of spatial order in a cellular automaton. In their model, N cells exist on a lattice, and there is no cell growth or division, thus the effect of the distribution of fate-determining proteins in the mother cell on the fates of the daughters (component #3 of the NDD model) does not factor in. When they include noise in their model, it consists of a Gaussian noise added to the activation threshold of the switch. In the NDD model presented here, we consider different sources of noise affecting biological systems separately, and in this work we looked at the effects of some of these sources on the emergence of spatiotemporal heterogeneity. In the NDD model, the intrinsic noise alone generates heterogeneity without including cell signaling.
Perez-Carasco et al. [20], use a bistable genetic switch to translate the analog morphogen gradient into a spatial pattern. Similar to Olimpio et al. [44], cells do not divide and the noise, when introduced, is in the form of transcriptional bursts. The point is not to downplay the relevance of these works, far from it; the NDD model takes the idea of a toggle switch and the role of noise and attempts to (i) tease apart the different sources of noise and (ii) organize the observations and the assumptions needed to explain the emergence of spatial heterogeneity in living organisms.
The results from our cell-aggregation model shows that the NDD model can provide an explanation for the emergence of spatial heterogeneity, simply based on the intrinsic noise in living systems. Without cell-cell communication (component #8), the mere presence of noise (components #1-3), coupled with a genetic toggle switch (components #4-7), can generate spatial heterogeneity. The addition of signaling generates the spatial pattern. Our results further cement the notion that a noise-driven system can generate cell differentiation. In their criticism of a noise-driven alternative to their model, Suzuki [39] considered it unlikely for a noise-driven model to maintain the exact levels of stochasticity needed to produce the desired proportion of differentiated cells to stem cells. In our view, if a switch is robust (component #7 of the NDD model), then it will be able to maintain its bias in the face of new mutations.
It has not escaped our notice that a bistable switch, as described here, can only create reversible cell differentiation, as seen in prokaryotes. In fact, Suzuki [39] points out this problem with a model of cell differentiation based on a bistable toggle switch. Can the NDD model be used to simulate the irreversible eukaryotic differentiation from stem cell to differentiated cell? Replacing the bistable switch with a multi-stable variant might solve this problem [45]. In fact, in a follow-up work by Khorasani, Sadeghi, and Nowzari-Dalini [46], replacing the bistable switch by a tri-stable one was enough to generate irreversible differentiation.
The ability of the cells to differentiate into different types was the crucial step that enabled the ancient solitary cells to leave the primordial soup behind and evolve into the vast array of specialized cells we see today. As Queller [47] point out, there are different shades of organismality -i.e., the ability for components to work together with little conflict among them-, each shade resulting from the affinity of the members of the system to cooperate versus the temptation to cheat. We can sidestep the problem of conflict since in prokaryotic multicellularity, e.g., biofilm, and in most truly multicellular eukaryotes, the cells are highly related, thus lowering the probability of cheating [48]. Without tangible levels of conflict, multicellularity as a trait becomes patently advantageous. In their seminal work, Maynard Smith and Szathmáry [49] considered two possible mechanisms to account for the emergence of cell differentiation: one relies on the presence of determinants that prohibit the stem cell to differentiate, and the other postulates the cell-cell contact as a mechanism that determines cell fate. While these suggestions account for how the multicellularity might be sustained, they do not explain how this major evolutionary transition could have occurred in the first place.
One of the quintessential aspects of the discussed model is its population-level perspective. Population-level thinking is one of the main points of the evolutionary theory, and bringing it to explain a cellular phenomenon can lead us to reap valuable insights. While a population of cells has, on average, certain properties relevant to differentiation, e.g., the mean number of key proteins, the average position of cell division plane, and etc., these average values do not tell the whole story. Instead, the variance in these values, i.e., the non-genetic variation present amongst individuals, is the key to understand differentiation (as observed in studies such as in [12,50]). This noise in the population is essentially the fuel that propels cellular differentiation, be it in the reversible differentiation in prokaryotes or the more complicated irreversible ones in higher organisms. We believe that this population-level vantage point is the necessary tool to understand this otherwise mind-boggling biological process. Without this perspective, the task of explaining such a seemingly fine-tuned process devolves into an attempt to come up with complex cellular interactions that would make climbing this improbable biological mountain feasible.
The fact that NDD model breaks down the components required for spatiotemporal heterogeneity, makes it easy to muse about the events that led to the different levels of multicellularity during the evolution of life on Earth. It is easier for cell differentiation to evolve via the emergence of a switch, rather than the less plausible path that involves the evolution of a clockwork mechanism. According to the NDD model, the emergence of early stages of multicellularity only requires the evolution of a suitable switch-the rest of the necessary ingredients needed for the transition into self-organization is provided by the stochastic elements affecting the switch. The major transition from unicellularity to multicellularity -i.e., from phenotypic diversity in a population to from an ordered and stable spatial heterogeneity-only requires one more step: the evolved switch should be simply affected by the signal(s) released by its neighbors (components #8). The spatial information received in this way would bias the switch such that the population-level organization is retained. The NDD model is the logical extension of earlier ideas describing the role of stochasticity in phenotypic variation and the switchlike behavior of genetic circuits vis-à-vis differentiation and multicellularity (e.g., see [51]). It is tempting to postulate a connection between the cell-differentiation switch, postulated in the NDD model, and the toggle switch used in quorum sensing in bacteria [52]. Quorum sensing enables bacteria to regulate their phenotypes apropos of their neighbors and is more robust in a dense community [53]. It seems plausible to consider this type community-based phenotypic regulation as a precursor to similar switch-based mechanisms for cell differentiation in multicellular organisms.

Conclusions
The NDD model is an attempt to construct a model of cell differentiation which takes into account the mounting evidence concerning the role of noise in biological systems. The components of this model consists of biological observations and a few assumptions. Using a cell aggregation model, we show that these components are sufficient to generate spatiotemporal heterogeneity in population of cells. The bistable switch allows for the simulation of the reversible cell differentiation. In a followup work, we used a tri-stable switch to explore the irreversible differentiation using the NDD framework. We believe that breaking down the minimum biological components necessary for cell differentiation provides a framework to reconstruct the chain of events that were required for different levels of multicellularity to emerge during the history of life on Earth.
Supporting information S1 Movie. The change in the distribution of TFs within cells just before they divide. Parameters used are the same as Fig 1. S2 and S3 Movies show 3-dimensional simulations of a community of cells in a layer. Simulation performed in a L × L × h cube and starts with one cell at the centre, the red cell with type A. There is a single layer of cells with height h, corresponding to the diameter of a single cell. The cells grow in volume; after reaching a critical volume they divide-the same as the two dimensional case-their cytoplasmic content distributes between the two daughter cells.