Madelynn McElroy,Kaylie Green,Nikolaos K. Voulgarakis
Self-Regulated Symmetry Breaking Model for Stem Cell Differentiation
cell differentiation,phase transitions,symmetry breaking,self-tuned criticality,mean field theory,bifurcation theory
In conventional disorder–order phase transitions, a system shifts from a highly symmetric state, where all states are equally accessible (disorder) to a less symmetric state with a limited number of available states (order). This transition may occur by varying a control parameter that represents the intrinsic noise of the system. It has been suggested that stem cell differentiation can be considered as a sequence of such symmetry-breaking events. Pluripotent stem cells, with their capacity to develop into any specialized cell type, are considered highly symmetric systems. In contrast, differentiated cells have lower symmetry, as they can only carry out a limited number of functions. For this hypothesis to be valid, differentiation should emerge collectively in stem cell populations. Additionally, such populations must have the ability to self-regulate intrinsic noise and navigate through a critical point where spontaneous symmetry breaking (differentiation) occurs. This study presents a mean-field model for stem cell populations that considers the interplay of cell–cell cooperativity, cell-to-cell variability, and finite-size effects. By introducing a feedback mechanism to control intrinsic noise, the model can self-tune through different bifurcation points, facilitating spontaneous symmetry breaking. Standard stability analysis showed that the system can potentially differentiate into several cell types mathematically expressed as stable nodes and limit cycles. The existence of a Hopf bifurcation in our model is discussed in light of stem cell differentiation.
Self-Regulated Symmetry Breaking Model for Stem Cell Differentiation

In conventional disorder–order phase transitions, a system shifts from a highly symmetric state, where all states are equally accessible (disorder) to a less symmetric state with a limited number of available states (order). This transition may occur by varying a control parameter that represents the intrinsic noise of the system. It has been suggested that stem cell differentiation can be considered as a sequence of such symmetry-breaking events. Pluripotent stem cells, with their capacity to develop into any specialized cell type, are considered highly symmetric systems. In contrast, differentiated cells have lower symmetry, as they can only carry out a limited number of functions. For this hypothesis to be valid, differentiation should emerge collectively in stem cell populations. Additionally, such populations must have the ability to self-regulate intrinsic noise and navigate through a critical point where spontaneous symmetry breaking (differentiation) occurs. This study presents a mean-field model for stem cell populations that considers the interplay of cell–cell cooperativity, cell-to-cell variability, and finite-size effects. By introducing a feedback mechanism to control intrinsic noise, the model can self-tune through different bifurcation points, facilitating spontaneous symmetry breaking. Standard stability analysis showed that the system can potentially differentiate into several cell types mathematically expressed as stable nodes and limit cycles. The existence of a Hopf bifurcation in our model is discussed in light of stem cell differentiation.

