This study used clinical diagnosis data from the Medical Information Mart for Intensive Care IV (MIMIC-IV; version 2.2) [25-27]. MIMIC-IV is a comprehensive, publicly available electronic health record (EHR) database developed by the Massachusetts Institute of Technology (MIT) and Beth Israel Deaconess Medical Center (BIDMC), encompassing inpatient medical records of patients with critical illness admitted to BIDMC between 2008 and 2019. Using individual hospital admissions as the primary unit of analysis, we collected hospitalization records for patients aged 20 years and older, organized chronologically by admission time. To construct a structured dataset capturing both cross-sectional and longitudinal features, historical hospitalization information from the preceding 3-year window was integrated for each index admission. Consequently, each record was treated as an independent sample, comprising a sequence of ICD codes from hospitalizations within the prior 3 years, and the corresponding study label.

The final cohort consisted of 71,212 patients and 226,992 hospital admissions, covering 1663 unique International Classification of Diseases, 10th Revision, Clinical Modification (ICD-10-CM) diagnosis codes. Each admission contained a mean of 12.36 (SD 6.99) diagnosis codes, providing a heterogeneous representation of patients' comorbidity profiles. Baseline demographic and clinical characteristics were summarized in Table 1, including age, sex, length of stay, mortality, and admission statistics. To characterize the availability of longitudinal information for each index admission, we summarized the distribution of historical hospital visits within the 3-year observation window. Each index admission was associated with a mean of 4.82 (SD 8.14) historical hospital visits (median 2, IQR 1-5). Among all index admissions, 33.8% had only one historical visit available, whereas 66.2%, 28.3%, and 12.1% had at least 2, 5, and 10 historical visits, respectively. These statistics demonstrate that the cohort included both sparse and dense longitudinal trajectories, enabling evaluation of the proposed framework across varying amounts of historical information.

Since the data collection period of MIMIC-IV coincided with a transition in the US disease coding system, both the ICD, Ninth Revision, Clinical Modification (ICD-9-CM) and ICD-10-CM coding systems coexist within the dataset. To standardize the feature format, the primary preprocessing step involved integrating these 2 coding systems. We downloaded the General Equivalence Mappings (GEMs) from the official website of the Centers for Medicare & Medicaid Services (CMS) [28] to convert all mappable ICD-9-CM codes into their corresponding ICD-10-CM codes. Specifically, the mapping logic was applied as follows:

Finally, to balance feature granularity and interpretability, all resulting ICD-10-CM codes were truncated to their first 3 characters. Furthermore, we excluded Z codes, which represent health care encounters unrelated to active diseases or injuries. Unless specified otherwise, all subsequent mentions of ICD codes in this study refer to the ICD-10-CM system.

Target-specific analyses were conducted for HF (ICD-10-CM code I50). For outcome labeling, hospitalizations were categorized as positive (label=1) if the target ICD code was present in either the primary or secondary hospital discharge diagnosis list; otherwise, they were labeled as negative (label=0). These discharge diagnoses were assigned by certified professional medical coders based on comprehensive clinical documentation using standardized ICD-10-CM coding practices. Although administrative diagnosis codes may be subject to coding inaccuracies, ICD-based outcome definitions have been widely adopted in large-scale EHR studies, thereby providing a standardized and reproducible approach for outcome ascertainment.

To evaluate the representational capacity of the model and prevent information leakage, a robust feature masking protocol was implemented. Specifically, the target code (I50) was removed from all diagnosis sequences prior to feature construction. Furthermore, to avoid reliance on direct clinical manifestations or established etiologies, we excluded highly correlated codes, including hypertensive heart disease (I11, I13), cardiomyopathies (I42, I43), complications and ill-defined descriptions of heart disease (I51), and pulmonary edema (J81). This approach encourages the model to learn latent risk representations from the underlying comorbidity structure and longitudinal disease history. The masked diagnosis sequences were used for all predictive feature construction, model training, and performance evaluation, including the risk stratification and architectural ablation analyses. In contrast, unmasked diagnostic records were used solely for post hoc interpretation and visualization, including the IG and positive pointwise mutual information (PPMI) comorbidity analysis and risk-aware pseudotime visualization. These unmasked records were analytically separated from the predictive modeling pipeline and not used for feature construction, model training, or performance evaluation.

Only historical ICD records preceding the index admission were used as input features, whereas diagnosis codes from the index admission were reserved strictly for outcome definition. This design ensured that risk stratification relied solely on historical medical information available prior to the index admission while reflecting a clinically realistic setting. Patients with at least one prior hospitalization within the 3-year look-back period were included in the trajectory construction. Patients without any historical hospitalizations during this period lacked a longitudinal trajectory and were therefore excluded from the trajectory-based classification analysis. For patients with a single admission, trajectory-derived features were set to zero, while pooled representations remained available. Consequently, these admissions were represented by a single-state embedding rather than being excluded from the analysis.

