Emergent Criticality in Coupled Boolean Networks

Early embryonic development involves forming all specialized cells from a fluid-like mass of identical stem cells. The differentiation process consists of a series of symmetry-breaking events, starting from a high-symmetry state (stem cells) to a low-symmetry state (specialized cells). This scenario closely resembles phase transitions in statistical mechanics. To theoretically study this hypothesis, we model embryonic stem cell (ESC) populations through a coupled Boolean network (BN) model. The interaction is applied using a multilayer Ising model that considers paracrine and autocrine signaling, along with external interventions. It is demonstrated that cell-to-cell variability can be interpreted as a mixture of steady-state probability distributions. Simulations have revealed that such models can undergo a series of first- and second-order phase transitions as a function of the system parameters that describe gene expression noise and interaction strengths. These phase transitions result in spontaneous symmetry-breaking events that generate new types of cells characterized by various steady-state distributions. Coupled BNs have also been shown to self-organize in states that allow spontaneous cell differentiation.


Introduction
The development of multicellular organisms depends largely on the ability of stem cells to self-replicate indefinitely (proliferation) and differentiate into specialized cells (pluripotency). Macroscopically, the outcome of this process is predictable in an almost deterministic fashion, i.e., they are almost always primed to functional cells in their lineage. Microscopically, stem cell dynamics appear to be stochastic due to the molecular nature of gene expression. The intrinsic noise causes noticeable cell-to-cell variability in stem cell populations [1]. Although the functional role of gene expression noise and cell variability remains elusive, many have suggested that gene expression noise and cell variability are integral to stem cell pluripotency. Huang et al. proposed that pluripotent (multipotent) stem cells are a "balanced, undecided state" of multiple gene expression patterns [2,3]. Noise and environmental signaling simply destabilize this state and force stem cells to overexpress only in certain genes, thus acquiring a new type of cell with specialized functions. Others have argued that this transition may only emerge or, at least, is significantly amplified in stem cell populations due to the interplay between cell-cell interaction, gene expression noise, and environmental signals [4][5][6]. Irrespective of the detailed mechanism, the general theory states that pluripotency is a state of higher symmetry (i.e., a higher amount of coexpressed genes), while differentiation leads to states of lower symmetry (lower amount of similarly expressed genes).
Symmetry breaking is the hallmark of disorder-order phase transitions in large physical systems of interacting agents [7]. The system transitions from a state of higher symmetry (disorder) to a state of lower symmetry (order) by reducing the intrinsic noise (control parameter) below a critical point. Examples of this type of phase transition include A Boolean network is defined on a set of n interacting nodes N = {1, . . . , n} representing regulatory genes. Each node i is assigned a Boolean value x i ∈ {0, 1} (inactive and active, respectively). Each gene i is regulated by k i ≤ n genes J i = {j 1 , . . . , j k i } ⊆ N. A synchronous update of the state vector x(t) = {x i (t)} i∈N at discrete time t is given by where y i (t − 1) = {x j (t − 1)} j∈J i is the state vector of the interacting nodes at the previous time and f i : {0, 1} k i → {0, 1} is the Boolean function. At any given time t, the state of BN is defined by the state vector x(t) or the decimal-encoded state, s(t) ∈ S, of x(t).
Here, S = {1, . . . , 2 n } is the decimal representation of the states indexed from 1 to N = 2 n possible configurations of the state space.
In reality, it is uncommon to have full knowledge of the specific gene connectivity and Boolean functions of a GRN. For this purpose, Kaufmann proposed an ensemble approach to study GRNs, where statistical averages of random Boolean networks (RBNs) are taken. In an RBN, for all i ∈ N (a), the number of connections is the same (that is, k i = k), (b) the interacting nodes are randomly chosen, and (c) for any input y i , f i (y i ) = 1 with probability bias of p and f i (y i ) = 0, otherwise. By varying the control parameters k and p, RBNs undergo a phase transition between order and chaotic network dynamics. The edge of chaos exists at the critical connectivity k c = 1/[2p(1 − p)] [19,20].
As the dynamics of a BN are deterministic in finite space, the network eventually entails repetition of gene states. These dynamical attractors are an important feature of complex systems and have long been established as gene expression patterns that characterize cell types [21][22][23]. One way to observe attractors of networks with a reasonably small number of genes is to numerically approximate the steady-state distribution through model simulations [24][25][26]. In a BN, Monte Carlo simulations can approximate the steadystate probability distribution of gene states by introducing a random perturbation to the synchronous update of the network [13]. A possible BN realization with perturbation is to update gene states according to Equation (1) with probability 1 − (1 − q) n and update with x i (t) = γ i , otherwise, where i ∈ N, γ i ∈ {0, 1} and q = P{γ i = 1}. We note that the original formulation of random perturbation in [13] is a binary "flip" of a randomly selected gene. Here, we perturb by updating the entire gene set x with a random vector of binary values {γ 1 , . . . , γ n }. With perturbation, there is a nonzero probability of arriving at any state, and thus the Boolean network is an ergodic Markov chain where the states converge given sufficient time [27]. We note that, although random perturbation is used here to require ergodicity of the dynamics, it also mimics the gene expression noise observed in biological systems. For example, in Escherichia coli, the gene expression noise modeled as stochastic dynamics has been well studied to be necessary for its regulation, fluctuation of transcription rate, and cell division [28][29][30][31]. Figure 1 shows the wiring diagram, truth table, state transition diagram, and steadystate distribution of a 6-gene Boolean network with connectivity of k = 2 and bias p = 0.5 with the perturbation probability of q = 0.1. Throughout this work, we use this model 6-gene network for demonstration.

