Multi-omic enriched blood-derived digital signatures reveal mechanistic and confounding disease clusters for differential diagnosis

Introduction

Blood serves as an unparalleled medium for monitoring systemic health, reflecting both localized pathology and organism-wide physiological changes. The vital role of blood is also evidenced by its critical function in the prevention, diagnosis, and management of chronic diseases, providing critical insights into ongoing physiological and pathological processes (Cabalar et al., 2024). Routine blood tests, such as complete blood counts and biochemical panels, remain among the most frequently ordered investigations worldwide due to their accessibility, standardization, and diagnostic relevance (Hoffbrand, 2024). The hematological and biochemical composition of blood captures dynamic signatures of immune activity, metabolic status, and organ dysfunction, thereby offering a high-resolution snapshot of disease states across multiple physiological systems (Liotta et al., 2003). In recent years, the concept of the digital twin has emerged in healthcare, referring to computational representations of biological entities that are continuously updated with real-world data (Sun et al., 2023; Vall´ee, 2023; Erol et al., 2020). Applications have ranged from cardiovascular modeling to surgery planning and medical imaging analysis (Sun et al., 2023; Li et al., 2025), highlighting the potential of computational models to refine disease stratification, predict progression, and guide therapeutic decisions. Conventional disease classification systems, such as the International Classification of Diseases (ICD), have long provided a standardized framework for diagnosis, billing, and epidemiology (World Health Organization, 2015; WHO, 1992). While indispensable for healthcare communication, these systems are grounded in clinical presentation and anatomical categorization, prompting investigation into whether blood-derived digital profiles can offer complementary perspectives that may align with, extend, or challenge existing taxonomies.

Despite the ubiquity of blood testing in clinical practice, large-scale computational studies using routine blood analytes to map disease similarity remain underdeveloped. Network medicine has demonstrated the value of constructing disease–disease associations from genes, proteins, or clinical comorbidities, most prominently in the seminal “human disease network,” which revealed how shared genetic associations drive non-random clustering of disorders (Erol et al., 2020; Goh et al., 2007), and the Phenotypic Disease Network (PDN), which mapped epidemiological comorbidity patterns across populations (Hidalgo et al., 2009). However, these frameworks typically bypass standard laboratory biomarkers, leaving an unexplored opportunity for systematic analysis of blood-derived disease signatures. Although digital twin models have been successfully developed in cardiology, surgery, and medical imaging (Sun et al., 2023; Li et al., 2025), very few studies operationalize routine blood analytes as the basis for disease-level digital twins (Chu et al., 2023), despite their accessibility, cost-effectiveness, and interpretability relative to omics or imaging-based data. At the same time, ICD-style taxonomies are based on symptoms, anatomy, or pathology rather than biological mechanisms (World Health Organization, 2015; WHO, 1992), potentially obscuring latent immuno-inflammatory connections or conflating clinically similar but mechanistically distinct conditions. These gaps highlight the need to determine whether blood-derived digital signatures can recover established clinical categories, reveal mechanistic overlaps across systems, and refine precision medicine frameworks.

The objectives of this study are threefold: (i) to quantify inter-disease similarity by constructing a standardized matrix of 103 disease signatures from blood analytes and applying correlation and clustering methods; (ii) to evaluate whether clinically meaningful categories such as hematopoietic, metabolic, endocrine, or respiratory disorders exhibit coherent clustering or cross-system overlap; and (iii) to identify key analytes and shared mechanisms driving disease stratification using enrichment analysis and dimensionality reduction, thereby uncovering biological links that transcend conventional classification boundaries.

Our findings show that blood-derived, disease-level digital signatures reliably recover hematopoietic coherence while revealing biologically unified yet clinically heterogeneous disorders via shared immune and inflammatory pathways. These findings demonstrate the potential of analyte-based similarity to complement traditional taxonomies and inform a more mechanistic, data-driven disease ontology.

Methods

Aims and objectives

The methodological design of this study was structured to translate the overarching research objectives into a reproducible analytical framework. Specifically, the aims were to integrate heterogeneous disease signatures into a standardized matrix that enables direct comparison across 103 conditions, quantify pairwise similarity between diseases using correlation-based measures and a phylogenetic framework, interrogate the biological and clinical relevance of resulting clusters through enrichment analyses and categorical grouping, and apply machine learning, feature selection, and dimensionality reduction techniques (PCA, UMAP, and Random Forest) to validate clustering patterns and identify key analytes driving disease stratification. Together, these steps operationalize the research objectives by moving from raw disease profiles to standardized data, quantifying relationships, evaluating biological meaning, and refining insights with advanced computational methods.

Data source and preprocessing

The dataset consisted of 103 synthetic disease signatures that we created based on medical knowledge and physiological reference ranges. We created profiles with hematological and biochemical trends across multiple timepoints for each disease with five consecutive analyte timepoints. An illustrative generic example is in Table 1, which contains measurements for white blood cells (WBC), lymphocytes, monocytes, segmented neutrophils, eosinophils, red blood cells (RBC), hemoglobin, mean corpuscular volume (MCV), platelet count, and mean platelet volume (MPV) across five sequential timepoints.

TABLE 1. Representative set of longitudinal time series from blood marker values for a disease profile illustrating the structure of five sequential clinical timepoints modelled as random variables.

Synthetic data generation and statistical structure

To enable cross-disease comparison, disease profiles were constructed as synthetic hematologic signatures derived from curated clinical knowledge of expected laboratory abnormalities. More precisely, the data points were extracted from literature and textbooks and curated by a junior clinician (final-year medical student with formal training in healthcare technologies) and a senior clinician (board-certified physician in oncology). Each disease was represented as a time-series vector of analyte values across five sequential timepoints (t1–t5) for each of the 11 analyte markers: WBC, lymphocytes, monocytes, segmented neutrophils, eosinophils, basophils, RBC, haemoglobin, MCV, platelet count, and MPV, capturing progression or variation in haematological changes.

The data were preprocessed into a standardized disease–analyte matrix using Python (Pandas, Regex, and OS libraries). A canonical analyte list was defined (including WBC, RBC, hemoglobin, MCV, platelet count, MPV, among others), and analyte names in each TXT file were normalized (e.g., trimming whitespace and resolving naming inconsistencies) using regular expressions (Friedl, 2006). This ensured that the same analyte recorded under slightly different labels would be consistently aligned across profiles.

Analyte values were assigned using clinically informed ranges rather than purely stochastic sampling, scaled to reflect expected deviations (e.g., low, normal, elevated) relative to physiological reference ranges, and propagated across timepoints to reflect plausible temporal trajectories (such as progressive anemia, inflammatory escalation, or recovery). Values were parsed into numerical format; when a particular analyte was absent, placeholder values (None) were assigned to preserve dimensional consistency.

For each disease profile, a time-series feature vector was constructed by concatenating all standardized analytes across five timepoints (Analyte t1–t5). Profile identifiers (e.g., Profile 01–103) were extracted and numerically sorted using regular expressions. The final output was a matrix (.csv) containing 103 rows (disease profiles) and 55 columns (11 haematological markers × 5 timepoints), providing a uniform representation of heterogeneous disease signatures for downstream correlation analysis, clustering, enrichment testing, and dimensionality reduction. This matrix follows an 11 × 5 analyte-by-time feature structure, where 11 refers to the standardized analyte set and 5 refers to the five simulated progression stages used in the analysis. Where an analyte was unavailable for a disease, the corresponding five-timepoint block was left empty.

Biological variability and noise

To approximate real-world variability, small perturbations were introduced across timepoints within each disease profile for noise-based feature space perturbation analysis in principal component analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) space. However, the dataset primarily represents structured, knowledge-driven signals rather than fully stochastic patient-level variability.

Correlation structure between analytes

Physiological relationships between analytes, such as RBC–hemoglobin–hematocrit coupling and inflammatory co-variation, were implicitly preserved through coordinated assignment of values across related variables. While these dependencies were not explicitly modeled using covariance matrices, they arise from shared clinical constraints embedded in the construction process. These correlations were therefore not explicitly simulated using statistical structures but instead emerge from clinically informed co-variation patterns embedded during profile construction.

Interpretation of timepoints