This study was a secondary analysis of a publicly available, deidentified MIMIC-IV database. The creation and maintenance of the MIMIC-IV database were approved by the Institutional Review Boards (IRBs) at MIT and BIDMC. Access to the database was granted to the author (JLH) after successful completion of the “Data or Specimens Only Research” training (CITI Program, Record ID: 61470359) and the execution of a formal data use agreement. All patient records contained in the MIMIC-IV database were rigorously deidentified prior to release in accordance with the Health Insurance Portability and Accountability Act (HIPAA) Safe Harbor regulations. Consequently, the requirement for individual patient informed consent was waived by the IRBs. There was no direct contact with human subjects, and participants received no compensation. Furthermore, no identifying information or clinical images of individual participants were included in this manuscript or its supplementary materials.

In machine learning, information gain (IG) is a fundamental metric for quantifying a feature’s contribution to class discrimination and is based on information-theoretic entropy [29]. In this study, IG was employed to assess the discriminative power of individual ICD codes for target disease identification, measuring the reduction in uncertainty regarding the outcome given the presence or absence of a specific code. For each ICD code in the training dataset, the IG value was calculated using equation 1:

where H(Y) denotes the marginal entropy of the target outcome and H(Y|ICD) represents the conditional entropy. This metric captures the extent to which each ICD code facilitates effective classification, where higher IG values signify superior discriminative capacity. Since the target outcome was HF (ICD code I50), the resulting IG scores specifically quantified the contribution of each ICD code to the identification of HF. These disease-specific IG scores were subsequently used in a weighted pooling mechanism to emphasize ICD diagnosis codes with higher clinical relevance for identifying the target disease.

To capture the underlying comorbidity structures within hospitalization records, we constructed an ICD code co-occurrence network from the training data using pointwise mutual information (PMI) [30]. PMI is widely used in machine learning and natural language processing (NLP) to quantify co-occurrence strength between discrete entities. In our framework, each hospitalization was treated as a “document” comprising all assigned ICD codes. A co-occurrence event was defined by the concurrent appearance of 2 ICD codes within the same admission. Based on the marginal probabilities of individual ICD codes and their joint probabilities estimated from the training dataset, the PMI between a pair of ICD codes (i, j) was defined as equation 2:

where P(i) denotes the occurrence probability of code i, and P(i,j) represents their joint probability. To reduce noise from negative associations, which may lack clinical interpretability and be mathematically unstable, we excluded negative PMI values. Instead, we adopted PPMI, as defined in equation 3, to retain only diagnostically meaningful positive associations between ICD codes. Finally, ICD codes were modeled as nodes, with PPMI values serving as edge weights, thereby constructing a weighted, undirected ICD co-occurrence graph.

PPMI(i,j)=max (PMI(i,j),0)

Following the construction of the ICD co-occurrence graph, the Node2Vec algorithm was implemented to learn complex relational patterns among ICD nodes [16]. Node2Vec employs second-order random walks that interpolate between breadth-first search (BFS) and depth-first search (DFS) through the return parameter (p) and in-out parameter (q). In this study, the default values (p=q=1) were adopted, providing an unbiased exploration of the disease network. Accordingly, no explicit bias toward either BFS- or DFS-like exploration was imposed during random walk generation. Through this approach, the high-dimensional comorbidity graph was mapped into a low-dimensional continuous vector space. This ensured that ICD codes with strong co-occurrence patterns, or those sharing similar structural roles within the network despite lacking direct connections, were positioned in proximity within the embedding space. Node2Vec was implemented using an embedding dimension of 64, a walk length of 10, and 30 random walks per node. The complete hyperparameter configuration was summarized in Multimedia Appendix 1.

For each index admission, all geometric features derived from the corresponding longitudinal patient history were concatenated into a 389-dimensional feature vector that served as the input to the downstream classifiers. This architecture preserves both the longitudinal disease context and the immediate diagnostic signals of the index admission, thereby enhancing the model’s ability to identify clinical risk states. To prevent data leakage and overestimation of model performance resulting from recurrent patient admissions, we implemented a patient-stratified data splitting strategy during 10-fold cross-validation. Under this framework, all admissions associated with a single patient were assigned to the same fold, ensuring complete independence between training and testing sets. We employed an extreme gradient boosting (XGBoost) classifier for risk identification and used logistic regression (LR) as a baseline to compare XGBoost with a conventional linear approach. The detailed hyperparameter settings for the XGBoost classifier and the LR baseline were summarized in Multimedia Appendix 2. These hyperparameters were determined through preliminary experiments and remained fixed throughout all cross-validation folds, except for the fold-specific scale_pos_weight parameter, which was recalculated according to the class distribution within each training fold. To mitigate class imbalance (18.83% positive prevalence), a cost-sensitive learning strategy was adopted instead of resampling, thereby preserving the original class distribution during model training. Specifically, for the XGBoost model, the scale_pos_weight parameter was determined dynamically for each training fold according to the ratio of negative to positive samples in the corresponding training data. For the LR model, the class_weight parameter was set to balanced. This approach penalized the misclassification of the minority class during loss optimization. Finally, model performance was evaluated using the area under the curve (AUC) and the area under the precision-recall curve (AUPRC) [31]. In addition, 95% CIs for the AUC and AUPRC were estimated using patient-level bootstrap resampling with replacement on the independent test set over 1000 bootstrap iterations. The resulting empirical bootstrap distributions were used to derive the corresponding 95% CIs.

