Attractor Landscapes and Network Dynamics of iPSC-Derived Neural Stem Cell Plasticity: A Complexity-Theoretic Benchmark of Cell-Fate Transitions

1. Introduction

Waddington's epigenetic landscape, also referred to as developmental landscape or cell fate decision landscape, provided an intuitive geometric metaphor and complex systems model for developmental plasticity dynamics: cells move through a multi-dimensional state-space orchestrated by complex interactions and feedback loops of underlying gene regulatory (or multiomic) networks, whose valleys canalize progressively constrained fates [1]. Kauffman's Random Boolean Networks (RBNs) subsequently supplied a mechanistic complex systems framework for identifying stable and transient gene-regulatory network configurations, and their dynamical network features including phenotypic attractors, transition paths, and perturbational switching [2,3]. In modern systems medicine and complex systems theory, these ideas have been distilled into the proposition that cell phenotypes can correspond to high-dimensional attractor states of gene regulatory networks [4], and that fate commitment may involve destabilization of one attractor before entry into another, or “pathological attractors” in disease-context [5]. By pathological attractors, we imply the cell fates can be stuck along the trajectories/canalizations prior to terminal differentiation or be trapped in a deep basin in an alternative ridge, valley, or transition region. One of the central goals of systems medicine is to algorithmically decode the re-steerability or reprogrammable of these ‘stuck’ attractor pathologies and help identify candidate network biomarkers for “differentiation therapies”. Stem-cell self-renewal and differentiation are therefore complex systems processes for the application of nonlinear dynamics, nonequilibrium physics, network theory, and multi-stability rather than only clustering with statistical inference pseudotemporal methods and differential expression analysis seen in traditional genomics data science [6].

Single-cell RNA sequencing now resolves the phenotypic plasticity and heterogeneity of differentiating cultures at high dimensionality, but most datasets remain destructive snapshots: each cell is measured once, and cells collected at different times are not observed lineages nor cross-matched. Cell fate trajectory inference methods can order cells or reconstruct branching structures using some a priori assumed statistical models, yet their outputs depend on assumptions about network topology, sampling, and progression [7]. This motivates complementary systems medicine approaches that ask not only where cells lie on a differentiation manifold, but how the observed representation learning, including dynamical features such as population geometry, local vector field(s), network (re)organization, and candidate in silico perturbations jointly steer/reprogram or constrain future cell (attractor) states.

To advance complex systems approaches in the study of stem cell fate decision-making, here we treat neural differentiation as a complex dynamical system and apply a systems medicine workflow to the GSE138121 human iPSC-derived motor-neuron dataset [8], whose differentiation protocol follows established developmental patterning principles [9]. The primary analysis was intentionally restricted to the Batch A control line at day 12 and day 18 so that the computational problem was a focused developmental-state transition rather than an ALS case-control analysis from this primary, publicly available dataset. We combine biologically informed neural-state signatures with Non-Negative Matrix Factorization (NNMF) (to identify meaningful gene metaprograms/modules) and Leiden community discovery (for modularity detection) [10], continuous quasi-potential landscapes, forward-time transport, network centrality, Block Decomposition Method (BDM)-based algorithmic perturbation [11,12,13], Kauffman’s Boolean attractors, and model-based in silico steering. This complexity-informed computational design is conceptually related to prior attractor-landscape and Boolean-control frameworks such as ATLANTIS [14], analyses of Boolean differentiation landscapes [15], and BENEIN [16], but integrates these elements with neural transition forecasting and explainable AI (XAI)/deep learning for identifying predictive cell fate transition signatures.

Our goals were to: (i) quantify the redistribution of plastic stem/progenitor and differentiated neural states; (ii) characterize dynamical systems features such as the underlying basin structure, local fixed-point behavior, dwelling time, and transition accessibility/phenotypic switching; (iii) identify network and algorithmic complexity-informed perturbation biomarkers of the differentiation program; (iv) test whether Boolean perturbations reveal candidate steerability and control variables; and (v) benchmark diverse predictive models of barycentric future state, including Echo-state network (ESN)/reservoir computing [17], Neural ODEs [18], Transformers [19], physics-informed hybrids [20], XGBoost [21], and the Continuous Thought Machine [22]. Integrated Gradients (IF) [23], implemented through Captum-compatible attribution tooling [24], was used alongside network and Boolean evidence to form consensus steering in silico candidates/hypotheses.

Figure 1. Complexity-theoretic workflow for neural stem-cell plasticity modeling. The primary GSE138121 analysis proceeds from QC and hierarchical neural-state annotation to metacell construction, forward-time barycentric transition targets, continuous attractor landscapes, network/BDM analysis, Kauffman-style Boolean dynamics, and a multi-model transition benchmark with explainable-AI steering. The workflow explicitly treats cross-sectional edges as inferred transport rather than lineage-traced parent-daughter relationships.

Figure 1. Complexity-theoretic workflow for neural stem-cell plasticity modeling. The primary GSE138121 analysis proceeds from QC and hierarchical neural-state annotation to metacell construction, forward-time barycentric transition targets, continuous attractor landscapes, network/BDM analysis, Kauffman-style Boolean dynamics, and a multi-model transition benchmark with explainable-AI steering. The workflow explicitly treats cross-sectional edges as inferred transport rather than lineage-traced parent-daughter relationships.

2. Methods

2.1. Dataset and Primary Analysis Scope

GSE138121 contains single-cell transcriptomes from human iPSC-derived motor-neuron differentiation cultures generated across several experimental batches and cell lines [8]. The public GEO record includes Batch A samples at multiple differentiation stages. The present primary analysis used the Batch A control line encoded in the analyzed cell identifiers as CTR_0083 and retained only day 12 (D12) and day 18 (D18). The archived result set contained 950 D12 cells and 1,010 D18 cells after quality control, for a total of 1,960 cells. Because this is a secondary computational analysis of publicly available, de-identified data, no new participant recruitment or intervention was performed. The analysis is therefore a two-snapshot differentiation study of a single control iPSC line, not a longitudinal single-cell lineage experiment. All transition probabilities, dwell times, and inferred vector fields are population-level model estimates.

