Cell-Fate Determination from Embryo to Cancer Development: Genomic Mechanism Elucidated

Elucidation of the genomic mechanism that guides the cell-fate change is one of the fundamental issues of biology. We previously demonstrated that whole genome expression is coordinated by the emergence of a critical point at both the cell-population and single-cell levels through the physical principle of self-organized criticality. In this paper, we further examine the genomic mechanism that determines the cell-fate changes from embryo to cancer development. The state of the critical point, acting as the organizing center of the cell fate, determines whether the genome resides in a super- or sub-critical state. In the super-critical state, a specific stochastic perturbation can spread over the entire system through the “genome engine”, an autonomous critical-control genomic system, whereas in the sub-critical state, the perturbation remains at a local level. The cell-fate changes when the genome becomes super-critical. We provide a consistent framework to develop a time-evolutional transition theory for the biological regulation of the cell-fate change.


Introduction
A mature mammalian somatic cell can reprogram its state and consequently acquire a very different gene-expression profile through exposure to a few reprogramming stimuli [1]. Such a drastic state change involves the coherent on/off switching of thousands of functionally heterogeneous genes [2]. However, there are fundamental physical difficulties in achieving such large-scale coordinated control on a gene-by-gene basis. These difficulties become more evident in a situation where there is (1) a lack of a sufficient number of molecules to reach a stable thermodynamic state (i.e., breakdown of the central limit theorem) and (2) a consequent stochastic noise due to the low copy numbers of specific gene mRNAs, which induces a substantial instability in genetic product concentrations [3,4].
In our previous studies [5][6][7][8][9][10], we demonstrated that the self-organization of whole-genome expression constitutes a "physically motivated" alternative to gene-specific regulation at both the population and single-cell levels. The mechanism of self-organization, through global genome reprogramming, eliminates the need for physically unfeasible gene-by-gene control of expression.
The core of the self-organization mechanism is the presence of massive system changes elicited by "apparently minor" external causes. To address this problem, Bak and colleagues [11] proposed self-organized criticality (SOC; the Bak-Tang-Wiesenfeld sandpile model). SOC is a general theory of complexity that describes self-organization and emergent order in non-equilibrium systems Our previous studies (see details in [7]) from embryo to cancer development have demonstrated that whole-genome expression is dynamically self-organized through the emergence of a critical point, with the co-existence of three distinct response domains (critical states). Furthermore, dynamic pictures of between-state flux provide a potential universal mechanism of self-organization interpreted in terms of a "genome engine". An autonomous critical-control genomic system is developed through the highly coherent behavior of low-variance genes (sub-critical state), which in turn, generates a dominant cyclic expression flux with high-variance genes (super-critical state) through the nuclear environment of the cell. This is evident given that gene expression is sorted and grouped according to the temporal variance of expression (nrmsf : normalized root mean square fluctuation), which acts as an order parameter for the self-organization of whole-genome expression [6,7]. On the contrary, randomly shuffled gene expression exhibits no evidence of cooperative behavior (no emergence of coordinated activity).
In this report, we further investigate time-series of whole-genome expression data for the cell-fate change: cell differentiation (heregulin (HRG)-stimulated MCF-7 human breast cancer cells compared with non-differentiated epidermal growth factor (EGF)-stimulated MCF-7 cells [19,20]; all-trans retinoic acid (atRA)-and dimethyl sulfoxide (DMSO)-stimulated HL-60 human leukemia cells differentiated to neutrophil cells [21]; differentiation of T helper 17 (Th17) cells from naïve T helper (Th0) cells [22]) and reprogramming (mouse and human early embryo development, [23] and [24], respectively). It is worth noting that these cases must be considered examples of a general mechanism and have been selected based on the availability of a sufficient number of experimental time points for the nrmsf value. . . . Figure 1. Schematic representation of the genomic mechanism for the cell-fate change. (1)(2)(3)(4) The critical point (CP) corresponds to the center of mass (CM) and represents a specific set of critical genes acting as a genome attractor. The CP has activated (ON) and deactivated (OFF) states. ON/OFF switching of the CP state occurs through switching of its singular behaviors, i.e., change in the critical transition of a specific critical gene set. Changes in the state of the CP, such as ON-OFF or OFF-ON switching, guide the genome into a super-critical state; thus, these perturbations can spread over the entire system in a highly cooperative manner. Chromatin remodeling plays a role as the material basis of the self-organized critical (SOC) control of genome expression [10]. (5) Due to the CP acting as a genome attractor, the self-organization of gene expression develops an autonomous critical-control genomic system (genome engine) through the formation of dominant cyclic flux between local critical states (distinct expression domains according to the degree of nrmsf : Section 4.2.1), where the local sub-critical state is the generator (see details in [7]). Coherent perturbation of the genome engine through changes in the CP (ON to OFF, OFF to ON, etc.) drives the cell-fate change. Before the cell-fate change, the genome (expression) system passes through a non-equilibrium fixed point (stable point of the thermodynamically open system). These five points support the development of a time-evolutional transition theory of biological regulation. Throughout this report, the whole-genome expression vector c(t j ) at t = t j represents the CM expression vector, in which the center of mass of whole-genome expression, CM(t j ) is subtracted from each expression (Section 4.2.3).
(1) The CP corresponds to the center of mass (CM(t j ): the average value of whole-genome expression at t = t j ) (Section 2.1). (2) The dynamics of the center of mass of any stochastic expression converge to that of whole-genome expression (i.e., CM(t j )). This shows that the CP is the genome attractor (Section 2.2). (3) The switching singular behaviors of the CP transform the genome into a super-critical state (i.e., super-critical genome). This in turn induces a "global expression avalanche", which is revealed by (4) the probability density function (PDF) of whole-genome expression. (5) Lastly, a cell-fate change occurs after the genome becomes a "super-critical genome" (Sections 2.3-2.6) and after the genome passes over a stable point (non-equilibrium fixed point) of the thermodynamically open system. This passing indicates symmetry-breaking which induces coherent perturbation of the genome engine (Section 3.1).
These five points constitute a framework for developing a time-evolutional transition theory of biological regulation.

