Self-Organizing Global Gene Expression Regulated through Criticality: Mechanism of the Cell-Fate Change

Abstract

Background: A fundamental issue in bioscience is to understand the mechanism that underlies the dynamic control of genome-wide expression through the complex temporal-spatial self-organization of the genome to regulate the change in cell fate. We address this issue by elucidating a physically motivated mechanism of self-organization.

Principal findings: Building upon transcriptome experimental data for seven distinct cell fates, including early embryonic development, we demonstrate that self-organized criticality (SOC) plays an essential role in the dynamic control of global gene expression regulation at both the population and single-cell levels. The novel findings are as follows: i) Mechanism of cell-fate changes: A sandpile-type critical transition self-organizes overall expression into a few transcription response domains (critical states). A cell-fate change occurs by means of a dissipative pulse-like global perturbation in self-organization through the erasure of initial-state critical behaviors (criticality). Most notably, the reprogramming of early embryo cells destroys the zygote SOC control to initiate self-organization in the new embryonal genome, which passes through a stochastic overall expression pattern. ii) Mechanism of perturbation of SOC controls: Global perturbations in self-organization involve the temporal regulation of critical states. Quantitative evaluation of this perturbation in terminal cell fates reveals that dynamic interactions between critical states determine the critical-state coherent regulation. The occurrence of a temporal change in criticality perturbs this between-states interaction, which directly affects the entire genomic system. Surprisingly, a sub-critical state, corresponding to an ensemble of genes that shows only marginal changes in expression and consequently are considered to be devoid of any interest, plays an essential role in generating a global perturbation in self-organization directed toward the cell-fate change.

Conclusion and significance: 'Whole-genome' regulation of gene expression through self-regulatory SOC control complements gene-by-gene fine tuning and represents a still largely unexplored non-equilibrium statistical mechanism that is responsible for the massive reprogramming of genome expression.

Fig 1. Correlation of gene expression profiles.

A) The same cell types: Left panels: a near-unity Pearson correlation, r, in whole gene expression (N: total number of mRNAs) within the same cell type is shown for different types of molecular stimulation (first row: HRG- vs. EGF-stimulated MCF-7 cells; second row: DMSO- vs. atRA-stimulated HL-60 cells). Right panels: 313 (n) gene expressions, which have a common probe ID among four transcriptome microarray expression data sets (see Methods) also show a near-unity Pearson correlation within the same cell type. B) Different cell types (n = 313 gene expressions): The near-unity Pearson correlation between independent samples of the same cell type breaks down when the gene expression profiles come from different (HL-60 and MCF-7) cell types. ε(t) represents the ensemble of expression at time point t (N: the whole set; n: an ensemble set) and ln(ε(t)) represents its natural logarithm, where the natural log of an individual expression value is taken. Plots show t = 90min for MCF-7-stimulated cells and t = 18h for HL-60-stimulated cells.

Fig 2. Transition of gene expression from a stochastic to a genome-wide attractor profile.

A) Plot shows the whole expression profiles at 10min (x-axis) and 15min (y-axis) for the HRG response in MCF-7 cells. A box is constructed from the center of mass (i.e., average of whole expression), (CM(10min), CM(15min)) (black dot); a box contains gene expression within the range from CM(t_j)—d to CM(t_j) + d with a variable box size, d. To highlight the scaling of the Pearson correlation with box size d, r_d for d = 0.05, 0.1 and 0.2 are reported. The plot in the upper left corner shows that, between gene expression profiles, the Pearson correlation r_d follows a tangent hyperbolic function: r_d = 0.97 ⋅ tanh(6.79 ⋅ d − 0.039) (p<10⁻⁴), which reveals a critical transition in the correlation development. B) Stream plots for the box sizes in Panel A. These plots are generated from the vector field values {Δx_i, Δy_i} at expression points {x_i(10min), x_i(15min)}, where Δx_i = x_i(15min)—x_i(10min), Δy_i = x_i(20min)—x_i(15min), and x_i(t_j) is the natural log of the i^th expression: x_i(t_j) = ln(ε_i(t_j)) at t = t_j (t_j = 10min or 15min; i = 1,2,.., N = 22,277). Blue lines represent streamlines and red arrows represent vectors at a specified expression point (plot every 2^nd, 6^th 10^th and 20^th point for d = 0.05, 0.1, 0.2, and the whole set, respectively). When we move from a small number of genes to the whole set, gene expression shifts from a stochastic to a genome-wide attractor profile.