Dynamic Control Kernel
The control of gene expression is a highly desirable action in many areas of molecular biology. For example, targeted drug delivery or the knockdown effect in gene therapy requires persistent external influence on a cell to achieve desired gene expression. In the context of BNs, there are numerous different control strategies. Shmulevich et al. initially proposed the optimal intervention strategy of probabilistic Boolean networks by altering Boolean functions [13]. For this type of intervention, one is often interested in finding the best candidate genes to intervene by minimizing the mean first-passage time. However, in recent work, direct modification of gene states to drive the dynamics to the desired attractors has been proposed [32]. Here, we focus on a specific notion of control, pinning, also known as node-state override, where the gene state is permanently fixed. Pinning differs from the original formulation of intervention in that Boolean functions are preserved for all genes except those that are kept static. This concept was first introduced by Serra et al. and was described as the knock-out of genes with an avalanche effect [33]. Kim et al. formalized the control kernel (CK) of a network as the minimal set of genes whose pinning reshapes the dynamics such that the basin of the attractor becomes the entire configuration space [34]. In a pinned network, only a small fraction of the total number of network nodes depends on the topological and logical characteristics of the network. Experimentally, Joo et al. have explored the control of a single gene as an external input for a 5-gene network that captures the molecular mechanisms and the cell-state transition from an epithelial to a mesenchymal stem cell [35].
Let r ≤ n be the cardinality of a CK set. Without loss of generality, we can assume that the CK nodes are the first r nodes of a BN. We define the state vector of the CK and the remainder network as x c = {x 1 , . . . , x r } and x o = x − x c , respectively. The set of decimalencoded 2 r states of CK is W r = {1, . . . , 2 r }. We denote the ordered sets of the pinned values of CK as X w , where w ∈ W r . For example, X 1 = {x 1 = 0, x 2 = 0, . . . , x r−1 = 0, x r = 0}, . , x r−1 = 0, x r = 1}, etc. Then, the dynamics of an RBN are as follows: We reserve the index w = 0 to describe the unpinned dynamics, that is, In other words, when w = 0, Equation (3) is Equation (1). For each w ∈ W r , the pinning procedure generates a new steady-state distribution g w (s)} s∈S that is attributed to a new type of cell. Thus, each CK set can generate up to 2 r new types of cells. An interesting property of this control approach is that it naturally partitions the state space, S, into 2 r equivalent sets. As a consequence, for i = j, we have the orthogonal condition where g j (s) is the scalar product of g (r) i and g (r) j . This property is illustrated in Figure 2, where we applied the pining procedure to the 6-gene model BN.
In the case of r = 1, pinning the CK to either X 1 or X 2 generates two nonoverlapping steady-state distributions g 3 (s), and g 2 (s), g 3 (s), and g Here, we define dynamic CK as x c (t) = X η t , where η t is a discrete-time deterministic or stochastic process with state space W ⊆ {{0}, W r }. Let g(s) be the new steady-state distribution as a result of X η t . If the mean waiting time in a given state of η t is finite and long enough, the distribution of the cell states will be a mixture of all g w s, i.e., where c w s are the scalar weights of the respective g w s.
To illustrate this point, we set r = 1 and assume that η t ∈ W 1 = {1, 2} is a discretetime Markov chain process with the transition diagram shown in Figure 3a. Figure 3b shows the time evolution of η t for the first 100 time steps. It is clear that η t = 2 is sampled more than η t = 1. The resulting steady-state distribution g(s), shown in Figure 3c, is a mixture of g  2 (s) with the corresponding weights c 2 > c 1 . Here, we interpret this as that the new cell type g(s) resulting from stochastic pinning is a weighted mixture of cell types g (1) 1 and g (1) 2 . This is true as long as the transition rates are not "too fast." In Section 2.4, we take advantage of the disjoint property of the g w s to develop a spectral decomposition method that provides these weights. Other methods, such as linear optimization and non-negative matrix factorization (NNMF), are discussed in the case where g