2.2. Preprocessing, Feature Selection, and State Annotation

Count data were log-normalized. The modeling space was formed from the union of highly variable genes and genes represented in curated neural-state signatures, yielding 815 genes in the archived run. Signature programs covered pluripotent, neural stem/progenitor, cycling progenitor, neuroblast/immature neuron, motor neuron, interneuron, excitatory neuron, astrocyte, oligodendrocyte lineage, and neural crest/mesenchymal states. For each cell, signature scores were used to assign a fine-state label and to derive hierarchy_class, coarse_program, branch_class, developmental_tier, plasticity_class, and a binary plasticity endpoint. The binary plasticity endpoint separated stem/progenitor-like from non-stem future outcomes for an auxiliary prediction task. The hierarchy was deliberately multi-resolution: the fine state preserved neurodevelopmental heterogeneity, the branch label represented broader developmental direction, and the plasticity tier distinguished plastic stem/progenitor states, intermediate states, and committed differentiated states.

2.3. Metacells and Unsupervised Communities

Cells were aggregated into 62 metacells (target size approximately 35 cells) to reduce single-cell sampling noise and stabilize state-space geometry. Non-negative matrix factorization (NNMF) was evaluated across candidate ranks, with best k=5 in the archived primary run. Leiden community detection was evaluated across resolutions 0.40, 0.70, 1.00, 1.30, 1.60, and 2.00 [10]. Agreement between unsupervised partitions and curated fine-state labels was quantified using adjusted Rand index, normalized mutual information, homogeneity, and completeness.

2.4. Forward-Time Transition Graph and Barycentric Future Targets

A directed forward-time graph was constructed between metacells using only edges consistent with increasing experimental time, yielding 300 temporal edges. Because multiple future metacells can be plausible neighbors of a source metacell, the supervised transition target was not an arbitrary single edge. Instead, each source metacell was assigned a weighted barycentric future latent state and a soft empirical future-state distribution aggregated over its forward neighbors. This converts destructive cross-sectional sampling into a probabilistic transport target while avoiding the claim that any edge is an observed lineage.

For prediction, the latent representation was reduced to up to 24 principal components and standardized. The benchmark jointly predicted latent displacement, fine future-state distribution, coarse branch, plasticity tier, and future non-stem probability. A training-only state/time Markov prior was estimated within each fold and, where used, its blend weight with the learned model was selected only on the validation subset of that fold. Both raw model metrics and prior-blended metrics were retained.

2.5. Continuous Landscape, Basins, Dwell Time, and Fixed Points

A two-dimensional landscape was constructed from the leading state-space coordinates. Relative cell/metacell density was converted to a quasi-potential, U=-log(P/Pmax), and an empirical drift field was estimated from forward-time transitions. Streamlines, three-dimensional quasi-potential surfaces, and directed state trajectories were rendered from the observed metacell geometry. Candidate fixed points were identified from local speed minima and local Jacobian eigenstructure, while high-density minima were additionally retained as attractor-like density seeds. Basin assignments were summarized by occupancy, a depth proxy, minimum and median U, two-dimensional fractal-dimension estimates when defined, and mean drift speed.

State transition matrices were used to calculate self-transition probability, expected dwell steps, transition entropy, effective outgoing-state number, and approximate dwell time in days. Approximate mean first-passage time (MFPT) matrices were also generated. Because the implementation caps effectively inaccessible transitions at a very large value (6×10^6 days in the archived table), MFPT was interpreted qualitatively rather than as a literal biological time for unreachable state pairs.

2.6. Gene-Transition and Gene-Landscape Associations

Genes were scored for association with source-to-target transition edges and for correlation with quasi-potential, drift speed, attractor score, and repeller score. Benjamini-Hochberg false-discovery-rate correction was applied where available. In the primary run, no robust repeller-specific fixed-point population was resolved and repeller-score correlations were undefined for many genes given the constraints and limitations of the dataset. The lack of time-resolved data sufficient to explore such complexity and dynamical features remained a major limitation. Thus, manifold geometry-based gene interpretation focused on quasi-potential, drift, attractor-associated, or mixed/inverse associations supported by finite statistical inference methods presented herein.

2.7. Spearman Networks and Algorithmic-Complexity Perturbation

A global and state-conditioned gene network was constructed using Spearman correlations among selected modeling genes. The primary global network used an 180-gene analysis set. Node descriptors included degree, absolute correlation strength, signed mean neighbor correlation, betweenness, weighted closeness, eigenvector centrality, PageRank, coreness, local clustering, and a consensus centrality percentile. Network visualizations were exported.

Algorithmic information contribution was estimated using the Block Decomposition Method (BDM) [11,12] through pyBDM [13]. For node perturbation, the BDM of the binary adjacency representation was recomputed after removal of a node. For edge perturbation, BDM perturbation signature was recomputed after deletion of an individual edge. The absolute change in BDM was used as a descriptive measure of how strongly the graph's algorithmic description changed under that perturbation. BDM approximates algorithmic (Kolmogorov) complexity, and hence, serves as a complexity estimator of the minimum description length generating the complex network/system. High-information nodes at binarization thresholds were therefore treated as candidate structural control points or steerable intervention points (algorithmic perturbation signatures) requiring experimental validation.

2.8. Kauffman Boolean Dynamics and in Silico Steering

A signed Boolean abstraction was constructed from 180 genes and 1,080 regulatory edges. Multiple seed states were simulated until convergence to fixed or cyclic attractors. The resulting landscape was summarized by observed attractor count, basin fractions, basin entropy, transient length, cycle length, convergence fraction, and one-bit Derrida sensitivity [2,3]. A sensitivity near one is conventionally associated with the ordered/chaotic boundary in random Boolean networks; values appreciably above one indicates damage amplification. Accordingly, the observed value was described as supercritical or damage-amplifying relative to this heuristic rather than as proof of criticality.

Candidate steering genes were evaluated through knockout, overexpression, node deletion, and influential edge deletion. Perturbation influences integrated changes in basin distribution, attractor number, basin entropy, Derrida sensitivity, transient length, cycle length, convergence, final Hamming damage, maximum damage, and where applicable delta-BDM. This extends attractor-control logic used in prior Boolean reprogramming frameworks [14,15,16] to the present neural differentiation setting.