Fig 3. Self-organized criticality (SOC) in the DMSO-stimulated HL-60 cell fate.

A) The grouping of whole expression at 18-24h generates 25 groups with an equal number of 505 elements (mRNAs) according to the fold change in expression. A plot of the average value for each group in log-log space (x: fold change at 18-24h vs. y: expression at 24h) reveals a sandpile-type critical behavior at the critical point (CP), where the CP in terms of the ensemble (group) average (< >) occurs at the near-zero-fold change (x = 0; null expression change; x-axis) with <nrmsf(CP)> = 0.0756. Orange dots represent single mRNA expression in the background. B) The probability density function for the ensemble of expression (coherent expression state: CES) is shown in the regulatory space (x: expression at 24h vs. y: fold change in expression at18h-24h). Plots show that around <nrmsf(CP)> = 0.0756, a CES (highest density: x< 2.0) is annihilated and a new CES (highest density: x> 2.0) is bifurcated. The left to right panels show the sequence of the bifurcation-annihilation event: before (I: 0.060 <nrmsf< 0.070: <nrmsf> = 0.0647), onset (II: 0.071 <nrmsf< 0.081: <nrmsf> = 0.0756)) and after the event (III: 0.078 <nrmsf< 0.088: <nrmsf> = 0.0828). Colored bars represent the probability density. C) The corresponding plot (B) reveals that a step function-like critical transition occurs at <nrmsf(CP)> = 0.0756 in the space (x: <nrmsf> of coherent state vs. y: its expression of the highest density at 18h). Blue and black arrows represent average values of the density trends before and after the transition, respectively. The plot shows, in the vicinity of the CP, the occurrence of a self-similar bifurcation (symmetry breaking) in the expression profile to that of the overall expression (see section III). D) Random mRNA DMSO expression matrix reveals that, in this case, a CP does not exist. This is confirmed by anomalous features of the corresponding SOC (Methods): non-scaling-divergent (3 different time points are shown by colors) and non-sandpile critical behavior of random expression between two different time points (orange dots: single random expression). The random matrix is made by randomly selecting each matrix component (i,j) from the original DMSO expression matrix (12625 expression (i) times 13 time points (j)). We observed similar linear correlative behaviors for other cells in both microarray and RNA-Seq data.

Fig 4. Time-development of the characteristic behaviors of SOC.

A) MCF-7 cells and B) HL-60 cells. At each experimental time point, t_j (A: 18 time points and B: 13 points; Methods), the (ensemble) average nrmsf value of the CP at t = t_j, <nrmsf(CP(t_j))> is evaluated at the sandpile-type critical point (top of the sandpile: right panels). In the top center panels, <nrmsf(CP(t_j))>, is plotted against the natural log of t_j (A: min and B: hr). Error bar represents the sensitivity of <nrmsf(CP(t_j))> around the CP(t_j), where the bar length corresponds to the change in nrmsf in the x-coordinate (fold-change in expression) from x(CP(t_j))—d to x(CP(t_j)) + d for a given d (A: d = 0.005; B: d = 0.01; due to much more mRNAs in MCF-7 cells). Temporal averages of <nrmsf(CP(t_j))> are A) $\overline{<nrmsf\left(\mathrm{C}\mathrm{P}\right)>}$_HRG = 0.094 and $\overline{<nrmsf\left(\mathrm{C}\mathrm{P}\right)>}$_EGF = 0.081, and B) $\overline{<nrmsf\left(\mathrm{C}\mathrm{P}\right)>}$_{DMSO, atRA} = 0.078, where an overbar represents temporal average. Note: An overbar for a temporal average is used when ensemble and temporal averages are needed to distinguish. A) MCF-7 cells: The temporal trends of <nrmsf(CP(t_j))> are different for HRG and EGF. The onset of scaling divergence (left panels: second and third rows) occurs at around $\overline{<nrmsf\left(\mathrm{C}\mathrm{P}\right)>}$ (black dashed line), and reflect the onset of a ‘genome avalanche’ (Methods). B) HL-60 cells: The trends of <nrmsf(CP(t_j))> for the responses to both DMSO (black line) and atRA (blue) seem to be similar after 18h (i.e., global perturbation; see section VI). The scaling-divergent behaviors for both DMSO and atRA reveal the collapse of autonomous bistable switch (ABS [10,11]) exhibited by the mass of groups in the scaling region (black solid cycles) for both DMSO and atRA. The onset of divergent behavior does not occur around the CP ($\overline{<nrmsf\left(\mathrm{C}\mathrm{P}\right)>}$ = 0.078), but rather is extended from the CP (see section III). The power law of scaling behavior in the form of 1- <nrmsf> = α<ε>^-β is: A) α = 1.29 & β = 0.172 (p< 10⁻¹⁰) for HRG, and α = 1.26 & β = 0.157 (p<10⁻⁹) for EGF; B) α = 1.60 & β = 0.301 (p< 10⁻⁶) for DMSO, and α = 1.42 & β = 0.232 (p< 10⁻⁶) for atRA. Each dot (different time points are shown by colors) represents an average value of A) n = 742 mRNAs for MCF-7 cells, and B) n = 505 mRNAs for HL-60 cells.