The five timepoints do not correspond to fixed clinical intervals or repeated measurements from a single real patient. Rather, they represent stages of disease progression or variation in silico, allowing a dynamical systems trajectory-like behavior to be encoded in analyte space rather than treated as static snapshots. These progression states were generated by assigning analyte values across timepoints based on clinically expected trajectories (e.g., recurrence or recovery patterns) derived from medical knowledge and reference ranges, thereby approximating disease evolution rather than longitudinal patient-specific measurements.

Handling of missing values

Missing analytes were initially encoded as null entries and subsequently handled prior to analysis by numerical coercion and imputation. Median imputation was used in the primary analysis because of its robustness for sparse and non-normally distributed clinical-style matrices, and alternative imputation strategies were evaluated in supplementary sensitivity analyses. Imputation was performed prior to all downstream analyses (including Pearson correlation, PCA, and machine learning models), ensuring full numerical compatibility of the dataset.

Pairwise correlation

To quantify disease–disease similarity, pairwise correlations were computed across all 103 disease profiles. Each profile was represented by its standardized analyte vector encompassing measurements across five clinical timepoints. Similarity was assessed using the Pearson correlation coefficient ρ_ij, defined as:where ρᵢ_j is the Pearson correlation between disease profiles i and j, ranging from −1 (perfect inverse association) to +1 (perfect concordance). Here, x_ik represents the value of analyte k in disease profile i, x̄_i is the mean analyte value in profile i, y_jk is the value of analyte k in disease profile j, and ȳ_*j* is the corresponding mean. The index k iterates over all analytes in the standardized disease–analyte matrix (Pearson, 1896).

ρ ᵢ j = Σ k x ᵢ k − x ¯ ᵢ y j k − ȳ j √ Σ k x ᵢ k − x ¯ ᵢ 2 × Σ k y j k − ȳ j 2

Clustering

The resulting 103 × 103 symmetric correlation matrix provided a full mapping of similarity among disease signatures. To convert similarity to a distance representation, each correlation coefficient was transformed using d_ij = 1 − ρ_ij. This transformation ensures that highly similar diseases appear closer together, while weakly correlated or divergent profiles appear more distant. The distance matrix was then subjected to hierarchical agglomerative clustering using the Unweighted Pair Group Method with Arithmetic Mean (UPGMA), also known as average linkage clustering (Sokal et al., 1958). This procedure iteratively merges the closest clusters based on average inter-cluster distance, producing a binary phylogenetic tree. Two dendrograms were generated: a full-range tree showing the complete hierarchy and a zoomed-range tree (0–0.07) focusing on fine-grained structure at high correlation thresholds. Both were annotated with disease identifiers to allow direct inspection of relationships. The y-axis of the dendrogram encodes clustering distance (1 − ρ_ij), enabling interpretation of relative proximity among profiles. Together, this combination of correlation analysis and hierarchical clustering provides a reproducible, interpretable framework for mapping disease similarity and comparing it with conventional taxonomies.

Statistical Validation of Hierarchical Clustering (SVHC)

To formally evaluate the statistical robustness of the hierarchical clustering structure, we applied Statistical Validation of Hierarchical Clustering (SVHC), a bootstrap-based framework that assigns statistical significance to clusters within a dendrogram. We implemented a Python-based svhc package developed by Bongiorno et al. (2022). The SVHC algorithm evaluates cluster stability by generating resampled datasets and computing the frequency with which specific clusters are recovered across bootstrap iterations, followed by false discovery rate (FDR) correction to account for multiple testing. The input matrix consisted of the fully numeric 103 × 55 disease-by-feature time-series matrix defined above (i.e., 103 rows corresponding to disease profiles and 55 columns corresponding to 11 haematological markers each measured at 5 timepoints). Clustering validation was performed using 1000 bootstrap resamples with a significance threshold of α = 0.05. Clusters were retained only if they satisfied FDR-adjusted significance criteria, thereby providing a statistically grounded alternative to arbitrary dendrogram cut thresholds.

Enrichment and cluster analysis

Cluster extraction

Hierarchical clustering produced a dendrogram in which diseases were grouped according to their pairwise correlation distances. To define discrete clusters, the tree was cut at a threshold of h = 0.02 on the distance scale d = 1 − ρ, corresponding to a minimum within-cluster similarity of ρ ≥ 0.98. This cut yielded 16 clusters in total, and the cluster assignments were exported in tabular form.

Disease ontology annotation