To benchmark the proposed maximum IG-weighted deviation geometric (mIGwDG) framework, we implemented 3 baseline approaches under the same training, validation, and evaluation protocol. First, to model longitudinal disease trajectories, we implemented a standard long short-term memory (LSTM) network using the same Node2Vec-based admission embeddings as sequential inputs. LSTM is a widely adopted deep learning architecture for modeling longitudinal EHR data and disease progression [32]. Second, we implemented LR and XGBoost using sparse ICD one-hot encoded features representing all historical diagnosis codes. These models are widely used baselines for structured EHR-based risk prediction and allow evaluation of the proposed representation without feature compression [33]. Finally, sparse ICD one-hot encoded features were concatenated with the proposed mIGwDG representation and used as input to both LR and XGBoost to evaluate whether the learned trajectory representation provided complementary predictive information beyond the original sparse ICD one-hot encoded features by incorporating longitudinal and geometric disease trajectory information.

This study presented a risk-stratified visualization framework to characterize disease progression patterns based on model-inferred risk estimates. Throughout this study, the term “pseudotime” denoted a risk-ordered representation of index admissions derived from calibrated model-predicted probabilities, rather than an optimization-based trajectory inference method, as commonly used in single-cell transcriptomics. Out-of-fold (OOF) risk probabilities were obtained through 10-fold cross-validation to ensure unbiased risk estimation and mitigate data leakage. Because the raw probability outputs of tree-based models such as XGBoost were not necessarily well calibrated, probability calibration was incorporated into the predictive pipeline before risk stratification. Specifically, Platt scaling was applied within each cross-validation iteration by fitting an LR model to the raw OOF probabilities from the training folds. This fold-wise calibration strategy ensured that probability calibration remained independent of the validation data, thereby preventing information leakage. Calibration performance was evaluated using reliability diagrams (calibration curves) across 10 decile bins, as shown in Multimedia Appendix 3.

The calibrated probabilities were subsequently used for pseudotime risk stratification. To avoid reliance on arbitrary probability thresholds, we employed a quantile-based binning to divide index admissions into 5 risk strata (Very low, Low, Medium, High, and Very high risk) according to the calibrated probabilities. Formally, the risk stage of the i-th index admission was defined using equation 17,

Where

denotes the calibrated OOF probability, and Q(⋅) represents the quintile mapping function that partitions the cohort into 5 equal-sized risk strata.

To characterize the accumulation of comorbidities across the inferred risk spectrum, we computed the prevalence of each ICD code within every risk stratum. For ICD code j in risk stage k, the prevalence was calculated using equation 18,

where n_jk is the number of index admissions in risk stage k containing ICD code j, with each ICD counted at most once per admission, and N_k is the total number of index admissions assigned to that risk stage.

To reveal the sequential emergence of comorbidities along the risk spectrum, ICD codes were ordered according to their weighted average risk stage (evolution score) and shown in equation 19,

Where w_k ∈ {0,1,2,3,4} denotes the ordinal index of the 5 risk strata from Very low to Very high. ICD codes were subsequently sorted in ascending order of S_j, thereby arranging diseases according to their predominant stage of occurrence.

We prioritized the top 25-IG×PPMI-weighted comorbid ICD codes for visualization. Finally, this prevalence matrix was visualized as a heatmap highlighting the gradual and structured accumulation of clinically relevant comorbidities along the latent risk spectrum. The resulting heatmap organized comorbidities along a pseudo-temporal risk axis that reflects increasing model-inferred risk rather than chronological time. This framework was intended solely as a descriptive, hypothesis-generating representation of disease evolution based on model-inferred risk and does not imply chronological progression, direct temporal causality, or validated clinical disease staging.

We evaluated our framework on 226,992 hospital admissions with HF (ICD-10: I50) as the target condition. Notably, the positive class (Label=1) accounted for 18.83% of the data, representing a significant class imbalance. To ensure the integrity of our findings and prevent data leakage, we implemented a patient-stratified 10-fold cross-validation strategy. This approach guarantees that all encounters belonging to a specific patient remain in either the training or the test set, never both. Furthermore, within each fold, both the disease co-occurrence graphs and the IG-based importance weights were derived exclusively from the training partition. This protocol ensures that the learned representations are validated on entirely independent data, preserving the methodological integrity of our risk stratification evaluation.

We first evaluated admission-level risk stratification performance using 10-fold cross-validation across different ICD embedding aggregation strategies. As shown in Table 2, XGBoost consistently outperformed LR in both AUC and AUPRC. For conventional static baselines, XGBoost achieved AUPRCs of 0.627 (SD 0.011; mean pooling) and 0.625 (SD 0.017; maximum pooling), compared with a mean of 0.538 (SD 0.013) and a mean of 0.527 (SD 0.029) for LR, respectively. In addition to the proposed framework and its internal ablation variants, 3 commonly used baseline approaches (LR, XGBoost with ICD one-hot encoding, and LSTM) were evaluated under the same experimental protocol.