In the course of early embryonic development, a group of identical stem cells actively self-replicates while simultaneously differentiating into all the specialized cells necessary for the formation of a multicellular organism. Although development progresses with remarkable precision, small populations of stem cells display significant cell-to-cell variability due to the molecular nature of gene expression [1,2]. It is truly remarkable how this seemingly fragile and heterogeneous cellular mass transforms into a multicellular organism with such exceptional accuracy. Waddington envisioned the differentiation process as analogous to a ball rolling down a hillside of a complex landscape containing various slopes and valleys [3]. In this picture, each valley represents a possible state of the cell during development. Stem cells start at the top of the hill, where they are in a relatively undifferentiated state. The ball, which represents the stem cell, starts rolling down the hill, passing through different slopes and valleys. The path the ball takes is determined by a combination of intrinsic and extrinsic factors that push it in a particular direction. As the ball rolls, it becomes increasingly specialized (differentiated) until it reaches its final, stable state at the bottom of the hill, where it becomes a committed cell. Even though Waddington’s epigenetic landscape was initially proposed as a metaphor, it has since become a valuable framework for understanding how the gene regulatory network of cells and environmental factors interact to shape the development of organisms. The present work is built on three major hypotheses that have been proposed over the years and are based on or inspired by Waddington’s epigenetic landscape. First, Huang and collaborators suggested that pluripotent stem cells can simultaneously express multiple genes that are specific to the progenitors [4,5]. Intrinsic and environmental noise can disrupt this state and cause the cell to express a limited number of those genes that correspond to a certain progenitor. In this view, stem cells are considered highly symmetrical cellular states as they can co-express a higher number of genes. In contrast, progenitor cells are considered to operate in a lower symmetry cellular state, as they express a smaller number of genes. This process of symmetry breaking may also occur during the differentiation of progenitor cells into specialized cells. The view of differentiation as a series of symmetry-breaking events is supported by experimental findings in hematopoietic stem cells, for example, which express genes common to all the cell types they can differentiate into [6,7]. This hypothesis indicates that the process of reaching a terminal cell fate is characterized by a decrease in entropy [8]. However, this decrease is preceded by an initial burst-like increase in entropy, which occurs before the differentiation events [8,9]. The crucial question that needs to be addressed is how cells regulate transcriptional noise to drive the system toward a state of lower entropy. This brings us to the second hypothesis, which suggests the existence of intrinsic molecular devices that regulate transcriptional noise [10,11]. Experimental evidence supports this viewpoint, showing that noise is indeed regulated in embryonic stem cells via molecular signaling and feedback mechanisms [10,12]. The third major hypothesis is that pluripotency may arise collectively within populations of stem cells through intercellular interactions, such as paracrine and autocrine signaling [13,14,15,16,17]. This hypothesis is primarily supported by experimental evidence indicating that paracrine signaling plays a crucial role in cell differentiation [18,19]. It is also supported by the emergence of a dynamic equilibrium of cell states that is observed only in stem cell populations [20]. The above hypotheses can be interpreted in standard statistical mechanics language as follows. The process of differentiation involves a collection of interacting agents (cells) that undergoes a transition from a highly symmetrical state (pluripotency) to a state of low symmetry (committed cells). This transition is achieved by decreasing the inherent noise within the system. This scenario exhibits similarities to symmetry breaking in phase transitions. Indeed, in traditional disorder–order phase transitions, a system shifts from a state of high entropy or disorder to a state of lower entropy or order by adjusting a control parameter such as temperature, pressure, or an external field [21]. The transition from a paramagnetic to a ferromagnetic state in a magnetic material is a classic example of a disorder–order phase transition. At high temperatures, strong thermal fluctuations force the system into a disordered state of zero magnetization, where each spin points in any possible direction (high symmetry). As the temperature of the magnet is lowered, all molecular spins start pointing in the same direction (low symmetry), providing the system with a nonzero net magnetization. Other examples of disorder–order phase transitions include the transition from a liquid to a solid. Various methods can be employed to test the validity of these hypotheses. For instance, if the differentiation of stem cells into progenitor cells represents a first-order phase transition, we would expect to observe the coexistence of both cell types at a critical intrinsic noise level. Alternatively, if this process is a second-order phase transition, significant spatiotemporal fluctuations and correlations may be detected when passing through the tipping point. Although the experimental evidence is currently challenging to directly relate to phase transitions, a recent study showed that an increase in both cell heterogeneity and coordination of gene expression occurred prior to induced pluripotent stem cell differentiation [22]. It has to be noted that the analytical theory of phase transitions is applicable only in the thermodynamic limit (infinitely large systems). In this limit, the system experiences strong ergodicity breaking, which stabilizes a long-range ordered state. However, in the case of small systems of interacting agents, such as stem cell populations, finite-size effects cause significant fluctuations in the system’s state. This type of fluctuation decreases with increasing system size. Thus, the observed cell-to-cell variability in stem cell populations could be attributed to a combination of intrinsic noise, finite-size effects, and the possibility of the population undergoing a phase transition. At the single-cell level, stem cell differentiation has primarily been studied through the lens of ordinary differential equations (ODEs), which display various types of bifurcations, such as supercritical pitchfork and saddle-node bifurcation [23,24]. In this framework, each stable branch in the bifurcation diagram represents a different type of cell, and transitions from one stable branch to another signify cell differentiation. It is possible to model the collective dynamics of stem cell populations by coupling such dynamical systems together through appropriate cell–cell interactions. A prominent example is the work by Furusawa and Kaneko, where they emphasized the significance of oscillatory expression dynamics for preserving an undifferentiated state [25]. In a different study, Huang and colleagues showed that gene regulatory motifs with mutual inhibitory interactions and self-activation can generate multiple attractors corresponding to stem and differentiated states [26]. Other studies have demonstrated that coupled systems of ODEs undergo bifurcations as a function of population size and cell–cell interactions, which occur solely at the population level [27,28,29]. These bifurcations lead to cell differentiation through spontaneous symmetry-breaking events. Similar collective dynamical states have also been observed in experimental settings [30,31,32,33,34]. An alternative method of modeling the dynamics of individual cells is through the use of Boolean networks (BNs) [35,36,37,38]. In BNs, molecular regulators are represented as binary nodes that can exist in an “on” (active) or “off” (inactive) state. These nodes interact with neighboring nodes according to predefined deterministic rules, resulting in various attractors that represent possible functions of the cell. Recently, we implemented coupled BNs to describe the collective dynamics of stem cell populations [39]. In this study, individual cell dynamics were described by isogenic BNs, while a multilayer Ising Hamiltonian was used to model cell–cell interaction through paracrine and autocrine signaling. Numerical analysis in two dimensions revealed that changes in intrinsic noise can induce a series of spontaneous symmetry-breaking events, resulting in various types of cells. In our model, pluripotency is defined as the state of a balanced mixture of next-generation cells (progenitors). When the intrinsic noise drops below a critical threshold, the system undergoes symmetry breaking, leading the population to assume a specific cell state (differentiation). Similar to experimental observations [22], prior to cell differentiation, our model exhibited significant cell-to-cell variability. This type of symmetry breaking resembles the disorder–order phase transition in statistical mechanics. While coupled ODEs and coupled BNs have distinct technical and theoretical characteristics, both can effectively model the collective behavior of stem cell populations. However, coupled BNs that use an Ising-type Hamiltonian offer a unique opportunity to investigate the applicability of phase transitions in stem cell populations. For instance, recent experimental observations of increased cell heterogeneity and the spatiotemporal correlation of gene expression [22] can be qualitatively described by our coupled BN model [39]. Here, we aimed to simplify our previous work by implementing Landau’s mean-field theory of phase transitions (see, for example, [40]). Although this approach is a crude approximation of the collective dynamics of coupled BNs, it provides a simpler way to understand the underlying phenomenology of noise-regulated cell differentiation. Here, we hypothesize that the noise regulation depends on the current state of the population and is governed by a negative feedback equation that considers the interplay between cell cooperativity and heterogeneity. This methodology led to several compelling outcomes, including the existence of diverse types of attractors (cell types). As also predicted by Landau’s original theory, our dynamical system exhibits hysteresis loops that could potentially model phenotypic memory. To address the impact of finite-size effects in small stem cell populations, we incorporated a stochastic forcing term that decreases proportionally with system size. Notably, finite-size effects can cause intriguing phenomena, such as noise-induced bistability and overall phenotype switching [41,42,43,44]. The paper’s structure is straightforward. In Section 2, we examine Landau’s potential function and briefly overview its bifurcation and phase transition analysis. We subsequently introduce the feedback mechanism that self-regulates the system’s intrinsic noise. A standard way of incorporating finite-size effects is also discussed. In Section 3, we analyze the stability of our model and examine the impact of self-regulation and finite-size effects on the attractors’ properties. Finally, Section 4 summarizes our results.