2.9. Predictive Model Benchmark

The archived V1.4 benchmark produced evaluable aggregate results for 12 models: a persistence baseline, a training-only annotated-state/time Markov prior, ExtraTrees, XGBoost [21], ESN/reservoir computing [17], residual MLP, Neural ODE [18], graph neural dynamics, context-aware Transformer [19], a Fokker-Planck 2D PINN hybrid, a reaction-diffusion graph PINN, and CTM [22].

Five-fold cross-validation was used. The transition model used latent dimension 24, up to 360 epochs, batch size 80, learning rate: 2×10^-4, AdamW with weight decay 10^-4, ReduceLROnPlateau (factor 0.5; patience 8; minimum learning rate 2×10^-6), early-stopping patience 55, and gradient-norm clipping at 3.0. The joint training objective combined robust latent-displacement loss, directional cosine loss, soft fine-state cross entropy, coarse-branch cross entropy, plasticity-tier cross entropy, and binary future non-stem loss. The reaction-diffusion and Fokker-Planck hybrids included additional derivative/PDE consistency penalties. Transformer and CTM inputs received small Gaussian noise during training as regularization.

The ESN used a fixed recurrent reservoir with reservoir size 512, spectral radius 0.92, leak 0.35, sparsity 0.92, and ridge readout penalty 1.0. ExtraTrees used 450 regression trees and 500 classification trees with minimum leaf size 2. XGBoost used 350 estimators, depth 4, learning rate 0.035, subsampling 0.85, column subsampling 0.80, minimum child weight 2, L2 regularization 2.0, and L1 regularization 0.02. The context-aware Transformer used gene patches of size 16, model dimension 128, depth 3, 8 attention heads, and dropout 0.15. Neural ODE and graph-dynamics hidden widths scaled with latent dimension, with minimum widths of 96. CTM used the Sakana AI Continuous Thought Machine core with a context-fusion head for current latent state, predecessor context, and elapsed time [22].

Metrics included endpoint MSE/MAE/R2, delta-state errors and cosine similarity, fine/coarse/plasticity classification metrics, future non-stem classification, branch-switch metrics, soft cross entropy, Brier score, Jensen-Shannon divergence, and transition-kernel RMSE/row-wise JS divergence. AUROC/AUPRC were reported only when held-out folds contained the required positive and negative classes (undefined values were not replaced by imputed scores).

2.10. Explainability and Consensus Steering

Integrated Gradients [23] were calculated across deep models/folds and aggregated to gene-level mean absolute attribution. The final consensus steering score combined six evidence streams when available: landscape geometry strength, transition-gene score, Boolean perturbation influence, model Integrated Gradients, classical-model feature importance, and network centrality. Percentile-normalized evidence streams were averaged over available components, and the number of supporting evidence streams were retained for each gene.

2.11. Gene Set Enrichment and Pathways

To test whether the prioritized network/steering genes mapped onto the predefined neural-state programs and broader biological pathways, we analyzed the nominal top-30 and top-50 Boolean node-removal |Delta BDM|, Boolean intervention, model-integrated IG, attractor-network consensus and archived six-stream consensus rankings. Fine-state and coarse-program enrichment used the same curated neural marker panels as the annotation workflow within the 180-gene Boolean universe. Functional over-representation testing used frozen GO Biological Process 2023, Reactome 2022 and KEGG 2021 Human libraries. One-sided hypergeometric tests were BH-adjusted across all terms within each signature/library/top-k family. The displayed pathway analysis used the library annotation universe, while the exact 180-gene network was evaluated as a conservative background sensitivity analysis. These tests provide functional annotation and cross-method concordance, without implying independent lineage or causal validation.

2.12. Software and Reproducibility

The pipeline was implemented in Python using NumPy, pandas, SciPy, scikit-learn, PyTorch, matplotlib, NetworkX, python-igraph, leidenalg, XGBoost, pyBDM, Captum, einops, and related utilities. See Data and Code Availability section.

3. Results

3.1. Differentiation Strongly Redistributes Plastic Stem/Progenitor and Neuronal States

The primary analysis contained 950 D12 cells and 1,010 D18 cells. D12 was dominated by cycling progenitor (428 cells) and neural stem/progenitor (231 cells) states, with additional pluripotent (44), astrocyte (59), oligodendrocyte-lineage (75), excitatory-neuron (37), neural-crest/mesenchymal (36), neuroblast/immature-neuron (18), motor-neuron (17), and interneuron (5) annotations. By D18, neuroblast/immature-neuron cells dominated (723), accompanied by motor neurons (102), interneurons (73), excitatory neurons (36), astrocytes (29), and smaller residual progenitor/pluripotent populations.

At the coarse level, the stem/progenitor branch fell from 703 cells at D12 to 16 at D18, whereas the neuronal branch increased from 77 to 934. Plasticity-tier annotations similarly shifted from 703 plastic cells at D12 to only 16 at D18, while the intermediate category increased from 18 to 723. These changes establish a strong developmental reorganization over the six-day interval and motivate treating the system as a state-transition problem rather than a static cluster comparison.

Unsupervised analyses captured related but non-identical structure. NNMF selected k=5 and showed the strongest agreement with the curated state labels (adjusted Rand index 0.301; normalized mutual information 0.424; completeness 0.694). Multi-resolution Leiden solutions produced normalized mutual information values of approximately 0.33-0.39. Thus, unsupervised communities were informative about latent organization but did not trivially reproduce the ten biologically defined states.

Figure 2. Redistribution of neural states and developmental hierarchy between day 12 and day 18. (A) Observed fine-state composition from the archived analysis. (B) Coarse developmental-branch composition. (C) Plasticity-tier composition. (D) Three-dimensional hierarchy-class landscape. The dominant shift is from plastic stem/progenitor states at day 12 toward neuronal and neuroblast/immature-neuron states at day 18.

Figure 2. Redistribution of neural states and developmental hierarchy between day 12 and day 18. (A) Observed fine-state composition from the archived analysis. (B) Coarse developmental-branch composition. (C) Plasticity-tier composition. (D) Three-dimensional hierarchy-class landscape. The dominant shift is from plastic stem/progenitor states at day 12 toward neuronal and neuroblast/immature-neuron states at day 18.

