Integrating molecular dynamics and machine learning to identify potential apo-state conformational and solvent-exposure signatures associated with resistant KRAS mutants

1 Introduction

Mutation-induced drug resistance presents a major challenge in cancer therapy, particularly for targeted treatments. Such resistance can arise through various mechanisms (Hou et al., 2014), including activation of alternative signaling pathways, e.g., BRAF (V600E) inhibition by vemurafenib in colon cancer (Prahallad et al., 2012); impaired apoptosis, as seen in p53 (R273W) resistance to 5-fluorouracil (Huang et al., 2019); and enhanced DNA repair, e.g., BRCA1/2 reversion mutations conferring resistance to polyadenosine diphosphate ribose polymerase inhibitors (Collet et al., 2024). This study focuses on a resistance driven by direct mutations in the drug target itself. Such mutations frequently alter the chemical and conformational properties of the target, often within the binding site, disrupting drug-receptor interactions. This behavior, among others, is especially concerning for one of the most frequently mutated oncogenes, Kirsten rat sarcoma (KRAS) viral oncogene homolog (Yang et al., 2023; Huang et al., 2021).

The KRAS oncogene encodes a 21 kDa guanosine triphosphate hydrolase (GTPase) K-Ras (KRAS henceforward) that regulates key signaling pathways involved in cell-cycle control (Huang et al., 2021; Pantsar, 2020). Structurally, KRAS consists of a G-domain (residues 1–166), which binds guanine nucleotides, and a hypervariable region (residues 167–189), which mediates membrane association (Pantsar, 2020). Its secondary structure includes 6 β-sheets and 5 α-helices connected by flexible loops. Important regions (shown in Figure 1 top row) include the switch I and switch II loops, which interact with effector and regulatory proteins, and the P-loop, which binds nucleotide phosphate groups (Pantsar, 2020). KRAS homolog alternates between an inactive guanosine diphosphate– (GDP) bound and an active guanosine triphosphate– (GTP-) bound state. Guanine nucleotide exchange factors (GEFs) promote the release of GDP, allowing GTP, which is more abundant in the cell, to bind and stabilize the active conformation (Bos et al., 2007). Binding of GTP activates KRAS by inducing conformational changes in the switch I and switch II loops. The activated protein then engages downstream signaling pathways, including the mitogen-activated protein kinase (MAPK) pathway, the Ras-like small GTPases pathway, and the phosphoinositide-3-kinase pathway, which are critical regulators of cell proliferation, survival, growth, and differentiation (Li et al., 2023). Termination of signaling requires GTP hydrolysis, facilitated by GTPase-activating proteins, returning KRAS to its inactive state (Bos et al., 2007).

Mutations in signaling pathways such as Ras/Raf/MAPK are frequently implicated in cancer, with KRAS being one of the most commonly mutated oncogenes (Huang et al., 2021) accounting for 11.6% of all carcinomas (Yang et al., 2023). In pancreatic ductal adenocarcinoma, KRAS mutations are present in 81.7% of cases, with G12D as the predominant subtype (Yang et al., 2023). In colorectal cancer, KRAS mutations occur in 38% of patients, primarily G12D and G12V (Yang et al., 2023), while in non–small-cell lung cancer, the prevalence is 21.2%, with G12C representing 45% of cases (Yang et al., 2023; Huang et al., 2021; Tan and Tan, 2022; El Osta et al., 2019). Point mutations at KRAS residues G12, G13, and Q61 impair GAP-mediated GTP hydrolysis, leading to constitutive KRAS activation and tumorigenesis (Parker and Mattos, 2018).

Until recently, KRAS-driven tumors were considered difficult to treat because of KRAS’s picomolar affinity for GTP and the absence of other druggable pockets on its surface. As a result, therapeutic strategies focused on targeting KRAS indirectly (Blasco et al., 2019; Gi et al., 1997). A breakthrough came in 2013, when Ostrem et al. (2013) reported the first covalent inhibitor targeting the KRAS G12C mutant. The compound bound specifically to a cysteine residue within an allosteric pocket located beneath the switch II region in the inactive (GDP-bound) KRAS conformation. This breakthrough led to the development of the first KRAS G12C inhibitors approved for clinical use. In 2021, AMG510 (Canon et al., 2019) (sotorasib, Lumakras) was approved by the FDA for treating non–small-cell lung cancer, followed by MRTX849 (Fell et al., 2020) (adagrasib, Krazati) in 2022. Both drugs irreversibly bind to the cysteine residue at position 12 and form noncovalent interactions with other amino acids in the switch II binding site, including tyrosine (Y96), glutamine (Q99), and histidine (H95). This interaction network stabilizes the inactive, GDP-bound conformation of KRAS, thereby inhibiting oncogenic signaling (Zhang and Nagasaka, 2021; Mahran et al., 2024).

Despite initial clinical success, resistance to sotorasib and adagrasib has emerged through various mechanisms (Liu et al., 2022; Mohanty et al., 2023). One involves RTK-mediated activation of wild-type RAS (Ryan et al., 2020). Another features secondary mutations in KRAS G12C, such as Y96D (Tanaka et al., 2021; Rosen et al., 2023; Koga et al., 2021), Y96S (Rosen et al., 2023; Koga et al., 2021), and Y96C (Rosen et al., 2023; Awad et al., 2021), which confer resistance to both approved drugs. Other mutations, such as G12C/G13D, R68M, A59S, and A59T, are resistant to sotorasib alone, while Q99L is associated with resistance to adagrasib (Rosen et al., 2023; Koga et al., 2021). These secondary mutations reactivate ERK signaling and are likely to promote broader resistance across tumor-cell populations (Zhao et al., 2021).

Current efforts to target KRAS-driven tumors include inhibitors of GDP-bound KRAS G12C and G12D, as well as proteolysis-targeting chimeras for KRAS G12D and G12V, all in clinical trials (Isermann et al., 2025). Alternative strategies focus on the active, GTP-bound form of KRAS and binding regions outside of switch II to circumvent resistance from secondary mutations (Holderfield et al., 2024). However, the molecular basis of resistance arising in the inactive-state binding site, particularly in KRAS secondary mutants, remains poorly understood, underscoring the potential benefit of targeted approaches against these specific variants.

This study moves beyond prevailing studies focused on conformational analyses (Tu et al., 2022), single KRAS mutants (Pandey and Roy, 2025; Mir et al., 2024; Mehmood et al., 2021) or drug design (Srisongkram and Weerapreeyakul, 2022). Instead, we investigate the underlying factors of mutation-induced resistance caused by secondary mutations in the KRAS switch II binding site: KRAS G12C/Y96D, G12C/Y96C, and G12C/Y96S, which exhibit resistance to both sotorasib and adagrasib and are located on the same residue in the binding site for all cases (unlike, e.g., R68S or H95D, not studied here). As both drugs interact with residue Y96 (Mahran et al., 2024), substitution at this site represents a direct structural change relative to the G12C single mutant, which may influence drug binding. However, other molecular features may also contribute to resistance. To address this, we developed a framework combining molecular dynamics– (MD-) derived descriptors with machine learning (ML) algorithms to identify features that differentiate treatment-resistant KRAS mutants from treatment-sensitive ones.