In our recent study [39], we utilized a system of coupled BNs to model the collective dynamics of stem cell populations. Through numerical analysis, we discovered that, by adjusting the level of intrinsic noise, the system undergoes a series of symmetry-breaking events that can be interpreted as cell differentiation. Interestingly, these symmetry-breaking events exhibit characteristics similar to order–disorder phase transitions in statistical mechanics. Here, we aimed to simplify our model by implementing Landau’s phenomenology of phase transitions. In the mean-field approximation, the order parameter m evolves under the influence of the following dimensionless potential function: where is the intrinsic noise of the system and , and . The time evolution of m is dictated by . This approach was originally developed to explain order–disorder phase transitions that occur by varying the temperature () of systems. In such systems, signifies a disordered state, while describes ordered states. As the temperature is lowered below a critical point, the system undergoes a phase transition from a state of high symmetry (disorder) to a state of lower symmetry (order). Note that this approach is applicable only in the thermodynamic limit, i.e., infinitely large systems. In the case of finite systems, it is possible to assume that the system is driven by the same potential function, but is subject to random fluctuations that are inversely proportional to the size of the system, N. In this approximation, the time evolution of the order parameter is given by the following stochastic differential equation: where t represents time, is the standard deviation of the random fluctuations, and is the standard Brownian motion. In the thermodynamic limit, where , the differential Equation (2) can be studied through the lens of bifurcation or phase transition theory. A bifurcation is described as a sudden change in the stability of the system; either new equilibrium states emerge, or the stability of existing states changes. On the other hand, a phase transition occurs when the global minimum of the potential function (ground state) shifts from one state to another. This shift may be continuous (second-order phase transition) or discontinuous (first-order phase transition) [21]. As we will see below, the bifurcation and phase transition diagrams may not necessarily coincide. In this work, we drew the analogy from traditional phase transitions and allowed to represent the strength of gene expression noise and m to define the average cell type of stem cell populations. Specifically, we considered the case where corresponds to pluripotent stem cells, while describes differentiated cells. In what follows, we used to denote the pluripotent state. Let us briefly discuss the results of three distinct models through the prism of this analogy [40]. The bifurcation and phase transition analysis presented below is based on four critical values of the intrinsic noise, , , , and . Model 1: Here, we set , , and . For , the pluripotent cell state () is both stable and the ground state of the potential function (solid blue line in Figure 1a). Below the critical noise , the pluripotent state becomes unstable, and two new stable cell states, , with the same energy emerge (solid yellow line in Figure 1a). Since both new cell types are energetically equally favorable, the ultimate cell fate is decided spontaneously. The change of m as a function of is continuous, and as a result, new cell states are reversible under moderate changes in intrinsic noise [23]. This bifurcation class is called supercritical pitchfork bifurcation, and the phase transition is of second order. In this case, the bifurcation and phase transition diagrams coincide (see Figure 1d). Model 2: In this model, , , and . The bifurcation analysis of this model is based on and . The behavior of the potential function is presented in Figure 1b. For , only the pluripotent cell state, , exists (solid blue line). At , the system undergoes a double-saddle-node bifurcation (black dashed line). For , two new stable and two unstable cell types emerge (solid yellow line). By reducing the noise below , disappear, and the pluripotent state becomes unstable (red solid line). The complete bifurcation diagram is presented in Figure 1e. Interestingly, this type of model exhibits a hysteresis loop that may be attributed to phenotypic memory (see blue arrows in Figure 1e). In other words, new cell type formation is usually irreversible under moderate changes of intrinsic noise [23]. In this model, the ground state of the system shifts from to at (dotted black line in Figure 1b). This shifting occurs discontinuously, and thus, the phase transition is of the first order. At this critical point, the pluripotent and differentiated cell states co-exist. Since both are energetically equally favorable, the differentiation process is spontaneous in this case as well. This model belongs to the saddle-node bifurcation class and is accompanied by a first-order phase transition. As seen in Figure 1d, the bifurcation and phase transition diagrams are notably different. One may say that first-order phase transitions are preceded by saddle-node bifurcations. Model 3: This model utilizes the values of , , and . The bifurcation analysis is based on the same critical values of intrinsic noise as in Model 2. Similar to the previous two models, several plots of the potential energy are presented in Figure 1c. When exceeds , only the pluripotent state, , exists (solid blue line). At , the system undergoes a single-saddle-node bifurcation, resulting in the emergence of two equilibrium points (dashed black line). Among the two, only is stable, indicating the emergence of a new cell type. Decreasing the intrinsic noise at triggers a transcritical bifurcation, causing the stability of and to switch (dashed–dotted black line). When falls below , both become stable, leading to the emergence of two new cell types (solid red line). The complete bifurcation diagram is presented in Figure 1f. In this model, the ground state of the system shifts discontinuously from to at a critical point of (first-order phase transition). Even after becomes stable, is always the energetically most-favorable. Notably, this type of model exhibits a hysteresis loop only for (see blue arrows in Figure 1f). Thus, the differentiation of into may be considered irreversible under moderate intrinsic noise changes. However, the transition to a stable is continuous, and as in Model 1, it can be reversed relatively easily. This model also belongs to the saddle-node bifurcation class and is also accompanied by a first-order phase transition. Figure 1f outlines the main difference between the bifurcation and phase transition diagrams.