3.2. Continuous Landscape Analysis Resolves Four Basins and Directional Funneling

The quasi-potential analysis resolved four basins containing 32, 15, 6, and 9 metacells, corresponding to occupancies of 51.6%, 24.2%, 9.7%, and 14.5%. Basins 0 and 2 showed the largest depth proxies (0.954 and 0.966), whereas basin 1 was comparatively shallow (0.066). Mean drift speed varied markedly across basins, from approximately 2.10 to 19.78 in the latent-space units used by the pipeline.

Local fixed-point analysis returned multiple saddle-like points and four attractor-like density seeds concentrated in high-density regions of the landscape. No robust repeller-like fixed-point population was resolved in the primary run. Consistent with this, gene-repeller correlation fields were undefined in the exported association table. We therefore interpret the continuous landscape as evidence for multi-basin funneling and local saddle structure.

Gene-landscape associations nevertheless separated genes tracking stable/high-density regions from genes tracking movement or inverse geometry. ENO1, HMGA1, PLP1, MDK, HES5, HNRNPD, NPM1, ZFP36L1, SLC2A1, and MGST1 showed strong positive relationships with drift/attractor-oriented geometry, whereas TUBA1A, DCX, L1CAM, STMN2, CDKN1C, MAFB, SYT4, and RGS4 displayed strong inverse or mixed associations. These patterns provide a developmental landscape-level distinction between candidate state-maintenance and state-transition programs.

Figure 3. Continuous attractor-landscape architecture. (A) Three-dimensional quasi-potential surface. (B) Quasi-potential contours with empirical transition streamlines and candidate fixed points. (C) Three-dimensional transcriptomic state landscape with directed transition edges. (D) Basin occupancy, depth proxy, and scaled mean drift speed. Four basins were recovered; the fixed-point table contained saddle-like points and attractor-like density seeds but no robust repeller-like fixed-point class.

Figure 3. Continuous attractor-landscape architecture. (A) Three-dimensional quasi-potential surface. (B) Quasi-potential contours with empirical transition streamlines and candidate fixed points. (C) Three-dimensional transcriptomic state landscape with directed transition edges. (D) Basin occupancy, depth proxy, and scaled mean drift speed. Four basins were recovered; the fixed-point table contained saddle-like points and attractor-like density seeds but no robust repeller-like fixed-point class.

3.3. Inferred Transport Converges Strongly on a Neuroblast/Immature-Neuron Destination

The forward-time transition kernel was highly asymmetric. Neuroblast/immature-neuron was the dominant non-self destination from several source states, with probabilities of 0.667 from excitatory neuron, 0.527 from neural crest/mesenchymal, 0.515 from pluripotent, 0.501 from neural stem/progenitor, 0.494 from motor neuron, 0.494 from astrocyte, 0.492 from cycling progenitor, and 0.486 from oligodendrocyte lineage. The inferred transport therefore formed a broad funnel into an immature neuronal region rather than a diffuse all-to-all transition pattern.

Most non-absorbing states had expected dwell times of approximately 6.0-6.75 days under the inferred kernel. Interneuron and neuroblast/immature-neuron had self-transition probability 1.0 in the estimated transition matrix, making their expected dwell time infinite under the discrete Markov approximation. This should be understood as sink-like behavior in the fitted kernel. MFPT values further supported rapid relative accessibility of the neuroblast/immature-neuron state, whereas many unsupported transitions reached the algorithmic cap and were treated as effectively inaccessible.

Transition-gene analyses repeatedly recovered neuronal maturation, cytoskeletal, progenitor, and patterning genes across multiple source-target edges. Among genes appearing most frequently in the top transition sets were STMN2, NEFL, NEFM, VIM, MDK, HMGA1, SLC2A1, HOXB-AS3, MAFB, STMN4, DCX, HOXB2, L1CAM, TUBA1A, CRABP2, ZFP36L1, PRTG, and MLLT11. For the astrocyte-to-neuroblast/immature-neuron transition, high-scoring genes included DPPA4, PFN1, CRABP2, NASP, KIF1A, HOXA5, HOXA3, HOXB2, RUFY3, MEIS1, DUSP4, NOVA1, and AUTS2.

Figure 4. Inferred state-transition dynamics. (A) Transition-probability matrix. (B) Approximate expected dwell time for states with finite estimates. (C) Log10-transformed approximate MFPT matrix; capped values represent effectively inaccessible transitions in the implementation and should not be read as literal biological times. (D) Genes recurrently ranked among the top 20 transition genes across source-target edges.

Figure 4. Inferred state-transition dynamics. (A) Transition-probability matrix. (B) Approximate expected dwell time for states with finite estimates. (C) Log10-transformed approximate MFPT matrix; capped values represent effectively inaccessible transitions in the implementation and should not be read as literal biological times. (D) Genes recurrently ranked among the top 20 transition genes across source-target edges.

3.4. Network Topology and BDM Perturbation Identify Candidate Structural Control Points

The global Spearman network identified SEMA6D, SERINC1, TERF2IP, MLLT11, MAP2, ONECUT2, MAFB, MAP1B, RPS6, RPS3, ELAVL4, RPS27A, RUFY3, BACH2, DCX, ATP1B1, PRTG, MAPK8, and ZFP36L1 among the highest consensus-centrality nodes. This list spans neuronal morphogenesis, cytoskeletal remodeling, transcriptional control, translation, and stress-response processes, suggesting that the inferred differentiation network is organized around both lineage-specific and general cellular-control modules.

BDM perturbation provided a complementary ranking based on the algorithmic description of network structure rather than graph centrality alone. Removal of SERINC1 produced the largest absolute delta-BDM (679.9), followed by TERF2IP (502.8), MAP1B (489.9), MLLT11 (477.0), TUBA1A (468.3), NEFL (432.4), SEMA6D (431.2), CRMP1 (411.6), NPM1 (392.5), and ZFP36L1 (369.9). High-information edges included CKS2-PTTG1, LIN28A-RPSAP58, NEFL-BEX2, MLLT11-RTN1, TUBA1A-STMN1, MAP2-NEUROD1, NCAM1-PGM2L1, TUBA1A-TUBB3, MLLT11-INA, and TUBA1A-HOXB2. The partial overlap between BDM and centrality rankings indicates that topological prominence and algorithmic-information contribution are related.