To further quantify the contribution of each feature module, we conducted an additional feature-module ablation study by progressively incorporating risk deviation, intervisit trajectory transition, and geometric variance/dispersion features (Multimedia Appendix 4). Starting from the risk deviation representation, trajectory change and geometric features resulted in only marginal differences in predictive performance, with changes generally within the range of cross-validation variability. In contrast, IG weighting produced a more pronounced improvement compared with maximum-based aggregation. The complete mIGwDG representation achieved the best overall performance, although the incremental contributions of the individual trajectory-change and geometric modules were not significantly obtained.

The complete mIGwDG framework achieved an AUC of 0.888 (SD 0.003) and an AUPRC of 0.667 (SD 0.016) under XGBoost in 10-fold cross-validation while maintaining comparable predictive performance to sparse ICD one-hot encoding with substantially lower feature dimensionality. Under LR, mIGwDG achieved an AUC of 0.881 (SD 0.004) and an AUPRC of 0.647 (SD 0.013). Importantly, the proposed hierarchical representation substantially narrowed the performance gap between the linear and nonlinear classifiers, suggesting that it transforms longitudinal disease trajectories into a more discriminative feature space that can be effectively exploited even by a linear classifier. Figure 1 presents precision-recall curves across the evaluated models. Solid lines indicate the mean precision at each recall level across the 10 folds, and the shaded area represents the standard deviation of the proposed XGBoost-mIGwDG model. Mean AUPRC values and cross-validation variability are reported in the figure legend.

To ensure an unbiased evaluation, we generated OOF risk predictions using 10-fold cross-validation and calculated AUPRC based on aggregated OOF predictions. The proposed mIGwDG representation combined with XGBoost achieved an aggregated OOF AUPRC of 0.677 (95% CI 0.657-0.677). The ICD one-hot XGBoost baseline achieved a comparable cross-validated AUPRC of 0.666 (SD 0.016). Although the one-hot representation retained the full ICD feature space, it resulted in a substantially higher-dimensional feature representation. In contrast, mIGwDG achieved similar predictive performance using a compact trajectory-based representation.

To evaluate patient-level risk stratification performance, we aggregated admission-level OOF predictions into patient-level risk scores using a max-pooling strategy across each patient’s longitudinal sequence. This resulted in 71,212 patient-level risk profiles derived from 226,992 admissions (19.52% positive prevalence).

As shown in Table 3, patient-level aggregation improved performance across all evaluated models compared with admission-level evaluation. The proposed XGBoost-mIGwDG model achieved the best performance, with an OOF-AUC of 0.902 (95% CI 0.899-0.905) and an OOF-AUPRC of 0.730 (95% CI 0.722-0.737). The LSTM baseline achieved an OOF-AUPRC of 0.702 (95% CI 0.694-0.710), whereas LR using the full feature set achieved an OOF-AUPRC of 0.705 (95% CI 0.697-0.713).

To evaluate performance within the most clinically vulnerable population, we performed a subgroup analysis using predefined clinical criteria that were independent of model predictions. The Charlson Comorbidity Index (CCI) was calculated for each index admission using all diagnosis codes recorded during the preceding 3-year observation window, thereby reflecting the cumulative burden of comorbidities prior to the index admission. Patients aged ≥75 years and with a CCI ≥4 were classified as the high-risk subgroup, isolating 24,574 high-acuity admissions (10.83%) and 7871 unique patients (11.05%). The predictive performance of all evaluated models within this subgroup is summarized in Table 4.

We further evaluated model performance within the high-risk subgroup at both the admission level (event rate: 54.98%) and the patient level (event rate: 61.34%). At the admission level, XGBoost models incorporating the proposed trajectory representations consistently outperformed conventional aggregation strategies. The XGBoost-mIGwDG model achieved the best performance, with an AUC of 0.776 and an AUPRC of 0.796, compared with an AUPRC of 0.787 for the LSTM baseline. LR also showed consistent performance improvements after incorporating the proposed trajectory representation. After patient-level aggregation, performance further improved across all evaluated models. The XGBoost-mIGwDG model achieved an AUPRC of 0.884, comparable to the LSTM baseline (AUPRC=0.878). LR also improved, with AUPRC increasing from 0.815 using mean aggregation to 0.873 with the proposed trajectory representation.

These results demonstrate consistent performance of the proposed representation across admission-level and patient-level high-risk stratification settings. Because the high-risk subgroup has a higher baseline event rate than the overall cohort, these performance metrics should be interpreted within the subgroup rather than directly compared with those of the full population.