Cluster 9—the largest and most heterogeneous group—was selected for deeper analysis. Each disease in this cluster was manually annotated using Disease Ontology (DO) categories (Schriml et al., 2012), with the underlying mapping retrieved from the Disease Ontology database in July 2025 from the website (https://disease-ontology.org/) (Project TDO, 2025). Each condition was assigned both a clinical category (e.g., hematopoietic, metabolic, endocrine, respiratory, renal, cardiovascular, gastrointestinal/hepatobiliary, rheumatologic/autoimmune, or infectious) and a broader super-category to support higher-level interpretation. These labels were later used to contextualize PCA/UMAP projections and Random Forest outputs.

Gene set construction from DisGeNET

To explore biological mechanisms associated with Cluster 9, gene–disease associations were retrieved from DisGeNET (Pin˜ero et al., 2020) (downloaded July 2025) (Project TDO, 2025). For each disease, records included disease name, gene symbol, score, and NCBI identifier. Gene entries were merged into a non-redundant list based on Entrez IDs and treated as a binary gene set for downstream enrichment, providing a knowledge-based mechanistic context.

Pathway over-representation analysis

Exploratory over-representation analysis (ORA) was performed using both KEGG (Kanehisa and Goto, 2000) and Reactome (Gillespie et al., 2022) databases (accessed July 2025), implemented in R via clusterProfiler (Xu et al., 2024) and ReactomePA (Yu and He, 2016). Both packages use the hypergeometric test with Benjamini–Hochberg (BH) correction (Benjamini and Hochberg, 1995), and q < 0.05 was considered significant. The background universe was defined as the union of all DisGeNET-retrievable genes for Cluster 9 diseases (project TD, 2025). Output included enrichment ratios, gene counts, and adjusted p-values.

Reproducibility

Cluster assignments, enrichment outputs, and analysis parameters were exported as structured tables. All R scripts documented database versions, species (“Homo sapiens”), and thresholds to ensure reproducibility. Interpretation emphasized recurrent pathway themes, given regular database updates.

Dimensionality reduction and feature selection

The standardized disease–analyte matrix was merged with category and super-category labels to form a unified dataset of 103 profiles. Missing values were imputed using median substitution, a robust technique for sparse clinical matrices (Little and Rubin, 2019). For linear structure, PCA was applied (Jolliffe and Cadima, 2016) and the first two components were visualized with profiles coloured by super-category. To capture non-linear structure, UMAP was applied (McInnes et al., 2025), generating 2-D embeddings preserving both local and global similarity profiles. Related nonlinear manifold-learning approaches, such as t-distributed stochastic neighbor embedding (t-SNE) can also be used, to enable low-dimensional representation of complex high-dimensional structure (van der Maaten and Hinton, 2008). For supervised feature selection, a Random Forest classifier (Breiman, 2001) with 500 trees was trained on labelled profiles, and feature importance was computed via mean decrease in Gini impurity. The classifier was also used to generate predicted super-category labels for “Unknown” profiles. This step was performed independently of the downstream machine learning–based validation analyses, which employed multiple classifiers and cross-validation as described below.

Robustness, effect-size, and machine-learning validation analyses

To assess clustering robustness, we compared multiple unsupervised configurations across imputation strategies, distance metrics, linkage rules, and cluster numbers. Specifically, median, most-frequent, and k-nearest-neighbour imputation were evaluated in combination with Euclidean, cosine, and correlation distances and average or Ward linkage. For each candidate configuration, silhouette scores, Davies -Bouldin and Calinski -Harabasz indices were calculated as clustering quality metrics, and the best-performing configuration was retained for downstream validation. Bootstrap branch support was estimated for the dendrogram, and the fully imputed matrix was exported for statistical dendrogram validation.

To quantify pairwise disease similarity more rigorously, we performed correlation analyses with uncertainty estimation. Pearson correlations were accompanied by Fisher-transformed confidence intervals and effective sample sizes, and sensitivity analyses were also conducted using Spearman and Kendall rank-based statistics. Feature-level contrasts were evaluated by comparing hematological disease profiles against non-hematological profiles, and each disease system against the remainder of the dataset. Standardized effect sizes were calculated using Hedges’ g, with statistical significance assessed using Welch’s t-test where relevant, after a Shaprio-Wilk test and Levene test assessed normality and homogeneity of variance. Multiple testing across features within each contrast was controlled using the Benjamini–Hochberg procedure (See Supplementary Table S2).

Post-hoc machine learning classifiers: multiclass validation

For post hoc machine learning (ML) validation and Algorithmic Information Dynamics (AID) analyses, a complementary feature matrix was derived from a structured clinical knowledge resource used to obtain our disease-feature matrix, which encoded the 103 disease conditions across 10 organ systems (haematological, nephrology, endocrine/metabolic, cardiovascular, infective, gastrointestinal/hepatobiliary, rheumatological/vasculitis, neurological, respiratory, and oncology). Each condition was represented by binary and categorical features reflecting the utility of the full blood count (FBC) under three sampling paradigms: longitudinal sampling without a predefined signature, single sampling with a defined signature, and longitudinal sampling with a defined signature, across four diagnostic contexts: screening, possible diagnosis, diagnosis with additional tests, and disease monitoring. Analyte names were normalized using regular expressions, and the resulting structured matrix was used as the input feature space for all downstream multiclass supervised classification and AID perturbation analyses, providing a richer, ontology-informed representation of disease categories than the analyte-only profile matrix used for primary clustering. The system-level labels derived from this resource served as ground-truth class labels for classifier training and evaluation. To evaluate generalizability of the learned structure, repeated stratified multiclass ML classification was performed using Random Forest (RF), XGBoost, support vector machine (SVM) with radial basis function (RBF) kernel, and multilayer perceptron (MLP) models. ML algorithms were implemented in Python using scikit-learn and gradient boosting classifier.

The ML analysis used repeated stratified 5-fold cross-validation with 5 repeats (25 folds total), and performance was summarized using balanced accuracy, macro-F1 score, and one-vs-rest macro- Area Under the Receiver Operating Characteristic (AUROC) as ML performance metrics. Hyperparameters were fine-tuned a priori: Random Forest used 800 trees with balanced class weights and minimum leaf size of 2; SVM-RBF used standardization followed by an RBF-kernel classifier with C = 1.0, gamma = scale, and balanced class weights; XGBoost used 500 trees, maximum depth 4, learning rate 0.03, subsample 0.9, colsample_bytree 0.9, and regularization λ = 1.0; and the MLP used two hidden layers (128, 64), ReLU activation, Adam optimization, early stopping, and a maximum of 600 iterations. In addition, leave-group-out sensitivity analyses were performed by removing hematology-related feature groups and reassessing silhouette scores, and explainable feature robustness/importance analysis was further evaluated using Gini importance, permutation importance, and SHAP (SHapley Additive exPlanations) (Lundberg and Lee, 2017).

Sensitivity analysis: imputation, noise, and embedding stability

To further evaluate the statistical robustness of the feature space and clustering structure, we performed a multi-stage sensitivity analysis incorporating imputation strategy, stochastic embedding variability, and noise perturbation analysis.

Missing data were handled using three imputation strategies implemented via scikit-learn: k-nearest neighbors (KNN), median imputation, and most-frequent imputation. Following imputation, all features were standardized and subjected to PCA using NumPy and SciPy-backed routines to assess global variance structure and separability across clinical categories. To characterize nonlinear manifold stability, UMAP was performed using the umap-learn Python library under multiple random initializations (seeds from 0 to 20). Default parameters were retained (n_neighbors and min_dist consistent across runs), isolating stochastic initialization as the primary source of variation (See Supplementary Information).

Noise robustness was assessed by injecting Gaussian noise into the feature matrix with σ (sigma) defined relative to feature standard deviation (i.e., sigma scaling applied proportionally across features). For each noise level, embeddings and variance structure were recomputed, and stability was quantified via preservation of cluster topology and variance explained. All visual analyses were conducted in Python using pandas, matplotlib, and seaborn for visualization.

Algorithmic information dynamics

Algorithmic Information Dynamics (AID) was implemented as an in silico causal inference post hoc network analysis using the pyBDM package to estimate changes in graph algorithmic complexity under feature perturbations. Features were binarized across multiple thresholds and evaluated using block decomposition method (BDM), Shannon entropy, and compression metrics (gzip, bz2, lzma, zlib) implemented via Python’s standard libraries. Perturbation-based sensitivity was quantified as the change in information content (bits) following feature-wise perturbations, and features were ranked based on mean absolute shifts in information content across binarization thresholds and complexity measures.

Results

Hierarchical clustering of the 103 disease profiles generated a dendrogram summarizing pairwise similarities across the dataset. Two complementary visualizations were produced: a full-scale tree covering the complete distance range (0.0–0.6; Figure 1) and a zoomed tree restricted to the high-similarity region (0.0–0.07; Figure 2). Together, these representations allow inspection of both global and fine-grained structure within the disease space. At the full scale, the tree demonstrates the broad heterogeneity of disease signatures, with branches spanning a wide spectrum of distances. Certain profiles appear as long, isolated branches, reflecting weak correlations with the majority of other conditions, in contrast to several groups of profiles that form tightly linked subtrees with minimal branch lengths, indicating high correlation and internal consistency. The zoomed view highlights the fine resolution of clustering at short distances, where closely related profiles merge into compact groups, emphasizing the presence of disease pairs or small clusters with very high similarity (ρ > 0.93), which are not as apparent in the full-scale visualization. This level of detail reveals that while some categories exhibit highly cohesive clustering, many profiles are distributed more diffusely across the tree, pointing to variable degrees of internal relatedness.

Cluster analysis outcomes

Cutting the hierarchical tree at a distance threshold of 0.02 produced 16 discrete clusters encompassing all 103 disease profiles. The distribution of cluster sizes was highly uneven, ranging from large, heterogeneous groupings to single-profile clusters. Most profiles (n = 80, 77.7%) aggregated within a single cluster (Cluster 9), which encompassed a wide spectrum of conditions spanning hematopoietic, metabolic, endocrine, cardiovascular, respiratory, gastrointestinal, and immune-related categories. This concentration indicates a 3.5-fold overrepresentation relative to all other clusters combined. The breadth of Cluster 9 indicates a high level of cross-system similarity among many disease signatures, despite their diverse clinical classifications. This heterogeneity motivated further analysis of its composition and potential shared mechanisms.

Outside of Cluster 9, most clusters were small. Several comprised a single disease profile, such as Cluster 1 (acute myeloid leukemia, AML), Cluster 2 (pancytopenia), Cluster 7 (myelodysplastic syndrome), and Cluster 16 (chronic liver disease). These singleton clusters reflect profiles that were relatively distinct from the remainder of the dataset at the chosen threshold. A few clusters contained two to six profiles, often representing conditions with moderate similarity. For example, Cluster 5 grouped viral lower respiratory tract infection with bacterial pneumonia, suggesting partial convergence among respiratory infections. Cluster 6 combined six conditions spanning hematopoietic, renal, cardiovascular, and infectious categories, underscoring cross-system associations. Cluster 10 grouped together chronic hematologic malignancies (CML2, CLL2, lymphoma), which aligns with known clinical relatedness.

Thus, cluster partitioning revealed a spectrum of outcomes: a single large, heterogeneous cluster (Cluster 9), several small clusters with modestly related diseases, and multiple singleton clusters representing distinct outliers. These results confirm that while certain categories, such as hematopoietic malignancies, can form coherent groups, the majority of disease profiles cluster into a broad, mixed set. This variability sets the stage for subsequent enrichment analyses to evaluate whether the observed clusters correspond to biologically meaningful groupings. A complete list of disease profiles within each of the 16 clusters is provided in Supplementary Appendix Table S1.

Enrichment analysis findings

To investigate whether large heterogeneous clusters captured biologically meaningful relationships, enrichment analyses were performed on Cluster 9, which contained 80 disease profiles spanning hematopoietic, metabolic, endocrine, cardiovascular, respiratory, gastrointestinal, and immune-related conditions. Gene–disease associations for these profiles were retrieved from DisGeNET, and the consolidated gene set was subjected to KEGG and Reactome over-representation analysis. Results are presented in Figures 3, 4.

KEGG pathway analysis

To evaluate shared mechanisms within the largest heterogeneous cluster (Cluster 9), KEGG over-representation analysis was performed and the top results summarized in a dotplot (Figure 3A). In this figure, p-values indicate raw significance from the hypergeometric test, whereas adjusted p-values (hereafter referred to as q-values) represent values corrected for multiple testing using the Benjamini–Hochberg false discovery rate procedure. The dot colour encodes q-values (more significant results appear toward the red end of the scale), dot size reflects the number of genes mapped to each pathway, and the x-axis shows the gene ratio.

The most highly enriched KEGG pathways included Cytokine–cytokine receptor interaction, Lipid and atherosclerosis, AGE–RAGE signaling pathway in diabetic complications, Rheumatoid arthritis, Th17 cell differentiation, Leishmaniasis, Inflammatory bowel disease, and Malaria. As indicated by the colour legend in Figure 3A, q-values for these top pathways span from approximately 4.3 × 10⁻¹⁵ (e.g., Th17 cell differentiation) down to approximately 9.4 × 10⁻²³ (e.g., Cytokine–cytokine receptor interaction). This gradient confirms that all reported pathways are highly significant, with several reaching extremely small q-values. Such extreme q-values are consistent with exceptionally large enrichment effects, indicating strong biological convergence across immune and inflammatory signaling pathways. Collectively, the KEGG results capture both broad immunological signaling themes and clinically recognizable disease processes, providing a coherent mechanistic context for the cross-system clustering observed in Cluster 9.

KEGG cnetplot

Two complementary KEGG cnetplots were generated: (i) a network restricted to the top 8 enriched pathways (Figure 3A), and (ii) a broader network containing all pathways meeting a stringent significance threshold (q < 1 × 10⁻⁶; Figure 3C). In both panels, circular nodes represent pathways, smaller nodes represent genes, node size (for pathways) reflects the number of mapped genes, node colour encodes adjusted p-values (q-values), and edges indicate gene–pathway membership. Because no expression fold changes were supplied, gene node colours do not indicate differential regulation. In the top-8 cnetplot (Figure 3B), the eight pathways, such as Cytokine–cytokine receptor interaction, Lipid and atherosclerosis, AGE–RAGE signaling pathway in diabetic complications, Rheumatoid arthritis, Th17 cell differentiation, Leishmaniasis, Inflammatory bowel disease, and Malaria, form a compact, highly interconnected network. Multiple genes are shared across several of these terms, visually apparent as gene nodes linked to more than one pathway, which creates a bridge between the immune/inflammatory module (e.g., cytokine signaling, Th17 differentiation, rheumatoid arthritis, inflammatory bowel disease) and infection-related terms (e.g., leishmaniasis, malaria). The metabolic–vascular pathways (Lipid and atherosclerosis, AGE–RAGE in diabetic complications) are not isolated; instead, they connect into the same mesh via shared inflammatory mediators, consistent with the broader immune–inflammation theme observed in the dotplot (see Figure 3). The q < 1 × 10⁻⁶ cnetplot (Figure 3C) expands this view, retaining the relationships seen in the top-8 network while introducing additional significant pathways. The resulting graph is larger and denser, revealing sub-modules in which immune signaling and infection pathways cluster together through overlapping gene sets. This broader network confirms that the core structure observed in the top-8 panel is not an artifact of pathway selection; instead, it persists and becomes more granular when all highly significant terms are included. Practically, the top-8 panel offers a readable summary of the core mechanism, whereas the q-filtered panel provides a comprehensive map of pathway–gene connectivity at stringent significance.

KEGG emapplot

The emapplot (Figure 3D) provides a network physiology-based view of complex relationships among significantly enriched KEGG pathways, with edges indicating gene overlap greater than 40%. Node colour reflects adjusted p-values (q-values), node size corresponds to the number of mapped genes, and edges denote pathways sharing a substantial proportion of genes. Several functionally coherent modules emerged. First, an infection-related cluster was observed, linking Leishmaniasis, Tuberculosis, and Toxoplasmosis, all of which rely on overlapping immune defense genes. Second, an inflammation and autoimmunity module connected Th17 cell differentiation, Inflammatory bowel disease, and Rheumatoid arthritis, reflecting the central role of Th17-driven inflammatory processes in both systemic autoimmunity and gastrointestinal pathology. Third, a cytokine signaling module was evident, in which Cytokine–cytokine receptor interaction was directly linked with Viral protein interaction with cytokine and cytokine receptor, illustrating convergence between endogenous immune signaling and pathogen-driven modulation of cytokine networks. The majority of these nodes appeared in red or orange hues, consistent with extremely low q-values, indicating robust enrichment. Together, the emapplot highlights how the enriched KEGG pathways in Cluster 9 are not isolated entities but form interconnected modules centered on infection response, cytokine signaling, and chronic inflammatory processes.

Reactome pathway analysis

Reactome dotplot

The Reactome dotplot (Figure 4A) summarizes the top enriched pathways for Cluster 9. Node colour encodes adjusted p-values (q-values), node size reflects the number of mapped genes, and the x-axis denotes the gene ratio. The leading terms comprise an interleukin-centered immune axis—Signaling by interleukins, Interleukin-4 and Interleukin-13 signaling, and Interleukin-10 signaling—and a hemostasis/ECM axis—Platelet activation, signaling and aggregation, Platelet degranulation, Response to elevated platelet cytosolic Ca²⁺, Extracellular matrix organization, and Degradation of the extracellular matrix. The predominance of warm hues (red/orange) indicates very low q-values across these pathways, consistent with robust enrichment. Together, the dotplot highlights two coherent Reactome themes in Cluster 9: cytokine/interleukin signaling and platelet/ECM biology.

Collectively, these Reactome enrichment results demonstrate that the Cluster 9 gene set is organized around interconnected cytokine/interleukin signaling, platelet activation and degranulation, and extracellular-matrix remodeling pathways. Gene sharing across these Reactome pathways helps explain the mechanistic convergence of clinically heterogeneous disease profiles within Cluster 9.

Reactome cnetplot

The Reactome cnetplot restricted to the top eight pathways (Figure 4B) visualizes gene–pathway membership as a bipartite network (pathways as larger circular nodes; genes as smaller nodes; edges indicate inclusion of a gene in a pathway; pathway-node colour encodes q-value; size reflects mapped-gene count). Three subgroups are apparent. First, an interleukin module is defined by Signaling by interleukins and Interleukin-4 and Interleukin-13 signaling, which lie in close proximity and share extensive gene overlap. Interleukin-10 signaling appears as a more isolated node: although positioned near the interleukin cluster, it is linked by fewer genes and thus forms a partially distinct subnode. Second, an extracellular matrix module is represented by Extracellular matrix organization and Degradation of the extracellular matrix, which cluster tightly together. Third, a platelet module encompasses Platelet activation, signaling and aggregation, Platelet degranulation, and Response to elevated platelet cytosolic Ca²⁺, which are positioned in close proximity and indicate a coherent theme related to platelet function. This layout highlights that the Cluster 9 gene set maps onto discrete Reactome submodules, namely, interleukin signaling, platelet activation, and matrix remodeling, with partial connectivity among them rather than a single unified hub.

The Reactome emapplot (Figure 4C) present a pathway–pathway similarity network in which edges indicate substantial gene overlap (similarity cutoff ≈0.35), node colour encodes q-value, and node size reflects mapped-gene count. Three features are notable. First, a platelet triad forms a tight triangle among Platelet activation, signaling and aggregation, Platelet degranulation, and Response to elevated platelet cytosolic Ca²⁺, all rendered in red, indicating very strong significance and high mutual overlap. Second, an interleukin cluster links Signaling by interleukins with Interleukin-4 and Interleukin-13 signaling, also in red, consistent with shared upstream/downstream mediators across interleukin pathways. Third, Degradation of the extracellular matrix (blue) connects to Extracellular matrix organization (red), indicating a functional pair with asymmetric significance: both participate in matrix remodeling, but the organization term is more strongly enriched.

Dimensionality reduction insights

Principal component analysis (PCA)

The PCA projection (Figure 5A) shows that PC1 explains 92.7% of the variance, with PC2 explaining 3.3%. Accordingly, the distribution of profiles is stretched predominantly along PC1, reflecting a single dominant axis of variation in the analyte matrix. Despite this, the two-component projection does not yield clean global separation between disease systems. Hematopoietic diseases (blue diamonds) exhibit a relatively narrow vertical spread along PC2, indicating partial convergence on secondary variance dimensions, while spanning an intermediate range across PC1 rather than aligning strictly with one side of the axis. By contrast, endocrine, metabolic, cardiovascular, respiratory, gastrointestinal, and other system categories remain broadly intermixed, clustering near the center of the projection with extensive overlap. These patterns suggest that although PCA captures dominant trends in the data, the variance structure is largely one-dimensional and insufficient on its own to resolve distinct system-level groupings.

Uniform Manifold Approximation and Projection (UMAP)

The UMAP embedding (Figure 5B) yields clearer structure than PCA, revealing three visually coherent zones along a curved manifold. On the rightmost island (UMAP1 ≈ 9–10), hematopoietic system diseases (blue diamonds) form a compact, internally consistent cluster with limited dispersion along UMAP2. A central zone (UMAP1 ≈ 6–8) contains a mixture of cardiovascular, urinary, gastrointestinal, endocrine, and other profiles, indicating cross-system proximity in blood-analyte space without a single dominant class. A leftward chain (UMAP1 ≈ −2–3) shows a graded progression that includes multiple non-hematopoietic super-categories together with several “Unknown” profiles, suggesting a continuous, non-linear spectrum rather than discrete partitions. Compared with PCA, UMAP thus enhances class coherence for hematopoietic diseases and reduces overlap among non-hematopoietic groups by capturing non-linear relationships.

Random forest feature importances

A Random Forest classifier (500 trees; Gini impurity; feature importances normalized to sum to 1) was trained on profiles with known super-category labels. The top 15 features are summarized in the bar plot. Importances are distributed across a small set of hematological indices repeated over multiple timepoints, rather than concentrated in a single marker, indicating that several related analytes jointly drive discrimination. Segmented neutrophils (1000 cells/µL) appears four times among the top features (t5, t4, t1, t3), with t5 ranked first, indicating that absolute neutrophil counts carry strong class-discriminative signal. MCV (fL) appears at all measured timepoints (t1–t5) within the top 15, and RBC (million cells/µL) appears at t5, t1, and t2. Haemoglobin (g/dL) also enters the top 15 (t5). This recurrence across multiple timepoints indicates a dominant discriminative contribution of neutrophils relative to other analytes. Together, these features point to red-cell morphology and mass as major contributors to class separation. Platelet count (1000 cells/µL) is present at t4 and t5 among the highest-ranked features, complementing neutrophil and erythrocyte measures. Several late timepoints (t5) recur among the highest-ranked features, while earlier timepoints also appear, suggesting that the discriminative signal is stable across time rather than confined to a single measurement. The Random Forest highlights a hematology-centric signature: neutrophils, MCV, RBC, haemoglobin, and platelets repeatedly rank at the top across timepoints. This is concordant with earlier findings from clustering and non-linear embedding, in which hematopoietic disorders showed the most coherent grouping. We note that tree-based importances reflect model-specific attributions and can be diluted by collinearity among related features; hence, rankings should be interpreted as relative contributions rather than causal effects.

Statistical Validation of Hierarchical Clustering (SVHC)

To address concerns regarding the arbitrariness of dendrogram-based clustering, we applied Statistical Validation of Hierarchical Clustering (SVHC), an unsupervised clustering algorithm, which provides formal statistical testing of cluster significance through bootstrap resampling with FDR control on the intrinsic data structure. Importantly, clustering was performed on disease profiles (n = 103) represented as time-series hematologic feature vectors, rather than on analytes. Thus, the identified clusters reflect groups of disease states sharing modularity (hierarchical organization) among similar blood-based signatures as dynamical systems. Using 1000 bootstrap iterations (α = 0.05), SVHC identified 25 statistically significant clusters (p < 0.01).

As shown in Figure 7, the validated clustering structure reveals both fine-grained subclusters and higher-order groupings, consistent with a multi-nested, multiscale organization of hematologic phenotypes. The correlation heatmap demonstrates pronounced diagonal block structures, indicating statistically robust intra-cluster similarity, alongside heterogeneous off-diagonal patterns reflecting inter-cluster relationships. These findings support the presence of modular organization within disease space, where subsets of conditions exhibit coherent hematologic profiles.

Clinically related disease categories, including inflammatory conditions, malignancies, and metabolic disorders, tended to co-localize within statistically validated modules, suggesting shared underlying biological signatures detectable from routine blood data. These results demonstrate that the observed clustering is not arbitrary, but reflects robust, statistically significant structure, reinforcing the presence of digital blood signatures as predictive biomarkers for modeling disease state-dynamics within a complex systems-level, data-driven framework and identifying candidate targets for precision therapies in personalized medicine.

Top: Hierarchical clustering dendrogram based on correlation distance between 103 disease profiles derived from time-series hematologic feature vectors. Bottom: Corresponding correlation matrix (Pearson), with red indicating positive correlation and blue negative correlation. Green boxes enriched along the diagonal denote statistically validated clusters identified by SVHC (α = 0.05, 1000 bootstrap samples, FDR-controlled). It shows 25 statistically validated clusters with distinct colour-coded branches at varying distance thresholds, and the heatmap reveals genuine heterogeneity including both strong positive (red) and negative (blue) correlations, indicating that the bootstrap resampling is discriminating real structure from noise. This modular organization supports the concept of digital blood signatures, suggesting that disease state-dynamics occupy structured regions in a biological state-space, akin to attractors in dynamical systems theory, enabling digital twin models to map, predict, and track patient-specific trajectories across these modules.

Multiclass classification and feature ablation analysis

To further quantify the robustness of disease groupings, we performed multiclass ML validation of disease separability and feature contribution analysis using simple classifiers. Across 25-fold cross-validation, ensemble methods demonstrated strong discriminative performance, with tree-based algorithms, namely, Random Forest (balanced accuracy: 0.80 ± 0.13; AUROC: 0.97 ± 0.02) and XGBoost (balanced accuracy: 0.80 ± 0.13; AUROC: 0.97 ± 0.03) outperforming SVM (balanced accuracy: 0.71 ± 0.14; AUROC: 0.96 ± 0.03) and the basic neural network, MLP (balanced accuracy: 0.40 ± 0.12; AUROC: 0.78 ± 0.13). Class-wise ROC curves (Figure 8) showed near-perfect separability for hematological and infectious categories (AUC ≈1.00), with lower performance observed in nephrology and cardiovascular classes, reflecting greater heterogeneity and potential overlap in their hematologic profiles.

Feature ablation analysis further demonstrated the contribution of hematologic features to clustering structure. Removing hematology-related features increased the silhouette score from 0.22 to 0.27, indicating that while these features strongly drive ML classification performance, their removal enhances global cluster separation by reducing dominance of hematopoietic signals. This highlights a dual role of hematologic variables: they provide high discriminative power while simultaneously compressing broader inter-system variability.

Feature importance and robustness analyses further supported these findings, demonstrating consistent identification of key discriminative features across Gini importance, permutation importance, SHAP values, and bootstrap stability (Supplementary Table S5). The convergence of multiple ML-derived feature importance metrics indicates that classification performance is driven by stable, reproducible feature contributions rather than model-specific artifacts or chance effects. Together, with our statistical structure validation, these results reinforce the interpretation of disease phenotypes as modular regions within a shared hematologic state-space, with both high discriminability and biologically coherent variability across biosystems. These findings are consistent with a complex systems view in which disease states form separable yet interconnected attractor-like regions within a phenotypic landscape (i.e., routine blood feature state-space).

Feature space sensitivity and stability analysis

Across all imputation strategies (KNN, median, most frequent), PCA demonstrated highly consistent variance structure, with PC1 explaining ∼37–38% and PC2 ∼12–13% of total variance. Visual inspection further confirmed that cluster organization and relative separation between disease groups were preserved across imputation methods, indicating minimal sensitivity to missing data handling (Figures 9A–C).

UMAP embeddings also exhibited stable global organization across random seeds. While minor local rearrangements were observed as expected due to stochastic initialization of the UMAP algorithm, the global pattern/topology, including distinct clustering of haematological and endocrine/metabolic groups and transitional positioning of gastrointestinal and nephrology samples, remained highly conserved as seen by the colocalization of distinct clusters (colors). This indicates strong manifold stability and robustness of the learned feature space in nonlinear representation space (Figures 9D–F).

Noise perturbation analysis further demonstrated resilience of the embedding structure. Increasing Gaussian noise resulted in gradual dispersion of clusters. However, core group separation persisted across moderate noise levels, suggesting that the signal captured by the feature space exceeds noise-induced variability. These findings confirm that both linear (PCA) and non-linear (UMAP) representations are robust to imputation choice, stochastic initialization, and moderate perturbations. Results in Tables 2, 3 complement Figure 9 by demonstrating that the observed pattern is robust not only to noise perturbation but also to preprocessing choices.

TABLE 2. Noise sensitivity analysis of time-series and learned global feature spaces. Gaussian noise (σ = 0.01–0.2) was applied proportionally to standardized features in both the raw 55-feature longitudinal time-series matrix and the learned global feature space derived from machine learning/AID representations. Stability was assessed using variance explained by the first principal component (PC1; and PC2 for time-series) and structural preservation of UMAP embeddings via procrustes similarity relative to the baseline (σ = 0). The time-series feature space demonstrates extremely high variance preservation (PC1 ≈ 0.93) and near-complete embedding stability across all noise levels, indicating a highly dominant low-dimensional structure. In contrast, the learned global feature space exhibits lower initial variance concentration but maintains strong embedding stability across noise levels, particularly at moderate perturbations (σ = 0.1), reflecting a more distributed yet robust representation.

TABLE 3. Effect of imputation strategy on variance structure in time-series and learned global feature spaces. Comparison of principal component variance across imputation strategies (median, most frequent, k-nearest neighbors) for both the raw 55-feature longitudinal time-series matrix and the learned global feature space derived from machine learning/AID representations. The proportion of variance explained by the first two principal components (PC1, PC2), along with cumulative variance, remains highly consistent across imputation methods in both spaces.

The time-series feature space is dominated by a single principal component (PC1 ≈ 0.93), with over 99% of variance captured by the first five components, indicating a highly compressible, low-dimensional structure in temporal ordering. This suggests that longitudinal trajectories are governed by a limited set of dominant transitions, consistent with a coarse-grained, two-state organization and attractor-like dynamics.

In contrast, the learned global feature space distributes variance more evenly (PC1 ≈ 0.38), reflecting a more complex and heterogeneous spatial organization of disease profiles. Together, these findings indicate that while temporal trajectories collapse onto a simplified, low-dimensional manifold, the global feature space captures a richer multiscale structure with greater separability between disease states. The consistency of these patterns across imputation strategies further indicates that the observed variance structure is preserved irrespective of missing data handling.

To further characterize structural feature importance beyond statistical and machine learning–based metrics, we applied an AID framework incorporating BDM, Shannon entropy, and compression-based complexity measures (gzip, bz2, lzma, zlib). Features were systematically perturbed across multiple binarization thresholds, and the resulting change in algorithmic information content was quantified. BDM perturbation analysis identified a small subset of predictive features, particularly hematology_platelet, neuro, hematology_splenomegaly, renal_urinary, and endocrine_metabolic, as exhibiting the largest BDM shifts.

From a causal inference perspective, these predictive features can be interpreted as structurally influential components of the system, or disease signatures: perturbations to these signatures induce disproportionately large changes in the global information structure of the disease-feature matrix, or graph network. This suggests that they act as potential control or regulatory elements within the hematologic state-space, consistent with an attractor-like dynamics in which specific network features govern critical transitions and stability across disease phenotypes. These findings complement our supervised feature importance and clustering results, reinforcing that the observed disease structure is driven by a subset of biologically meaningful patterns rather than random signals (See Supplementary Figure S3).

Algorithmic Information Dynamics (AID) perturbation analysis reveals a small subset of hematologic features, primarily WBC, haemoglobin, and segmented neutrophils across multiple timepoints, that exert disproportionate causal influence on the global information structure. BDM-based complexity shifts consistently dominate over entropy and compression-based metrics, indicating superior sensitivity to structural perturbations (Figures 10A–C). The recurrence of these hematological features across thresholds and timepoints among the diverse disease groups suggests temporally stable, complexity signatures robust to noise perturbations or randomness. Collectively, these results support an attractor-like self-organization of the routine blood data-derived disease feature space, where a limited set of predictive features governs transitions within a structured manifold (See Supplementary Figure S4).

Discussion

Hematopoietic clustering

A central and consistent observation across all analytical approaches was the robust clustering of hematopoietic diseases. In the phylogenetic tree, disorders such as acute and chronic leukemias, lymphomas, myelodysplastic syndromes, and anemias either formed discrete compact subclusters or appeared as isolated but internally coherent branches. In correlation-based clustering, these profiles repeatedly grouped together, in stark contrast to the dispersion of most endocrine, metabolic, or cardiovascular diseases within a large heterogeneous cluster. Dimensionality reduction reinforced this distinction: in PCA, hematopoietic conditions exhibited a relatively narrow vertical spread along PC2, indicating convergence on secondary variance dimensions, while spanning an intermediate range across PC1 without forming complete separation from other systems. In UMAP, by contrast, they formed a particularly tight and distinct cluster that was clearly segregated from the broader manifold of non-hematopoietic categories. Random Forest classification further underscored this pattern by identifying neutrophil counts, red blood cell indices (MCV, RBC), and platelet counts as the top discriminative features, pointing to classical hematological markers as the most informative drivers of disease separation.

The consistency of hematopoietic clustering across correlation metrics, enrichment, dimensionality reduction, and feature selection methods, as well as the algorithmic information dynamics-driven perturbation analyses, underscores the stability of blood-based phenotyping for these conditions. It also highlights a contrast with other complex systems: whereas hematopoietic disorders form coherent entities in analyte pattern space, endocrine, metabolic, and cardiovascular profiles appear more diffuse, reflecting greater heterogeneity and overlap. This divergence also suggests that blood-based disease similarity measures may be most directly applicable to disorders of the hematopoietic system, and that expanding their interpretability to other systems will require integrating additional modalities ranging from subjective data enrichment (signs and symptoms) to emerging multi-omics frameworks, including genetic (e.g., spatial transcriptomics), or tissue-specific precision biomarkers.

However, it should also be noted that the prominence of hematopoietic clustering observed across methods is, in part, expected given that the input data consists of blood-derived analytes, which are encompassed within hematologic processes. Yet, this structure may also be modality-dependent: blood-based features preferentially capture systemic and hematologic perturbations, while other disease classes are reflected indirectly through secondary effects, and whose biomarkers may not be as routinely accessible. Consequently, alternative data modalities, such as urinalysis, microbiological assays (e.g., swabs), electrophysiological signals (e.g., ECG), or imaging, would be expected to yield distinct but complementary clustering structures, emphasizing organ-specific or modular and multi-scale, functional relationships.

Importantly, this suggests that the observed hematopoietic cluster should be interpreted not as a universal disease taxonomy, but as a projection of disease similarity within a blood-based feature space. Extending this framework to circulating biomarkers (e.g., liquid biopsies), multi-omics profiles, and next-generation sequencing data could enable a more integrated, multimodal disease landscape, capturing both systemic and tissue-specific dynamics. Thus, our results highlight both the strength and the inherent limitations of single-modality representations, motivating future work in multimodal integration for more comprehensive disease relationship mapping.

This interpretation is further supported by algorithmic information dynamics (AID)–based perturbation analysis, in which BDM-derived perturbation shifts identify a small subset of hematologic disease signatures whose perturbation in silico induces disproportionately large changes in the global information structure, consistent with causal influence and emergent self-organization (Figure 10; Supplementary Figure S3; Supplementary Table S6). By approximating generative mechanisms via minimal description length (Kolmogorov complexity) under network perturbation, AID identifies structurally emergent features beyond statistical correlation. Complementary machine learning analyses, including Random Forest feature importance (Figures 6, 8), independently converge on the same hematologic signatures (e.g., neutrophils, RBC indices, platelets), reinforcing that the observed clustering reflects reproducible, biologically meaningful patterns rather than arbitrary feature selection. Thus, complexity metrics further demonstrate that these findings extend from descriptive clustering to causal, generative, multi-method–validated disease signatures and emergent patterns in silico.

Further, as a sensitivity analysis, we SVHC alongside effect size–based feature validation, both of which confirmed the robustness and reproducibility of the observed clustering structure. SVHC demonstrated that clusters are statistically stable under bootstrap resampling with FDR control, supporting a non-arbitrary, biologically grounded organization of disease phenotypes. These findings are consistent with an attractor-like network physiology, in which hematologic states self-organize into hierarchical, modular groupings within a shared phenotypic space.

The validated clusters correspond to coherent blocks along the diagonal of the correlation matrix, indicating strong intra-cluster similarity and a structured, multiscale organization comprising both fine-grained subclusters and higher-order groupings. Complementary effect size analysis (Supplementary Tables S2, S3) further supported this structure, revealing moderate-to-large and statistically significant differences across disease systems. Hematological profiles exhibited consistent effect size deviations across multiple feature domains (Hedges’ g up to ∼0.8, q < 0.05 after correction), whereas endocrine/metabolic and cardiovascular profiles showed comparatively attenuated or opposing shifts.

The recurrence of the same features across multiple timepoints indicates temporally stable disease signatures rather than isolated single-timepoint deviations. These effects were driven by core hematologic indices, including RBC, neutrophils, and MCV, supporting a modular and biologically meaningful organization driven by stable longitudinal trajectories. Together, these complementary results indicate that clustering is driven by reproducible system-specific hematologic patterns rather than noise fluctuations, reinforcing the interpretation of disease states as structured, attractor-like regions within a shared state-space.

These results support the notion that hematopoietic clustering reflects both methodological robustness and biological plausibility, establishing a benchmark against which the more complex patterns of non-hematopoietic diseases can be interpreted using complex systems methods.

Heterogeneity of non-hematopoietic systems

In contrast to the coherence of hematopoietic disorders, non-hematopoietic disease categories exhibited marked heterogeneity across all analyses. In the phylogenetic tree and correlation-based clustering, most endocrine, metabolic, cardiovascular, gastrointestinal, respiratory, and nervous system diseases were distributed across a large composite cluster, lacking the compact substructures observed in hematopoietic profiles. Dimensionality reduction reinforced this dispersion: PCA captured global variance predominantly along a single axis without resolving system-level structure, while UMAP positioned non-hematopoietic conditions along a continuous manifold with extensive category overlap. Although enrichment analysis showed that these diseases contributed to immune and inflammatory pathways, their analyte profiles converged on broad systemic mechanisms rather than system-specific signatures. Several biological and methodological factors likely explain this pattern. Many non-hematopoietic conditions exert only indirect effects on circulating analytes, with endocrine disorders such as diabetes and thyroid disease frequently coexisting and disrupting shared metabolic pathways with highly variable manifestations (Hage et al., 2011), and cardiovascular or gastrointestinal diseases presenting secondary hematological changes (e.g., anemia of chronic disease, inflammatory responses, electrolyte shifts) that overlap extensively across systems (Weiss and Goodnough, 2005). These shared downstream phenotypes, compounded by syndromic complexity and comorbidities, dilute the specificity of blood-based signatures. Methodologically, Pearson correlation and PCA emphasize global variance structures, masking finer-grained distinctions, while UMAP reveals gradients rather than discrete partitions, and Random Forest prioritizes hematology-related features, further reflecting the lack of dominant analyte drivers in non-hematopoietic systems. This dispersion indicates that blood-analyte signatures alone are insufficient for stable categorization outside hematopoietic conditions, highlighting both the interpretive limitations of analyte-based clustering and the need for integrative multimodal frameworks to resolve cross-system biological overlaps.

Inflammation

Cluster 9, which encompassed most non-hematopoietic diseases alongside a subset of hematopoietic and “Unknown” profiles, emerged as the most heterogeneous grouping in the correlation-based clustering. Despite its size and diversity, enrichment analyses provided clear evidence that this cluster is unified by a series of shared immune and inflammatory mechanisms. KEGG pathway analysis identified highly significant enrichment for cytokine–cytokine receptor interaction, Th17 cell differentiation, AGE–RAGE signaling in diabetic complications, rheumatoid arthritis, inflammatory bowel disease, and infection-related pathways including leishmaniasis and malaria. The q-values for these pathways were extremely low (spanning approximately 4.3 × 10⁻¹⁵ to 9.4 × 10⁻²³), underscoring the robustness of the enrichment. Notably, these processes encompassed both broad immunological signaling axes and clinically recognized disease categories, indicating that Cluster 9 is not simply a residual grouping of dissimilar conditions, but one that reflects convergence on fundamental immune–inflammatory biology.

Network physiology-based visualizations further refined this interpretation. The KEGG cnetplot revealed that multiple pathways within Cluster 9 were connected through overlapping sets of hub genes, particularly pro-inflammatory cytokines and transcriptional regulators, which bridged autoimmune diseases, infection-driven processes, and metabolic complications. The emapplot corroborated these findings by highlighting functional modules: one linking leishmaniasis, tuberculosis, and toxoplasmosis, another connecting Th17 differentiation, rheumatoid arthritis, and inflammatory bowel disease, and a third coupling cytokine–cytokine receptor interaction with viral protein interaction with cytokine signaling. These modules illustrate how disparate clinical diagnoses can converge on shared gene regularity networks, explaining why such diverse entities are grouped together in correlation space.

Reactome analysis provided complementary evidence, emphasizing interleukin signaling and platelet/ECM biology as key bridging processes. Interleukin pathways (Interleukin-4/13 and Interleukin-10 signaling) clustered tightly, reflecting central cytokine cascades in autoimmune and inflammatory contexts, while platelet activation, degranulation, and calcium-mediated responses formed a parallel submodule linked to extracellular matrix organization and degradation. These findings indicate that although Cluster 9 is clinically heterogeneous, it is biologically coherent at the mechanistic level, converging on shared immune, inflammatory, and hemostatic pathways. This demonstrates that analyte-based similarity can uncover cross-cutting biological processes that transcend traditional system-based classifications and may inform refinements to disease taxonomy.

Limitations and system-level interpretation

The following subsections are retained verbatim, now grouped under a single interpretive heading to reduce structural fragmentation while preserving all original content exactly as provided. While the analyses presented here provide a coherent account of disease (dis)similarity based on blood-analyte signatures, several methodological considerations and limitations warrant careful discussion.

Another important limitation of this study is the use of synthetic disease signatures derived from curated medical knowledge rather than real patient-level data. While this approach enables controlled exploration of computational frameworks and mechanistic structure, it does not fully capture the complexity, measurement noise, and population heterogeneity inherent to patient-derived clinical datasets. Consequently, the present findings should be interpreted as demonstrating the feasibility and internal consistency of analyte-based disease modeling in silico using network physiology methods, rather than definitive biological conclusions. Future validation using real-world clinical cohorts will be essential to confirm the robustness and translational applicability of these network patterns.

Choice of similarity metric

Disease–disease similarity was quantified using Pearson correlation, which captures linear associations between analyte profiles. Although widely used and computationally efficient, this approach may underestimate relationships that are non-linear or driven by rare analytes. Alternative metrics, such as Spearman correlation (De Winter et al., 2016), cosine similarity (Xia et al., 2015), or distance measures more robust to outliers, might reveal additional structure. Moreover, converting correlation to distance through 1 − ρ assumes symmetry and equal weight across all analytes, potentially obscuring clinically meaningful heterogeneity.

Clustering and dimensionality reduction

The hierarchical clustering with UPGMA provided an interpretable dendrogram and was useful for partitioning into clusters. However, cluster membership is sensitive to the choice of cut threshold, as demonstrated by the large and heterogeneous Cluster 9. Similarly, PCA captured global variance but failed to separate systemic categories beyond hematopoietic diseases, reflecting the limitations of linear projections (Jolliffe and Cadima, 2016). UMAP recovered more nuanced local structures but is inherently stochastic and parameter-dependent (McInnes et al., 2025); alternative embeddings or repeated runs might yield slightly different layouts. These issues highlight that observed separations are robust only when consistent across multiple methods, as was the case for hematopoietic clustering.

Feature attribution

Random Forest classification identified neutrophils, red blood cell indices, haemoglobin, and platelets as key features. While these results align with biological plausibility, tree-based feature importance is influenced by collinearity among analytes and repeated measurements across timepoints. Importances therefore indicate relative contribution within the model, rather than causal effects.

The findings parallel, yet diverge from organ-based ICD taxonomies, suggesting that inflammatory and immune-network topology may provide a complementary axis of biological disease classification.

Sensitivity and robustness of feature space organization

Given that manifold and representation learning approaches in network physiology are inherently sensitive to methodological variability and noise, we performed post hoc complex systems theory-driven analyses to assess robustness and reproducibility of the observed network patterns. PCA structure remained consistent across imputation strategies, indicating that missing data handling does not bias global variance patterns. Similarly, UMAP embeddings were stable across random seeds, with high Procrustes similarity and preserved distance relationships (Supplementary Tables S5, S6), confirming that clustering is not driven by stochastic initialization but rather reflects self-organized structure.

Noise perturbation analysis further demonstrated that both the time-series and learned global feature spaces are resilient to low-to-moderate Gaussian noise, with minimal disruption to variance structure and embedding topology (Figure 9; Tables 2–4). Importantly, imputation strategy did not alter variance structure in either feature space representations (both, global and time-series). However, a key distinction emerges: the time-series feature space is highly compressible and low-dimensional, whereas the learned global feature space captures more distributed and informative variation in the clustering patterns by disease states. Complementary SVHC validation, ML classification, and AID-based feature perturbation analyses provide convergent evidence that the observed clustering pattern is reproducible and not an artifact of clustering methodology or data construction.

TABLE 4. Clustering performance across imputation strategies and distance metrics (55-feature time-series space). Clustering quality metrics computed on the full 55-feature longitudinal time-series matrix (103 disease profiles) across multiple imputation strategies (median, most frequent, k-nearest neighbors) and distance metrics (Euclidean, correlation, cosine). Performance was evaluated using silhouette score, Davies–Bouldin index, and Calinski–Harabasz index across cluster numbers (k = 2–7). Higher silhouette and Calinski–Harabasz scores and lower Davies–Bouldin values indicate improved clustering quality.

Attractor inference and state-space geometry

Analysis of the 55-feature longitudinal time-series space revealed a self-organized clustering pattern, with clustering consistently optimal at k = 2 and reduced cohesion at higher k, indicating a dominant two-state structure. PCA further showed that variance is overwhelmingly captured by a single component (PC1 ≈ 0.93), with >99% explained by the first five PCA components. These findings suggest that temporal trajectories are highly compressible and governed by constrained transitions among the disease states’ hematological features.

This supports an attractor-like interpretation in routine blood data, where disease trajectories converge onto a limited number of stable hematological states. In contrast, the learned global feature space exhibits a more distributed variance profile (PC1 ≈ 0.38), revealing a higher-dimensional manifold geometry with richer multiscale separation between disease spatial profiles. Thus, temporal dynamics are compressible and attractor-like, whereas the global (spatial) feature space encodes complex geometric structure.

Furthermore, nonlinear manifold analysis reinforces this attractor-like interpretation. UMAP stability across seeds demonstrates that both spaces preserve highly conserved relational structure, indicating that the topology of the disease state space is intrinsic and robust (Supplementary Info: Supplementary Tables S5, S6; Supplementary Figure S2). AID and ML analyses further support this interpretation, identifying consistent feature patterns that align with this manifold geometry. As shown in Figure 10 and Table 5, AID-based perturbation analysis supports an attractor-like organization of the disease state-space, where a small subset of hematologic features exerts control over global system structure. The dominance of BDM over entropy and compression metrics highlights its sensitivity to algorithmic, rather than purely statistical, complexity. The recurrence of these features across timepoints further reinforces the presence of network regulatory/control points. Within the broader context of network physiology and dynamical systems interpretation, these findings suggest that disease states mapped by blood-derived features occupy structured regions of a shared state-space, where stable, low-dimensional temporal trajectories are embedded within a more complex global geometric landscape.

TABLE 5. Top global features identified by algorithmic information dynamics (BDM and compression metrics perturbation analysis).

Future directions

The current analysis was based on a curated repository of 103 disease profiles, each derived from knowledge-based synthetic blood-analyte signatures. The disease profiles were curated as part of a synthetic digital signature repository, each containing selected analytes measured at multiple timepoints. While this dataset enabled proof-of-concept analyses, the size and diversity remain limited. Future studies should incorporate clinically derived data beyond synthetic, larger and more heterogeneous cohorts ideally spanning multiple populations, age groups, and clinical settings. Such expansion would increase statistical power, reduce biases, and enable stratified analyses that account for demographic and comorbidity effects. Linking profiles to longitudinal data would further allow investigation of disease progression and treatment response.

Blood-analyte profiles capture only one dimension of disease biology. Integrative frameworks that combine hematology with symptoms, and imaging should be explored. Multi-omic integration would help disentangle heterogeneous systems such as endocrine or gastrointestinal diseases, where blood indices provide only indirect readouts. Incorporating electronic health records and patient-reported outcomes could bridge molecular signatures with clinical trajectories, enhancing translational relevance.

Ultimately, future work should explore how analyte-based clustering can contribute to a hybrid taxonomy that integrates ICD-based classification with cross-system mechanistic axes. This may enable more biologically meaningful groupings, support early disease interception, and guide mechanism-based treatment allocation. Moving from proof-of-concept toward clinical-grade classification will require iterative refinement between computational analyses, mechanistic interpretation, and translational implementation. Future studies should also incorporate causal-mechanistic frameworks capable of revealing how minimal perturbations to blood-analyte signatures propagate through disease clusters and shared inflammatory axes. Approaches grounded in algorithmic information dynamics (AID) offer such a way to infer generative structure, identify causal drivers of cluster formation, and model counterfactual intervention strategies directly from digital blood twins (Zenil et al., 2019; Uthamacumaran and Zenil, 2022; Zenil et al., 2023; Abrahão and Zenil, 2022). Our findings support the AID framework within network physiology, where perturbation analysis reveals a small subset of structurally influential features that disproportionately shape the global information structure, consistent with control-like behavior in an attractor landscape. Integrating these causal tools with larger, real-world datasets may help refine mechanistic disease taxonomy, expose latent confounding pathways, and enhance the interpretability of analyte-based disease stratification.

Conclusion

This study demonstrates that blood-derived digital signatures can recover highly coherent disease clusters for hematopoietic conditions while simultaneously revealing a large, mechanistically convergent cluster of systemic disorders driven by immunological and inflammatory pathways. Rather than replicating ICD-style taxonomies, hematology-based similarity exposes confounding factors and a deeper mechanistic substrate, one in which cytokine signaling, interleukin networks, platelet activation, and extracellular matrix remodeling form shared axes that unify clinically disparate diseases.

In translational and systems medicine, this suggests that routine laboratory data, already embedded in healthcare workflows, could serve as the basis for a complementary mechanistic disease ontology. Such a hybrid taxonomy, grounded partly in analyte biology and partly in clinical categorization, would support earlier risk stratification, more granular comorbidity mapping, and biologically informed patient grouping. Importantly, this network physiology framework does not require new tests, technologies, or invasive sampling; it reuses existing clinical data to surface latent biological structure.

By linking hematologic coherence with immune–inflammatory convergence across complex systems, this work outlines a pathway toward precision medicine that is computationally accessible and clinically scalable. As future validation studies expand the dataset, integrate multi-modal signals, and evaluate prognostic value, analyte-based digital twins may evolve into real-time risk stratification tools capable of informing continuous disease monitoring, mechanistic (causal) discovery, and precision treatment allocation in patient-centered care.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2026.1747187/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Plamen Ch. Ivanov, Boston University, United States

Reviewed by: Salvatore Micciche’, University of Palermo, Italy

Alexey Zaikin, University College London, United Kingdom