Integrating dynamic properties of biological systems with ML offers a powerful approach to tackling complex challenges in drug discovery (Stitou et al., 2024). While MD combined with ML has been applied to characterize mutation-induced drug resistance using ligand-bound complexes (Xie et al., 2025; Intravaia et al., 2025), the apo (unbound) state, capturing intrinsic mutant dynamics prior to drug engagement is less explored in this context. Here we present an apo-state MD with ML classification supported by clustering-aware statistical validation, and we report both naive frame-level (MD-derived conformers treated as independent observations) and group-aware (mixed held-out, system-level) performance estimates to explicitly quantify the inflation that pseudoreplication introduces.

2 Results

This study investigated mutation-induced drug resistance using the KRAS oncogene as a proof-of concept study. Based on experimental treatment responses, KRAS structures were classified either as treatment-sensitive (G12C) or treatment-resistant (G12C/Y96C, G12C/Y96S, and G12C/Y96D). Because the secondary mutations occur within the switch II binding site, we first assessed structural differences in this region using SiteMap (Schrödinger, LLC, 2024) (Figure 1 bottom row). One key output, SiteScore, reflects the likelihood of small-molecule binding, with values near one indicating promising binding sites. KRAS WT exhibited a shallow switch II pocket with a low SiteScore (0.754) and small volume (123 ų), suggesting unlikely small-molecule binding, consistent with the absence of endogenous ligands in native KRAS. In contrast, the G12C single mutant and all secondary mutants had a SiteScore near 1, but their binding pocket volumes varied: G12C had a smaller volume (256 ų), while the secondary mutants ranged from 325 to 353 ų. This volume increase is likely attributable to substitution of bulky tyrosine at position 96 with smaller residues (cysteine, serine, or aspartic acid).

To capture dynamic structural changes not evident from static analysis, we performed MD simulations of each KRAS mutant in the absence of inhibitors. From the MD trajectories, we extracted a comprehensive set of molecular descriptors reflecting global structural behavior and localized changes within key regions of the protein. These features were categorized into four groups: structural, energetic, thermodynamic, and intra- and intermolecular contacts. Descriptors strongly associated with resistance, including those related to the site of the secondary mutation (residue 96) and the average solvent-accessible surface area (SASA (AV)) of residues G10 and A11, were excluded from the subsequent analysis to avoid biasing the models.

2.1 Univariate statistical analysis with correction for conformer non-independence

Prior to formal statistical analysis of individual molecular descriptors, we first examined overall simulation behavior to confirm that treatment-sensitive and treatment-resistant systems sampled distinct conformational regimes. Protein backbone RMSD relative to the reference structure was used to assess global structural stability throughout the simulations (Figure 2). Treatment-sensitive systems showed comparatively lower and more consistent deviation (RMSD range 0.13–0.32 nm, maximum 0.41 nm, Figure 2A), whereas treatment-resistant systems showed greater overall fluctuation and a wider range (RMSD range 0.13–0.39 nm, maximum 0.50 nm, Figure 2B), consistent with the resistant mutants sampling a broader conformational landscape. Residue-level flexibility was further examined using root mean square fluctuation (RMSF; Figures 3A,B), which showed consistently high peaks (>0.4 nm) at residues T35, E63, and Y64 across both groups, indicating these as inherently flexible regions independent of resistance status. Notably, P34 and E62 showed visually greater fluctuations specifically in treatment-sensitive systems, whereas E37 and Q61 showed elevated fluctuation specifically in treatment-resistant systems, a pattern broadly consistent with the subsequent residue-level statistical analysis.

To identify molecular descriptors that differed significantly between treatment-sensitive and treatment-resistant systems, we performed univariate statistical testing on each of the 101 descriptors spanning structural, energetic, thermodynamic, and intra- and intermolecular contact categories. Treating each of the 126 MD conformers as an independent observation inflates the effective sample size, because conformers within a single trajectory are temporally autocorrelated and should not be treated as statistically independent (Grossfield and Zuckerman, 2009; Bolker et al., 2009). Consistent with this concern, the naive conformer-level analysis, which ignores this clustering, identified 38 of 101 features as significant after Benjamini-Hochberg FDR correction (Benjamini and Hochberg, 1995), we present these results in the Supplementary Table S1.

To correct for this pseudoreplication, we applied two complementary clustering-aware approaches. First, following the block-averaging strategy recommended for MD-derived statistical comparisons (Knapp et al., 2018), each of the 18 independent trajectories (6 systems × 3 replicas) was summarized by its mean descriptor value prior to statistical testing, as a result collapsing the seven conformers per trajectory into a single observation. Mann-Whitney U tests (Mann and Whitney, 1947; Virtanen et al., 2020) were then performed on these 18 trajectory-level means (9 sensitive, nine resistant), with Benjamini-Hochberg FDR correction applied across all 101 descriptors. This substantially reduced the number of significant features from 38 to seven (q = 0.024–0.044), comprising LJ 1,4 energy, RMSF H95, RMSF M72, SASA (AV) H95, and SASA (SD) for G10, E62, and Q99 (Table 1).

TABLE 1. Univariate statistics summary: features statistically significant in Replica-level Mann-Whitney U and LMM after Benjamini-Hochberg FDR correction across all 101 tested features.

r: rank-biserial correlation (range −1, to +1); P (replica MWU): Mann-Whitney U on 18 replica-means (nine sensitive, nine resistant); LMM: Linear Mixed-Effects Model; q: Benjamini-Hochberg FDR, correction across all 101 tested features; β: standardized LMM, coefficient; ICC: intraclass correlation (system level); Marginal R2: variance explained by Resistance fixed effect. Grey cells–feature not significant in either model (q > 0.05). Abbreviations: AV, average; SD, standard deviation; SASA, solvent accessible surface area; LJ, Lennard-Jones; RMSF, root mean square fluctuation.

Second, as an independent check operating directly on all 126 conformers, we fitted a linear mixed-effects model (LMM) (Seabold, 2010) for each descriptor (feature ∼ Resistance + (1|System)), with System (n = 6) included as a random intercept to absorb system-level clustering. Six of the seven replica-level significant features were also significant in the LMM after FDR correction (q = 3.04 × 10^-6–1.17 × 10^-3). Because the LMM relies on only six random-effect levels, its fixed-effect p-values may be anti-conservative (Bolker et al., 2009); we therefore treat LMM significance as corroborating rather than independent confirmatory evidence, interpreted alongside the non-parametric replica-averaging result. Notably, the two approaches address the same non-independence concern at different levels of the data hierarchy, block-averaging collapses conformers before testing, whereas the LMM retains all 126 conformers but partitions the variance through the random intercept, making their convergence stronger evidence than either method alone. The LMM additionally flagged SASA (SD) V9, SASA (AV) K16, RMSF Q99, RMSF G60, and total protein SASA as significant, which were not significant in the replica-level analysis and are therefore considered weaker candidates pending further validation.