Fixed Critical Point: A Specific Group of Genes Corresponding to the Center of Mass of Whole-Genome Expression
The existence of a critical point (CP) is essential for determining distinct response domains (critical states) [7]. To generate a unified model of biological regulation, we must go in depth into specific features of the CP in terms of "sandpile criticality" (Figure 2A). Sandpile criticality emerges when whole-genome expression is sorted and grouped according to the fold-change in expression between     of gene expression. Note: In (C), the removal of all genes with zero expression value at a specific cell state in RNA-Seq data does not change their correlation behavior (similar behavior is also observed in the mouse embryo case). Furthermore, the higher correlation in the embryo case makes the convergence to the CM, shown in (D), slower than in the cancer case (B). Panels (B,D) illustrate, in terms of Euclidean distance (y-axis; refer to the convergency of critical-state attractors: Section 2.6), the convergence of the dynamics of the CM of a randomly selected group of gene expression to that of the CM of whole-genome expression as the number of selected genes increases (average of 200 repetitions). The CP, the CM of whole-genome expression according to nrmsf, acts as the genome attractor (see Section 2.6 for heteroclinic attractors). Note: For the rest of the biological regulations throughout this report, similar behaviors are observed. The dashed lines in (A,B) indicate that coherent behaviors emerge at n = 50 single-gene expression values. Refer to experimental time points in Section 4.1.
To give a proof-of-concept of CSB, we performed a bootstrap simulation to capture two basic signatures of CSB: (1) The stochastic behavior of gene expression shows (relatively low) correlation convergence for randomly selected gene ensembles as the number of elements (n) is increased ( Figure 3A). (2) The CM of randomly selected gene ensembles (with hundreds of repetitions for the gene ensembles) must dynamically converge to that of whole-genome expression as the number of elements (n) is increased ( Figure 3B).
This happens at both the population and single-cell levels ( Figure 3C,D). The existence of a threshold n at around 50 randomly picked genes [6,25] allows us to reproduce CSB with a random choice of N genes with N > n, which is further proof of the reliability of CSB. The constancy of the minimum number of randomly selected genes to reach convergence ( Figure 3) is remarkable and suggests the presence of a sort of "percolation threshold" reminiscent of the size of the genetic networks operating in the system. This convergence clearly reveals that the dynamics of the CM of genome expression describe an attractor of the dynamics of stochastic expression (see also Section 2.6). Therefore, as shown in Figure 3, the CP, the CM according to the degree of nrmsf, acts as the genome attractor. Therefore, a change in the CP provokes a genome-wide avalanche over the entire genome expression. Whole-genome expression follows the change in the CP, and this is the origin of coherent gene expression behaviors [26].
Here, it is important to note that CSB does not stem from the law of large numbers under an equilibrium condition, but rather is deduced from the general specific property of an emergent property of an open (non-equilibrium) system. We consider the cell-fate change equivalent to a genome-state change (refer to biological discussions in [7]). The genome-state change occurs in such a way that the initial-state SOC control of whole-genome expression, which is in charge of maintaining a largely invariant global gene-expression profile, is destroyed through the erasure of initial-state sandpile criticality [7,9]. Given that the CP acts as a genome attractor, it is essential to understand how changes in the CP lead to the super-critical genome. To elucidate the timing of the change in the CP, we focus on CM correlation dynamics (see Section 4.2.2), i.e., the CM correlation between the initial and other time points:ĉ k (t 0 ).ĉ k t j (k = 1,2, . . . , K) over the experimental point, t j , whereĉ k t j is the unit vector (unit length) of the k th group vector, c k t j . CM grouping is used to examine singular behavior at the CP.
In MCF-7 breast cancer cells, activation of the ErbB receptor by HRG and EGF elicits two different biological responses (Figure 4). HRG stimulation induces cell differentiation, whereas EGF stimulation provokes cell proliferation [19,20]. The temporal CM correlation in HRG-stimulated MCF-7 cells reveals a divergent behavior at t j = 15 min ( Figure 4A, left panel), whereas EGF-stimulated MCF-7 cells (right panel) do not show any divergent behavior. Both responses exhibit sandpile CPs ( Figure 4D). These results suggest that the CP possesses both activated and inactivated states, i.e., ON/OFF expression states for a set of genes (critical gene set) corresponding to the CP. In EGF-stimulated MCF-7 cells, the CP is in the inactivated (OFF mode) state, whereas in HRG-stimulated MCF-7 cells, the CP is ON at 10-15 min and thereafter turns OFF.  2): 1 −ĉ k t j .ĉ k (t 0 ) (y-axis) is evaluated. At 15 min, the correlation response for HRG stimulation diverges from other responses, indicating that the CP is activated between 10 and 15 min for HRG, whereas for EGF stimulation, such a distinctly diverging correlation response does not occur (i.e., the CP is inactivated). Activation or inactivation of the CP marks the occurrence of cell differentiation for HRG and only proliferation (no differentiation) for EGF. Panels Direct evidence of the ON/OFF state of the CP is as follows: (1) In Figure 4B, for HRG stimulation, the switching of singular behaviors at 15-20 min occurs at the CP. At the boundary of the CP (ln<nrmsf >~−2.5), the singular behavior exhibits bimodal behavior in the fold-change on the CM grouping. At 15 min, a dominant positive (i.e., fold-change greater than one) singular behavior (ln<nrmsf > > −2.5) suggests that the swollen coil state of DNA occurs at the CP, while a negative singular behavior (ln<nrmsf > < −2.5) suggests the compact coil state of DNA. At 20 min, this bimodal behavior switches to a dominant negative singular behavior (> −2.5) with a positive singular behavior (< −2.5). With regard to EGF stimulation, such a switching transition does not occur during the early time points. Note: A negative fold-change occurs given that the reference point is the CM of whole-genome expression. Subtraction of the CM(t j ) from each gene expression gives negative expression. (2) In Figure 4C, the probability density function (PDF) of whole-genome expression [27] shows that at 15 min, around the CP region, the maximum probability density occurs in a positive area, whereas it becomes negative at 15-20 min. This validates that the CP is in the ON state at 15 min and in the OFF state at 20 min. The PDF clearly shows that ON-OFF switching of the CP induces a global avalanche in whole-genome expression.
The fold-changes within expression groups <c k (t j+1 )>/<c k (t j )> exhibit a clear first-order phase transition involving genome-sized DNA molecules (see more in Section 3.2). Through our studies, it became evident that the transition occurs as coherent behavior (at a mega-bp scale on the chromosome) and emerges from stochastic expression (coherent-stochastic behavior [8]). In other words, the averaging behavior (mean-field) of group expression corresponds to a coherent transition. In contrast, the ensemble average of the fold-change in individual expressions between two temporal groups, <c k (t j+1 )/c k (t j )>, does not reveal such characteristics in the transition. Interestingly, the ensemble average of the time difference in the expression group <c k (t j+1 )c k (t j )> supports the coherent scenario. This indicates that fluctuation (noise) of coherent dynamics is eliminated (see attached Supplementary Figure S1).
Regarding the cell-fate change in HRG-stimulated MCF-7 cells, the erasure of initial sandpile criticality occurs after 2 h ( Figure 4D). While the cell fate-guiding critical transition occurs at an earlier time point (15-20 min), the cell-fate change happens after 2 h. This time lag is needed to develop a new cell-fate attractor through coordinated local chromatin interaction (see details in [10]). Regarding the cell-fate change, the timing of the state-change of the CP through the switching of singular behaviors at the CP coincides with the timing of the erasure of the initial-state sandpile: at 24-48 h for atRA stimulation ( Figure 5D) and at 12-18 h for DMSO stimulation ( Figure 6D), suggesting that the cell-fate change occurs at 24-48 h for atRA and at 12-18 h for DMSO. The divergent behavior of the temporal CM correlation points to the occurrence of the cell-fate change: the change after the divergent behavior for atRA ( Figure 5A) and the change before it for DMSO ( Figure 6A). As for the DMSO response, it is interesting to observe a multi-step process of erasure of initial-state sandpile criticality ( Figure 6D): erasure at 8-12 h; recovery of criticality at 12-18 h; and then erasure again at 18-24 h. Multiple erasures suggest that the cell-population response passes over two SOC landscapes [7] at 8-12 h and 18-24 h.
The results obtained in cancer cells suggest that activation-deactivation of the critical gene set (CP) as the genome attractor plays an important role in the cell-fate change. Embryo development starts from a single cell (zygote) and thus corresponds to a completely different situation than in the cases discussed above, where the data refer to averages over millions of cells in a plate (i.e., cell-population level). The CP acts as a genome attractor at the single-cell level and in cell populations ( Figure 3). In temporal CM correlation, the CP is a point with no differential (Figures 7A, 8A and 9A), while it appears as a divergent point in cell populations. This feature reveals distinct response domains (critical states) in single-cell genome expression.
In a reprogramming event, the temporal CM correlation for the CP traverses a value of zero (corresponding to random-like behavior) at the 8-cell state for a human embryo cell ( Figure 7A). This implies that, after the 8-cell state, the human embryo cell completely erases the initial zygote criticality. Furthermore, this coincides with the erasure of the sandpile-type zygote CP ( Figure 7D). In biological terms, this corresponds to erasure of the initial stage of embryogenesis (driven by maternal heredity; see "SOC Control in Human and Mouse Related to the Developmental Oocyte-to-Embryo Transition" in the Discussion [7]).
It is important to note that groups of low-nrmsf presenting flattened CM correlations in time (ln<nrmsf> < −8.0) do not point to a no-response situation; on the contrary, they behave in a highly coherent manner to generate the autonomous SOC mechanism (see, e.g., Figure 6 in [6]). higher-order structure at the CP (K = 35: n = 360 mRNAs) is suggested to occur at 24-48 h, when initialstate criticality is erased (D). The switching of singular behaviors at the CP occurs with bimodal behaviors of the fold-change around the CP. The dominant negative fold-change indicates that the OFF state (compact globule DNA) occurs at 18-24 h, while the dominant positive fold change (swollen coil DNA) occurs at 24-48 h. (C) The temporal change in the PDF shows that a global avalanche occurs from 18-24 h (PDF: contracted) to 24-48 h (swollen) to 48-74 h (contracted), which supports the switching behaviors of the CP (B). (D) The full erasure of initial-state criticality occurs at 24-48 h, which indicates a cell-fate change. Underlined numbers represent the timing of the cell-fate change.