To evaluate the robustness of the proposed disease trajectory representation to incomplete diagnostic information, we conducted an additional analysis by retaining admissions that lacked disease-specific diagnostic information after preprocessing. This expanded cohort comprised 71,323 patient-level risk profiles across 227,328 admissions, with positive prevalence rates of 18.80% and 19.49%, respectively. In the primary analysis, these admissions were excluded because they did not contribute meaningful disease signals to the trajectory representation. In the robustness analysis, these admissions were retained to represent increased structural uncertainty arising from incomplete diagnostic information, a common characteristic of real-world EHR data. Compared with the primary analysis, predictive performance showed only marginal reductions, suggesting that the learned trajectory representation is relatively robust to moderate levels of incomplete diagnostic information. The detailed predictive performance of this robustness analysis is presented in Table 5.

To interpret model risk identification specifically for HF (ICD: I50), we constructed an IG-weighted comorbidity subgraph that integrates statistical co-occurrence structure with IG-based feature relevance. For this post hoc interpretability analysis, the PPMI-based co-occurrence structure and IG values were calculated separately using unmasked diagnostic records. The IG values used to determine node weights in this visualization were derived specifically for the interpretability analysis, and they were not used for predictive feature construction, model training, or performance evaluation. In contrast, the predictive analyses were conducted using the masked diagnosis sequences described above. To synthesize structural relevance with feature relevance, node weights were assigned based on the separately calculated IG values, while edge weights represented the co-occurrence strength measured by PPMI. The resulting IG×PPMI comorbidity subgraph highlights diagnostic codes that were strongly represented in the observed comorbidity structure and informative in the IG-based interpretability analysis. It also provided an interpretable lens into the local comorbidity architecture surrounding the target diagnosis. Notably, this visualization serves as a post hoc interpretability tool, designed to elucidate informative comorbidity configurations rather than to infer direct causal pathways among diseases.

To investigate how diagnostic roles influence comorbidity architectures, we compared I50-centered subgraphs across different contexts (Figure 2). In general diagnostic contexts (Figure 2A), the neighborhood is characterized by a broad spectrum of systemic conditions, including hypertension-related disorders (I11, I13) and metabolic syndromes (E11, E78, and M10), alongside acute critical events such as myocardial infarction (I21). These findings suggest that, in broader clinical contexts, I50 often manifests as a comorbid component within a milieu of chronic systemic burden and acute dysregulation. Conversely, when I50 is designated as the primary diagnosis (Figure 2B), the comorbidity architecture becomes more focused on intrinsic structural and functional pathologies. Notably, conditions such as renal anemia (D63), conduction disorders (I44), and venous congestion (I87) emerge as top-20 components exclusively in this primary diagnostic role. This divergence demonstrates that the IG×PPMI framework effectively captures nuanced contextual semantics, identifying distinct clinical profiles that transcend simple co-occurrence frequencies.

Furthermore, the IG×PPMI framework revealed 2 prominent structural modules centered on E11 (type 2 diabetes) and N18 (CKD), respectively (Figure 3). The E11-centered module was primarily associated with cardiovascular and metabolic diseases, including dyslipidemia (E78) and chronic ischemic heart disease (I25), whereas the N18-centered module was closely associated with cardiorenal complications, including hypertensive heart and CKD (I13), electrolyte imbalance (E87), and hypertensive nephropathy (I12). These findings reveal 2 distinct comorbidity modules surrounding HF, providing a structural basis for subsequent interpretation of heterogeneous disease progression patterns.

To assess the clinical face validity of our framework, we visualized the disease histories of 3 low-risk and 3 high-risk patients, selected based on their calibrated OOF risk probabilities, using principal component analysis (PCA; Figure 4A). To provide a common visualization space, the PCA basis was fitted on the mean-pooled admission embeddings derived from all admissions in the cohort, and the visit-level maximum information gain-weighted (mIGw) embeddings were subsequently projected onto this shared coordinate system. Here, mIGw denotes the visit-level information-gain-weighted embedding (Equation 6), obtained by applying the IG weights to the embedded ICD codes within each hospitalization. These visit-level representations are subsequently aggregated along each patient’s longitudinal hospitalization history to construct the mIGwDG path-level representation used for model training and prediction. The first 2 principal components accounted for 41.3% of the total variance (25.2% and 16.1%, respectively). Each point represents an individual admission embedding, and consecutive admissions were connected chronologically with directed arrows to form a patient trajectory, illustrating the evolution of the patient's clinical state over time.

As illustrated in Figure 4A, representative high-risk (calibrated OOF>60%) and low-risk (calibrated OOF<0.1%) patients occupy visually distinct regions in the learned disease geometry. The representative high-risk trajectories remain localized near the learned high-risk centroid (marked by the black star), whereas the representative low-risk trajectories consistently remain distal to the high-risk centroid. Although this visualization is illustrative rather than quantitative, it qualitatively shows that patient trajectories are organized into a coherent geometric structure within the learned representation. This alignment between geometric movement and calibrated model-estimated risk probabilities provides evidence for the interpretability of our latent representation.