The primary purpose of this work was to incorporate a feedback mechanism in Landau’s phenomenology that can effectively buffer the intrinsic noise of the system and ultimately guide stem cell populations toward a critical state. In this critical state, the system is empowered to determine its own fate through spontaneous symmetry-breaking events. To mathematically represent this process, we assumed that the evolution of is governed by the following equation: where G and M are the control parameters of our model. Without loss of generality, we assumed . It is worth noting that this represents a negative feedback mechanism, since an increase in m beyond M results in an increase in , which in turn reduces m (see Figure 1a,d). The system of Equations (2) and (3) is the dynamical system of interest in this work. As we will demonstrate below, this negative feedback mechanism gives rise to novel cell types.

In this section, we demonstrate the stability of the equilibrium points of the two-dimensional system described by Equations (2) and (3). Since represents the strength of intrinsic noise, we will restrict our analysis to . Each asymptotically stable equilibrium point (SEP) corresponds to a different cell type, with pluripotent stem cells represented by SEPs at () and differentiated cells represented by SEPs at . The eigenvalues of the Jacobian matrix at the equilibrium points are , where and are the trace and determinant of J, respectively. The stability of each point can then be determined using the Poincare diagram (refer to [45] for details). Moreover, if there exists an M such that and (non-hyperbolicity condition), the system may exhibit Hopf bifurcation, provided that the transversality and genericity conditions are also satisfied (see, for example, [46]). As we will see below, this model predicts the existence of stable limit cycles (LCs), which may be considered another possible cell type of our model. To demonstrate spontaneous symmetry-breaking events and cell fate decisions, stability analysis is complemented by representative simulations. The numerical solution of Equations (2) and (3) was obtained using the standard Euler–Maruyama method, with an integration step of . Distributions of the order parameter m were obtained over 1000 independent stochastic trajectories.