C) Cell Fate-Guiding Genomic Avalanche
Global Avalanche

Human Embryo Development: Genome Avalanche Along Low-Expressed Genes
Embryo development starts from a single cell (zygote) and thus corresponds to a completely different situation than in the cases discussed above, where the data refer to averages over millions of cells in a plate (i.e., cell-population level). The CP acts as a genome attractor at the single-cell level and in cell populations ( Figure 3). In temporal CM correlation, the CP is a point with no differential (Figures 7a,8a,9a), while it appears as a divergent point in cell populations. This feature reveals distinct response domains (critical states) in single-cell genome expression.
In a reprogramming event, the temporal CM correlation for the CP traverses a value of zero (corresponding to random-like behavior) at the 8-cell state for a human embryo cell (Figure 7a). This of collective behaviors that emerge in low-expressed genes (refer to local sub-critical state as the generator of the genome engine in Section 3.1). It is intriguing to observe that scaling behaviors (in a log-log plot) emerge, and the most vivid scaling behaviors develop at the 8-cell-morula state. Therefore, with the results regarding temporal CM correlation, we conclude that a major the cell-fate change (embryonic reprogramming) occurs after the 8-cell state. This coincides with the timing of inverse coherent perturbation of the genome engine (cyclic flux flow), where the genome system passes over the non-equilibrium fixed point after the 8-cell state (the onset of symmetrybreaking occurs; see Section 3.1).  . The bimodal foldingunfolding feature at the CP suggests the occurrence of intra-chain segregation in the transition for genome-sized DNA molecules (see more in Section 3.2). (C) A genome avalanche occurs as a traveling density wave along low-expressed genes (around non-fold change: y = 0; the yellow arrow indicates the next direction of travel). Furthermore, three distinct scaling behaviors (linear behaviors in the density profiles in the log-log plot) emerge, which reveals that coordinated chromatin dynamics emerge to guide the reprogramming process. (D) The erasure of zygote criticality after the 8-cell state coincides with the timing of (A) the zero temporal CM correlation through (B) a switching singular transition at the CP. This timing further matches the timing of coherent perturbation of the genome engine passing through the non-equilibrium fixed point (see Section 2.6). Note: A linear behavior emerges in the grouping of randomly shuffled whole-genome expression, as expected; therefore, higher nrmsf is associated with higher expression (see more in Figure 3D in [7]).