To qualitatively examine how the proposed representation influences the learned disease geometry, representative patient trajectories and the corresponding high-risk centroid were projected into a 2D space using PCA. Figure 4B compares the trajectory representations generated by the conventional mean pooling strategy and the proposed mIGw representation using the same representative for patients with high-risk and low-risk. Both representations were projected onto the same PCA coordinate system to enable direct visual comparison. In both representations, the representative high-risk trajectory remained in the vicinity of the learned high-risk centroid, whereas the representative low-risk trajectory occupied a distinct region of the embedding space.

Compared with the conventional mean pooling strategy, the proposed mIGw representation produced more compact trajectory distributions with clearer separation between the representatives of patients with high-risk and low-risk. The proposed representation also appeared to produce more localized high-risk trajectories, with the representative patient with high-risk remaining closer to the learned high-risk centroid while exhibiting reduced geometric variability around the centroid. Although these visualizations are qualitative rather than statistically conclusive, they provide a geometric perspective that complements the quantitative predictive performance reported for mIGwDG. These qualitative observations suggest that the proposed representation better preserves the local geometric organization of patients with similar disease risk, resulting in more coherent trajectory organization within the embedding space.

In clinical practice, disease progression is highly heterogeneous; some patients remain stable over a decade of follow-up, while others experience rapid functional decline within only 2 years. Traditional temporal analyses often struggle to align these disparate trajectories when relying solely on chronological time or visit counts. To overcome this limitation, we visualized the top 25 comorbid ICD codes identified through the IG×PPMI framework using unmasked diagnostic records along a model-derived, risk-aware pseudotime axis derived from the masked predictive representation (Figure 5). Individual patient trajectories and their corresponding longitudinal comorbidity patterns were visualized using unmasked diagnostic records, whereas patients were ordered along the pseudotime axis according to mIGwDG risk scores rather than absolute time, thereby treating risk magnitude as the primary metric of progression. This ensures that comorbidity trajectories reflect the pathological evolution of the disease rather than mere temporal accumulation, providing a standardized scale to observe how comorbidity profiles transition across 5 risk strata (from Very low to Very high).

The analysis reveals a clear monotonic accrual of chronic systemic conditions rather than stochastic acute episodes. Foundational metabolic disorders, such as hyperlipidemia (E78) and type 2 diabetes (E11), are markedly elevated even in low-risk stages, with E78 reaching a prevalence of 0.78 at the very high-risk stage. As patients progress, a cardiorenal-metabolic signature emerges with increasing synchrony; CKD (N18), hypertensive kidney disease (I12), and electrolyte disorders (E87) exhibit a pronounced gradient. This trajectory mirrors the pathophysiological shift from localized injury to systemic homeostatic failure.

Distinct temporal patterns were observed across disease categories. Cardiovascular comorbidities, including cardiomyopathy (I42) and atrial fibrillation (I48), show concentrated gains in the moderate-to-high risk range. In contrast, acute conditions like shock (R57) and myocardial infarction (I21) remain infrequent until the very high-risk stage, suggesting they serve as indicators of terminal decompensation rather than early risk markers.

Collectively, the model-inferred risk stratification identifies a progressive accrual of chronic systemic comorbidities rather than stochastic acute episodes. This suggests that the model effectively prioritizes structural disease evolution. Consequently, as patients transition into high-risk strata, our findings advocate for a clinical shift toward integrated multiorgan management, addressing the complex interplay between cardiac, renal, and pulmonary systems.

The y-axis represents ICD codes, and the x-axis represents 5 risk strata along the model-inferred pseudotime axis. Color intensity indicates the normalized prevalence of each ICD code within each risk stratum. The continuous pseudotime was divided into 5 ordinal risk strata (Very low, Low, Medium, High, and Very high) using quantile-based binning based on the empirical 20th, 40th, 60th, and 80th percentiles of the calibrated pseudotime distribution. This pseudotime representation is intended to summarize model-inferred risk progression and should not be interpreted as a direct temporal or causal ordering of diseases. Instead, it illustrates how clinically relevant comorbidities vary across increasing levels of model-estimated risk.

In this study, we developed a longitudinal ICD trajectory representation that incorporates geometric structural information for high-risk patient identification from EHR data. By integrating geometric features, including historical means and maxima, Euclidean distances, directional shifts, and cosine dispersion, with IG-based weighting, the proposed approach summarizes disease histories in a structurally informed manner. Rather than relying solely on conventional pooling or sequential modeling, the proposed approach is designed to preserve structural characteristics of disease trajectories that may reflect progression toward higher model-estimated risk. In addition, we introduced a risk-aware pseudotime framework that organizes heterogeneous patient trajectories along a latent risk continuum inferred from model predictions. Rather than relying on chronological alignment, this approach enables characterization of disease trajectory patterns according to increasing risk profiles, providing a complementary perspective for longitudinal risk stratification.

The proposed framework consistently outperformed conventional embedding aggregation strategies while achieving predictive performance comparable to the sparse ICD one-hot representation despite reducing the feature space by approximately 75%, from approximately 1600 ICD variables to a 389-dimensional trajectory representation. Beyond prediction, the learned representation facilitates visualization of disease trajectories, risk-aware pseudotime analysis, and geometric characterization of disease evolution, providing additional interpretability through visualization and trajectory characterization that high-dimensional sparse one-hot features alone do not readily provide. Furthermore, combining mIGwDG with the original ICD one-hot representation consistently improved the performance of both LR and XGBoost, suggesting that the proposed representation captures complementary information beyond conventional sparse ICD encoding.