Figure 5. Network biology and algorithmic-complexity perturbation. (A) Global Spearman gene network with leader-line callouts identifying the archived top-five centrality (C), node |ΔBDM| (B), model Integrated Gradients (I), and edge |ΔBDM| (E1-E5) signatures. (B) Top consensus-centrality genes. (C) Top pyBDM node information contributions measured as absolute change in network BDM after node removal. Concentric node rings retain overlapping C/B/I memberships; edge callouts terminate at edge midpoints. Network IG highlights use the archived maximum attribution across models. Complete callout values are provided with the figure data.

Figure 5. Network biology and algorithmic-complexity perturbation. (A) Global Spearman gene network with leader-line callouts identifying the archived top-five centrality (C), node |ΔBDM| (B), model Integrated Gradients (I), and edge |ΔBDM| (E1-E5) signatures. (B) Top consensus-centrality genes. (C) Top pyBDM node information contributions measured as absolute change in network BDM after node removal. Concentric node rings retain overlapping C/B/I memberships; edge callouts terminate at edge midpoints. Network IG highlights use the archived maximum attribution across models. Complete callout values are provided with the figure data.

3.5. Boolean Dynamics Reveal Multistability and Damage Amplifying Sensitivity

The Kauffman-style Boolean model contained 180 nodes and 1,080 edges and produced 18 observed attractors across the simulated seed ensemble. Basin entropy was 2.75 bits, mean transient length was 3.5 steps, median transient length was 2 steps, mean cycle length was 1.15, and all simulated seeds converged. The largest Boolean basin mapped to cycling progenitor and contained 41.9% of seeds; the second largest mapped to neuroblast/immature neuron and contained 27.4%. The remaining attractors had substantially smaller basins.

One-bit Derrida sensitivity was 1.69 (SD 1.71). In the classical Kauffman heuristic, sensitivity near 1 marks the boundary between ordered and damage-amplifying regimes. A value of 1.69 is therefore better described as supercritical or damage-amplifying/propagating (consistent with one-step damage propagation) than simply 'edge of chaos'. It indicates that one-bit perturbations expanded on average over one update in this fitted Boolean abstraction.

Such one-step damage amplification is compatible with chaotic-like or dynamically unstable propagation, but it is not by itself evidence of “deterministic chaos”, a strange attractor, or self-organized criticality: demonstrating chaos would require persistent divergence across successive updates rather than expansion after only one perturbation [25,26]. In classical Kauffman Random Boolean Networks, K denotes the number of regulatory inputs influencing each node; increasing connectivity generally promotes perturbation propagation, whereas lower K favors ordered dynamics, with the critical boundary depending on Boolean-function bias and update scheme [2,25,27]. For the canonical unbiased random Boolean case, the annealed approximation predicts K_C ≈2 separating ordered from damage-amplifying regimes [25]. Our empirically weighted threshold networks therefore cannot be mapped directly onto this classical K_C criterion, and because a finite deterministic Boolean system ultimately enters fixed points or periodic cycles, we do not infer classical continuous-state “strange attractors.” Rather, the observed supercritical Derrida sensitivity suggests a damage-amplifying regime that could permit rapid phenotypic switching and transitional adaptive behaviors among network states or attractor basins, motivating future time-resolved perturbation experiments. Time-resolved multi-omics, multi-step Derrida analyses, asynchronous-update tests, and longitudinal measurements to determine whether this computational instability corresponds to chaotic-like transitions, enhanced phenotypic plasticity, or an experimentally supported critical regime [28,29].

Boolean perturbation influence prioritized ONECUT2, ZFP36L1, TMEM123, ENO1, TYMS, IPO5, RPL13A, CENPF, CCNB1, MAP1B, DCX, TUBA1A, and ELAVL4. For example, ONECUT2 node deletion changed the basin distribution (total variation 0.516), increased attractor count by one, increased Derrida sensitivity by 0.687, and produced a large negative delta-BDM, yielding the highest influence score in the archived perturbation table. These computational results nominate genes whose intervention may reorganize attractor occupancy rather than merely correlate with a particular state.

Figure 6. Boolean attractor dynamics and in silico steering. (A) Signed Boolean gene network with the same C/B/I node and E1-E5 edge callout convention as Figure 5; arrowheads show the directions of highlighted Boolean edges. (B) Largest observed attractor basins and mapped signature labels. (C) Highest perturbation influence scores across knockout, overexpression, and node-deletion experiments. (D) Summary of Boolean dynamical statistics. Derrida sensitivity above one indicates damage amplification in this model; instead of evidence of criticality.

Figure 6. Boolean attractor dynamics and in silico steering. (A) Signed Boolean gene network with the same C/B/I node and E1-E5 edge callout convention as Figure 5; arrowheads show the directions of highlighted Boolean edges. (B) Largest observed attractor basins and mapped signature labels. (C) Highest perturbation influence scores across knockout, overexpression, and node-deletion experiments. (D) Summary of Boolean dynamical statistics. Derrida sensitivity above one indicates damage amplification in this model; instead of evidence of criticality.

3.6. Predictive Models Capture Complementary Aspects of Future-State Transition Dynamics

The archived benchmark produced aggregate results for 12 models across five folds, with six held-out source examples per fold in the barycentric transition benchmark. Fine-state classification was saturated for the Markov prior, CTM, ESN, ExtraTrees, Transformer, and XGBoost after validation-selected prior blending, and reached 1.0 for most remaining learned models after blending. This saturation, together with the small held-out sample size and undefined AUROC/AUPRC for folds lacking valid one-vs-rest classes, means that fine-state accuracy is not a reliable basis for ranking model quality in this experiment.