Mouse Embryo Development: Switching Scaling Behaviors in the Cell-Fate Change
In mouse embryo development, the scenario of the cell-fate change described above is confirmed: 1) Complete erasure of the memory of the zygote CP occurs after the late 2-cell state (Figure 8a  . The bimodal folding-unfolding feature at the CP suggests the occurrence of intra-chain segregation in the transition for genome-sized DNA molecules (see more in Section 3.2). (C) A genome avalanche occurs as a traveling density wave along low-expressed genes (around non-fold change: y = 0; the yellow arrow indicates the next direction of travel). Furthermore, three distinct scaling behaviors (linear behaviors in the density profiles in the log-log plot) emerge, which reveals that coordinated chromatin dynamics emerge to guide the reprogramming process. (D) The erasure of zygote criticality after the 8-cell state coincides with the timing of (A) the zero temporal CM correlation through (B) a switching singular transition at the CP. This timing further matches the timing of coherent perturbation of the genome engine passing through the non-equilibrium fixed point (see Section 2.6). Note: A linear behavior emerges in the grouping of randomly shuffled whole-genome expression, as expected; therefore, higher nrmsf is associated with higher expression (see more in Figure 3D in [7]).  which is confirmed by the temporal change in the PDF of whole-genome expression. However, there is no distinct scaling behavior as in embryo development. The timing of the cell-fate change from Th0 to Th17 cells, as in other biological regulations, is determined by the erasure of initial-state sandpile criticality (Figure 9d). This appears at 6 h, where inversion of the singular behaviors of the CP takes place, which is further confirmed by the timing of the inversion of perturbation of the genome engine (see Section 3.1). The OFF-ON switching at the CP guides the cell-fate change.

Systematic Determination of Local Critical States in Genome Expression
We demonstrate that the CP is a fixed point relative to a given biological regulation. Next, based on this fact, we show that critical states in genome expression can be determined systematically for both single-cell and cell population genome expression: 1) Single-cell level: Temporal CM correlation (Figures 7a,8a,9a) manifests distinct response The switching transition of the CP ( Figure 7B) occurs at every cell state change from zygote to blastocyst (only results from the 4-cell state to blastocyst are shown). This suggests that in early embryogenesis, changes in the cell state involve a global change (genome avalanche) in whole-genome expression. Notably, the temporal change in the PDF profile ( Figure 7C) shows that a genome avalanche occurs along low-expressed genes (around zero on the y-axis), shown as a traveling density wave (higher to lower nrmsf ) from the 8-cell state to the morula state. The reverse traveling wave (lower to higher nrmsf ) occurs after the morula state. This traveling wave points to the important role of collective behaviors that emerge in low-expressed genes (refer to local sub-critical state as the generator of the genome engine in Section 3.1). It is intriguing to observe that scaling behaviors (in a log-log plot) emerge, and the most vivid scaling behaviors develop at the 8-cell-morula state.
Therefore, with the results regarding temporal CM correlation, we conclude that a major the cell-fate change (embryonic reprogramming) occurs after the 8-cell state. This coincides with the timing of inverse coherent perturbation of the genome engine (cyclic flux flow), where the genome system passes over the non-equilibrium fixed point after the 8-cell state (the onset of symmetry-breaking occurs; see Section 3.1).

Mouse Embryo Development: Switching Scaling Behaviors in the Cell-Fate Change
In mouse embryo development, the scenario of the cell-fate change described above is confirmed: (1) Complete erasure of the memory of the zygote CP occurs after the late 2-cell state ( Figure 8A There are also three distinct scaling behaviors, as in human embryo development. Notably, during the middle 2-cell to 4-cell state through the late 2-cell state, whole-genome expression shifts from upto downregulation, i.e., marking the occurrence of a global avalanche, and furthermore, there is a reflection of the scaling behavior in the local sub-critical state (ln(nrmsf ) < −11) along the zero-fold change axis (y = 0). This reflection illustrates that a coherent switching of folding and unfolding chromatin dynamics occurs, which supports the important role in the local sub-critical state in early embryo development [8]. (3) The timing of the cell-fate change is further confirmed by the erasure of zygote sandpile criticality after the late 2-cell state.

Differentiation of Th17 Immune Cells from Th0 Cells: Partial Avalanche Guides the Cell-Fate Change
The cell-fate change in Th17 immune single-cell development again follows the same general scenario. Compared with embryo development, there are several distinct differences worth noting: (1) The CP does not pass over zero temporal correlation with the initial state (t = 0 h; Figure 9A), although the correlation response increases in time, which indicates that cell differentiation induces a partial-scale (specific set) change in whole-genome expression (vs. whole genome-scale change in embryonic reprogramming). (2) The switching singular transition at the CP (genome attractor) occurs sequentially from 3 h to 12 h ( Figure 9B). This manifests as the OFF state of the CP at 6-9 h and as the ON state at 9-12 h. Figure 9C illustrates that the change in expression for switching is still on a large scale, which is confirmed by the temporal change in the PDF of whole-genome expression. However, there is no distinct scaling behavior as in embryo development.
The timing of the cell-fate change from Th0 to Th17 cells, as in other biological regulations, is determined by the erasure of initial-state sandpile criticality ( Figure 9D). This appears at 6 h, where inversion of the singular behaviors of the CP takes place, which is further confirmed by the timing of the inversion of perturbation of the genome engine (see Section 3.1). The OFF-ON switching at the CP guides the cell-fate change.

Systematic Determination of Local Critical States in Genome Expression
We demonstrate that the CP is a fixed point relative to a given biological regulation. Next, based on this fact, we show that critical states in genome expression can be determined systematically for both single-cell and cell population genome expression:  Table 2), where the CP exists at the boundary between near-and sub-critical states (vs. between super-and sub-critical states in a single cell). This occurs for both MCF-7 and HL-60 cancer cell populations. In our previous studies on microarray data (cell-population level), critical states were determined by a change in expression profile by means of Sarle's bimodality coefficient ( Figure 10A; see Figure 9 for MCF-7 and HL-60 cells in [7]). This was accomplished by adding evidence of distinct response domains (super-, near-and sub-critical domains according to nrmsf ).   The existence of distinct critical states at both the single-cell and cell-population levels (Figures 7-10) is further demonstrated through the different convergence behaviors of stochastic expression, where the center of mass (CM) of a randomly selected expression of a critical state converges to the CM of the critical state ( Figure 11). This further explains the co-existence of distinct local attractors within the genome attractor (Figure 3), i.e., the formation of a heteroclinic critical-state attractor system. The co-existence of distinct local critical states points to the existence of an autonomous critical-control genomic system, i.e., the genome engine. Distinct coherent dynamics in a local critical state emerge from stochastic expression (refer to Figure 10 in [7]) due to the critical-state attractor self-organizing the dynamics of stochastic expression within the local critical state. distance (second column) between two ensemble averages, < ( = 0) and , (k = 1, 2,..,K = 25) from higher nrmsf, and confirmed three distinct behaviors with a boundary indicated by dashed vertical lines. Similarly, for HL-60 human leukemia cells, the corresponding results are shown for (D) atRA stimulation and (E) DMSO stimulation. In the third column (B-E), the corresponding region of critical states is shown in a histogram of gene expression according to ln(nrmsf). Note: The CP lies at the boundary between near-and sub-critical states, different from the case with a single cell (Figures  7-9).
Thus, dynamic expression flux analysis (Section 4.2.4) for the CM of critical states can be applied to reveal the genome engine mechanism for describing how autonomous SOC control of genome expression arises. Figure 12 shows that the sub-critical state acts as an internal source of expression flux and the super-critical state acts as a sink. The cyclic flux forms a dominant flux flow that generates a strong coupling between the super-and sub-critical states accompanied by their antiphase expression dynamics ( Figure 13). This results in a change in oscillatory feedback to sustain autonomous SOC control of overall gene expression. The average between-state flux ( Figure 12) represents a stable manifold for the thermodynamically open system. Formation of the dominant cyclic flux provides a genome engine metaphor for SOC control mechanisms pointing to a universal mechanism in the gene-expression regulation of mammalian cells for both the population and singlecell levels. As demonstrated in Section 3.1, global perturbation stems from the change in CP, which either enhances or suppresses the genome engine, and thus transforms the cell-fate change.   Table 1 and Table 2). The x-axis and y-axis correspond to experimental time points (Section 4.1) and the natural logarithm of the average of group expression, respectively. These confirm the existence of local critical states, and furthermore, their CM corresponds to local attractors (also confirmed in other biological regulations). Note: A mixture of critical states does not converge to the CM (see Figure 3 in [8]), showing that the genome attractor develops to form a heteroclinic critical-state attractor system. Thus, dynamic expression flux analysis (Section 4.2.4) for the CM of critical states can be applied to reveal the genome engine mechanism for describing how autonomous SOC control of genome expression arises. Figure 12 shows that the sub-critical state acts as an internal source of expression flux and the super-critical state acts as a sink. The cyclic flux forms a dominant flux flow that generates a strong coupling between the super-and sub-critical states accompanied by their anti-phase expression dynamics ( Figure 13). This results in a change in oscillatory feedback to sustain autonomous SOC control of overall gene expression. The average between-state flux ( Figure 12) represents a stable manifold for the thermodynamically open system. Formation of the dominant cyclic flux provides a genome engine metaphor for SOC control mechanisms pointing to a universal mechanism in the gene-expression regulation of mammalian cells for both the population and single-cell levels. As demonstrated in Section 3.1, global perturbation stems from the change in CP, which either enhances or suppresses the genome engine, and thus transforms the cell-fate change. The sub-critical state acting as a "large piston" for short moves (low-variance expression) and the super-critical state acting as a "small piston" for large moves (high-variance expression) with an "ignition switch" (a critical point: the genome attractor) are connected through a dominant cyclic state flux as a "camshaft", resulting in the anti-phase dynamics of two pistons (see the next Figure). Numerical values represent average between-state expression flux, whereas in the cell population, the values are based on a 10 −1 scale.   3), where whole-genome expression is ordered by the degree of nrmsf (see Figure 1).