Six descriptors, LJ 1,4 energy, RMSF H95, RMSF M72, SASA (AV) H95, SASA (SD) G10, and SASA (SD) E62, were significant in both clustering-aware analyses (Table 1, in bold; Figure 4) and are reported as the primary hypothesis-generating candidates. All six showed large effect sizes (replica-level rank-biserial correlation |r| = 0.852–1.000; standardized LMM coefficients |β| = 0.776–1.277), with consistent direction of effect between the two methods: RMSF H95, SASA (AV) H95, SASA (SD) G10, and RMSF M72 were elevated in resistant conformers, while LJ 1,4 energy and SASA (SD) E62 were reduced (Figure 4). LJ 1,4 energy showed perfect rank-order separation between classes (r = −1.000), with every sensitive replica-mean exceeding every resistant replica-mean; RMSF H95 approached perfect separation (r = 0.988). These two features also had the lowest BH-FDR corrected p-values in the replica-level analysis (q = 0.024 each), consistent with their near-complete class separation.

To characterize the extent to which system-level clustering affected each feature’s LMM estimate, we report the intraclass correlation coefficient (ICC) and marginal R ² for all six convergent features. Four features showed low ICC (<0.05: LJ 1,4 energy = 0.014, RMSF H95 ≈ 0.000, SASA (SD) E62 = 0.009, RMSF M72 = 0.004), indicating that system identity contributed minimally to their variance and that the between-class (sensitive vs. resistant) signal dominated; their marginal R ² ranged from 0.150 (RMSF H95) to 0.242 (LJ 1,4 energy), indicating that the fixed effect of resistance status explained 15–24% of the total variance after accounting for system-level random effects. SASA (AV) H95 and SASA (SD) G10 showed the largest standardized fixed effects coefficients (β = 1.235 and 1.277, respectively) but also measurable system-level clustering (ICC = 0.112 and 0.125, respectively), meaning a larger proportion of their variance was attributable to system identity; their LMM derived p-values should therefore be interpreted with somewhat greater caution, though both remained significant in the independent replica-level analysis.

2.1.1 Generalizability within resistant mutants

Beyond statistical reliability evaluated using ICC and marginal R², a separate and equally important concern is whether these six convergent features reflect a generalizable resistance phenotype or instead capture Y96 substitution-specific effects. Because the resistant class is composed of three chemically distinct mutations at residue 96 (G12C/Y96C, G12C/Y96D, G12C/Y96S), and because cysteine, aspartate, and serine differ in size, charge, and side chain chemistry, nearby conformational and solvent-exposure descriptors could in principle reflect the specific identity of the substitution rather than a shared resistance mechanism. This possibility remains even after excluding descriptors that directly measure residue 96. We reasoned that if a convergent feature tracks the Y96 substitution rather than resistance, its value should differ systematically among Y96C, Y96D, and Y96S; conversely, if the feature reflects a shared resistance phenotype, its values should remain similar across all three resistant mutants despite differences in the physicochemical properties of the substituted residues.

To test this, we compared the distribution of each convergent feature across the three resistant mutants using the Kruskal-Wallis test (Kruskal and Wallis, 1952), complementing the system-level variance assessment provided by ICC and marginal R ². Where ICC and marginal R² characterize how much of a feature’s overall variance is attributable to system identity (sensitive and resistant systems combined), the within-resistant Kruskal-Wallis test addresses this question specifically within the resistant class, directly testing whether mutant identity, rather than resistance status, explains the observed variation. Five of six features showed no significant difference in central tendency across the three resistant mutants (Kruskal-Wallis p = 0.196–0.870), indicating that Y96C, Y96D, and Y96S, despite their distinct physicochemical properties at residue 96, share a consistent typical value for these descriptors. This finding is consistent with the low ICC values observed for four of these features (LJ 1,4 energy, RMSF H95, RMSF M72, SASA (SD) E62; ICC <0.05) and supports interpreting the corresponding sensitive-versus-resistant differences as reflecting a shared resistance-associated effect rather than a Y96 substitution-specific effect (Table 2; Supplementary Figure S1).

TABLE 2. Kruskal-Wallis test across three resistant mutants (Y96C vs. Y96D vs. Y96S) for six convergent features identified by univariate models.

P (KW), Kruskal-Wallis; η2, Eta-squared. Abbreviations: AV, average; SD, standard deviation; SASA, solvent accessible surface area; LJ, Lennard-Jones; RMSF, root mean square fluctuation.

The exception was SASA (SD) G10 (Kruskal-Wallis p = 0.005, η² = 0.144), which differed significantly among the three resistant mutants, driven primarily by a difference between Y96D and Y96S (Dunn’s post-hoc test (Dunn, 1964), BH-FDR corrected q = 0.0009; Table 3). This result is consistent with SASA (SD) G10 also showing one of the highest ICC values among the convergent features (ICC = 0.125) and the largest standardized LMM effect (β = 1.277), together suggesting that a meaningful portion of this feature’s apparent association with resistance may instead reflect substitution-specific structural variation rather than a shared resistance-associated mechanism. We therefore treat SASA (SD) G10 as a less reliable candidate for a generalizable resistance marker, despite its strong performance in the primary convergence analysis. SASA (AV) H95, by contrast, showed both moderate system-level clustering (ICC = 0.112) and no significant within-resistant heterogeneity (Kruskal-Wallis p = 0.247), indicating that its moderate ICC reflects broader sensitive-versus-resistant variation rather than substitution-specific effects among the resistant mutants themselves.

TABLE 3. Post-hoc Dunn’s test for features significant in Kruskal-Wallis test across three resistant mutants (Y96C vs. Y96D vs. Y96S).

P, Post-hoc Dunn’s test; q, Benjamini-Hochberg FDR, correction across all 101 tested features. Abbreviations: SD, standard deviation; SASA, solvent accessible surface area.

2.2 Machine-learning models

We utilized two complementary validation schemes for machine learning models (Logistic Regression [LR], Random Forest [RF], and Support Vector Machine [SVM]) (Pedregosa et al., 2011) to identify the most informative features capable of distinguishing treatment-resistant mutants from treatment-sensitive ones. The first validation scheme consisted of a random 80/20 split, in which each model was trained on 80% of the conformers and tested on the remaining 20%. More rigorous validation strategy consisted of mixed held-out validation, in which all nine possible combinations of one held-out sensitive system and one held-out resistant mutant system were used as test sets, with the model trained on the remaining four systems in each case.