Fig 5. Genome-state change revealed through erasure of the initial-state criticality in overall expression (cell population level): The grouping of overall expression at t = t_j (j ≠ 0) according to the fold change in expression from the initial overall expression (t = 0) shows how the initial-state critical behavior is erased over time, i.e., a sandpile profile in overall expression is destroyed at t = t_j from t = 0.

This event points to the time when the genome-state change occurs. On the x-axis, ln(<ε(t)/ε(0h)>) (t: min or hr) represents the natural log of the ensemble average (< >) of the fold change in expression, ε(t)/ε(t = 0h), and on the y-axis, ln(<ε(t)>) represents the natural log of the ensemble average of expression, <ε(t)>. A) MCF-7 cells: In HRG-stimulated cells (black), the divergent behavior in up-regulation is no longer observed, and the same shape is apparent after 3h. This suggests that the genome-state change occurs at 3h (red) through erasure of the initial-state up-regulation process (partial erasure). In contrast, in EGF-stimulated cells (blue), almost the same sandpile profile remains for up to 36h, which suggests that no genome-state change occurs (see the local perturbation in Fig 14). B) HL-60 cells: A CP is erased at 24h (red) and 48h (red) in DMSO- (black) and atRA-stimulated (blue) cells through the disappearance of divergent behaviors in both up- and down-regulation of the initial state (full erasure), respectively. Furthermore, these plots suggest that the epigenomic states in DMSO- and atRA-stimulated cells become the same after 48h. Note: The Pearson correlations of overall expression between different time points are near-unity (Fig 1A). Each dot represents the average value of A) n = 742 mRNAs for MCF-7 cells, and B) n = 505 mRNAs for HL-60 cells.

Fig 6. The self-similar bifurcation or annihilation of a characteristic coherent expression state in the vicinity of a critical point in different cell types.

A coherent expression state (CES), which contains approximately 1000 expression points, is bifurcated or annihilated around the CP through nrmsf grouping. The average nrmsf of an expression group, <nrmsf>, is evaluated for a variable range: x + 0.01^.(m-1) <nrmsf< x + 0.01^.m (integer, m≤10 for each x = 0.7, 0.8 and 0.9). The probability density function (PDF) in the regulatory space (z-axis: probability density) around <nrmsf(CP)> shows. A) MCF-7 cells: In the HRG response at 15-20min, around <nrmsf(CP)> (<nrmsf(CP)>_HRG = 0.0948 and <nrmsf(CP)>_EGF = 0.0809; Fig 4A), the PDF exhibits a bimodal-flattened-bimodal transition, in which a bimodal profile points to the existence of two CESs: one represents a low-expression state (LES) and the other represents a high-expression state (HES); the valley defines the boundary between low and high expression [10]. In the EGF response (second row), above <nrmsf(CP)>, a low-expression state (LES) is annihilated and only a high-expression state (HES) exists, i.e., a unimodal-bimodal transition is present. In the HRG response, at a time period other than 15-20min, a unimodal-bimodal transition occurs [10]. This result shows the self-similar symmetry-breaking event to the overall expression for MCF-7 cells, even at a transient change in the HRG response: a bimodal-flattened-bimodal transition at 15–20 min. B) HL-60 cells: A unimodal-flattened-unimodal transition occurs at 18-24h (pseudo-3-dimensional PDF plots of Fig 3B) for DMSO and at 0-2h for atRA, which again reveals the self-similarity to the overall expression for HL-60 cells (section III).