Genomic Dynamics Determining the Cell-Fate Change from Embryo to Cancer: Cell-Fate Change Passing through a Non-Equilibrium Fixed Point
The different examples of cell-fate determination (cell differentiation and reprogramming in embryo development) have highlighted two crucial common elements: (1) the critical point (CP) acts as the organizing center of the cell fate, and (2) a sort of genome engine (Figure 12) fuels the transition

Sub-Critical State
Antiphase

Genomic Dynamics Determining the Cell-Fate Change from Embryo to Cancer: Cell-Fate Change Passing through a Non-Equilibrium Fixed Point
The different examples of cell-fate determination (cell differentiation and reprogramming in embryo development) have highlighted two crucial common elements: (1) the critical point (CP) acts as the organizing center of the cell fate, and (2) a sort of genome engine (Figure 12) fuels the transition and constitutes the basic mechanism of the cell-fate change. This genome engine works by canalizing the stochastic expression of local critical states into coherent behavior, provoking the motion of the center of mass (CM) of genome expression (i.e., change in the genome attractor). This is exactly what happens when an artificial engine transforms combustion (or any other energy source) into motion.
The CM acts as an attractor at both the whole-genome ( Figure 2) and critical-state levels ( Figure 11). Based on this fact, the expression flux (effective force) approach (Section 4.2.4) was developed to reveal dynamic interaction flux between critical-state attractors. Figure 13 demonstrates that the expression flux approach is instrumental for revealing a heteroclinic critical transition at the occurrence of a global avalanche that guides the coherent behaviors that emerge in the critical states. Interaction flux of between-state attractors is the underlying basic mechanism of epigenetic self-regulation, incorporating a rich variety of transcriptional factors and non-coding RNA regulation to determine the coherent oscillatory behaviors of critical states.
The flux dynamics approach has been further developed to analyze the quantitative evaluation of the degree of non-harmonicity and time-reversal symmetry-breaking in nonlinear coupled oscillator systems [28].
Thus, we summarize how and when the cell-fate change occurs as follows: (1) How the cell-fate change occurs: The intersection of interaction fluxes ( Figure 14 for cell population and Figure 15 for single cell; see the dashed ovals in interaction expression flux) occurs just before a cell-fate change. This intersection means that, thermodynamically, the genomic system lies at a non-equilibrium fixed point (e.g., see Figure 14C: atRA stimulation at 24 h with small interaction fluxes; see also Section 4.2.4), which suggests that before the cell-fate change, time-reversal symmetry is broken through passing over a non-equilibrium fixed condition. In terms of the genome engine, around the cell-fate change, a global perturbation induces either enhancement or suppression of the genome engine, where there is a dominant cyclic flux flow between the super-and sub-critical states (Table 3: single  cell; Table 4: cell population). In HL-60 cells (cell population), the genome engine is enhanced before the cell-fate change and suppressed (enhancement-suppression) thereafter. On the contrary, a reverse process of suppression-enhancement takes place in MCF-7 cancer cells ( Figure 14A; see more in [7]). In single-cell cases (embryo development and Th17 immune cells), a suppression-enhancement of the genome engine occurs (Figure 15). The varying sequences of perturbation of the genome engine may stem from different stages of the suppressive pressure on cell differentiation against cell proliferation.      (2) When the cell-fate change occurs: Cancer cell differentiation is accompanied by a change from a compact globule (OFF) to a swollen coil (ON) or vice versa at the CP (see . As shown in Figure 5C and Figure 6C (HL-60 cells), a global genome avalanche occurs at the timing of the cell-fate change. For MCF-7 cells (HRG stimulation), the cell-fate change (after 2 h: Figure 4D) occurs after the genome avalanche (see details of the attractor mechanism in [10]). The global impact provoked by the CP is further supported by the fact that (1) EGF-stimulated MCF-7 cells, where the CP is OFF with no cell differentiation, induce only local perturbation ( Figure 16C) and (2) HL-60 cells coincide with the timing of global perturbation ( Figure 16D).
Regarding the cell-fate change in single-cell dynamics, a global perturbation ( Figure 16A,B) occurs for the genome engine (after the late 2-cell and 8-cell stages and after 6 h for mouse, human embryo, and Th17 cells, respectively), where the timing coincides with that of the cell-fate change (i.e., erasure of the initial-state CP memory). This suggests that the cell fate-guiding change in the CP corresponds to the erasure of initial-state sandpile criticality. For embryo reprogramming, this picture is supported by the fact that temporal CM correlation of the CP from the initial state (zygote) passes zero-correlation (i.e., erasure of the initial-state CP memory: Figures 7A and 8A). Due to the genome attractor, the change in the CP (the edge of criticality) has a global impact on the entire genome expression through perturbation of the genome engine: global perturbation leads to a genome avalanche at the timing of the cell-fate change (Figures 7C-9C; also refer to Figure 13 and Biological Discussion in [7]). Therefore, elucidation of the activation-deactivation mechanism for the CP (involving more than hundreds of genes) through a coordinating change in chromatin structure [10] uncovers how the genome engine is either enhanced or suppressed around the timing of the cell-fate change. Furthermore, our SOC hypothesis is expected to predict when and how the cell-fate change occurs in different situations (e.g., cancer, induced pluripotent stem cells (iPS cells), etc.). The universality of the proposed model does not stem from our experimental results, but rather from the existence of very basic physical constraints (e.g., chromatin dimension and general organization, and phase transition phenomenology) independent of biological specificities.

Synergetic Behaviors of Mega Base Pair-Size DNAs through a Phase Transition
Up to this point, we have focused on the genomic dynamics that convert the signatures of the cell-fate change into self-organized complexity. Here, we should recall that our first objective was to recognize the unfeasibility of precise gene-by-gene regulation in the presence of the extremely compact folding of a 2-meter-long molecule that is compressed within a few microns of space. This calls for a search for a structural counterpart of the above-discussed model. In this respect, it is worth noting that the CP shows a clear bimodal expression distribution reminiscent of the swollen coil (ON) and compact globule (OFF) DNA transition, corresponding to inversion of the intra-chain phase segregation of the coil and globule states. This behavior closely resembles the intrinsic characteristics of the first-order phase transition inherent in genome-sized DNA molecules. The activation-deactivation of the CP is essential for the occurrence of the cell-fate change.
As for the phase transition of mega base pair-size DNA, until the late 20th century, it was believed that single polymer chains, including DNA molecules, always exhibited a cooperative but mild transition between the elongated coil and compact globule states, which was considered neither a first-order nor a second-order phase transition [30,31]. More recently, it has become evident that long DNA molecules above the size of several tens of kilo base pairs exhibit characteristics that help them undergo a large discrete transition, i.e., first-order properties related to the coil-globule transition [32][33][34]. Such first-order characteristics are rather general for semi-flexible polymers, especially polyelectrolyte chains such as giant DNA molecules. Additionally, it is important to note that insufficient charge neutralization in the globule state causes instability and leads to the generation of intra-chain segregation of individual single genomic DNA [35,36] (see, e.g., Figure 8B). When such instability is enhanced with long-range Coulombic interaction, the characteristic correlation length tends to become shorter, corresponding to generation of the critical state in the transition. Such transitional behavior that accompanies the folding transition of DNA is also observed for reconstructed chromatin [37][38][39].
Thus, the state change in the CP is suggested to guide the synergistic effect on mega base pair-size DNAs through the phase transition in the cell-fate change. Our recent study based on non-linear signal analysis techniques (recurrence quantification analysis and principal component analysis) revealed that chromatin remodeling played a role as the material basis of SOC-controlled genome expression [10]. A genome-wide expression avalanche stems from the coordinated activity of positional-based (not biological-functional) local chromatin interaction. The cell fate-guiding critical transition occurs based on the non-linear resonance of two major expression fluctuation modes, which establish an autonomous genome engine through activation of the CP. Elucidation of the link between the spatial position on the chromosome and co-regulation together with the identification of specific locations on the genome devoted to the generalization of perturbation stimuli provides a molecular basis for the self-organization dynamics of genome expression and the cell-fate decision.

CP Potentially Acting as the Center of Genome Computing
The singular behavior at the CP exhibits bimodal folding and unfolding chromatin dynamics through synergistic or cooperative behavior in higher-order structural transition of genomic DNA, which suggests the existence of a complex network of mega-sized ON/OFF DNA phase transitions in the cell-fate change. This interpretation was recently confirmed by the evidence of a massive re-arrangement of peri-centromeric nuclear regions correspondent to gene expression transitions [40]. The necessary link between chromatin remodeling and gene expression regulation was experimentally assessed [41]. These further imply that the CP may act as the center of genome computing, where chromatin remodeling is the material basis of such computation.
The time-dependent behavior of the genomic DNA transition follows a kinetic equation to exhibit cubic nonlinearity. A simultaneous change in the translational and conformational entropy of giant DNA together with surrounding counter ions causes bimodality in the free energy, which in turn speeds the transition under cubic nonlinearity in the kinetic equation [42,43]. The incorporation of negative feedback in cubic nonlinearity leads to fundamental characteristic nerve firings. Therefore, the occurrence of SOC in the brain's complex neural networks suggests that the state change in the CP is a computing process through coordinated DNA structural transitions to guide the cell-fate change. Elucidation of the computing process at the CP could reveal the control mechanism for the cell-fate change in a desirable manner, i.e., the existence of "genome intelligence". Here, note that cubic nonlinearity in the time-differentiation equation corresponds to a kinetic representation of the bimodal free energy, and criticality is generally interpreted based on the symmetry argument under this type of energy landscape [44].

Biological Meaning of the CP
As noted above, we investigated a new interpretation of the specificity of biological regulation without relying on an improbable "Maxwell's demon" ability to select the "right genes" interspersed in a two-meter long molecule constrained within the few cubic microns of the cell nucleus. The genes exert specific biological functions and an analysis of such functions could give us further hints regarding the mechanism of the cell-fate change. The strict link between chromatin remodeling and the phase transition is in fact mirrored by the gene composition of the super-critical (genome) state. Table S1 based on principal component analysis shows the 200 genes that are most strongly affected by the second principal component (PC2) of gene expression variance in HRG-treated MCF7 cells, as explained in [10]. PC2 corresponds to the super-critical set referring to the component with genes that show higher temporal variation. The continuous fluctuation in time of this component (orthogonal by construction to the first principal component (PC1), which corresponds to the equilibrium position, i.e., the attractor) allows for both continuous adaptation to the changing micro-environment (small avalanches) and fueling of the state change (phase transition). The biological role of these genes (endowed with statistically significant z-scores on PC2 and thus being the "top players" in the supercritical state, see Table S1) confirms that chromatin remodeling is the basic structural driver of the genome engine. The bolded genes in Table S1 are all consistent with structural chromatin remodeling obtained by both histone modifications and DNA breaks that allow for global nuclear re-organization. This nuclear re-organization requires a cross-talk with the perinuclear environment (note the prevalence of actin, myosin, cadherin and integrin genes) and activates known signaling hubs such as proto-oncogenes c-FOS and c-Myc, Mitogen-activated protein (MAP) kinase, and transformation-related protein 53 (TP53) ( Table S1). This is a stringent, albeit indirect, proof of the proposed model in terms of mainstream biology.

Final Remarks
We highlighted the basic invariance of the model in systems from embryo development to cancer. The different model cases are considered as examples of a general mechanism, where the need for a relevant time course to be analyzed is the only inclusion criterion, together with the presence of a relevant perturbation of the studied system. We can anticipate that other kinds of cells will behave according to the SOC model with differences linked to the presence (or absence) of a stable phenotype change during the time course considered. As aptly stated in [45], criticality is a general feature of biological regulation.
Finally, we summarize our findings in the following points: (1) The critical point (CP) acts as the center of the cell fate.
(2) How the cell-fate change occurs: Before the cell-fate change, the genomic system passes over the stable point (non-equilibrium fixed point) of the thermodynamically open system. The genome engine, an emergent dynamic property of an autonomous interaction flux between local critical states (distinct expression domains in whole-genome expression), is enhanced or suppressed to induce coherent perturbation of the dominant cyclic flux between local super-and sub-critical states.
(3) When the cell-fate change occurs: The change occurs at the timing of the erasure of the initial-state sandpile CP. At this time, a global genome avalanche occurs, except in HRG-stimulated MCF-7 cell differentiation, with a time lag between the genome avalanche and the cell-fate change (for further details of this mechanism, see [10]).
The genome engine suggests that the activation-deactivation mechanism of the CP elucidates how global perturbation occurs on self-organization through a change in signaling due to external or internal stimuli. Coordinated chromatin dynamics emerge to guide reprogramming through SOC control. A recent study found that the dynamics of the high-order structure of chromatin exhibit liquid-like behavior [46], which could be a crucial characteristic that enables the genome to exhibit SOC gene expression control for determination of the cell fate.
Further studies on these matters are needed to clarify the underlying fundamental molecular mechanism(s). The development of a theoretical foundation for the autonomous critical control mechanism in genome expression as revealed in our findings is expected to lead to new insights regarding for a general control mechanism for determination of the cell fate and genome computing.
For now, we can safely affirm that the strong interactions among genes with very different expression variance and physiological roles push for a complete reshaping of the current molecular-reductionist view of biological regulation focused on a single "significantly affected" gene to explain these regulation processes. The view that the genome acts as an integrated dynamic system is here to stay.

Biological Data Sets
We analyzed mammalian transcriptome experimental data for seven distinct the cell-fates in different tissues.
The colors used in the various plots throughout this report are based on the experimental events and were assigned as follows: black as the initial event, purple as the second event, and subsequent events as blue, dark cyan, dark green, dark yellow, brown, orange, red, dark pink, and pink.
For microarray data, the robust multichip average (RMA) was used to normalize expression data for further background adjustment and to reduce false positives [47][48][49].
For RNA-Seq data, RNAs with RPKM and FPKM values of zero over all of the cell states and experimental time points were excluded. In the analysis of sandpile criticality, random real numbers in the interval [0, a] generated from a uniform distribution were added to all expression values (only in Figures 7D-9D). This procedure avoids the divergence of zero values in the logarithm. The robust sandpile-type criticality through the grouping of expression was checked by changing a positive constant: a (0 < a < 10); we set a = 0.001. Note: The addition of large random noise (a >> 10) destroys the sandpile CP.

Normalized Root Mean Square Fluctuation (nrmsf )
Nrmsf (see also Methods in [6]) is defined by dividing rmsf (root mean square fluctuation) by the overall maximum {rmsf i } (Equation (1)): where rmsf i is the rmsf value of the i th RNA expression, which is expressed as ε i (s j ) at a specific cell state s j or experimental time (e.g., in mouse embryo development, S = 10: s 1 = zygote, early 2-cell, middle 2-cell, late 2-cell, 4-cell, 8-cell, morula, early blastocyst, middle blastocyst, and s 10 = late blastocyst), and ε i is its expression average over the number of cell states. Note: Nrmsf is a time-independent variable, where nrmsf acts as an order parameter of self-organization [6,7].

CM Correlation Analysis
To investigate transition dynamics, the correlation metrics based on the center of mass (CM) grouping, the CM correlation, was built upon the following basic statistical formalization: (1) CM grouping: genome expression is considered as an N-dimensional vector, where the average value (CM) of whole-genome expression at t = t j is subtracted from each expression (refer to Figure 1). Next, the whole-genome expression is sorted and grouped according to the degree of nrmsf, where CM grouping has K groups and within each group there are n number of expressions (N = K . n): N-dimensional CM grouping vector,C t j = c 1 t j , c 2 t j , . . . , c k t j , . . . , c K t j ; c 1 t j and c K t j are the highest and lowest group vectors of nrmsf, respectively. Here, the unit vector of the k th vector c k t j is . Note that the elements less than n in the last group (the lowest nrmsf ) have been removed from the analysis.
(2) If we keep in mind that correlation corresponds to the cosine of the angle between unit vectors, i.e., the inner product of unit vectors (cos θ =â.b, θ : angle;â.b: dot product (scalar) of unit vectors:â,b), two different CM correlations can be considered: (i) Spatial CM correlation: for a given time point (t = t j ), the development of CM correlation between the first group (highest nrmsf group) and other vectors:ĉ 1 t j .ĉ k t j ; (k = 2,3, . . . K).
(ii) Temporal CM correlation: for a given group (k), the development of CM correlation between the initial and other experimental points:ĉ k (t 1 ).ĉ k t j (k = 1,2,3, . . . K) over experimental time points, t j (see Biological Data Sets).
Next, we consider these vectors with regard to their CM (called the CM expression vector), the CM natural log of the whole-genome expression vector and its fold-change vector. The CM natural log of the whole-genome expression vector is defined as ln(c(t j )) = ln(ε(t j )) − CM(ln(ε(t j ))) and CM(ln(ε(t j ))) = 1 N N i=1 ln i t j and its fold-change vector is ln(c(t j )/c(t k )) = ln(ε(t j )/ε(t k )) − CM(ln(ε(t j )/ε(t k ))). This vector representation again satisfies the logarithm of a quotient: ln(c(t j )/c(t k )) = ln(c(t j )) − ln(c(t k )).
From embryo to cancer development, global avalanches in changes in the CP are captured, and furthermore, this presents scaling behaviors in addition to the traveling expression wave in human and mouse embryo development (RNA-Seq data) (see Figures 7C and 8C), which suggests coordinated DNA folding-unfolding dynamics.

Expression Flux Analysis: Expression Flux between Critical-State Attractors
We developed an expression flux analysis to describe the genome engine mechanism for both single-cell and population genome expression (see Figure 12). The key fact is that the dynamics of coherent behavior emerging from stochastic expression in distinct critical states (coherent-stochastic behavior: CSB) follow the dynamics of the CM of a critical state. This convergence shows that the CM of a critical state acts as a critical-state attractor for stochastic expression within the critical state ( Figure 11). We developed the expression flux approach to reveal the dynamic interaction flux between critical-state attractors [6][7][8].
The CSB in a critical state corresponds to the scalar dynamics of its CM. The numerical value of a specific critical state (i.e., super-, near-or sub-critical state) is represented by X(s j ) at a specific experimental event (s j ), where an experimental event (s j ) corresponds to a cell state or an experimental time point. The expression flux between critical states is interpreted as a non-equilibrium system and evaluated in terms of a dynamic network of effective forces, where interaction flux is driven by effective forces between different critical states and can be described by a second-order time difference. It is important to note that the oscillatory phenomenon is interpreted using a second-order difference equation with a single variable. This is equivalent to the inhibitor-activator dynamics given by a couple of first-order difference equations with two variables.
The flux dynamics approach has been further developed to analyze the quantitative evaluation of the degree of non-harmonicity and time-reversal symmetry-breaking in nonlinear coupled oscillator systems [28].
The basic formulas of expression flux dynamics are given as follows: Net self-flux of a critical state: The net self-flux, the difference between IN flux and OUT flux, describes the effective force on a critical state. This net self-flux represents the difference between the positive sign for incoming force (net IN self-flux) and the negative sign for outgoing force (net OUT self-flux); the CM from its average over all cell states represents up-(or down-) regulated expression for the corresponding net IN (or OUT) flux.
The effective force is a combination of incoming flux from the past to the present and outgoing flux from the present to the future cell states (Equation (2) where ∆P is the change in momentum with a unit mass (i.e., the impulse: F∆s = ∆P) and the natural log of the average (< . . . >) of a critical state, X s j = ln 1 N C N C i=1 ε i s j with the i th expression ε i s j at the j th experimental event, s = s j (N C = the number of RNAs in a critical state; refer to Tables 1 and 2); the average of net self-flux over the number of critical states, <f (X)> = <INflux> − <OUTflux>.
Here, scaling and critical behaviors occur in log-log plots of group expression, where the natural log of an average value associated with group expression such as ln<nrmsf > or ln<expression> is taken. Thus, in defining expression flux, the natural log of the average expression (CM) of a critical state is considered.
It is important to note that each embryo cell state is considered as a statistical event (a statistical event does not necessarily coincide with a biological event) and its development is considered as a time arrow (time-development) when evaluating the average of group expression: the fold-change in expression and temporal expression variance (nrmsf ). This implies that an interval in the dynamic system (Equation (2)) is evaluated as the difference in events, i.e., ∆s j = s j − s j−1 = 1 and ∆s = s j+1 − s j−1 = 2 in embryo development, as well as the difference in experimental times such as in cell differentiation (the actual time difference can be considered as scaling in time). We then evaluate a force-like action in expression flux.
The interaction flux of a critical state: The interaction flux, representing the flux of a critical state X(s j ) with respect to another critical state (super-, near-, sub-) or the environment (E: milieu), Y j , can be defined as (Equation (3)): where, again, the first and second terms represent IN flux and OUT flux, respectively, and the net value (i.e., IN flux − OUT flux) represents incoming (IN) interaction flux from Y for a positive sign and outgoing (OUT) interaction flux to Y for a negative sign. Y represents the numerical value of either a specific critical state or the environment, where a state represented by Y is different from one represented by X.
With regard to the global perturbation event, the net kinetic energy flux [7] clearly reveals the timing of the global perturbation event ( Figure 16) (Equation (4)): where the kinetic energy of the CM for the critical state with unit mass at s = s j is defined as 1/2 . v s j 2 with average velocity: v s j ≡ X(s j )−X(s j−1 ) ∆s j . Net self-flux as summation of interaction fluxes: Due to the law of force, the net self-flux of a critical state is the sum of the interaction fluxes with other critical states and the environment (Equation (5)): where state A i ∈ Super, Near, Sub with A i X, and M is the number of internal interactions (M = 2), i.e., for a given critical state, there are two internal interactions with other critical states. Equation (5) tells us that the sign of the difference between the net self-flux and the overall contribution from internal critical states, f X s j − M=2 i=1 f (X(s); A i ), reveals incoming flux (positive) from the environment to a critical state or outgoing flux (negative) from a critical state to the environment. When the difference in all critical states is zero, the genome system itself lies in a stable point (non-equilibrium fixed point) of a thermodynamically open system (no average flux flow from the environment). Average between-state flux ( Figure 12) represents a stable manifold of the thermodynamically open system. Order of Averaging: ln 1 n n i=1 ε i s j vs. 1 n n i=1 ln ε i s j Here, we need to address the previous result of expression flux dynamics in mouse single-cell genome expression [8], where the expression of a critical state was taken as X s j = 1 N C N C i=1 ln ε i s j , which has a different order of operations: first take the natural log of the expression and then calculate the average. Hence, in flux dynamics, we examine whether or not the mathematical operations of averaging and the natural log, i.e., between ln 1 n n i=1 ε i s j and 1 n n i=1 ln ε i s j , can be exchanged (mathematically commuted). In microarray data, flux behaviors do not change much between different orders of actions, whereas in RNA-Seq data, they are not commuted because its data structure has many zero values. While the addition of small random noise to the log of expression, ln ε i s j (used for the previous result [8], a noise-sensitive case) has a good effect, the addition of such noise to ln 1 n n i=1 ε i s j ), used for noise-insensitive cases (such as in this report), does not have a good effect. Therefore, due to this sensitivity, micro-array data may be better for expression flux analysis (see Figure 13). Although the detailed dynamics of interaction flux changes by ordering actions differently for RNA-Seq data (e.g., Figure 6 in [8]), two important characteristics in the genome engine, the formation of dominant cyclic flux between super-and sub-critical states and the generator role of the sub-critical state, do not change (invariant features). Thus, we conclude that the concept of the genome engine is quite robust. Acknowledgments: M.T. sincerely thanks the following institution and individuals who helped complete this research project: the SEIKO Life Science Laboratory, Osaka, Japan, his family (particularly, his daughters, Kimiko and Kazumi Tsuchiya, and Harry Taylor with any editing), Andrzej Kasperski for fruitful discussions through the proof-reading and Jekaterina Erenpreisa for biological discussions.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Normalized root mean square fluctuation PC1, PC2

Abbreviations
The first principal component, the second principal component PDF Probability density function SOC Self-organized criticality Th0 Naïve T helper Th17 T helper 17