While deep learning architectures (such as RNNs and Transformers) offer substantial flexibility for temporal modeling, their performance in clinical risk prediction may depend on large datasets, adequate representation of positive cases, and extensive hyperparameter tuning. In the context of HF, hospitalization data are characterized by high heterogeneity, sparse events, and pronounced class imbalance, conditions under which deep models may struggle with stability and generalizability.

Recent GNN frameworks have demonstrated excellent predictive performance by jointly modeling longitudinal EHRs and comorbidity connectivity through heterogeneous graph attention mechanisms [18]. These approaches represent an important direction in medical representation learning but generally require more sophisticated model architectures and greater computational resources. As a complementary alternative, the proposed framework employs Node2Vec embeddings coupled with geometric trajectory descriptors, providing a compact yet explicit representation of longitudinal comorbidity patterns. This design ensures robust learning and consistent performance even under moderate sample sizes. A key advantage of this approach is its intrinsic interpretability; it decomposes historical hospitalizations into IG-weighted representations and geometric measures (displacement and dispersion), and each component can be independently linked to clinically meaningful patterns.

Furthermore, the computational efficiency of the Node2Vec and XGBoost pipeline facilitates deployment in standard information systems without specialized hardware. Notably, the proposed representation remains highly effective even when coupled with a simple LR classifier. This suggests that the proposed representation reshapes longitudinal clinical representations into a more separable risk space, reducing the reliance on high-capacity, noninterpretable nonlinear models. Recent studies have likewise demonstrated that constructing informative representations from high-dimensional clinical data prior to supervised learning can substantially improve predictive performance [34]. Although our framework does not employ explicit feature selection through statistical correlation analysis or optimization-based search algorithms, it follows a similar methodological motivation by transforming sparse longitudinal disease histories into compact graph-based geometric representations before classification. Consequently, predictive performance is enhanced through more informative representations rather than by relying solely on more complex predictive models.

Traditional pooling-based embedding approaches, such as mean pooling or max pooling, summarize diagnostic information either by representing an overall central tendency of a patient’s disease space or by emphasizing a small number of extreme diagnoses that emphasize the most prominent diagnostic features. While effective for dimensionality reduction, these methods fail to distinguish between stable disease burden and abrupt pathological deviations along the disease course.

In contrast, our framework is designed to jointly capture chronic burden and acute shifts. By employing geometry-based features, we quantify the spatial relationships between diagnostic combinations and their proximity to a high-risk region within the embedding space. The average distance was designed to quantify diagnostic dispersion, indicating whether a patient’s profile is concentrated on a single clinical axis or spans multiple comorbidity modules. The variability of these distances characterizes systemic complexity; higher dispersion may reflect greater heterogeneity across comorbidity patterns.

A distinctive innovation is the incorporation of deviation vectors, providing a directional geometric perspective on comorbidity evolution. We observed that high-risk trajectories deviate from low-risk regions, gravitating toward the high-risk centroid. This displacement may reflect systematic changes in disease trajectory organization rather than the influence of isolated diagnostic events. Consequently, this trajectory-level representation provides additional information for risk stratification, enabling the model to characterize the internal architecture of diagnostic combinations. By capturing this structural information, our framework enhances risk discrimination without necessitating the computational overhead of complex deep learning architectures. The primary contribution of these geometric features is to characterize changes in disease trajectories rather than individual diagnostic states.

The IG×PPMI comorbidity analysis provides an additional perspective for interpreting the learned representations. By integrating model-derived attribution (IG) with population-level co-occurrence (PPMI), the framework highlights clinically relevant comorbidity patterns that may not be apparent from either source alone. An advantage of this approach is its ability to identify comorbidities that receive relatively high IG attribution despite only moderate population co-occurrence. These findings suggest that the proposed representation may capture diagnostic patterns that contribute disproportionately to model-estimated risk beyond highly prevalent comorbidities such as hypertension. The combined IG×PPMI framework therefore complements conventional feature importance analyses by jointly considering model-derived predictive attribution and population-level disease co-occurrence.

In the top 20 comorbidities in the HF-centered subgraph (I50), anemia (D63), and electrolyte disorders (E87) exhibit high IG weights, often emerging before overt symptoms like orthopnea or marked edema. Anemia, with a reported prevalence of 35%-70% in patients with HF [35,36], is a well-established marker of disease severity and is strongly associated with increased hospitalization and mortality [37]. Similarly, electrolyte disturbances (E87), including hyponatremia and potassium imbalances, reflect both intrinsic neurohormonal dysregulation and the impact of pharmacological interventions. These imbalances serve as critical indicators of pathophysiological shifts and are robustly linked to adverse outcomes [38,39].