Fig 7

Genome-state change revealed through erasure of the initial-state critical point in overall expression (single cell level): ‬‬‬‬‬‬‬Transcriptome RNA-Seq data (RPKM) analysis for A-C) human and mouse embryo development, and D) T helper 17 cell differentiation. A) In both human and mouse embryos (red: human; blue; mouse; refer to development stages in Methods), a critical transition is seen in the development of the overall expression correlation between the cell state and the zygote single-cell stage, with a change from perfect to low (stochastic) correlation: the Pearson correlation for the cell state with the zygote follows a tangent hyperbolic function: a − b ⋅ tanh(c ⋅ x − d), where x represents the cell state with a = 0.59, b = 0.44 c = 0.78 and d = 2.5 (p< 10⁻³) for human (red dashed line), and a = 0.66, b = 0.34 c = 0.90 and d = 3.1 (p< 10⁻²) for mouse (blue). The (negative) first derivative of the tangent hyperbolic function, -dr/dx, exhibits an inflection point (zero second derivative), indicating that there is a phase difference between the 4-cell and 8-cell states for human, and between the middle and late 2-cell states for mouse (inset); a phase transition occurs at the inflection point. Notably, the development of the SOC control in the development of the sandpile-type transitional behaviors from the zygote stage (30 groups; n: number of RNAs in a group: Methods) is consistent with this correlation transition: B) In human, a drastic change in criticality occurs after the 8-cell state; a sandpile-type CP (at the top of the sandpile) disappears, and thereafter there are no critical points. This shows that the zygote SOC control in overall expression (i.e., zygote self-organization through criticality) is destroyed after the 8-cell state, which indicates that the memory of the initial stage of embryogenesis is lost through a stochastic pattern as in the linear correlation trend (refer to the random expression matrix in S1 Fig or to Fig 3D). The results suggest that reprogramming (massive change in expression) of the genome occurs after the 8-cell state. C) In mouse, a sandpile-type CP disappears right after the middle stage of the 2-cell state and thereafter a stochastic linear pattern occurs, which suggests that reprogramming of the genome after the middle stage of the 2-cell state destroys the SOC zygote control. D) In Th17 cell differentiation, a sandpile-type CP disappears after 3h through a stochastic linear pattern. Therefore, the plot suggests that the genome-state change occurs at around 3-6h in a single Th17 cell. ‬‬‬‬‬‬‬‬‬‬‬‬‬‬‬‬‬‬‬

Fig 8. The SOC control landscape as revealed by a sandpile-to-sandpile transition in mouse embryo development.

The development of a sandpile transitional behavior between sequential cell states suggests the existence of an SOC control landscape (first row: schematic picture of a valley-ridge-valley transition; x-axis: cell state; y-axis: degree of SOC control). The second and third rows show that a sandpile-to-sandpile transition occurs in mouse embryo development from the zygote single-cell stage to the morula cell state: I: a sandpile (i.e., critical transition) develops from the zygote single-cell stage to the early 2-cell state = > II: a sandpile is destroyed from the middle to the late 2-cell state, which exhibits stochastic expression (i.e., no critical transition; refer to the random mouse expression matrix in S2 Fig) = > III: a sandpile again develops from the 8-cell state to the morula state. These results show that a significant perturbation (reprogramming) in self-organization occurs from the middle stage to the late 2-cell stage through a stochastic overall expression (refer to Fig 7). Note: Qualitatively, there is a high degree of SOC control for a well-developed shape of the sandpile-type transition (SOC control), an intermediate degree of SOC control for a weakened (broken) sandpile, and a low degree for non-SOC control. The latter is due to stochastic expression. The linear behavior (absence of a critical point) in mouse embryo development is also reflected in the low Pearson correlation (r ~0.21 after the 8-cell state) (Fig 7A and 7C).

Fig 9