The model exhibits two equilibrium points , which, due to symmetry, have the same stability. The positiveness of the intrinsic noise imposes the restriction . The Jacobian of the system is given by Note that is always negative and is always positive. Thus, both are asymptotically stable. Depending on the value of G, the equilibrium solutions can be a stable spiral (for ), a stable node (for ), or a stable degenerate node (for ). Hence, is always less than , and as shown in Figure 2a (stable spiral) and Figure 2b (stable node), the stem cells will inevitably differentiate into either of the cell types. To demonstrate self-tuned symmetry breaking in the system, we conducted 1000 independent simulations with initial conditions of (pluripotent state), (high noise limit), and . Figure 2c,d show the time evolution of the probability density function (pdf), . The results indicate that the system spends a significant amount of time in the pluripotent state before spontaneously differentiating and settling at . One may reconstruct the Waddington epigenetic landscape by computing the quasi-potential ∼ vs. time.

The equilibria for this model are , which are also symmetric and have the same stability. As in the previous model, we imposed the restriction to ensure the positivity of the intrinsic noise. In this case, the Jacobian reads Interestingly, this model exhibits a supercritical Hopf bifurcation. Note that the only instance where vanishes and is positive is when . Moreover, both the transversality and genericity conditions are satisfied in a manner that the system displays a stable limit cycle. To summarize, for , correspond to stable LCs. At , the system undergoes a supercritical Hopf bifurcation, and transform into two stable fixed points. Specifically, for and , represent a stable spiral and stable node, respectively. In our analogy, the supercritical Hopf bifurcation converts two types of cells (LCs) into two other distinctly different types of cells (SEPs). Figure 3a,b demonstrate the Hopf bifurcation for with . Both start from a pluripotent state, , and a high internal noise of . In Figure 3a, the system exhibits a stable counterclockwise limit cycle for , and it transitions to a stable spiral for (Figure 3b). If the simulation had started with , the clockwise stable LC or stable spiral for would have been observed. Figure 3c,d present the probability density function (pdf) for and , respectively, and . Both simulations start with and . Similar to Model 1, the system spends a considerable time in the pluripotent state before differentiating into either the clockwise or counterclockwise LC (Figure 3b) and the left or right SEP (Figure 3d) with equal probability. It should be noted that, if the amplitude of the LC is elevated (Figure 3e) and is strong enough, direct transitions from one LC cycle to another may occur (Figure 3g). This situation is known as noise-induced phenotype switching (see, for instance, [42]). When the amplitude of the LC is considerably increased, an interesting scenario arises. Figure 3f depicts an LC that periodically approaches the pluripotent state. This simulation was conducted with and . It is evident that there are two modes in this case. Firstly, we see two slow modes: one at the pluripotent state (highlighted in blue) and another at a differentiated state (highlighted in green). Secondly, there are two fast modes when the trajectory quickly transitions between the pluripotent and differentiated state. Practically, the system continuously resembles the hysteresis loop between and illustrated in Figure 1e. It can be argued that this LC is a combination of the pluripotent stem cell () and a differentiated cell (). With the appropriate initial conditions, it is possible to observe the LC oscillating between and , as well. This scenario becomes even more intriguing when finite-size effects are taken into consideration. In this case, as depicted in Figure 3h, the system stochastically transitions between and . As a result, a hybrid state of three cells is formed.