Furthermore, among the identified comorbidity modules, the E11- and N18-centered modules are of particular clinical interest because they correspond to 2 well-established systemic syndromes in HF. The E11-centered module is consistent with the well-established cardiometabolic continuum, whereas the N18-centered module is consistent with a common cardiorenal syndrome in patients with advanced HF. These findings indicate that the proposed IG×PPMI framework could identify clinically meaningful comorbidity structures beyond what can be inferred from diagnostic frequency alone. By integrating feature importance with disease co-occurrence, the framework may reveal potential early warning signals embedded within routine EHRs, thereby highlighting clinically relevant comorbidity patterns associated with HF and facilitating clinical risk stratification.

A second contribution of this study is the introduction of the risk-aware pseudotime framework, which provides a risk-oriented organization of longitudinal disease trajectories. Longitudinal EHR analyses are often challenged by temporal heterogeneity and irregular sampling, making direct alignment of clinical trajectories across individuals difficult. Our framework addresses this challenge by projecting discrete patient histories onto a “model-inferred risk continuum” derived from longitudinal disease representations. This continuum provides a common axis for characterizing heterogeneous patient trajectories according to increasing model-estimated risk rather than chronological timing alone.

The resulting risk-aware organization enables investigation of how comorbidity patterns vary across the inferred risk spectrum. The observed enrichment of cardiorenal-metabolic conditions, including type 2 diabetes (E11) and CKD (N18), along the pseudotime axis highlights their importance in the disease trajectories captured by the model. Furthermore, the striking contrast between the gradual accrual of these chronic comorbidities and the abrupt, infrequent appearance of terminal events like shock (R57) validates the framework’s capacity to prioritize structural disease evolution over episodic triggers.

Collectively, these findings suggest that the proposed geometric trajectory framework captures longitudinal patterns of multimorbidity associated with increasing HF risk. Beyond predictive modeling, the risk-aware pseudotime framework provides an interpretable approach for visualizing disease trajectory organization and identifying patterns of comorbidity evolution preceding high-risk states.

Real-world longitudinal EHR data are inherently sparse and heterogeneous, often containing irregular follow-up intervals, variable numbers of historical admissions, and incomplete diagnosis coding due to administrative encounters or heterogeneous documentation practices. Such uncertainty may affect the construction of longitudinal disease trajectories by introducing admissions with limited disease-specific information. In our robustness analysis, retaining additional diagnosis-free admissions resulted in only marginal reductions in predictive performance, suggesting that the proposed trajectory representation maintained stable performance under moderate levels of incomplete diagnostic information. One possible explanation is that the representation aggregates disease information across multiple longitudinal encounters rather than relying on individual admissions alone, thereby reducing the impact of isolated incomplete records on the overall trajectory representation.

In addition, the acute-care focus of the MIMIC-IV database may introduce cohort selection bias, as patients receiving frequent or severe hospital-based care may be overrepresented in the available longitudinal records. This may limit the generalizability of the observed disease trajectories and risk patterns to broader populations, particularly those managed predominantly in outpatient or community settings.

Nevertheless, several limitations should be acknowledged. First, although the proposed framework naturally accommodates variable-length trajectories, patients with only a single prior admission inevitably contribute less longitudinal information than patients with multiple readmissions. Consequently, the richness of the learned trajectory representation remains dependent on the amount of historical clinical information available. Second, the present robustness analysis evaluated only one practical form of structural uncertainty. Other sources, including varying degrees of missing diagnoses, systematic coding bias, and missing encounters, warrant further investigation. Finally, this work lacks external validation on independent multicenter cohorts. Although the framework demonstrates strong internal validity, its generalizability across heterogeneous health care systems remains to be established in future studies.

A further limitation of this study is the exclusion of diagnoses that are highly correlated with HF, including hypertensive heart disease and cardiomyopathies. Although this strategy is methodologically rigorous from a machine learning perspective, being adopted to reduce information leakage and ensure that the model learned broader longitudinal comorbidity patterns rather than relying on near-proxy indicators of the outcome, it also results in a simplified representation of the patient's clinical condition. In routine clinical practice, these diagnoses constitute important etiological and prognostic information and are frequently incorporated into clinical decision-making. Consequently, the present framework should be interpreted primarily as an evaluation of the predictive value of latent comorbidity structure rather than a fully comprehensive clinical risk assessment model. Future work should investigate the integration of both disease trajectory features and direct HF-related diagnoses to better reflect real-world clinical decision support scenarios.

Future work should focus on external validation across diverse cohorts, incorporation of multimodal clinical information, and further development of advanced trajectory modeling approaches to better characterize complex disease evolution across clinical domains.

In conclusion, we developed a longitudinal disease trajectory representation that integrates IG-based weighting, disease co-occurrence information, and geometric structural features to transform heterogeneous EHR data into compact and interpretable representations of clinical evolution. The proposed framework captures latent patterns associated with disease trajectories and provides a risk-aware perspective for longitudinal risk stratification in high-noise clinical data. Beyond predictive performance, this approach facilitates visualization and characterization of disease trajectories, providing a framework for further investigation of trajectory-based risk assessment in longitudinal clinical data.