Critical states revealed through distinct functional behaviors of the bimodality coefficient: A) MCF-7 cells, B) HL-60 cells and C) Random DMSO expression matrix (see Fig 3D). The left panels (A and B) show the frequency distribution of expression according to the degree of nrmsf. The center panels show the temporal average of Sarle's bimodality coefficient, <b_i>, of the i^th group over time; the value 5/9 represents the threshold between the unimodal (below 5/9) and bimodal or multimodal distributions (above 5/9). The grouping of expression is made at a specific sequence of discrete values of nrmsf (x_i: nrmsf_i = i/100; i: integers) with a fixed range: x_i—k^.d < x_i < x_i + k^.d. The values of k and d are set to be k = 150 and d = 0.0001 for MCF-7 and HL-60, and k = 100 and d = 0.001 for a random matrix based on the convergence of the bimodality coefficient. The convergence of the difference in bimodality coefficients at x_i with an increase in k (i.e., as the number of elements in a group increases) is shown between the next neighbors, <b_i(k;x_i)>—<b_i(k-1;x_i)> for HL- 60 cells in the right panels of B). The 6 colored dots represent the convergent behaviors of different nrmsf points. The behavior of the time average of the bimodality coefficient exhibits A) Tangent hyperbolic function, b_i = a − tanh(b + c⟨nrmsf⟩_i); a = 1.38 and 1.35; b = 0.123 and 0.301; c = 13.8 and 10.2 for HRG (p< 10⁻⁴) and EGF (p< 10⁻¹⁰), respectively, B) Heaviside step function-like transitions for HL-60 cell fates, and C) No transition for a random DMSO expression matrix, which importantly reveals that random noises through the formation of a Gaussian distribution destroy a sandpile critical behavior. Based on these distinct behaviors, we can determine the boundaries of averaged critical states (Table 1; see section III): A) Critical states are defined by two points: the average CP, $\overline{<nrmsf\left(\mathrm{C}\mathrm{P}\right)>}$, the onset of a genome avalanche (Fig 4A), for the upper boundary of the sub-critical state, and the point where the change in the bimodality coefficient, Δ<b_i>, reaches zero for the lower boundary of the super-critical state; the near-critical state is between them, and B) Step function-like transitions reveal the boundaries of averaged critical states, where the near-critical state corresponding to the transitional region separates the other states.

Fig 10

Coherent-stochastic behaviors in critical states: A) MCF-7 cells and B) HL-60 cells. 200 random-number ensemble sets are created, where each set has n (variable) sorted numbers, which are randomly selected from an integer series {1,2,..,N} (N: total number of mRNAs in a critical state: Table 2). These random sets are used to create random gene ensembles at t = t_j (j^th experimental time points; Methods) from each critical state. 1) Left panels (A and B): For each n at t = t_j, Pearson correlations are evaluated between different random gene ensembles in the critical states by averaging over the 200 ensemble sets, and then, their temporal average over experimental time points are evaluated. This gives a near-zero Pearson correlation, consistent with the global stochastic character of microscopic transcriptional expression regulation in critical states. 2) Right panels (A and B): For each random gene ensemble, the Euclidean distance of single time points from the center of mass CM(t_j) of the critical states is evaluated by averaging over the 200 ensemble sets. The sharp damping of variability confirms that the emergent coherent dynamics of the critical states correspond to the dynamics of the CM(t_j).

Fig 11. Dynamics of the center of mass (CM) of critical states compared with the net self-flux dynamics: Colored lines (red: super-critical; blue: near-critical; purple: sub-critical state) represent net self-fluxes of critical states from their temporal averages (effective force acting on the CM: Methods).

Black lines represent the dynamics of the CM of critical states from their temporal averages, which are increased three-fold for comparison to the corresponding net self-fluxes. The plots show that the net self-flux dynamics follow up- (down-) regulated CM dynamics, such that the sign of the net self-flux (i.e., IN and OUT) corresponds to activation (up-regulated flux) for positive responses and inactivation (down-regulated flux) for negative responses. The natural log of the experimental time points (MCF-7: minutes and HL-60: hours) is shown.

Fig 12. Genomic expression dynamics revealed through a flux analysis that includes crosstalk with the environment: A) Average values of interaction flux (colored arrows) for MCF-7 and HL-60 cells show that a sub-critical state acts as an internal ‘source’, where IN flux from the environment is distributed to other critical states.