Latent endpoint regression and transition-kernel fidelity were more discriminative. ESN/reservoir computing achieved the lowest endpoint MAE at 0.138; however, its endpoint R2 remained negative (-1.03), showing that even the best absolute-error model did not explain variance relative to the fold baseline under this stringent endpoint representation. Among prior-blended transition kernels, CTM achieved RMSE 0.0402 and Transformer 0.0404, followed by the Markov prior at 0.0420, Neural ODE at 0.0436, graph neural dynamics at 0.0445, and reaction-diffusion graph PINN at 0.0447.

Raw, unblended kernel metrics provide additional context. CTM and Transformer remained strong at 0.0432 and 0.0439, respectively, whereas raw ESN kernel RMSE was 0.169. Thus, the ESN's strength was direct latent endpoint regression rather than probability-kernel reconstruction. The validation-selected state/time prior materially improved several model probability outputs, which is why raw and blended metrics are both retained in the supplementary data.

Branch-switch accuracy ranged from 0.967 to 1.000 for most prior-blended learned models, while future non-stem accuracy was 1.000 for all non-persistence evaluable models. These simpler hierarchical endpoints were therefore substantially easier than precise latent-state forecasting. The model comparison supports a multi-metric interpretation: reservoir computing was strongest on endpoint MAE, whereas CTM and Transformer best reproduced the empirical transition distribution.

Figure 7. Predictive benchmarking and explainable-AI prioritization. (A) Endpoint MAE for the best-performing evaluable models. (B) Raw and validation-selected prior-blended transition-kernel RMSE; CTM and Transformer show the strongest kernel fidelity. (C) Top genes by mean absolute Integrated Gradients among their entries in the archived per-model top-50 lists. (D) Final consensus steering genes integrating geometry, transition evidence, Boolean perturbation, model attribution, classical feature importance, and network centrality.

Figure 7. Predictive benchmarking and explainable-AI prioritization. (A) Endpoint MAE for the best-performing evaluable models. (B) Raw and validation-selected prior-blended transition-kernel RMSE; CTM and Transformer show the strongest kernel fidelity. (C) Top genes by mean absolute Integrated Gradients among their entries in the archived per-model top-50 lists. (D) Final consensus steering genes integrating geometry, transition evidence, Boolean perturbation, model attribution, classical feature importance, and network centrality.

3.7. Explainability and Multi-Stream Consensus Converge on Candidate Steering Genes

The archived per-model top-50 Integrated Gradients lists prioritized NEFM, ONECUT2, NEFL, CRABP1, MLLT11, CPE, STMN2, HOXB-AS3, PBX3, SLC1A2, TYMS, DCX, ELAVL4, ANOS1, VIM, NEUROD1, and HOXB2 among the highest mean attributions among their available top-50 entries. These genes span neuronal structural maturation, transcriptional specification, retinoid biology, progenitor-to-neuron transition, and glial/neural lineage programs.

The final consensus steering score integrated up to six evidence streams. The top ten genes were ONECUT2 (0.890; six streams), ZFP36L1 (0.865; six), DCX (0.863; six), LOC400043 (0.838; five), HOXB2 (0.836; five), TMEM123 (0.830; six), TUBA1A (0.823; six), ENO1 (0.820; six), ELAVL4 (0.816; six), and MAP1B (0.815; six). PRTG, MLLT11, NCAM1, SLC2A1, TOP2A, and INA also ranked highly. Because these scores combine independent geometry, network, perturbation, and predictive evidence, they are better interpreted as prioritized intervention hypotheses than as conventional marker rankings.

3.8. Explainability and Multi-Stream Consensus Converge on Candidate Gene Set Enrichment Links Cell Fate Steering Signatures to Source and Destination Programs

The neuroblast/immature-neuron program was significantly enriched in the attractor-network consensus (DCX, ELAVL3, ELAVL4, STMN2 and TUBB3; BH q = 0.01353) and even more strongly in IG (q = 0.001548). At the coarser level, both the plastic stem/progenitor source program (q = 0.02721) and the intermediate-neuroblast destination program (q = 0.009017) were enriched, consistent with a transition-control rather than endpoint-only ranking. Exploratory pathway annotation further identified Boolean Delta BDM enrichment for axon guidance (Reactome q = 0.007243) and nervous-system development (q = 0.009507), IG enrichment for post-chaperonin tubulin folding (q = 0.0116), and attractor-consensus enrichment for regulation of stem-cell proliferation (GO q = 0.01297). None of these pathway-level signals survived the stricter exact 180-gene background, so they are retained as exploratory functional annotations, whereas the neural state panel enrichment itself was tested directly within the 180-gene network (Supplementary Figure S7-S8).

4. Discussion

This complex systems modelling study asks whether iPSC-derived neural stem-cell differentiation can be modeled as a structured dynamical system whose state geometry, transition kernel, network organization, and perturbational response carry complementary information. The primary result is a coherent shift from a plastic stem/progenitor-rich D12 population to a D18 population dominated by neuroblast/immature-neuron and other neuronal identities. This redistribution is mirrored by the inferred transition kernel, continuous landscape, Boolean basins, and several independent gene-prioritization layers.

The observed funnel toward neuroblast/immature neuron is consistent with the notion that developmental progression can canalize heterogeneous starting states toward a smaller set of accessible phenotypic basins [1,4,5,6]. The strongest statement supported by the data is not that a single transcriptomic attractor has been experimentally proven, but that multiple computational representations independently identify asymmetric state-space organization. In the continuous landscape, density minima and saddle structure describe local geometry. In the Boolean abstraction, discrete attractors and basin fractions provide a complementary coarse-grained dynamical picture. Their agreement around cycling/progenitor and immature neuronal configurations suggests that these regions are major organizing states of the sampled differentiation window.

Kauffman's Boolean network formulation provides a powerful systems medicine framework for cell fate dynamics inference because it connects explicit gene regulatory network update rules to recurrent states, their basins of attraction, and responses to perturbation analyses [2,3]. In this complex systems view, a phenotype is a coordinated configuration sustained by interacting genes, and a fate transition is a change in the system's dynamical organization. Hence, Kauffman characterizes cell fates as network attractors. Huang et al. provided experimental support by showing convergence of neutrophil differentiation trajectories initiated by distinct stimuli [4], while Huang, Ernberg, and Kauffman developed the broader interpretation of normal and abnormal phenotypes as network attractors [30]. Applied to neural differentiation dynamics, this complexity framework allows trajectory inference to ask which configurations are stable under a specified model, which initial states reach them, and which perturbations alter accessibility. Basin occupancy and perturbational switching thus provide testable dynamical hypotheses beyond a geometric ordering of transcriptomes.