This model exhibits two asymmetric equilibrium states . To maintain positive , . Here, the Jacobian is Let us first analyze the stability of . At , we have and , and both the transversality and genericity conditions are satisfied. Thus, this system also exhibits a stable LC for . However, at , the system experiences a supercritical Hopf bifurcation, causing the LC to turn into an SEP. More specifically, for , represents a stable spiral and for a stable node. Similar to Model 2, this bifurcation can be thought of as transforming one cell type (LCs) into another cell type (SEP). For the second equilibrium point , we note that for , both eigenvalues are negative, resulting in a stable equilibrium point. In general, according to this model, a pluripotent stem cell can spontaneously differentiate into two asymmetric cell types . The stability analysis of this model is illustrated in Figure 4. In particular, Figure 4b depicts the Hopf bifurcation for as we adjust the control parameter M from to , respectively. It is important to note that the stability of is not affected despite changes in M. The corresponding probability densities are displayed in Figure 4c,d with . The system was initialized with and . Similar to Models 1 and 2, the system remains in the pluripotent state for a period of time before eventually differentiating into . Increasing the amplitude of the LC yielded qualitatively different results compared to Model 2. If fluctuations due to finite-size effects are ignored, the system continues to switch quickly between two slow modes that represent and , as demonstrated in the hysteresis loop of Figure 1e. Therefore, for large systems, the system can function in a hybrid state of and . Note that the duration of time during which the system is in proximity to is notably shorter than the time it spends close to . In small systems where random fluctuations cannot be ignored, the system inevitably differentiates into cell type (see Figure 4f). The reason for this is straightforward: once the system approaches an asymptotic stable equilibrium point, it cannot return to the pluripotent state. This outcome is due to the asymmetry between the two equilibrium points presented in this model.

This work was based on three central hypotheses. Firstly, we assumed that pluripotent stem cells operate in a state of high symmetry that enables them to co-express a certain number of genes that are found in all of the subsequent progenitors. Conversely, progenitors can only express a limited number of those genes, and thus, they assume a state of lower symmetry. Similar symmetry breaking was hypothesized for the differentiation of the progenitors into committed cells. Therefore, cell differentiation can be viewed as a sequence of symmetry-breaking events. Secondly, we hypothesized that pluripotency is a collective property that emerges within stem cell populations. Thirdly, we postulated that the system possesses an internal mechanism that regulates transcriptional noise and drives the system through a critical state, allowing for spontaneous cell fate decisions. These three hypotheses allowed us to draw an analogy from statistical mechanics and assume that stem cell differentiation corresponds to symmetry-breaking events arising from a series of disorder–order phase transitions. Accordingly, disorder (or high-entropy states) represents pluripotent stem cells, while order (or low-entropy states) corresponds to differentiated stem cells. To effectively model these hypotheses, we expanded upon Landau’s mean-field theory of phase transitions, integrating a feedback mechanism to self-regulate the intrinsic noise of the system. Our methodology for investigating the dynamics of stem cell populations has the advantage of being simpler compared to other, more realistic models [25,26,27,28,29,39]. Moreover, as is the case with all mean-field models, it offers a straightforward way of understanding the underlying phenomenology of cell differentiation in terms of phase transitions. Finite-size effects were also incorporated by introducing a stochastic term, the strength of which is inversely proportional to the system size. Stability analysis of the two-dimensional dynamical system revealed the presence of several attractors (sinks and limit cycle) that correspond to different cell types. Of greater significance, the feedback mechanism has the capacity to navigate the system through various bifurcation points, thereby generating new attractors (differentiation) or switching attractors (phenotypic switch). Depending on the magnitude of the finite-size effects, stem cell populations may exist in hybrid states that encompass both pluripotent and differentiated cell types. Interestingly, our model exhibited a Hopf bifurcation, which can potentially describe the spontaneous emergence of collective oscillatory behavior in stem cell populations. Several future directions can be explored to expand our understanding of the suggested model. One of these is to conduct a more comprehensive analysis of the effects of random fluctuations on the stability of the attractors [47]. This project could involve investigating the existence of stochastic bistability and phenotypic switching [44]. A key question in this direction would be determining whether external noise bursts could trigger cell reprogramming into a pluripotent state [4]. Another promising avenue is to extend our model beyond the standard mean-field approximation by incorporating spatial dependence in the form of stochastic partial differential equations. This approach could potentially describe self-organized criticality, as demonstrated in [48,49]. It would be interesting to explore the applicability of this approach to self-organized criticality in gene expression dynamics [50,51,52,53].