Coupled Boolean Networks
The study of coupled BNs has gained increased interest in recent decades. Villani et al. and Serra et al. have shown that increasing interactions between RBNs results in more disordered states [36,37]. Damiani et al. revealed that short-distance interacting RBNs display robust generic properties [38]. The collective behaviors of coupled BNs have been investigated in [26,39,40]. In particular, [39] showed that the pattern formation in the tissues of BNs is the most information-rich in the near-critical complexity domain, while in [26], the authors showed that the pattern formation of long-term steady states is most often observed in networks of critical dynamics. In the most recent study by Kim et al., multilayer RBNs were investigated, where isogenic GRN were coupled according to a random selection of topology in silico with activation rules from [39]. They showed that a multilayer RBN structure facilitated the production of antifragile systems [40].
In this work, we model stem cell populations as L-coupled isogenic BNs (tissue). Let x i,m represent the m-th gene of the i-th cell. Each cell i is allowed to interact with a set of other cells denoted by Γ i . For each Boolean variable σ ∈ {0, 1}, we define a linear transformationσ = 2σ − 1 ∈ {−1, 1}. We assume that BNs interact with each other through their CK according to the multilayer Ising Hamiltonian: In terms of a standard spin model, the first term of Equation (5) is the interaction between cells i and j, where J ij is the strength of the interaction. It is simply the Hamiltonian of r-independent Ising models. The second and third terms represent local external fields acting on the control kernel, where h and h 0 are the corresponding strengths. In our coupled BN model, the first term describes paracrine signaling with neighboring cells. The second term describes the cell's tendency to follow the dynamics of its original cell (f i,m ). By default, it is time-dependent and indirectly couples the r Ising models. The third may describe autocrine signaling and/or external interventions (ψ m ) and can be static or time-dependent.
Cell-cell communication or cell interaction with external signals is facilitated by diffusing signaling molecules. To capture this stochastic process, we introduce an additional type of intrinsic noise, denoted as T, which acts at the population level. The analog of T in standard statistical mechanics is the temperature of the system. To simulate the system, we use the standard Metropolis algorithm for the coupled CKs and the Equation (2) for the synchronous BN update. During an MC step, the N MC ≤ L subset of nodes, x i,m , for each 1 ≤ m ≤ r, is chosen randomly. We sequentially flip their state. The new state is accepted with probability min(1, exp(−∆E/T), where ∆E is the energy difference between the current and the attempted state. After each MC step, all BN states x o are updated.

Detecting Cell Types
In Section 2.2, we illustrated that dynamic CKs can generate a new steady-state distribution as a mixture of the distributions. Let the vector G(t) be the instantaneous distribution of the states of all cells in a population and be the approximation given by Equation (4) In the case h 0 = 0, the population dynamics are dictated by the network of interacting CKs and external intervention. Thus, all CKs are expected to be in a state w ∈ W r . In other words, c 0 = 0. This assumption can significantly simplify the computation of coefficients. Taking advantage of the nonoverlapping (disjoint) property of g w , we define a set of orthonormal vectors asĝ w . Then, one can use a simple spectral decomposition method to obtain the coefficients: If the error (t) is not sufficiently small, the assumption of Equation (6) breaks down. In such cases, machine learning algorithms can be implemented to decompose G and gain insight into the dynamics of the system. Here, since the elements of G are non-negative, we use non-negative matrix factorization (NNMF) [41]. In short, NNMF can decompose a non-negative matrix G N ×T into a product of two non-negative matrices Ω N ×K and E K×T . Here, T is the simulation time, and K is the rank of the decomposition to be determined. The matrix Ω N ×K provides K-different steady-state distributions, which can also be interpreted as new types of cells. The corresponding time-dependent density numbers are stored in E K×T .

Results
This section presents the results of our numerical simulations based on the hybrid Monte Carlo and synchronous BN update method. We assume that we have a simple 2D tissue with periodic boundary conditions. The set Γ i consists of the four nearest neighbors of the ith cell. For simplicity, we choose J ij = J ≥ 0.
Some interesting limits of the multilayer Ising model include:

Case 3.
For h 0 ≈ J, the r Ising layers are coupled through the dynamics of each BN. This is a nontrivial, bidirectional, and time-dependent nonlinear coupling.
In this work, we set J = 1 and neglect the third term in Equation (5). Thus, the control parameters are intrinsic noise ("temperature") T and local coupling h 0 . For each set of parameters, we first equilibrate the system for t eq time steps and then compute the time average of an observable A(t) as A = ∑ T t=1 A(t)/T , where T is the sampling time. The time average of the numerical error is always less than 10 −2 . In Section 3.1, we use the linear optimization technique to compute the coefficients of Equation (6), while in Section 3.2, where we assume that h 0 = 0, we implement the spectral decomposition method from Section 2. 4.
In what follows, we present the results for r = 1 or r = 2. Unless otherwise stated, T = t eq = 10 4 . The parameters values are summarized in Appendix A.

Spontaneous Cell Differentiation
This section examines the behavior of the system as a function of the control parameters T and h 0 . Specifically, we show that varying the control parameters causes the system to undergo a series of spontaneous symmetry-breaking events. It should be noted that many biological phenomena, such as the stem cell differentiation process, are precisely explained by this phase transition. For example, the first step in the core embryonic stem cell cycle is the organization of the pluripotent state in cells [42]. Upon dissolution of pluripotency, stem cells reach a critical state, at which point they undergo symmetry breaking. The stem cell population spontaneously (and collectively) undergoes a cell fate specification for differentiated specialized cells.
In the first set of simulations, shown in Figure 4, we fix r = 1. In Figure 4a-c, we kept the temperature constant at T = 0.25, T = 1, and T = 1.25, respectively, and varied the parameter h 0 . Specifically, we started with a high value of h 0 = 3, where we expected the fraction of the original cell c 0 to be dominant and gradually reduced it to h 0 = 0. For h 0 = 3, the initial conditions of all nodes in the tissue were randomly chosen. The final state of this simulation is the initial condition for the next value of the local field h 0 , and so on. For high values of h 0 , we see that c 0 is the most dominant fraction in the colony. As h 0 decreases, c 0 drops to a critical value where spontaneous differentiation occurs. An arrow in Figure 4 approximately indicates the starting point of the differentiation process. Note that the system always differentiates to the g (1) 2 cell. This decision is related to the structure of the original BN. As seen in Figure 1d, the steady-state distribution is slightly higher in the second half of the state space, and this small difference always favors g (1) 2 . However, different BNs may exhibit alternative phenotypic behaviors. This transition is smooth for high temperatures (T = 1 and T = 1.25). However, the differentiation seems rather steep for low temperatures (T = 0.5). In Figure 4d-f, we fixed the local field constant at h 0 = 0, h 0 = 0.5, and h 0 = 1, respectively, and simulated with various temperature points from a high value of T = 3 to a low value of T = 0.1. Here, we also observe a differentiation at a critical temperature indicated by the arrow, which consistently rules in favor of g (1) 2 for nonzero h 0 . However, for h 0 = 0 (Figure 4d), the original cell has no influence on the dynamics of the system, and the cells can differentiate to g   1 (t) , and the brackets indicate the ensemble average, it can be shown that the system belongs to the same universality class as the 2D Ising model.
In Figure 5, we performed the same simulations as in Figure 4 but for r = 2. In this case, the BN has four potential steady states to differentiate. As observed in all the graphs in Figure 5, two consecutive cell differentiations exist; the arrow indicates the first sequence of differentiation, and the dashed arrow indicates the second. First, there is a differentiation between g (2) 1 , g (2) 2 and g (2) 3 , g (2) 4 , where the second group is always favored for the same reasons mentioned in Figure 4. Then, there is a second spontaneous differentiation in which the cell decides equally between the g (2) 3 and g (2) 4 cell types. When h 0 = 0, the colony can spontaneously differentiate into any of the four types of cells. Figure 5d illustrates a simulation in which the cell differentiates into g (2) 1 . Next, we illustrate the cell-to-cell variability in the system. As an example, we used the parameters of Figure 5e and plotted the time evolution of the cell fractions, as well as representative snapshots at three different temperature points: T = 3, T = 1.85, and T = 1 (Figure 6). At temperature T = 3, precisely before the first differentiation, all types of cells co-exist throughout the simulation time (Figure 6a). The snapshot in Figure 6b illustrates this high level of cell variability. At T = 1.85 (Figure 6c), where the second differentiation occurs, the system fluctuates significantly between g (2) 3 and g (2) 4 . This fluctuation resembles the features of the second-order phase transition observed in 2D Ising models. In support of this assumption, the snapshot of Figure 6d shows the characteristic fractal-like structure of g (2) 3 island sizes in a sea of g (2) 4 cells. At low temperature T = 1, the second differentiation has already taken place, and the decision to differentiate to cell g (2) particular simulation (Figure 6e,f). It must be underlined that, for this specific BN, the original cell g 0 can differentiate equally between the g (2) 3 and g (2) 4 .   ((a,b)), T = 1.85 for the second column ((c,d)), and T = 1 for the third column ((e,f)). The time unit is 10 2 MC steps.
A second source of cell variability is due to the size of the colonies. For colonies of small size, there is always a measurable probability of having the first two types of cells, although the structure of the original cell does not favor this decision. This probability decreases significantly as the size of the tissue increases and almost vanishes for sizes greater than 20 × 20. This is due to weak ergodicity breaking that makes it difficult for the system to visit the first half of the state space. This theoretical interpretation may contribute to the cell-to-cell variability observed in small stem cell colonies.
Phase transitions can be studied in terms of the order parameters c ww = c w − c w , where w, w ∈ {0, W r }. Note that for each r there are M = m!/(m!(m − 2)!), where m = 2 r + 1 unique fraction differences. Although this computation is manageable for r = 1, it becomes increasingly difficult to solve for r ≥ 2. In an ongoing project, we attempt to classify the phase transitions presented in Figures 4 and 5 and compute the corresponding critical exponents. Preliminary results indicate that partial construction of the Waddington epigenetic landscape [43] is possible. Specifically, one can derive the quasi-potential V({c ww }) ∝ −ln(P({c ww })), where P({c ww }) is the joint probability distribution of all {c ww } [44,45]. Note that this scalar function is defined in an M-dimensional space. This construction can be simplified considerably in cases where some transitions are not possible.
For example, in the model BN for r = 2, the transitions from g (1) 1 to g (2) 3 and g (2) 4 and from g (1) 2 to g (2) 3 and g (2) 4 were not observable in the limited number of simulations performed for this section.

Self-Tuned Cell Differentiation
In Section 3.1, we showed that by manually tuning T and h 0 through the critical values, the system undergoes a series of symmetry-breaking events. In this section, we demonstrate that cells can collectively self-tune through a critical state that allows them to decide their fate. For simplicity, we set r = 1 and neglect the local fields in Equation (5). This limit is equivalent to the standard Ising model. We recently showed that a mean-field approach of the Ising model with a negative feedback mechanism drives the system through a supercritical pitchfork bifurcation that can be interpreted as a cell fate decision [46]. Here, we apply this approach to the full Ising Hamiltonian Equation (5).
The heterogeneity of the system can be measured in terms of instantaneous "magneti- 1 and g (1) 2 (M t = 0) corresponds to the pluripotent state of the cell [2]. On the contrary, M t = −1 and M t = 1 correspond to homogeneous populations of cell types g (1) 1 and g (1) 2 , respectively. We consider an internal mechanism that allows the heterogeneity of a population, measured in terms of instantaneous magnetization of the tissue, M t , to regulate the intrinsic noise of the population (temperature T): where α is the relaxation coefficient. That is to say, Equation (8) captures the negative feedback response between cell-cell cooperativity and its intrinsic gene expression noise. The model Hamiltonian for coupled Boolean networks (Equation (5)) with h = 0 and h 0 = 0, combined with the internal temperature-magnetization feedback mechanism (Equation (8)), was simulated on tissues of size 32 × 32. The system was set to evolve for up to t = 8 × 10 5 MC steps, where N MC = 10. Equation (8) was solved using the Euler method with a step size of ∆τ = 5 × 10 −7 , and the parameter α = 0.8 was fixed. The model Boolean network from Figure 1 was instantiated for all cells with internal noise of q = 0.02.
The simulation of the model begins with random gene states for all cells in the tissue, except for the CK node, which, without loss of generality, is set: x c = {1}. We choose a high starting temperature of T = 2.8 for the system, because it generates a natural state of hypothesized heterogeneity in pluripotent cells in the early stage of the embryonic stem cell cycle. Through negative feedback on instantaneous magnetization, Equation (8) then self-tunes the system towards the critical and then subcritical temperatures, where symmetry breaking triggers spontaneous differentiation. Figure 7a shows the time evolution of magnetization (M t ) for two independent simulations. Here, it can be seen that both simulations begin with M t = 1, which quickly approaches 0 over time. As time passes, the simulations decide and split in their magnetization paths (M t ≈ 1 and M t ≈ −1), resulting in differentiation of population cell types (which correspond to g (1) 1 and g (1) trajectories of the two simulations. Here, the trajectories begin at T = 2.8, and with time, the temperature drops below the critical temperature and eventually equilibrates to a subcritical temperature. Combined, we see that as the temperature reaches the Ising critical temperature (T c ), the tissue magnetizations diverge, with an equal chance of the system choosing one of the two cell types. Observing tissue-level statistics provides additional insight into collective behaviors in pluripotent cells transitioning to two possible cell types. A total of 100 independent and identical tissue simulations of Figure 7 were carried out, and instantaneous state distributions G(t) and the average gene state (mean G(t) of a whole tissue) were collected at each time step. At the beginning of the simulation, t = 1, the system was initialized with a high temperature point (T = 2.8) and with the CK fixed to x c = {1}, as in Figure 7. Trivially, all cell states are in g (1) 2 , and the average gene state, which describes the tissuelevel distribution, is unimodal (Figure 8a). At time t = 1 × 10 5 , where the system reaches M t ≈ 0, two different cell types are probed, resulting in a mixture of g (1) 1 and g (1) 2 cell types, while the average state at the colony level remains unimodal, centered between g (1) 1 and g (1) 2 (Figure 8b). With time, the temperature self-tunes and reaches an equilibrium point below the critical temperature (T c ). At time t = 6 × 10 5 , approximately half of the tissues form a homogeneous cell type of g (1) 1 , and the other half form a cell type g (1) 2 . The tissue-level states reach a split bimodal distribution (Figure 8c). The system describes a full transition from a population of pluripotent tissue to two differentiated cell types. At the tissue level, this unimodal-bimodal transition at the critical junction of the phase transition occurs in several areas from mouse embryogenesis [12] to the development of cancer cell line [47].  As an application, self-tuned differentiation can replicate Okamoto et al.'s experimental work on the collective differentiation of mouse embryonic stem cells (mESCs) under strict conditions [12]. The authors have observed the gene expression levels of key transcription factors in mESC, Nanog, and Oct4 in the early stage of differentiation. According to the immunofluorescence markers of Venus and mKate2, which report Nanog and Oct4 gene expressions, respectively, colonies of mESCs exposed to leukemia inhibitory factor (LIF) demonstrated a high intensity of fluorescence, thus exhibiting single-state behavior in Nanog and Oct4. Here, LIF acts to enhance Nanog heterogeneity, in other words, to maintain the pluripotent state in the stem cell population. However, in the absence of LIF, high and low levels of Venus and mKate2 fluorescence were observed in cells, and cells are free to transition from pluripotent to differentiated state. This is termed unimodal-bimodal transition in the heterogeneity of gene expression levels at the cellular and colony level. This cellular and colony-level unimodal-bimodal transition in the mESC distributions of a pluripotent population [6,48,49] is characterized by the multilayer Ising Hamiltonian (Equation (5)) with the temperature-magnetization feedback mechanism (Equation (8)) when pluripotent cells and differentiated cells are assumed to be complementary in gene states (i.e., they are induced by a fixed CK of r = 1). Then, pluripotent and differentiated cell types form g (1) 1 and g (1) 2 , which can exhibit the unimodal-bimodal transition with the self-tuning mechanism, as seen in Figure 8.

Conclusions
In this paper, we developed a model that describes the interplay between stem cell cooperativity and gene expression noise at the population level. Our model uses isogenic BNs to represent individual cell dynamics and a multilayer Ising Hamiltonian to describe the cell-cell interaction. This approach captures various cell signaling effects (paracrine and autocrine) and external gene expression interventions. The time evolution of this model is obtained by a hybrid MC method and synchronous update of BN states. We showed that a mixture of steady-state distributions can accurately represent the collective dynamics of coupled BNs. By interpreting each steady-state distribution as different types of cells, we characterized the compositions of stem cell populations with the aid of linear optimization techniques, standard spectral decomposition methods, and unsupervised machine learning algorithms.
Numerical analysis of our model in two dimensions revealed that by varying the system parameters through some critical values, the system undergoes a series of symmetrybreaking events that alter the dynamics of the original cell. Before these events, the system exhibits considerable cell-to-cell variability. This property resembles signatures of disorderorder phase transitions. We hypothesize that such transitions can be interpreted as the differentiation of the original cells into specialized cells. This interpretation supports the hypothesis that stem cells operate in a highly symmetric state (disorder state), and differentiated cells have a reduced symmetry state (order state). Furthermore, we showed that by introducing a differential equation that describes the negative feedback between cell cooperativity and intrinsic noise, the system can self-tune through the critical points and spontaneously differentiate into various types of cells. The number of different types of cells that each BN can generate strongly depends on the complexity of the cell-cell interaction and the structure of the original BN.
It is known that BNs can exhibit ordered, critical, and chaotic dynamics. How do the different BN dynamics affect the number and type of symmetry-breaking events observed in our model? Is criticality essential at the population level, locally at a single-cell level, or both? [36]). In this work, we showed that our model generates a number of different cell states. Depending on the structure of the BN, some states are stable and others are metastable. Preliminary results show that metastable cell states can emerge in certain conditions in cell tissues. What is the biological interpretation of such metastabilities? Here, we only considered ferromagnetic cell-cell interaction (i.e., J > 0). Can a different type of interaction (e.g., antiferromagnetic or random field) lead to other types of phase transition? Is a spin-glass phase transition possible in our model, and how could this affect collective cell dynamics? These questions, along with the partial construction of the epigenetic landscape and characterization of the phase transitions, will be studied in detail in a future publication [46].

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Table of parameters, definitions, and parameter values utilized in the model development and simulations. There are two control parameters: the intercellular intrinsic noise (T) and the tendency of the system to retain its original BN dynamics (h 0 ). Both parameters are measured in J units.