The study by Zhou et al. illustrates the value of examining multiple routes through gene-expression state space [31]. Sixteen non-cytotoxic compounds induced MCF7 cells to leave a proliferative state, with transcriptomic trajectories that diverged before converging toward differentiated configurations. Although that experiment concerns cancer cells, its methodological lesson is relevant here: similar destination phenotypes/cell fates can be reached through distinct transient responses, so endpoint markers alone need not identify the processes controlling a transition. Our combination of inferred transport, Boolean basin perturbation, network centrality, and BDM perturbation signatures is intended to distinguish such candidate roles in neural differentiation and neural stem cell dynamics. The density-derived quasi-potential remains a descriptive projection, and attractors in an inferred Boolean model do not establish biological attractors or causal gene control without time-resolved perturbation experiments.

We used BDM perturbation analysis to identify information signatures, i.e., genes and network interactions, whose removal substantially changes the inferred network’s algorithmic information content, thereby supporting in silico causal hypothesis generation. Kolmogorov complexity, (K), is the length of the shortest program that reproduces an object on a fixed universal computer; BDM provides a computable estimate of this otherwise uncomputable quantity. In the spirit of minimum description length (MDL), changes in BDM after node or edge deletion probe a component’s contribution to a compact description of network structure, without constituting an exact MDL optimization. This rationale connects algorithmic intervention-based inference and reprogramming [12] with the broader computational concept of behavioural reprogrammability [32]. It also extends the network-based identification of candidate cell-fate regulators discussed by Uthamacumaran and Craig [33]. Combined with Boolean basin perturbations, these signatures prioritize candidate genes and interactions for testing whether interventions can steer differentiation toward alternative attractors; structural information changes alone do not establish biological causality or the direction of a fate transition. We did not resolve transition times through mixed-phenotype states, test partial reprogramming, or determine whether intermediate configurations represent stable spurious attractors. This distinction matters because Lang et al. [34] modeled partially reprogrammed cells as hybrid attractors combining existing cell-identity programs, while Mikkelsen et al. [35] documented incomplete repression of lineage programs in cells trapped in partially reprogrammed states. Such trapping should be distinguished from deliberate partial reprogramming through transient, cyclic expression of Oct4, Sox2, Klf4, and c-Myc (OSKM), which can modify cellular properties without requiring complete conversion to pluripotency [36]. These findings motivate future time-resolved perturbation experiments to assess reprogramming extent, cell fate reversibility, and unintended attractor occupanc, such as via CRISPR screens and pharmacological or RNA interference experiments.

The network and perturbation signatures results extend conventional marker analysis by separating several notions of importance. A gene can be central in a correlation network, highly informative under BDM removal, influential on Boolean basin structure, or strongly attributed by a predictive model. ONECUT2, ZFP36L1, DCX, MAP1B, TUBA1A, ELAVL4, MLLT11, TMEM123, and related genes repeatedly appear across these layers. Such convergence is conceptually aligned with control-oriented Boolean approaches such as ATLANTIS and BENEIN, which seek network interventions capable of moving a system toward desired attractors [14,16]. The present analysis does not claim validated reprogramming control; it identifies candidates for targeted perturbation experiments.

Complex-systems tools also connect these candidate control points to disease modelling. Uthamacumaran and Craig [33] combined attractor landscape reconstruction, network science and algorithmic-complexity analysis to investigate glioblastoma plasticity and glioma stem cell states. BENEIN [16] used Boolean-network modelling to identify differentiation-control targets; combined inhibition of MYB, HDAC2 and FOXA2 promoted enterocyte-like differentiation and reduced malignant properties in colorectal cancer models. Together, these studies motivate using attractor geometry, network dynamics and modelled interventions to nominate in silico systems medicine candidates for steering cell-fate decisions. Testing such candidates in donor-derived iPSC disease models could support reprogramming strategies and personalized genomic medicine, subject to experimental validation.

The enrichment analysis adds a functional layer to this convergence. The strongest state-level signal is the neuroblast/immature-neuron signature, matching the dominant D18 population and the sink-like destination of the inferred transition kernel, while simultaneous enrichment of the plastic stem/progenitor program indicates that the prioritized network also contains genes associated with the departing source state. The pathway annotations similarly emphasize nervous-system/axon-development and microtubule/tubulin-remodeling processes together with stem-cell proliferation control. These signatures suggest the identified network biomarkers’ cell fate steerability towards teleonomic endpoints (fate commitments) of neural stem cell dynamics. This pattern is compatible with a network transition from proliferative/progenitor organization toward immature neuronal structure, but it is not independent validation: the marker panels contributed to the original state labels, the pathway tests annotate genes selected from the same transcriptomic experiment, and pathway significance was background sensitive. Experimental and functional perturbation with denser longitudinal or lineage-resolved measurements is therefore required to determine whether these in silico candidates causally alter basin accessibility or neuronal differentiation.

The algorithmic-complexity layer is especially non-conventional. BDM perturbation evaluates how removal of a node or edge changes an approximation to the algorithmic description length of the network [11,12,13]. This is not equivalent to degree or betweenness: a structurally informative node can matter because it participates in a pattern that strongly changes the graph's generative regularity. The high BDM contribution of SERINC1, TERF2IP, MAP1B, MLLT11, TUBA1A, NEFL, and SEMA6D therefore supplies a distinct hypothesis about network reprogrammability. Related algorithmic-information approaches have been proposed for perturbational causal discovery and network reprogramming [12], while broader dynamical-systems perspectives have argued that attractor analysis can reveal otherwise hidden organization in biological networks [37,38].