Model performance differed substantially between the two validation schemes (Table 4 vs; Table 5). Under random 80/20 cross-validation, all seven model configurations achieved high and broadly similar performance (accuracy 0.900–0.931; AUC 0.945–0.984, Table 4), with minimal separation between methods or between default and optimized variants (run by run performance in Supplementary Table S2 and ROC-AUC curves in Supplementary Figure S2. For sensitivity check, we ran two additional ML models (XGBoost and PLS-DA), which remained in agreement with the results (Supplementary Tables S3, S4; Supplementary Figures S3 and S4). This uniformly strong performance is consistent with the conformer-level data leakage inherent to this scheme, in which highly correlated conformational states from the same system can appear in both training and test partitions, and should not be interpreted as evidence of generalization to new systems.

TABLE 4. Performance of the ML Models; random 80/20 validation scheme: Logistic Regression (LR), Random Forest (RF) and Support Vector Machine (SVM) with default and optimized parameters (summary accross runs).

Logistic Regression (LR):

LR (L2 default): C = 1.0, solver = lbfgs, penalty = L2 (default parameters).

LR (L1 optimized): penalty = L1, solver = saga; C selected via GridSearchCV (5-fold CV, re-run inside each of 10 runs): values ranged across runs: C ∈ [0.1, 1.0, 10.0, 100.0].

LR (Elastic Net optimized): penalty = Elastic Net, solver = saga; best parameters from GridSearchCV (5-fold CV, re-run inside each of 10 runs): C ∈ [0.1, 1.0, 10.0], l1_ratio ∈ [0.1, 0.3, 0.5, 0.7, 0.9].

Random Forest (RF):

RF (default): n_estimators = 100, max_depth = None, min_samples_split = 2, min_samples_leaf = 1, max_features = ‘sqrt’ (default parameters).

RF (optimized): best parameters from GridSearchCV (5-fold CV, re-run inside each of 10 runs): n_estimators ∈ [50, 100, 200], max_depth = None, max_features ∈ [sqrt, log2], min_samples_leaf ∈ [1, 2], min_samples_split ∈ [2, 5].

Support Vector Machine (SVM):

SVM (default): C = 1.0, gamma = ”scale”, kernel = RBF (default parameters).

SVM (optimized): best parameters from GridSearchCV (5-fold CV, re-run inside each of 10 runs): C = 10, gamma ∈ [0.001, 0.01], kernel = RBF.

TABLE 5. Performance of the ML Models; mixed held-out validation scheme: Logistic Regression (LR), Random Forest (RF) and Support Vector Machine (SVM) with default and optimized parameters (summary across folds and runs).

Logistic Regression (LR):

LR (L2 default): C = 1.0, solver = lbfgs, penalty = L2 (default parameters).

LR (L1 optimized): penalty = L1, solver = saga; C selected via GridSearchCV (5-fold CV, re-run inside each of 10 runs): values ranged across runs: C ∈ [1.0, 10.0, 100.0].

LR (Elastic Net optimized): penalty = Elastic Net, solver = saga; best parameters from GridSearchCV (5-fold CV, re-run inside each of 10 runs): C ∈ [0.1, 1.0, 10.0], l1_ratio ∈ [0.1, 0.3, 0.5, 0.7, 0.9].

Random Forest (RF):

RF (default): n_estimators = 100, max_depth = None, min_samples_split = 2, min_samples_leaf = 1, max_features = ‘sqrt’ (default parameters).

RF (optimized): best parameters from GridSearchCV (5-fold CV, re-run inside each of 10 runs): n_estimators ∈ [50, 100, 200], max_depth = ∈ [5.0, None], max_features ∈ [sqrt, log2], min_samples_leaf ∈ [1, 2], min_samples_split ∈ [2, 5].

Support Vector Machine (SVM):

SVM (default): C = 1.0, gamma = ‘scale’, kernel = RBF (default parameters).

SVM (optimized): best parameters from GridSearchCV (5-fold CV, re-run inside each of 10 runs): C = ∈ [1, 10, 100], gamma ∈ [0.001, 0.01, scale], kernel = RBF.

Under mixed held-out validation, performance dropped substantially and diverged sharply by method (Table 5). Random forest retained the strongest and most stable performance (AUC = 0.909 ± 0.098 default, 0.904 ± 0.104 optimized), followed by SVM default (AUC = 0.827 ± 0.137). Logistic regression generalized most poorly overall: the L2-regularized default model achieved AUC = 0.535 ± 0.388, close to chance, while L1-optimized and Elastic Net-optimized variants partially recovered performance (AUC = 0.747 ± 0.281 and 0.685 ± 0.292, respectively) but with substantially higher variance than RF or SVM. Specificity showed the largest gap between schemes across all models (e.g., LR L2 default: 0.869 ± 0.154 under 80/20 versus 0.286 ± 0.360 under mixed held-out), indicating that the leakage-inflated 80/20 scheme particularly overstated each model’s ability to correctly reject sensitive conformers. Run-by-run results and ROC-AUC curves of mixed held-out validation scheme are shown in Supplementary Table S5; Supplementary Figure S5.

The poor overall performance of LR under mixed held-out validation was driven almost entirely by a single PDB system. Examining AUC by individual held-out fold revealed that LR (L2, default) failed specifically and consistently when G12C (6mtb) was the held-out sensitive system, regardless of which resistant mutant was paired with it (AUC = 0.011, 0.000, and 0.000 for Y96C, Y96D, and Y96S, respectively; Figure 5A). In contrast, both RF and SVM (default) generalized successfully to 6mtb across all three resistant pairings (RF: AUC = 0.92–1.00; SVM: AUC = 0.76–1.00), including the Y96S pairing, which was the most difficult fold overall. This indicates that the poor performance of G12C from 6mtb may reflect a method-specific limitation of the linear decision boundary learned by L2-regularized LR, likely because the conformational feature space of this crystal structure with manual mutation differs from G12C (6oim) and G12C (6ut0) in a way that default L2 regularization cannot accommodate, rather than a fundamental property of the underlying data. Further support for this interpretation comes from L1-optimized LR, whose sparsity-inducing penalty selects a smaller, more discriminative feature subset and partially improved performance on the 6mtb + Y96C and 6mtb + Y96D folds (AUC = 0.968 and 0.660) yet still failed on 6mtb + Y96S (AUC = 0.000), suggesting that even feature selection cannot fully compensate for a insufficient linear decision boundary in the most difficult fold. Hyperparameter optimization did not consistently improve robustness: optimized SVM showed a collapse on the 6mtb + Y96C fold (AUC = 0.147, versus 0.952 for SVM default), illustrating that aggressive regularization tuning (C = 100) is prone to overfitting fold-specific structure in a small, system-correlated training set rather than improving generalization.

Across all nine held-out combinations, performance also varied systematically by which sensitive system was held-out, independent of model choice: folds withholding G12C (6ut0) were uniformly strong (AUC 0.81–0.99 across methods), G12C (6oim) folds were moderate (AUC 0.60–0.86), and G12C (6mtb) folds ranged from very poor (LR) to excellent (RF, SVM default), highlighting that model choice had the greatest impact for this held-out system. RF was the only method to perform consistently well across all three sensitive systems and all three resistant mutants, supporting that it is the most robust classifier among those evaluated for this dataset, despite the limited number of independent systems.

Taken together, the results are consistent with the methodological concern motivating this validation strategy: random conformer-level splitting substantially inflates apparent performance, whereas the mixed held-out validation strategy used here, which respects system-level independence, provides a credible estimate of generalization to unseen KRAS systems. Furthermore, the system- and method-specific nature of the 6mbt failure demonstrates that the performance limitations are more likely attributable to model choice and hyperparameter sensitivity than to a complete absence of generalizable signal in the underlying molecular descriptors.

Having established that mixed held-out validation provides the more credible estimate of model generalization, we next examined which molecular descriptors drive each model’s performance, and whether this discriminative importance reflects a shared resistance phenotype or instead dataset-specific structure exploited by the classifier. Feature importance rankings were highly consistent between the random 80/20 and mixed held-out validation schemes, with the top-ranked features overlapped substantially across both designs (12 of the top 15 features shared between schemes, Figure 5B vs. Supplementary Figure S2B). This stability indicates that the feature importance rankings is not an artifact of a particular train/test split. However, mixed held-out validation addresses only leakage between systems; within each training fold, the conformers remain correlated observations from four underlying systems, so feature importance estimates, even under this stricter validation scheme, remain susceptible to the same pseudo-replication concerns raised for the univariate analysis. This is illustrated by SASA (SD) G10, which ranked first in both validation schemes (normalized importance 0.88–0.94) yet was also the only convergent feature to show significant within-resistant-group heterogeneity (Kruskal-Wallis p = 0.005). This discrepancy reflects a fundamental distinction between what ML importance and univariate clustering-aware testing measure and is central to the interpretation of the following results. ML feature importance quantifies how useful a descriptor is for a classifier to discriminate sensitive from resistant conformers within this dataset, it does not indicate whether that discriminative power reflects a shared resistance-associated phenotype or a mutant-specific structural variation. A feature, could in principle, be highly useful for classification because it varies systematically with the identity of the resistant mutant, rather than with resistance status. Under mixed held-out validation, where the model is trained on two of the three resistant mutants and tested on the third, such a feature could still appear discriminative if the held-out mutant’s structural signature correlates with the resistant class overall. The within-resistant Kruskal-Wallis test addresses this directly by asking whether a feature’s value is consistent across G12C/Y96C, Y96D, and Y96S. A feature that varies substantially between these three chemically distinct substitutions is partly capturing mutant identity rather than a shared resistance phenotype, regardless of how well it performs in classification. We therefore treat ML importance and within-resistant univariate testing as answering distinct questions, discriminative utility versus cross-mutant generalizability, and treat only features satisfying both criteria as credible candidates for a generalizable resistance signal.

To formalize this distinction, we cross-referenced the top-ranked features from both validation schemes against the six statistically convergent descriptors identified by replica-level Mann-Whitney U testing and LMM (Figure 5C). Five of the six convergent features, SASA (SD) E62, SASA (AV) H95, LJ 1,4 energy, RMSF M72, and RMSF H95, showed no significant heterogeneity among resistant mutants. The same features were also highly ranked by the ML feature importance under mixed held-out validation: SASA (SD) E62 and SASA (AV) H95 each carried similarly high normalized importance (∼0.50), followed by LJ 1,4 energy (0.33), RMSF M72 (0.24), and RMSF H95 (0.15), all ranking within the top 16 features overall. This convergence across three independent lines of evidence, replica-level testing, mixed-effects modeling, and ML importance under system-independent validation, provides the strongest support for these five descriptors as shared resistance-associated signals (Figure 5C, Tier 1).

By contrast, SASA (SD) G10, despite ranking first in ML importance under both validation schemes (normalized importance 0.88–0.94), was flagged by significant within-resistant Kruskal-Wallis testing (p = 0.005) and is reported separately (Figure 5C, Tier 2) as likely reflecting Y96 substitution-specific structural variation rather than a shared resistance-associated phenotype. This feature illustrates precisely the distinction outlined above: its high discriminative power for the classifier does not establish it as a generalizable resistance marker, and we caution against interpreting ML importance rankings as direct evidence of biological relevance without independent univariate corroboration.

The remaining top-ranked ML features (Figure 5C, Tiers three and 4) were either significant only in the LLM and not in the more conservative replica-level test (Tier 3: SASA (SD) V9, SASA (AV) K16, RMSF Q99 and G60, or were not statistically supported by either clustering-aware univariate analysis at all (Tier 4). These features should be interpreted with particular caution: feature importance, even under the stricter mixed held-out validation scheme, remains susceptible to within-fold conformer correlation, since each training fold still comprises conformers drawn from only four underlying systems. High ML importance in the absence of independent univariate support therefore does not exclude the possibility that the feature reflects overfitting to system-specific structure rather than a genuine resistance-associated signal.

As the last analysis, we utilized Bayesian Network models (Gogoshin et al., 2017) to understand the underlying structure of relationships between the molecular features. First, BN inference was used to explore how treatment response status influences the probabilistic behavior of molecular descriptors. When the treatment resistance node was fixed to either the sensitive or resistant state, the resulting shifts in marginal probabilities revealed that treatment resistance exerts a structured, nonrandom influence on specific network regions. The most strongly affected descriptors included the SASA (SD) for G10 and E62, the SASA (AV) of H95, and LJ 1,4 energy (Figure 6A), aligning with the features previously identified by ML models.

To assess whether these descriptors could, in turn, influence treatment response, each was fixed to one of its discretized states, and the resulting marginal probabilities of the treatment resistance node were computed. For comparison, two additional descriptors, the MSD of the binding site backbone and the number of intraprotein hydrogen bonds, were also tested because they had little effect on resistance-state inference. This analysis revealed several nonlinear, descriptor-specific effects. Specifically, reduced values of SASA for G10 (SD) and H95 (AV), and increased values of SASA for E62 (SD) and LJ 1,4 energy, were associated with higher marginal probabilities of treatment sensitivity (Figure 6B). Conversely, the opposite descriptor states consistently shifted probabilities toward treatment resistance, consistent with their importance associated with the resistant phenotype (Figure 6B). These findings are consistent with the ML results identifying molecular features that distinguish the two groups.

To visualize the context of the top-ranked ML features, the most populated MD conformations from both treatment-sensitive and treatment-resistant KRAS systems were examined with the respect to the conformations of residues G10, E62, and H95 (Figure 7). Residues E62 and H95 exhibited distinct conformational dynamics across the most populated clusters, which may explain the differences in solvent exposure values and could influence ligand accessibility that warrant further investigation. Residue G10 was located within a highly flexible region, where neighboring residues may contribute to the variability of its solvent exposure.

3 Discussion

In this study, the unit of biological variation was the KRAS variant. We analyzed four structurally related KRAS variants (G12C; G12C/Y96C, G12C/Y96S, G12C/Y96D) using 126 conformers sampled from multiple PDB starting structures, independent replicas, and RMSD-based clusters. This design intentionally focused on a limited but clinically reported set of resistance-associated double mutants located within the inhibitor-binding site, with the explicit goal of generating hypotheses rather than establishing definitive mechanisms of resistance.

All molecular dynamics simulations were performed in the unbound (apo) state; no inhibitor-bound complexes were simulated. Conclusions are therefore restricted to apo-state conformational rearrangements and do not establish inhibitor-specific resistance mechanisms. Rather than identifying causal determinants of resistance, our results support the presence of apo-state conformational and solvent-exposure correlates associated with resistant KRAS mutants. These unbound-state findings provide a rationale for future testing with representative inhibitors but do not, by themselves, define the mechanisms of underlying resistance.

3.1 Feature importance for classification does not equal effect causality

We observed an important distinction between features that are useful for machine learning classification and features that may reflect a generalizable resistance phenotype across mutants. ML feature importance quantifies how useful a descriptor is for a classifier to discriminate sensitive from resistant conformers within this dataset; it does not indicate whether that discriminative power reflects a shared resistance-associated phenotype or a mutant-specific structural variation. A feature can therefore be important for classification without necessarily being a causal determinant of resistance, a distinction that is particularly important here given the small, structurally related set of mutants used to train the models.

Four complementary modeling approaches were used to interrogate this dataset: logistic regression, random forest, and support vector machine classifiers discriminated treatment-sensitive from treatment-resistant conformers based on molecular descriptors, while a Bayesian network (BN) was used to model probabilistic dependencies among features, offering a more interpretable, graph-based view of feature relationships. The BN identified largely consistent discriminative descriptors to those highlighted by the LR, RF, and SVM classifiers, providing an additional, methodologically distinct line of evidence for the discriminative relevance of these features and providing additional confidence that the ML importance rankings were not an artifact specific to any single modeling approach.

To assess whether this discriminative power reflects a generalizable resistance-associated phenotype, we combined these ML- and BN-derived importance rankings with clustering-aware univariate analysis (replica-level Mann-Whitney U testing and LMM) explicitly designed to account for the correlated nature of conformers sampled from the same MD trajectory. This integrated analysis revealed cases where the two forms of evidence diverged. SASA (SD) of G10, for example, was among the most discriminative features across the ML and BN models, yet also showed significant heterogeneity among resistant mutants (i.e., its distribution differed significantly among Y96C, Y96D, and Y96S), indicating that its classification utility likely reflects mutant-identity-specific variation rather than a shared property of the resistant mutants. By contrast, five features, SASA (SD) E62, SASA (AV) H95, LJ 1,4 energy, RMSF M72, and RMSF H95, showed no significant within-resistant-class heterogeneity while remaining prominently represented in both ML and BN importance rankings, making them the strongest hypothesis-generating candidates for a shared resistance-associated signature.

3.2 Structural plausibility of the identified features

Several of the descriptors highlight conformational and solvent-exposure changes of residues that have established roles in inhibitor binding in prior structural and biochemical studies (G10, K16, E62, and H95). Adagrasib forms a salt bridge with E62 and hydrogen bonds with H95 and K16, while its cyanomethyl group interacts with G10, displacing a conserved water molecule (Fell et al., 2020). Sotorasib similarly forms hydrogen bonds with K16, exploits a specific H95 orientation to enhance binding affinity to KRAS G12C, forms a hydrogen bond with E63, and engages in π–π interactions with Y96 (Canon et al., 2019). Solvation of G10 in particular has been discussed extensively in prior structural and MD studies (Fell et al., 2020; Tu et al., 2022; Leini and Pantsar, 2023). In our apo-state simulations, resistant systems showed increased solvent exposure at G10 and H95 and reduced at E62, suggesting altered solvation patterns and local rigidity around the binding site. We propose, as a hypothesis rather than a demonstrated mechanism, that these conformational features, even in the absence of ligand, could pre-organize resistant mutants into states less compatible with optimal inhibitor binding. This interpretation is consistent with independent ligand-bound MD studies reporting similar shifts in adagrasib and sotorasib contact frequencies at these same residues in related resistant mutants, including G12C/R68S, G12C/K16T, G12C/Y96C, and G12C/Y96D (Tu et al., 2022). Together, these convergent observations across independent apo- and ligand-bound studies are consistent with a model in which resistance-associated mutations are associated with subtle conformational changes that pre-organize the inhibitor-binding interface in the apo state, potentially altering drug contacts upon subsequent binding, though we emphasize that our data alone cannot establish this causal sequence.

3.3 Limitations

Several limitations should be considered when interpreting these findings. First, this analysis is restricted to a small number of structurally related variants (one sensitive background and three clinically reported Y96 resistant mutants), and should be regarded as a proof-of-concept comparison rather than a comprehensive survey of resistance mechanisms; the limited number of independent systems (n = 6) constrains the statistical power of all analyses and the generalizability of the identified features to other resistance-conferring mutations. Second, all simulations were performed in the apo state, without inhibitor-bound complexes; the conformational correlates identified here should therefore be interpreted as hypotheses about pre-existing structural tendencies rather than validated, inhibitor-specific resistance mechanisms, and require confirmation through ligand-bound simulations or experimental binding studies. Third, despite the use of clustering-aware statistical methods (replica-level averaging, LMM) to mitigate pseudoreplication, conformer-level non-independence remains an inherent feature of MD-derived datasets of this size, and residual correlation cannot be fully excluded. Finally, the machine learning models, trained on a small and structurally related set of mutants, should not be interpreted as having learned causal determinants of resistance; rather, they identify descriptors that are statistically associated with resistance status in this dataset, which we further prioritized using independent univariate evidence. Prospective computational and experimental testing on additional, structurally independent resistant variants, will be required to confirm whether the identified conformational and solvent-exposure signatures generalize beyond the mutants studied here.

3.4 Summary

In summary, we present a proof-of-concept analysis of KRAS mutation-associated conformational signatures, comparing one sensitive background with a small set of clinically reported Y96 resistant mutants. Comprehensive molecular descriptors, capturing structural, energetic, thermodynamic, and contact-based properties, were derived from MD simulations and analyzed using complementary methods. ML classifiers (LR, RF, SVM) discriminated treatment-sensitive from treatment-resistant conformers based on these descriptors, while a Bayesian network modeled probabilistic dependencies between features and identified largely concordant results, further supporting the discriminative relevance of the identified descriptors. Univariate statistical analysis provided independent insight into the generalizability of individual features across resistant mutants. Together, these approaches identified candidate apo-state conformational correlates associated with resistant mutants, particularly residue-level differences in flexibility and solvent exposure within the switch II binding region. While these findings do not establish causal or inhibitor-specific resistance mechanisms, they may inform the design of inhibitors targeting KRAS secondary mutants resistant to sotorasib and adagrasib and illustrate the value of integrating molecular dynamics with machine learning to address the complexity of mutation-induced drug resistance.

4 Methods

4.1 Computational chemistry

To investigate the dynamic behavior of KRAS mutants, we conducted all-atom MD simulations of the inactive KRAS conformations bound to Mg²⁺ and GDP. The initial protein structures were obtained from the Protein Data Bank (PDB) (Berman et al., 2000) under accession codes: 6OIM (Fell et al., 2020), 6UT0 (Canon et al., 2019), and 6MBT (Dharmaiah et al., 2019), and prepared before MD simulations. For the KRAS G12C mutant bound to the covalent inhibitor sotorasib (PDB ID: 6OIM), crystallographic water molecules were removed together with the ligand. Missing amino acid residues 105–107 were reconstructed using MODELLER 10.5 (Webb and Sali, 2016). Three engineered mutations, C51S, C80L, and C118S (Ostrem et al., 2013), present in the structure were reverted back to cysteines using the PyMOL 3.0 mutagenesis panel (Schrödinger, LLC, 2026). The protein was capped with an acetyl group at the N-terminus and with a methyl amide at the C-terminus using PyMOL’s builder module (Schrödinger, LLC, 2026). This structure served as the basis for building KRAS double mutant models, which have not yet been experimentally resolved. Subsequently, KRAS G12C (PDB ID: 6OIM) was manually mutated in PyMOL to construct models of the double mutants G12C/Y96C, G12C/Y96S, and G12C/Y96D. For the wild-type (WT) KRAS structure (PDB ID: 6MBT) and the KRAS G12C mutant bound to the covalent inhibitor adagrasib (PDB ID: 6UT0), the same preparation steps were applied. Additionally, the WT KRAS structure was mutated at position 12, from glycine to cysteine, to model an additional KRAS single-mutant system.

All protein models were processed using the CHARMM36 force field (Huang et al., 2017) in GROMACS (Van Der Spoel et al., 2005). Hydrogen atoms were added to each protein structure, and the protonation states of the residues were assigned at neutral pH. The GDP molecule was parametrized separately with the CHARMM General Force Field (Vanommeslaeghe et al., 2010). To prepare the systems for energy minimization, the complexes were placed in cubic boxes, extending at least 10 Å from the protein. The solvent was modeled using the CHARMM-modified TIP3P water model (Huang et al., 2017). The protein charges were neutralized by adding the appropriate amount of Na+ ions, and physiological NaCl salt concentration 0.1 M was ensured. Energy minimization was performed using the steepest descent algorithm (Press et al., 1986).

Next, NVT equilibration was conducted for 100 ps to stabilize the system’s temperature at 310 K with position restraints on protein and ligand heavy atoms to allow the solvent and ions to relax. Temperature was maintained at 310 K using the V-rescale thermostat (Bussi et al., 2007) with separate coupling groups for protein and GDP ligand, and solvent/ions. H-bonds were constrained using LINCS (Hess, 2008), enabling the 2 fs timestep. This was followed by 250 ps of NPT equilibration, continuing from the NVT step, with the same position restraints and temperature control, while pressure was maintained at 1 bar using isotropic Berendsen coupling (Berendsen et al., 1984). Production MD were performed for 200 ns with a 2 fs timestep. Temperature was maintained at 310 K using the V-rescale thermostat, and pressure was maintained at 1 bar using isotropic Parrinello-Rahman coupling. Coordinates and energies were saved every 10 ps. Each system underwent three replica runs, with different random seeds set before each simulation run. The water shell confirming octahedral arrangement of water molecules around Mg²⁺ ion is shown in Supplementary Figure S6.

To extract representative conformations from the trajectories, we performed trajectory clustering on the 50–200 ns portion of each simulation using the GROMOS (Daura et al., 1999) method. This method counts neighboring structures within a specified cutoff and identifies the representative frame as the one with the smallest average root mean square deviation (RMSD) to all its neighbors. Clustering was performed for six different Cα RMSD cutoffs, ranging from 0.10 nm to 0.15 nm, to find reasonable distribution of frames for each system. Representative structures from the top 7 clusters (∼90% of the trajectory) were extracted from each of the three replicas per system for molecular descriptors calculations. This resulted in 126 conformations across all six systems. Specific elements of each system were treated separately in the calculations, including extended switch II binding site residues (residues 9, 10, 11, 12, 16, 34, 35, 37, 58, 59, 60, 61, 62, 63, 64, 68, 69, 72, 95, 96, 99, 100, 102, and 103), which were defined as residues within 5Å of any site point identified in the KRAS G12C switch II binding site by the SiteMap module (Schrödinger, LLC, 2024; Halgren, 2009) of Schrödinger Suite and included motifs from P-loop, switch I region, and adjacent helices. A full list of molecular descriptors, grouped into categories, is provided in Table 6.

TABLE 6. The complete list of the molecular descriptors extracted from MD simulations.

First frame of the production run was used as the reference conformation.

4.2 Univariate analysis and ML

4.2.1 Dataset

A combined univariate and machine learning approach was used to identify molecular descriptors distinguishing treatment-sensitive from treatment-resistant KRAS mutants. The dataset consisted of three treatment-sensitive (Resistance = 0) and three treatment-resistant (Resistance = 1) KRAS systems, each contributing 63 conformers, for a total of 126 samples described by 101 molecular descriptors. Highly target-correlated features (correlation >75%) were removed, along with all descriptors directly related to residue 96, the site of the secondary mutation, to avoid bias toward known changes.

4.2.2 Univariate analysis

The Mann–Whitney U test (Mann and Whitney, 1947; Virtanen et al., 2020) was used to assess statistical significance between feature distributions in the treatment-sensitive and treatment-resistant groups, applied at two levels of resolution: (i) the full set of 126 conformers, and (ii) replica-level block averages. Feature distributions were first tested for normality using the Shapiro–Wilk test (Andrew and King, 2019), motivating the use of a non-parametric test.

To account for the non-independence of conformers sampled from the same MD trajectory, replica-level block averaging was performed: for each descriptor, the mean value was computed within each of the 18 replicas (6 systems × 3 replicas each), and the Mann–Whitney U test was applied to these 18 replica-level means rather than to individual conformers. This approach collapses correlated conformers into a single independent observation per replica. P-values were corrected for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995) across all 101 descriptors, applied separately to the conformer-level and replica-level analyses.

As a second, complementary approach to account for non-independence among conformers derived from the same molecular system, a LMM was fitted for each descriptor: feature ∼ Resistance + (1 | System), where Resistance (0 = sensitive, 1 = resistant) was the fixed effect of interest and System (n = 6) was included as a random intercept to capture system-level baseline differences. Models were fitted using restricted maximum likelihood (REML) via the MixedLM function in the Python statsmodels library (Seabold, 2010). For each feature, we report the standardized fixed-effect coefficient (β), its associated p-value, the intraclass correlation coefficient (ICC) quantifying the proportion of total variance attributable to system-level clustering, and the marginal R ² representing the variance explained by the Resistance fixed effect alone. P-values were corrected for multiple testing using the Benjamini-Hochberg FDR procedure across all 101 descriptors. Feature distributions were visualized using Violin plots (Waskom, 2021; Hunter, 2007), with overlaid boxes indicating the interquartile range (25th–75th percentile), median lines, and whiskers extending to 1.5× the interquartile range.

To assess whether descriptors identified as significant by both the replica-level Mann–Whitney U test and the LMM reflect a generalizable resistance phenotype rather than effects specific to the identity of the Y96 substitution, we performed within-resistant-class variability tests comparing the three resistant mutants (G12C/Y96C, G12C/Y96D, G12C/Y96S). For each feature, differences in central tendency among the three resistant mutants were tested using the Kruskal-Wallis H test (Kruskal and Wallis, 1952). Where the omnibus test was significant, pairwise post-hoc comparisons were performed using Mann–Whitney U tests with Benjamini-Hochberg FDR correction across the three pairwise comparisons, equivalent to Dunn’s procedure (Dunn, 1964). Effect size was quantified using eta-squared (η²), calculated as η² = (H − k + 1)/(n − k), where H is the Kruskal-Wallis statistic, k is the number of groups, and n is the total sample size (Tomczak and Tomczak, 2014).

4.2.3 Machine learning

Two complementary validation schemes were used to assess model performance and mitigate the risk of inflated performance estimates from conformer-level data leakage. The lenient scheme consisted of a random 80/20 split for training and testing sets. The rigorous scheme consisted of mixed held-out validation, in which all nine possible combinations of one held-out sensitive PDB system and one held-out resistant mutant system were used as test sets, with the model trained on the remaining four systems in each case.

Three ML methods were evaluated, each with default and/or optimized hyperparameter configurations, giving seven total model configurations: logistic regression (LR; L2-regularized default, L1-optimized, and Elastic Net-optimized), random forest (RF; default and optimized), and support vector machine (SVM; default and optimized). All scikit-learn estimators and utilities (StandardScaler, LogisticRegression, RandomForestClassifier, SVC, GridSearchCV, permutation_importance) were implemented in Python using scikit-learn (Pedregosa et al., 2011; Python Software Foundation, (2025)). Each configuration was evaluated over 10 independent runs per validation scheme.

All preprocessing was performed independently within each training fold to prevent data leakage. Features were standardized using StandardScaler (zero mean, unit variance), fitted on the training data only and applied to the corresponding test data. A correlation-based feature filter was then applied to the training data, retaining only features with pairwise Pearson correlation below r = 0.75. Hyperparameter optimization was performed using GridSearchCV with 5-fold cross-validation, repeated independently within each validation fold (and within each of the 10 runs) to prevent leakage from the optimization procedure itself.

For LR, L2 (Ridge) regularization was used with default parameters (C = 1.0, solver = lbfgs); L1 (Lasso) and Elastic Net regularization were optimized via grid search using the saga solver. The regularization strength C was varied over [0.001, 0.01, 0.1, 1, 10, 100] for L1, and [0.001, 0.01, 0.1, 1, 10] for Elastic Net, which additionally included tuning of the l1_ratio over [0.1, 0.3, 0.5, 0.7, 0.9]. For the SVM classifier, default parameters used C = 1.0, gamma = ‘scale’, and a radial basis function (RBF) kernel; the optimized grid searched C over [0.1, 1, 10, 100] and gamma over [0.001, 0.01, 0.1, ‘scale’], with the kernel fixed to RBF, which consistently outperformed linear and polynomial kernels in preliminary testing. For the RF model, default parameters used n_estimators = 100, max_depth = None, min_samples_split = 2, min_samples_leaf = 1, and max_features = “sqrt”; the optimized grid searched n_estimators over [50, 100, 200], max_depth over [None, 5, 10], max_features over [‘sqrt’, ‘log2’], min_samples_split over [2, 5], and min_samples_leaf over [1, 2].

Feature importance was determined using the absolute value of model coefficients for LR, Gini impurity-based feature importances for RF, and permutation importance (10 repeats, scored by ROC-AUC) for SVM, since the RBF kernel does not provide a native coefficient-based importance measure. Feature importance was computed within each validation fold using only the features retained after the within-fold correlation filter (r < 0.75) in at least 50% of fold-run iterations for a given model (≥5/10 runs for random 80/20; ≥45/90 fold-runs for mixed held-out validation), then averaged across runs (and folds, for the mixed held-out scheme).

4.2.4 Bayesian Network

To analyze probabilistic relationships among molecular descriptors, we used a Bayesian Network (BN), an interpretable ML model that represents joint probability distributions as a directed acyclic graph. A universal BN was constructed by combining treatment-sensitive and treatment-resistant KRAS datasets using BNOmics (Gogoshin et al., 2017), with 100 random restarts to ensure convergence and structural robustness. Molecular descriptors were discretized into three categories, low (0), medium (1), and high (2), using a maximum entropy algorithm (Gogoshin et al., 2017). The resulting network captured both direct associations with treatment response and higher-order dependencies among features. Two types of inference analyses were performed using a belief propagation algorithm. First, to assess how treatment response affects molecular descriptors, the treatment resistance node was fixed to either sensitive (0) or resistant (1), and the resulting posterior marginal probabilities were computed for each feature. These distributions quantify how each node’s probability shifts in response to the fixed phenotype. Second, to evaluate how selected molecular descriptors influence treatment response, each was fixed to one of its discretized states (0, 1, or 2), and the resulting marginal probability of the treatment resistance node was inferred. These probabilities were compiled into a heatmap to visualize the effect of each perturbation.

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/fchbi.2026.1900407/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Myungjin Lee, OncoC4 Inc., United States

Reviewed by: Marcin Gackowski, Nicolaus Copernicus University in Toruń, Poland

N. V. Chenchu Lakshmi.K, KLEF Deemed to be University, India