In contrast, a super-critical state acts as an internal ‘sink’ that receives IN fluxes from other critical states, and the same amount of expression flux is sent to the environment, due to the average flux balance (Methods). Furthermore, the formation of a dominant cyclic state flux is revealed between super- and sub-critical states through the environment. The average interaction flux is represented as i-j: interaction flux of the i^th critical state with the j^th critical state (i, j = 1: super- (Super; red), 2: near- (Near; blue), 3: sub-critical state (Sub; purple)), and a colored arrow for an internal i-j interaction points in the direction of interaction with the relative amount of flux (see details in Table 2; positive and negative values represent incoming and outgoing flux, respectively, at a critical state; outline and base colors are based on the i^th critical and j^th critical states, respectively). E represents the internal nucleus environment. B) An early flux dynamic event in the HRG response resulting from the HRG interaction flux dynamics (see C) are shown. Interaction flux dynamics i< = j (or i = >j; color based on j) represent the interaction flux from the j^th critical state to the i^th critical state or vice versa. The interaction fluxes (see HRG in C; the flux direction changes at y = 0) align to suppress the cyclic state flux at 10min (the first point in C), where the interaction flux shows 1< = E, 1 = >2, and 1 = >3 at the super-critical state, 2< = E, 2< = 1, 2 = >3 at the near-critical state, and 3 = >E, 3< = 1, 3< = 2 at the sub-critical state. They then align to enhance the cyclic state flux at 45 min (5^th point in C). This change in the dynamic flux structure is due to the global perturbation at 15-20min (see Fig 14; section VI). At each node, the net flux (Fig 11; 10min: first point; 45min: 5th point) is indicated as IN for net incoming flux (y >0), OUT for outgoing flux (y <0), or Balance (y~0). Note: The average flux balance at each node is maintained, but not at individual time points. C) Notably, for MCF-7 cell fates (HRG and EGF), near-synchronous interaction flux dynamics are seen at sub-critical states. The plots further show that the overall patterns are similar between the same cell types, which again supports cell-type-specific SOC control. Dashed lines represent interaction flux dynamics for i-j (color based on j).

Fig 13

SOC control mechanism of overall gene expression in terminal cell fates (MCF-7 and HL-60 cells): A) Flow regime shows how a dynamical change in criticality (critical behaviors) at the CP in terminal cell fates affects the entire genome (thick yellow arrow with a black dashed outline) as follows (Super: super-; Near: near-; Sub: sub-critical state): The dynamical change in criticality (critical behaviors) that originates from the CP at Near perturbs the net interaction flux (2–3 + 3–2; refer to Fig 12A) at Sub-Near, which determines the net self-flux of Sub (B: right panels), and thus directly perturbs the Sub-response. This induces a perturbation in the net interaction flux of Sub-Super (1–3+3–1), i.e., directly perturbs the Super-response (B: left panels). Furthermore, these perturbations on Sub and Super disturb the dominant cyclic flux between Sub and Super, which in turn has an impact on sustaining the critical dynamics (source: Sub and sink: Super); solid thick arrows represent average fluxes (Table 2). Note: the net self-flux of a critical state represents the effective driving force acting on CM, and thus dynamic interactions among states determine the coherent oscillatory dynamics of Sub and Super (Fig 11). This schematic picture also shows that the erasure of an initial-state criticality at the genome-state change reflects the destruction of the initial-state SOC control on the dynamical change in the entire genomic system (i.e., pruning procedure for regulating global gene expression at the initial state).

Fig 14. Local and global perturbations in self-organizing genome-wide expression: The kinetic energy self-flux dynamics (y-axis) for the CM of a critical state (see section VI) exhibit clear energy-dissipative behavior.

Notably, the results show the occurrence of global and local perturbations of self-organization. Pulse-like global perturbations show a transition from IN to OUT net kinetic energy flux or vice versa (IN-OUT switching) in more than one critical state: at 15-20min, HRG; at 12-18h, DMSO; and at 2-4h (significant) and 12-18h, atRA. In contrast, local perturbation is observed in the EGF response for up to 36h: there is only a marked response in the super-critical state, and almost no response in the other states (i.e., the dynamics of the CM of critical states are localized around their average: Fig 11). The results suggest that the global and local perturbations differentiate MCF-7 cell fates, whereas global perturbations drive the state change in HL-60 cells (see section I).