Two of the highest BDM candidates provide a specific link to chromatin regulation. TERF2IP (Spearman |ΔBDM| rank 2) encodes RAP1, which modulates TIP60/p400 histone acetyltransferase activity and pluripotent-state transitions in mouse embryonic stem cells [39]. MAP1B (rank 3) associates with the BRG1-containing chromatin-remodelling complex in neural stem cells. Its nuclear enrichment supports the stem-cell state, whereas its cytosolic function promotes neuronal differentiation [40]. These findings motivate chromatin-focused perturbation experiments, while the present transcriptomic analysis does not itself establish epigenetic mechanisms or causal control.

The Boolean analysis also illustrates why careful language around criticality is essential. Kauffman networks motivate an ordered-to-chaotic transition with sensitivity near one [2,3]. The observed Derrida sensitivity of 1.69 is not 'near one'; it indicates average one-step damage amplification. We therefore interpret it as a supercritical, flexible Boolean regime rather than direct evidence that the biological system sits precisely at the edge of chaos. Demonstrating biological criticality would require stronger finite-size, perturbational, and scaling analyses.

The predictive benchmark further shows why model evaluation should match the biological question. Fine-state classification became nearly trivial for several models because the barycentric targets are highly structured and the test folds are small. In contrast, latent endpoint error and transition-kernel fidelity remained informative. ESN/reservoir computing, which uses a fixed recurrent reservoir and learned readout, achieved the best endpoint MAE, while CTM and Transformer best reproduced the future-state distribution. CTM is particularly interesting in this context because it explicitly represents neuron-level temporal processing and synchronization [22], but its use here should be viewed as an architectural benchmark rather than a claim of mechanistic equivalence between artificial-neural synchronization and cellular differentiation.

Physics-informed models require similar precaution. The reaction-diffusion and Fokker-Planck hybrids impose dynamical regularization inspired by transport equations [20], but the transcriptomic state space is not a physical spatial medium with a known governing PDE. Their role here is to test whether structured dynamical inductive biases, including morphogenetic pattern formation systems, improve prediction. They should not be interpreted as establishing a literal reaction-diffusion law for gene expression dynamics nor the plasticity spectrum. The Neural ODE offers a more general continuous-time inductive bias [18] and performed competitively on transition-kernel fidelity.

Relative to static differential-expression workflows, the chief advantage of the present framework is its separation of state, transition, and control hypotheses. State markers describe occupancy; transition genes describe changes along inferred edges; landscape associations describe position relative to density and drift; centrality describes topological embedding; BDM captures algorithmic structural contribution; Boolean perturbation tests attractor reorganization; and IG identifies genes used by predictive models. Their intersection yields a more stringent candidate set for future perturbational validation.

4.1. Limitations

Several limitations materially constrain interpretation. A central limitation is the lack of sufficiently resolved, independently annotated neuronal iPSC time-course data in the present analysis to validate multistate trajectories across donors and differentiation protocols. GSE247600 samples a different neural lineage, but its unresolved annotation prevents it from closing this validation barrier. First, the primary analysis uses a single Batch A control iPSC line (CTR_0083) at only D12 and D18. The large cell count does not substitute for independent biological replicates. Second, the cells are destructively sampled; there is no observed cell-to-cell lineage continuity. The directed graph and barycentric targets are therefore inferred population transport. Third, the 62-metacell benchmark yields small held-out test sets, causing several classification metrics to saturate and some AUROC/AUPRC values to be undefined. Fourth, the continuous fixed-point analysis did not recover a robust repeller-specific class, so repeller-gene associations were not interpretable in the primary run. Fifth, MFPT contains capped values for inaccessible transitions and should be used qualitatively. Sixth, the Boolean network is an inferred coarse-graining whose update rules and perturbational effects require experimental testing. Finally, the external GSE247600 processing recovered only one signature state in the current annotation pipeline, preventing its use as a multistate validation cohort.

These limitations suggest clear next steps: analyze additional iPSC lines and batches; use datasets with more developmental time points; exploit lineage tracing or RNA-velocity-compatible measurements where available; evaluate donor/line-held-out generalization; remap GSE247600 Ensembl identifiers to gene symbols before signature annotation; and experimentally perturb top consensus candidates to test whether predicted steering effects alter state occupancy or differentiation efficiency.

Targeted CRISPRi/CRISPRa experiments in donor-derived iPSCs could initially test SERINC1, ONECUT2 and ZFP36L1, alongside TERF2IP and MAP1B. This panel combines the leading Spearman BDM candidate with the two highest-ranked Boolean perturbation and consensus-steering genes and chromatin-linked candidates. Pairing these interventions with time-series lineage tracing and epigenetic profiling across donors and disease genotypes could test effects on differentiation, maturation and reprogramming. MAP1B experiments should also assess subcellular localization [40]. This would evaluate whether the predicted control points can improve donor-specific iPSC disease modelling and support personalized genomic medicine.

5. Conclusion

Neural stem-cell differentiation in the analyzed GSE138121 control-line subset is organized by a strong redistribution from plastic progenitor states toward immature neuronal states, a multibasin continuous landscape, multistable Boolean dynamics, and a reproducible set of network and steering candidates. The analysis nominates ONECUT2, ZFP36L1, DCX, HOXB2, TMEM123, TUBA1A, ENO1, ELAVL4, MAP1B, MLLT11, and related genes for experimental follow-up. Predictive benchmarking showed complementary strengths: ESN minimized latent endpoint MAE, while CTM and Transformer most closely reproduced the transition kernel after appropriate prior-blending audit. More broadly, the framework demonstrates how complexity theory, network medicine, algorithmic information, Boolean control, and explainable machine learning can be combined to move from descriptive single-cell snapshots toward explicit, testable hypotheses about cell-fate dynamics and reprogrammability. With experimental validation across donors and disease genotypes, these candidate interventions could inform iPSC-based disease modelling and reprogramming strategies for personalized genomic medicine.

6. Data and Code Availability

The primary source dataset is publicly available through NCBI GEO accession GSE138121 [8]. The supplementary feasibility dataset is available through GEO accession GSE247600 [41]. A GitHub public code repository for our original codes and computational models will be made available post-publication.

7. Ethics Statement

This study is a secondary computational analysis of publicly available, de-identified transcriptomic datasets and involved no new recruitment or intervention.

8. Competing Interests

The authors declare no competing interests.