Interpretable Toxicity Prediction Using Fragment–Molecular Evidence Stacking: Explicit Integration of Local Fragment and Whole-Molecule Evidence

1. Introduction

Toxicity prediction is an important task in early-stage drug discovery for evaluating the safety of candidate compounds, assessing environmental risks of chemicals, and generating hypotheses about toxicity mechanisms. Public benchmarks such as Tox21 serve as a basis for systematically evaluating machine learning models across nuclear receptor (NR) and stress-response (SR) toxicity endpoints [1,2]. In recent years, quantitative structure–activity relationship models based on molecular descriptors and graph neural networks (GNNs) that directly handle molecular graphs have been widely used for molecular property prediction, including toxicity prediction [3,4]. Chemprop, which is based on the directed message passing neural network (D-MPNN), has been widely used as a representative molecular property prediction model that constructs learned molecular representations from molecular graphs [4,5]. Graph-based representations, such as molecular graph convolution and attention-based graph models, are powerful molecular representations; nevertheless, descriptor-based models also exhibit competitive performance depending on the dataset and evaluation setting [6,7,8]. Chemprop/D-MPNN has exhibited high performance across diverse molecular property prediction tasks compared with fixed descriptors and conventional GNNs, and it remains an important reference point in current molecular GNN research. However, in message passing models, including D-MPNN, directly extracting the contributions of local structure-derived predictive evidence and whole-molecule context from high-dimensional atom/bond-level representations into a final prediction equation is difficult. Therefore, this study does not focus on competing with general-purpose molecular GNNs solely on performance; instead, it separates fragment- and molecular-level evidence to explicitly express their final integration.

However, high predictive performance alone is insufficient for toxicity prediction [9]. Accurate chemical interpretation requires distinguishing which local structures contribute to the prediction, identifying how properties such as hydrophobicity, polarity, charge, and aromaticity are involved, and determining whether whole-molecule physicochemical properties reinforce or suppress these local structural effects. Molecular descriptor models can conveniently handle global physicochemical contexts, such as molecular partition coefficient (MolLogP), topological polar surface area (TPSA), molecular weight, and ring count; however, they do not explicitly indicate which local structural units within a molecule contribute to the prediction [3,8,9,10,11]. Conversely, atom-level GNNs can learn local structural patterns; however, their internal representations are high-dimensional and cannot be easily converted into rules or equations that chemists can interpret directly [4].

Post-hoc feature attribution and fragment contribution analyses are useful auxiliary tools for model interpretation [12,13,14]. However, they are not always directly associated with the prediction equation of the model. Merely analyzing fragment contributions as a post-hoc procedure cannot demonstrate that the prediction equation uses fragment-level evidence. In recent years, methods such as subgraph attribution and graph/molecular counterfactual explanations have been proposed to improve the interpretability of molecular GNNs [15,16,17,18]. These methods are useful for identifying atoms, bonds, substructures, or minimal structural changes that reverse a prediction. However, in many cases, they function as post-hoc explanations for trained models and do not provide, as part of the model architecture, which evidence logits are integrated with which coefficients. Therefore, improving the interpretability of molecular GNNs requires not only post-prediction attribution but also a framework that separates evidence generation and integration at the model-design stage. Fragment-level evidence must be explicitly incorporated into the model input and prediction structure and handled separately from molecular-level evidence.

Fragments are chemically meaningful structural units located between atoms and whole molecules [19,20,21]. Fragment-based representations connect local structure and chemical interpretation as an intermediate level between overly fine atom-level representations and overly coarse whole-molecule descriptors. In toxicity prediction, both the presence of a given fragment and its local hydrophobicity, polarity, charge, and aromaticity may be important. Fragment-level descriptors provide chemical information at a different hierarchy from whole-molecule descriptors. Motivated by similar concerns, recent work has increasingly incorporated substructural units such as motifs and fragments into molecular representations, highlighting the importance of using chemically meaningful substructures rather than treating molecules solely as global vectors [19,20,21]. Building on this direction, this study explicitly separates local evidence obtained from fragment-level graph representations and global evidence obtained from molecular descriptors and then integrates them at the final stage via two-logit linear stacking.

In this study, molecules are decomposed into chemically meaningful fragments, and a fragment graph is constructed by treating each fragment as a node. The fragment branch takes the fragment graph as input and outputs an evidence logit derived from local structural information. In contrast, the molecular branch takes compact molecular descriptors as input and outputs an evidence logit in the whole-molecule context. The final prediction uses two-logit linear stacking, where only these two evidence logits are used as inputs [22]. This study does not claim stacked generalization or a linear meta-learner as a new method. Its novelty lies in explicitly separating toxicity prediction evidence into local structural evidence derived from fragment graphs and whole-molecule evidence derived from molecular descriptors and in using only these two logits for final integration. This design ($f_{\mathrm{frag}}{f}_{\mathrm{mol}}$) enables direct analysis of not only the final predictive performance but also the weights assigned to fragment- and molecular-level evidence for each endpoint.

$$y=affrag\mathrm{}+bfmol\mathrm{}+c$$

Figure 1. Overview of proposed fragment–molecular evidence stacking framework. Molecules are decomposed into chemically meaningful fragments and represented as fragment graphs. The fragment branch generates fragment-level evidence logit f_frag, whereas the molecular branch generates molecular-level evidence logit f_mol from compact molecular descriptors. The final toxicity logit is obtained via two-logit linear stacking, y = a f_frag + b f_mol + c.

Figure 1. Overview of proposed fragment–molecular evidence stacking framework. Molecules are decomposed into chemically meaningful fragments and represented as fragment graphs. The fragment branch generates fragment-level evidence logit f_frag, whereas the molecular branch generates molecular-level evidence logit f_mol from compact molecular descriptors. The final toxicity logit is obtained via two-logit linear stacking, y = a f_frag + b f_mol + c.

With this design, the final prediction is expressed as an explicit linear combination of fragment and molecular evidence. The primary contributions of this study are as follows. First, we constructed a framework that represents molecules as chemically meaningful fragment graphs and explicitly learns fragment-level evidence.

Second, we introduced a separate molecular-level descriptor branch and treated the global physicochemical context as evidence independent of the fragment branch.

Third, two-logit linear stacking was used to express the final integration of fragment and molecular evidence as an explicit linear equation.

Fourth, a 10-seed production evaluation under the Tox21 scaffold splits demonstrated that linear_stack (proposed model) achieved a statistically comparable receiver operating characteristic-area under the curve (ROC-AUC) and a numerically higher mean precision-recall-area under the curve (PR-AUC) than current_gate, a feature-based gated fusion model. Furthermore, linear_stack significantly outperforms fragment_only and molecular_only, which use only the fragment and molecular descriptor branches, respectively.

Fifth, the structural removal ablation of fragment-level partition coefficient (LogP)/TPSA demonstrated that local physicochemical information belonging to the same chemical axes as global MolLogP/TPSA provided auxiliary predictive information when represented as distributions on the fragment graph.

2. Results

2.1 Production 10-Seed Evaluation

In the production 10-seed evaluation under Tox21 scaffold splits, linear_stack exhibited ROC-AUC statistically comparable to that of current_gate and a numerically higher mean PR-AUC. Furthermore, linear_stack exhibited superior ROC-AUC and PR-AUC compared with fragment_only and molecular_only [23,24]. current_gate, fragment_only, and molecular_only were the comparator models.

Table 1. Overview of model variants evaluated in this study.

Model	Fragment branch	Molecular branch	Fusion / final layer	Purpose
linear_stack	Yes	Yes	Two-logit linear stacking	Proposed explicit evidence integration model
current_gate	Yes	Yes	Feature-based gated fusion	Interpretable fusion baseline
fragment_only	Yes	No	Single-branch output	Fragment-evidence baseline
molecular_only	No	Yes	Single-branch output	Molecular-descriptor baseline

Table 1. Overview of model variants evaluated in this study.

Model	Fragment branch	Molecular branch	Fusion / final layer	Purpose
linear_stack	Yes	Yes	Two-logit linear stacking	Proposed explicit evidence integration model
current_gate	Yes	Yes	Feature-based gated fusion	Interpretable fusion baseline
fragment_only	Yes	No	Single-branch output	Fragment-evidence baseline
molecular_only	No	Yes	Single-branch output	Molecular-descriptor baseline

As shown in Figure 2, linear_stack achieved an ROC-AUC comparable to that of current_gate and a numerically higher mean PR-AUC. linear_stack exhibited higher ROC-AUC and PR-AUC than fragment_only and molecular_only. The production 10-seed performance is shown in Table 2.

As shown in Table 2, in the production 10-seed evaluation using the final documented universal scaler, linear_stack achieved an ROC-AUC of 0.8017 ± 0.0142 and a PR-AUC of 0.3470 ± 0.0263. current_gate achieved an ROC-AUC of 0.8028 ± 0.0155 and a PR-AUC of 0.3353 ± 0.0235. Although the mean ROC-AUC of current_gate was slightly higher than that of linear_stack, a paired comparison demonstrated no statistically significant difference. In contrast, linear_stack exhibited a numerically higher PR-AUC than current_gate even though the difference was not statistically significant.

The primary objective of this study is to clarify ranking/screening performance and evidence integration under scaffold-separated conditions rather than construct a threshold-optimized binary classifier. Therefore, ROC-AUC and PR-AUC were treated as primary metrics, whereas MCC and balanced accuracy were treated as auxiliary threshold-dependent metrics. For MCC and balanced accuracy, linear_stack did not outperform current_gate or fragment_only. Therefore, the primary advantage of linear_stack is its ability to maintain competitive ranking performance while explicitly expressing final evidence integration rather than its superiority in threshold-specific classification performance. Compared with fragment_only and molecular_only, linear_stack significantly improved both ROC-AUC and PR-AUC. Figure 3 shows the seed-wise paired differences. Table 3 summarizes the corresponding paired statistical tests.

As shown in Table 3, the ROC-AUC difference between linear_stack and current_gate was not statistically significant (paired t-test p = 0.6427; Wilcoxon p = 0.7695). For PR-AUC, linear_stack exhibited a higher mean value than current_gate; however, the difference was not statistically significant (paired t-test p = 0.1435; Wilcoxon p = 0.1309). In contrast, linear_stack significantly improved both ROC-AUC (p = 5.6 x 10^−5) and PR-AUC (p = 0.0027) compared with fragment_only. Compared with molecular_only, linear_stack also exhibited superior performance in terms of both ROC-AUC and PR-AUC. Therefore, the primary advantage of linear_stack is its ability to maintain explicit final evidence integration while significantly outperforming single-branch baselines rather than its statistical performance superiority over current_gate. The Benjamini–Hochberg-adjusted p-values and Cohen’s dz as a paired effect size are shown in Supplementary Table S8. The interpretation of the comparisons between linear_stack and both fragment_only and molecular_only remained unchanged after applying an auxiliary multiple-comparison correction. In contrast, no adjusted significant difference was observed between linear_stack and current_gate for either ROC-AUC or PR-AUC. Thus, even after applying a multiple-comparison correction, the conclusion remains unchanged: linear_stack does not exhibit statistical performance superiority over current_gate; however, it outperforms the single-branch baselines. This analysis relies on internal validation using held-out scaffold tests within Tox21 rather than external validation on an independent dataset.

2.2. Endpoint-Wise Standardized Evidence Contribution Analysis

To clarify how the explicit two-logit stacking layer integrates the two evidence streams, endpoint-wise standardized evidence contributions from the fragment and molecular branches were calculated.

Figure 4. Endpoint-wise standardized evidence contribution spectrum. For each endpoint, standardized fragment and molecular evidence contributions were calculated as C_frag = a · SD(f_frag) and C_mol = b · SD(f_mol), where a and b represent the learned linear stacking coefficients, and SD(f_frag) and SD(f_mol) denote the seed-specific standard deviations of the fragment- and molecular-branch logits computed from held-out scaffold test molecules, respectively. The fragment evidence share was defined as |C_frag| / (|C_frag| + |C_mol|). Points indicate mean values across production seeds, and error bars indicate standard deviations. The dashed vertical line represents equal standardized contributions from fragment and molecular evidence. This analysis visualizes endpoint-dependent reliance on fragment- and molecular-level evidence rather than raw coefficient magnitude alone.

Figure 4. Endpoint-wise standardized evidence contribution spectrum. For each endpoint, standardized fragment and molecular evidence contributions were calculated as C_frag = a · SD(f_frag) and C_mol = b · SD(f_mol), where a and b represent the learned linear stacking coefficients, and SD(f_frag) and SD(f_mol) denote the seed-specific standard deviations of the fragment- and molecular-branch logits computed from held-out scaffold test molecules, respectively. The fragment evidence share was defined as |C_frag| / (|C_frag| + |C_mol|). Points indicate mean values across production seeds, and error bars indicate standard deviations. The dashed vertical line represents equal standardized contributions from fragment and molecular evidence. This analysis visualizes endpoint-dependent reliance on fragment- and molecular-level evidence rather than raw coefficient magnitude alone.

In the standardized evidence contribution analysis, the mean fragment evidence share exceeded 0.5 for all Tox21 endpoints under the branch-logit-scale standardization defined in this study. This indicates that under the present analysis conditions, the final linear stacking layer generally assigns a larger standardized contribution to fragment-branch evidence than to molecular-branch evidence. However, the degree of fragment dominance differed across endpoints. NR peroxisome proliferator-activated receptor gamma (NR-PPAR-gamma) exhibited the highest fragment evidence share, followed by three SR endpoints: heat shock response element (SR-HSE), antioxidant response element (SR-ARE), and alpha thalassemia/mental retardation syndrome X-linked blocked (SR-ATAD5). In contrast, the molecular evidence contribution was more pronounced for three NRs: androgen receptor (NR-AR), estrogen receptor (NR-ER), and aryl hydrocarbon receptor (NR-AhR). Therefore, linear_stack should be interpreted as an endpoint-wise evidence integration model that combines fragment and molecular evidence in different proportions for each endpoint rather than a fragment-only classifier. The corresponding numerical summary is shown in Supplementary Table S7. A stacked-bar representation of the same evidence shares is shown in Supplementary Figure S5.

2.3. Interpretability via Explicit Evidence Integration

Linear stacking provides an explicit final evidence integration rule. The final integration rule of linear_stack can be expressed as follows:

y = a f_frag + b f_mol + c.

Here, f_frag denotes the fragment-level evidence logit output by the fragment branch, f_mol denotes the molecular-level evidence logit output by the molecular branch, and y represents the final toxicity logit.

This equation is the actual final integration layer of the proposed model rather than a post-hoc surrogate. Therefore, with respect to final evidence integration, the integration rule is explicit. The stacking law fidelity R² = 1.0 reflects the fact that this layer is originally designed as a linear combination of two evidence logits; it does not imply that the internal mappings of the branches have been fully symbolized.

However, this result does not imply that the internal mappings of the fragment or molecular branches are completely mechanistically transparent. The interpretability claim in this study is the explicit expression of the final integration of fragment and molecular evidence rather than the complete closed-form description of the entire internal neural computation.

2.4. Target-Level Analysis

Target-level ROC-AUC analysis demonstrated that linear_stack achieved performance comparable to or higher than that of both fragment_only and current_gate for many targets; however, it was not consistently superior for all targets, indicating endpoint-dependent behavior. The details of the target-level ROC-AUC analysis are presented in Table 4. linear_stack exhibited the highest or near-highest ROC-AUC for several endpoints, including NR-AR-ligand binding domain (NR-AR-LBD), NR-AhR, NR-Aromatase, NR-PPAR-gamma, SR-ARE, SR-ATAD5, SR-HSE, SR-mitochondrial membrane potential (SR-MMP), and SR-p53. However, it was not uniformly best across all targets: current_gate outperformed linear_stack for NR-ER and NR-ER-LBD, and molecular_only slightly outperformed linear_stack for NR-AR. These results indicate that the primary advantage of linear_stack is its ability to simultaneously maintain mean performance and provide explicit evidence integration rather than achieving uniform superiority across all endpoints.

The target-level ROC-AUC values shown in Table 4 and Supplementary Figure S3 are calculated from pooled predictions for each endpoint and thus differ from the seed-wise mean aggregation shown in Table 2. Because the unit of aggregation differs, the relative ordering of linear_stack and current_gate may differ in some cases. Although this auxiliary analysis visualizes endpoint-specific trends, the primary statistical evaluation is based on production 10-seed paired comparisons.

These results indicate that the significance of linear_stack lies in its structure, which maintains average performance while explicitly expressing final evidence integration rather than achieving uniform performance superiority across all targets. In particular, because the chemical factors involved in toxicity endpoints differ by endpoint, the relative importance of fragment and molecular evidence is also expected to vary by target.

Therefore, in this study, linear_stack is positioned as a fusion framework that separates fragment- and molecular-level evidence and enables analysis of endpoint-specific contribution structures rather than as an absolutely superior predictor for all endpoints.

2.5. Complementary Local Distribution Insights Using Fragment-Level LogP/TPSA Descriptor

Global descriptors represent the average or overall physicochemical properties of a whole molecule[10,11]. In contrast, fragment-level descriptors represent the distribution of hydrophobicity and polarity across local structural units. Thus, even within the same chemical axes, fragment- and molecular-level descriptors represent information at different hierarchical levels.

To evaluate whether the fragment-level LogP/TPSA descriptors were redundant with global MolLogP/TPSA, the local LogP/TPSA-related fragment descriptors were structurally removed from the input features, the input dimensions and scaler were updated, and the model was retrained. Structural removal decreased the ROC-AUC of linear_stack compared with the full-feature condition. This result indicates that although Frag_LogP_Sum and Frag_TPSA_Sum belong to the same chemical axes as MolLogP/TPSA and are additively related quantities, treating them as locally placed physicochemical information on the fragment graph can provide auxiliary predictive information that is not fully captured by global descriptors alone. The fragment feature ablation results are presented in Supplementary Figure S2 and Table S2.

2.6. Performance of KAN-Style Heads at Branch Level

The Kolmogorov–Arnold network (KAN)-style head was evaluated as a candidate branch-level nonlinear approximator [25]. In the head-architecture ablation for the fragment branch, the KAN-style head exhibited a slightly higher ROC-AUC than the multilayer perceptron (MLP) and linear heads.

In the fragment-branch ablation, fragment_only with the KAN-style readout exhibited an ROC-AUC of 0.778, whereas the MLP and linear heads exhibited ROC-AUC values of 0.769 and 0.768, respectively. These results indicate that the KAN-style head is at least competitive as a branch-level nonlinear function approximator. However, the KAN was not required for final evidence integration. In two-logit fusion, linear_stack was the simplest and was effective in production evaluation. Therefore, the KAN is reported as a branch-level ablation in the Supplementary Information rather than being treated as a central contribution of this study.

2.7. Compressibility Variations Across Model Levels via Symbolic Distillation

Table 5 presents an overview of symbolic distillation [26,27]. Because the final stacking law is the model structure, it achieved a fidelity of R² = 1.0 for all representative targets. Although this confirms the designed explicitness of the final integration layer, the nontrivial finding is that the symbolic compressibility of the fragment- and molecular-branch surrogates is endpoint-dependent. The fidelity of the fragment-tier symbolic surrogates varied substantially across targets, and some targets could not be sufficiently recovered by low-complexity symbolic regression.

In particular, NR-AR, SR-p53, and NR-ER-LBD exhibited low fidelity for the fragment-tier symbolic surrogate, indicating that the internal mapping of the fragment branch was difficult to completely explain using a simple closed-form expression. Thus, the primary value of symbolic interpretation in this study lies in the ability to separate and analyze fragment evidence, molecular evidence, and their final integration rather than in the complete formalization of the internal mappings of the fragment branch. The branch-level symbolic surrogate is a low-complexity approximation based on a diagnostic feature pool, and some equations include auxiliary global/gate descriptors recorded for symbolic distillation. Therefore, the fragment-tier surrogate should not be interpreted as a pure mechanistic equation composed only of fragment input features.

Table 5 shows the symbolic distillation results for the stacking, fragment, and molecular tiers for representative targets. Because the stacking tier is the final two-logit linear stacking layer designed as a linear combination of two evidence logits, it achieved a fidelity of R² = 1.000 for all representative targets. This value does not indicate the symbolic compressibility of the internal branch mappings.

In contrast, the fidelity values of the fragment and molecular tiers varied substantially by target, demonstrating that the low-complexity symbolic approximability of the branch-level outputs was endpoint-dependent. This result indicates that the interpretability of the proposed framework lies in the hierarchical separation of fragment, molecular, and final evidence integration rather than in obtaining a complete closed-form representation of the internal branch mappings. In representative examples, such as NR-AhR, NR-Aromatase, and SR-MMP, the symbolic surrogates contained descriptors related to aromaticity, hydrophobicity, polarity, molecular size, and charge. Although these terms are qualitatively consistent with established chemical knowledge for the corresponding assay endpoints, individual equations should be interpreted as summaries of predictive evidence rather than direct proof of toxicity mechanisms.

3. Discussion

3.1 Interpretability of Linear Stacking over Feature-Based Gating

Linear stacking provides a clearer fusion mechanism than feature-based gating. In this study, the performance of the feature-based current_gate and two-logit linear_stack models was compared for integrating fragment- and molecular-level evidence. current_gate integrates the two streams using a gate output based on whole-molecule features.

$$g\left(x_{\mathrm{gate}}\right)y=g\left(x_{\mathrm{gate}}\right)f_{\mathrm{frag}}+\left(1-g\left(x_{\mathrm{gate}}\right)\right)f_{\mathrm{mol}}.$$

In this formulation, fragment and molecular evidence are integrated using a mixing coefficient. This structure is reasonable because it dynamically changes the contributions of local and global evidence according to whole-molecule physicochemical properties [28]. The gate output can also be approximated with high fidelity via symbolic regression, and current_gate is not an uninterpretable model.

However, from the perspective of the explicitness of final evidence integration, linear_stack is simpler. In linear_stack, the final logit is expressed as follows:

y = a f_frag + b f_mol + c

Because the model is expressed as a linear combination with fixed coefficients, the integration of fragment and molecular evidence can be interpreted directly. This depicts the model structure rather than being a post-hoc explanation. In the production 10-seed evaluation, linear_stack exhibited no statistically significant difference in ROC-AUC compared with current_gate but achieved a numerically higher mean PR-AUC although the difference was not statistically significant.

Therefore, current_gate is not treated as a failed model in this study but as an interpretable alternative fusion baseline. Based on the simplicity of its final integration rule, its numerical increase in PR-AUC, an out-of-fold (OOF)-to-test gap nearly comparable to that of fragment_only, and its consistency with symbolic interpretation, linear_stack was selected as the primary model. The objective of this study is not to replace powerful general-purpose molecular GNNs, such as D-MPNN/Chemprop [5]. Instead, this study aims to separate local structure-derived evidence from whole-molecule context-derived evidence for toxicity prediction and explicitly express the integration of both evidence into the final prediction. Thus, the significance of the proposed framework lies in its ability to separately analyze fragment-level evidence generation and final evidence integration rather than in achieving superior performance over black-box GNNs. The OOF estimates were higher than held-out scaffold test performance, and a decline was particularly observed for PR-AUC on scaffold-separated test molecules. However, this decline was not specific to linear_stack; a similar decline was observed for fragment_only. Therefore, the robust performance of linear_stack can be attributed to the effective integration of fragment- and molecular-level evidence on the held-out scaffold test set rather than to the excessive fitting of the stacking layer to OOF logits. Although this result confirms the robustness of the model within the Tox21 scaffold split, it does not directly prove generalization to an independent external compound set. current_gate and linear_stack are not a perfectly controlled comparison; they differ only in their fusion mechanisms. current_gate is an end-to-end fusion baseline using a feature-based gate, whereas linear_stack learns the final stacker using OOF logits from independently trained fragment and molecular branches. Thus, this comparison evaluates differences in the design philosophy of final integration and is not a strict isolation study of the fusion module alone. The held-out test values shown in Table S6 were aggregated from archived OOF diagnostic outputs stored to evaluate the OOF-to-test gap. Thus, these values differ in aggregation target and analytical purpose from the final production 10-seed summary shown in Table 2. Although the absolute values do not match, the qualitative conclusion is maintained: the OOF-to-test gap of linear_stack is comparable to that of fragment_only.

3.2. Fragment Evidence as Major Driver of Toxicity Prediction with Complementary Molecular-Level Information

In this study, molecular_only exhibited lower ROC-AUC and PR-AUC than the other models. In contrast, fragment_only exhibited high performance, indicating that fragment-level evidence is a major information source for Tox21 toxicity prediction.

This result is consistent with the importance of local structural information in toxicity prediction. In many toxicity endpoints, specific aromaticity, hydrophobicity, polarity, charge distribution, functional groups, or fragment connectivity contributes to prediction performance [1,9]. Global descriptors such as MolLogP and TPSA alone cannot sufficiently represent which local structures contribute to toxicity prediction.

Nevertheless, the high performance of fragment_only does not imply that molecular-level evidence is unnecessary. linear_stack significantly improved both ROC-AUC and PR-AUC compared with fragment_only, suggesting that the molecular branch complements the global context that the fragment branch alone does not fully capture.

The proposed model should be interpreted as a framework that uses fragment evidence as the primary axis and integrates molecular evidence as an auxiliary source rather than as a framework that uses molecular descriptors as the primary predictor. Standardized evidence contribution analysis supports this interpretation. Under the coefficient-weighted branch-logit standard deviation standardization used in this study, the mean fragment evidence share exceeded 0.5 for all endpoints. However, for NR-AR, NR-AhR, and NR-ER, the fragment evidence share was close to 0.5, and the molecular evidence contribution was significant; therefore, fragment dominance should not be overinterpreted for these endpoints. The molecular evidence share was larger at these endpoints than at the other endpoints. Thus, molecular-level evidence cannot be regarded as unnecessary; its contribution was endpoint-dependent and nonnegligible for multiple targets. This result supports interpreting linear_stack as an endpoint-wise evidence integration model rather than a fragment-only classifier.

3.3. Local Physicochemical Descriptor Distributions as Complementary Features to Global Molecular Descriptors

In this study, a structural removal ablation study was conducted to examine whether the fragment-level LogP/TPSA descriptors are redundant with global MolLogP/TPSA [10,11]. The removal of the local LogP/TPSA-related fragment descriptors reduced the ROC-AUC of linear_stack by approximately 0.015 compared with the full-feature condition. This result indicates that fragment-level hydrophobicity and polarity are not redundant with global MolLogP and TPSA. MolLogP and TPSA represent whole-molecule physicochemical properties, whereas Frag_LogP_Sum and Frag_TPSA_Sum are local contributions obtained by decomposing the same chemical axes at the fragment level. Therefore, global MolLogP/TPSA and fragment-level LogP/TPSA are not completely independent information sources; instead, they are additively related. This nonredundancy can be attributed to the hierarchical and spatial representation of local physicochemical properties across fragments in the fragment graph rather than to the independence of the chemical axes.

This hierarchical difference is important in toxicity prediction. For example, even when the overall molecular hydrophobicity is similar, the local structural information available to the model differs depending on whether hydrophobicity is concentrated in a particular fragment or distributed throughout the molecule. Similarly, even when the whole-molecule TPSA is similar, the information represented as fragment-level evidence can differ depending on whether polar sites are locally concentrated or distributed across multiple fragments. Thus, the fragment-level physicochemical descriptors used in this study should be considered effective local physicochemical axes rather than redundant inputs duplicating global descriptors.

3.4. KAN as Supplementary Branch-Level Module

In the initial design of this study, the KAN-style nonlinear head was considered an important component [25]. However, additional experiments demonstrated that the KAN was not required for the final two-logit evidence integration. In the final integration, linear_stack was the simplest model and was effective in the production evaluation.

In addition, the branch-level ablation study demonstrated that the KAN-style head achieved a slightly higher ROC-AUC than the MLP and linear heads. Therefore, the KAN cannot be dismissed as meaningless. The KAN-style head may be useful as a candidate nonlinear function approximator within the fragment or molecular branch.

However, these results are insufficient to establish the KAN as a central contribution of this study. These results do not demonstrate that the KAN is essential for performance improvement or that it consistently outperforms MLP in symbolic distillation. Therefore, the KAN is not placed in the title or main claims of this study; instead, it is reported in the Supplementary Information as a candidate branch-level nonlinear head.

This organization allows the main claim of the study to focus on the proposed framework that separates fragment and molecular evidence and integrates them through explicit linear stacking rather than on the KAN.

3.5. Scope of Symbolic Interpretation

In this study, PySR-based symbolic distillation was used to separately analyze the fragment-tier symbolic surrogate, molecular-tier symbolic surrogate, and final stacking law. In particular, the stacking law is the architecture of linear_stack.

$$y=a{f}_{\mathrm{frag}}+b{f}_{\mathrm{mol}}+c.$$

Consequently, the final evidence integration layer exhibits high interpretability by design. However, the explicitness of the stacking law does not imply that the internal mapping of branch-level evidence generation is also simple.

This does not imply that the entire model is fully transparent. The fragment branch uses fragment graph connectivity, edge attributes, and local descriptor interactions to generate evidence logits. Therefore, it is difficult to fully recover the internal mapping of the fragment branch using only low-complexity equations based on pooled fragment descriptors. Furthermore, the symbolic compression fidelity was target-dependent.

Therefore, the symbolic interpretation in this study is not intended to derive causal equations for toxicity mechanisms. The obtained equations approximate and organize the model output predictive evidence in a human-readable form [26,27]. Even when a specific symbolic term appears, it should not be directly interpreted as a biological causal mechanism.

The interpretability of the proposed framework should be understood as the following hierarchical separation: First, the fragment branch generates fragment-level evidence; second, the molecular branch generates molecular-level evidence; third, linear_stack integrates them through an explicit linear equation.

This hierarchical organization allows the overall prediction of the proposed model to be discussed at least in terms of evidence generation and integration. Therefore, the symbolic interpretation in this study is complementary to subgraph attribution methods, such as GNNExplainer and SubgraphX, and to molecular counterfactual explanations [14,15,17]. Existing molecular explainable AI methods focus on identifying which atoms, bonds, substructures, or structural changes are important for prediction. In contrast, the symbolic distillation in this study explicitly describes the low-dimensional rule by which the evidence logits output by the fragment and molecular branches are finally integrated. Therefore, the interpretability of the proposed model provides a means to understand the hierarchical generation and integration of predictive evidence rather than directly identifying important sites in the input structure.

3.6. Relationship to Known Toxicity Mechanisms

Although the symbolic surrogate equations and stacking coefficients obtained in this study represent predictive evidence for toxicity endpoints, they do not directly prove toxicity mechanisms[9,18,27]. However, because fragment-level descriptors include local hydrophobicity, polarity, aromaticity, and charge-related features, the obtained evidence can be connected to known chemical interpretations. For example, at endpoints in which aromaticity or hydrophobicity may be involved, fragment-level aromaticity or local LogP-related descriptors may appear as predictive evidence. Such results can serve as clues for generating subsequent validation hypotheses by comparing them with known chemical properties.

The symbolic surrogates for representative targets included features consistent with known assay biology and chemical knowledge.

The following chemical mappings are qualitative comparisons between features included in low-complexity symbolic surrogates and known chemical knowledge; they do not imply complete reconstruction of branch outputs or proof of toxicity mechanisms. In particular, the fidelity of the fragment-tier surrogates was low for some targets; these equations should be interpreted only as approximations that summarize part of the predictive evidence.

For NR-AhR, the fragment-tier surrogate contained terms such as Frag_Aromatic, Frag_LogP_Sum, and Frag_TPSA_Sum. AhR is associated with hydrophobic and aromatic xenobiotic ligands, including polycyclic aromatic hydrocarbons and dioxin-like compounds. This result is qualitatively consistent with the importance of aromaticity and hydrophobicity in AhR ligand recognition [29]. However, because the equations in this study did not measure direct binding to AhR, these terms should be interpreted as predictive evidence for the NR-AhR assay rather than as receptor-binding mechanisms.

For NR-Aromatase, the fragment- and molecular-tier surrogates contained Frag_LogP_Sum, Frag_TPSA_Sum, Frag_Aromatic, MolLogP, LabuteASA, RingCount, and related descriptors. Aromatase is a cytochrome P450 enzyme that processes steroidal substrates, and hydrophobic scaffolds, aromaticity, molecular size, and polarity are typically involved in the activity of aromatase inhibitors. Therefore, the appearance of local LogP/TPSA, aromaticity, and molecular surface area-related descriptors in the NR-Aromatase surrogate is qualitatively consistent with known aromatase ligand/inhibitor chemistry [30,31].

For SR-MMP, the fragment-tier surrogate contained Frag_Aromatic, Frag_LogP_Sum, Frag_Q_Max, and related descriptors, whereas the molecular-tier surrogate contained MolLogP, TPSA, MolWt, FractionCSP3, and related descriptors. SR-MMP is an SR endpoint related to mitochondrial membrane potential, and membrane permeability, hydrophobicity, and charge distribution can be associated with mitochondrial localization and membrane potential-dependent accumulation. Therefore, the inclusion of hydrophobicity, polarity, and charge-related descriptors in the symbolic surrogate for SR-MMP is consistent with known physicochemical factors for membrane potential-related toxicity[24,32]. However, because this study did not perform docking, experimental binding assays, or causal intervention experiments, these correspondences should be treated as predictive interpretations for hypothesis generation rather than mechanistic proof. This point is important for the practical use of toxicity prediction models. The proposed framework does not simply predict positive or negative labels for unknown compounds; instead, it organizes whether the prediction depends on fragment- or molecular-level evidence. This approach can support subsequent experimental validation and hypothesis generation for chemical structure modification.

Representative molecule-level examples linking fragment-branch predictive evidence to structural alert-like motifs are provided in Supplementary Table S9 and Figure S6 as illustrative chemical readings of model-derived evidence.

3.7. Limitations

This study has several limitations. First, the evaluation was primarily conducted on the Tox21 dataset, and future studies must examine whether the proposed framework is also effective for other toxicity datasets, different assay systems, or external compound collections. In particular, because Tox21 differs in positive rate, number of compounds, and assay characteristics across targets, the obtained trends in performance and interpretability cannot be assumed to generalize directly to all toxicity prediction tasks [1].

Second, although the fragment representation used in this study is useful for organizing predictive evidence based on substructures within molecules, it does not represent all chemical and biological factors involved in toxicity expression. For example, three-dimensional conformations, conformational changes, metabolic transformations, reactive intermediates, concentration dependence, cell-type dependence, and exposure conditions were not explicitly considered. Therefore, the fragment-derived predictive evidence obtained in this study should be interpreted as prediction evidence based on two-dimensional structural information rather than as a complete explanation of toxicity mechanisms.

Third, although the linear stacking model has the key advantage of explicitly integrating fragment- and molecular-level descriptor-derived logits, this simplification may not fully represent higher-order nonlinear interactions between them. In this study, no statistically significant difference in ROC-AUC was observed between linear stacking and nonlinear gated integration; however, although the PR-AUC of linear stacking was numerically higher, the difference was not statistically significant. These results demonstrate that for Tox21, separating and integrating fragment and molecular descriptor information can provide sufficient predictive performance and interpretability without the need for an overly complex integration mechanism.

Fourth, the symbolic approximation equations obtained in this study are intended to concisely summarize the internal predictive behavior of the proposed model and do not directly identify the causal laws of toxicity expression. The features and coefficients included in the symbolic equations represent statistical associations used by the prediction model; they are not proof of biological or chemical reaction mechanisms. Therefore, these equations should be used for visualizing predictive evidence, summarizing model behavior, and generating hypotheses for experimental validation rather than for definitively establishing toxicity mechanisms. The same applies to standardized evidence contribution analysis: it should be interpreted as a model-level decomposition of the learned two-logit stacking layer rather than as causal attribution of toxicity mechanisms. This analysis depends on trained branch logits and their seed-specific distributions and quantifies how the trained model integrates fragment- and molecular-branch evidence. Therefore, it does not prove that a particular fragment causally determines endpoint activity.

Fifth, the validity of symbolic approximation differs across targets and model outputs. In this study, the final integration component could be explicitly expressed as a simple two-logit linear equation, whereas the internal representations of the two branches, particularly the fragment branch, could not be compressed into equally simple symbolic equations for all targets. Therefore, the interpretability of the proposed framework does not imply that the entire model is fully symbolic; rather, it lies in the separation of fragment-level descriptor-derived evidence from molecular-level descriptor-derived evidence and the explicit final integration of both evidence.

Sixth, this study used a 10-seed scaffold split evaluation, leakage-free OOF stacking, and ablation analysis and adopted a design intended to reduce additional overfitting introduced by the stacking layer. A comparison of the OOF estimates with held-out scaffold test performance demonstrated performance degradation on the scaffold-separated test molecules. However, the magnitude of this degradation was not specific to linear_stack and was comparable to that observed for fragment_only. This result supports the view that under the Tox21 scaffold split condition, the stacking layer introduced only a limited possibility of additional meta-model overfitting. However, this evaluation is not an external validation based on an independent toxicity dataset. Target-specific performance differences, positive rate differences, splitting method dependence, and retraining conditions during feature removal can affect result interpretation. Future work should examine the generalizability of the proposed method using more diverse data splits, independent external datasets, and detailed target-wise analyses.

Seventh, this study did not directly compare performance with extensively validated general-purpose molecular GNNs, such as D-MPNN/Chemprop [5]. Therefore, the results do not demonstrate comprehensive performance superiority over existing powerful molecular GNNs. The primary focus of this study is to separate fragment- and molecular-level evidence and explicitly express the rule of their final integration [7,8]. Future studies should systematically evaluate the trade-off between predictive performance and interpretability, including comparisons with standard GNN baselines such as Chemprop.

4. Materials and Methods

4.1 Dataset and Prediction Tasks

This study used the Tox21 benchmark dataset to evaluate binary classification tasks for 12 toxicity endpoints[1,2]. The analysis included 7,831 compounds that had at least one target endpoint label after the removal of metal-containing compounds and the application of standardization. The endpoints included 12 targets belonging to NR and SR systems: NR-AR, NR-AR-LBD, NR-AhR, NR-Aromatase, NR-ER, NR-ER-LBD, NR-PPAR-gamma, SR-ARE, SR-ATAD5, SR-HSE, SR-MMP, and SR-p53.

To avoid overestimation due to molecular structural similarity, scaffold splitting was used for evaluation [20]. For each seed, the molecules were divided into train, validation, and test sets at a ratio of 80:10:10 based on the scaffold. The production evaluation was performed using 10 seeds, namely, Seeds 0–9. The mean and standard deviation across seeds were reported.

ROC-AUC and PR-AUC were used as the primary metrics. Because Tox21 contains class imbalance, PR-AUC was treated as a particularly important metric. MCC and balanced accuracy were also recorded as auxiliary metrics.

4.2. Molecular Fragmentation and Fragment Graph Construction

Each molecule was decomposed into chemically meaningful fragment units. The fragmenter was designed to identify structural units, such as rings, linkers, substituents, and functional caps, and to represent each fragment as a graph node. Connections between fragments were represented as edges, and each molecule was treated as a fragment graph.

The fragmenter used in this study was designed with reference to existing fragment-space construction methods, such as BRICS; however, it generated fragments at a finer granularity to obtain fragment-level evidence for toxicity prediction [21]. In BRICS, some molecules may produce too few fragments, and the entire molecule may be treated as a single fragment. In this study, protective splitting procedures were introduced to ensure that large ring systems did not remain as overly large fragments, thereby enabling interpretable fragment-level evidence.

Fragmentation diagnostics across Tox21 demonstrated zero failed molecules. The mean number of fragments was approximately 5.76, median was 5.0, and one-fragment molecule ratio was limited to approximately 6%. These results indicate that the fragmenter generates fragment graphs for Tox21 molecules stably and without excessive coarseness.

4.3. Fragment-Level Node Features

Each fragment node was assigned 16-dimensional fragment-level descriptors. These features represent local physicochemical properties and structural characteristics at the fragment level.

The fragment features included local LogP contribution, local TPSA contribution, fragment size, heteroatom count, aromaticity, formal charge, halogen count, rotatable bond count, partial charge-related descriptors, double bond count, valence-related descriptors, hydrogen-bond acceptor/donor features, and aliphatic ring features [10,11].

Among these, Frag_LogP_Sum and Frag_TPSA_Sum are not simply redundant with the molecular-level descriptors MolLogP and TPSA. In this study, MolLogP/TPSA were treated as whole-molecule global properties, whereas Frag_LogP/TPSA were treated as physicochemical axes at the local fragment level. This distinction was tested through the following structural removal ablation. Frag_LogP_Sum and Frag_TPSA_Sum are local contributions belonging to the same chemical axes as whole-molecule MolLogP and TPSA; therefore, whole-molecule values are related to the additive aggregation of fragment-level contributions. This study focuses on whether local physicochemical information distributed over the fragment graph contributes to branch-level evidence generation rather than emphasizing the global descriptor as a total quantity.

4.4. Molecular-Level Descriptors

The molecular branch received 10-dimensional descriptors representing a compact global context for the whole molecule [10,11]. The descriptors used were MolWt, MolLogP, TPSA, RingCount, Log_Num_Fragments, fragment_count, MolMR, LabuteASA, FractionCSP3, and Has3DConformer [10,11].

This molecular descriptor set was not designed to build a high-dimensional descriptor model. Instead, it was designed as a compact representation of the global context, including whole-molecule size, hydrophobicity, polarity, ring structure, surface area, three-dimensionality, and fragment complexity, which are not directly represented by the fragment graph alone.

Has3DConformer was a constant feature in the reference distribution used in this analysis, and its standard deviation was determined as 1 via zero-variance handling. Although this feature was retained for compatibility with the descriptor set, it was not a central element of predictive performance or interpretation in this study.

4.5. Fragment Evidence Branch

The fragment branch takes the fragment graph as input and outputs a fragment-level toxicity evidence logit $\left(f_{\mathrm{frag}}\right)$. Fragment node features, fragment graph connectivity, and edge attributes were used as inputs.

The role of this branch is to learn evidence related to toxicity prediction from local structural units within molecules and from connections between fragments. In this study, the output of the fragment branch was treated as a fragment-level evidence logit passed to the final stacking layer rather than as the final classification probability.

4.6. Molecular Evidence Branch

The molecular branch takes the 10-dimensional molecular descriptors as input and outputs a molecular-level toxicity evidence logit $\left(f_{\mathrm{mol}}\right)$. This branch learns complementary information contributed by the whole-molecule physicochemical context to toxicity prediction$.$The molecular branch was designed as an auxiliary evidence generator rather than an alternative to the fragment branch. Therefore, the architecture of the proposed framework primarily uses fragment evidence as the central axis while integrating molecular evidence as a global context.

4.7. Two-Logit Linear Stacking

For the final prediction, two-logit linear stacking was used, with only the two evidence logits output by the fragment and molecular branches as inputs.

y = a f_frag + b f_mol + c

Here, y represents the final toxicity logit, f_frag denotes the fragment-level evidence logit, and f_mol denotes the molecular-level evidence logit. The coefficients a and b and intercept c are learned using the stacking model$\left(f_{\mathrm{frag}}{f}_{\mathrm{mol}}\right).$This stacking law is the final integration layer of the model rather than a post-hoc surrogate. Therefore, the final evidence integration can be directly interpreted as an explicit linear equation.

4.8. Leakage-Free OOF Stacking Protocol

When training a stacking model, using in-sample predictions on the train set for the meta-model can cause leakage. To avoid this, this study used a leakage-free OOF stacking protocol [33]. For each seed, the train split was divided into three folds. In each fold, the fragment and molecular branches were trained on training folds other than the held-out fold, and outputs were generated for the held-out fold. Repeating this process for all folds produced OOF expert logits for the entire train split. The stacking model was not trained separately within each fold; instead, the OOF logits obtained from the entire train split were concatenated, and a single two-logit linear stacker was optimized. During validation and test evaluation, the fragment and molecular branches were retrained on the full train split, and expert logits were generated for the validation/test molecules. These logits were then input into the trained stacking model to obtain the final prediction. This design ($f_{\mathrm{frag}}{f}_{\mathrm{mol}})$ prevented test set or in-sample train predictions from contaminating the training of the stacking model.

4.9. Feature-Based Current Gate Baseline

In this model, fragment and molecular evidence are integrated as follows:

y = g f_frag + (1 - g) f_mol

Here, the gate module ($g(x\_gate))$ outputs a mixing coefficient that is computed from a few molecular features.

current_gate is structured to change the mixing ratio of fragment and molecular evidence according to whole-molecule physicochemical properties. Therefore, this baseline was treated as an interpretable alternative fusion baseline for comparison with linear_stack rather than as an uninterpretable black-box fusion model.

4.10. KAN-Style Branch Head Ablation

The KAN-style nonlinear head was evaluated as a nonlinear function approximation module within each branch [25]. The KAN was not used for the final two-logit evidence integration; instead, it was compared with the MLP and linear heads as a candidate readout head for the fragment and molecular branches.

The purpose of this ablation was to evaluate whether the KAN-style head is competitive with MLP or linear heads for branch-level evidence generation rather than to evaluate whether the KAN is required for the final fusion. Therefore, the KAN results are reported in the Supplementary Information as an ablation study and are not part of the central claims of the proposed model.

4.11. ChEMBL-Derived External Reference Feature Standardization

In this study, a ChEMBL-derived precomputed universal scaler was used to scale descriptors related to the molecular branch and gate. This procedure involves descriptor standardization based on an external reference distribution rather than model pretraining or model weight initialization. The universal scaler was constructed to standardize descriptors independently of the Tox21 train/validation/test split.

The final documented universal scaler was constructed from a ChEMBL-derived source file. In the reproducibility metadata, this source file was recorded as <LOCAL_CHEMBL_SOURCE_PATH>/chembl_sparse.csv.gz. The MD5 hash of this source file was b409b867ae87ab70a361e799f288bb1f.

The ChEMBL-derived source contained 244,245 compounds. The molecular structures were converted to canonical SMILES using RDKit, and duplicate canonical SMILES were removed before sampling. There were zero RDKit sanitization failures. In addition, Tox21 overlapping molecules were explicitly excluded based on canonical SMILES overlap, removing 907 molecules from the reference set. Tox21 endpoint labels and scaffold split assignments were not used for scaler estimation.

The reference distribution used 50,000 ChEMBL-derived molecules sampled with a fixed sampling seed (42). For each molecule, a fragment graph was constructed using the same fragmenter configuration as in the Tox21 experiments, and the observed fragment_count was calculated. Log_Num_Fragments was defined as ln(fragment_count + 1). The molecular descriptor vector comprised 10 descriptors: MolWt, MolLogP, TPSA, RingCount, Log_Num_Fragments, fragment_count, MolMR, LabuteASA, FractionCSP3, and Has3DConformer. The means and standard deviations of these descriptors were estimated from the ChEMBL-derived reference distribution.

In the final documented universal scaler, fragment_count and Log_Num_Fragments for the ChEMBL reference molecules were calculated using the same fragmenter as for Tox21. In the final documented scaler, the mean +/ standard deviation values of fragment_count and Log_Num_Fragments were 8.367440 +/ 6.904302 and 2.147384 +/ 0.370096, respectively. The final molecular scaler was saved as final_universal_scaler_molecular_x.joblib, and the gate scaler was saved as final_universal_scaler_gate_x.joblib. The combined scaler MD5 hash was 3cbcf7fe9f4c3d9c7bfeeb966dcb9a95_9bc96d41dda8f567adecf5f543a02778. Scaler generation was performed on 2026-06-15 at 13:17:42 on the Linux platform Linux-6.17.0-35-generic-x86_64-with-glibc2.39. The execution environment used Python 3.10.19 and RDKit 2022.09.5.

Because the ChEMBL release version and Git commit hash could not be confirmed from the available files, they were explicitly listed as unresolved reproducibility metadata in Supplementary Table S1. In addition, no additional preprocessing steps, such as salt removal/parent selection, metal exclusion, organic-only filtering, or charge/tautomer/stereochemistry handling, were implemented in the final scaler generation script; molecular processing was based on standard RDKit sanitization. The detailed metadata, descriptor means/standard deviations, and unresolved items for the final documented universal scaler are shown in Supplementary Table S1.

4.12. Symbolic Distillation Using PySR

To organize the predictive behavior of the proposed model into human-readable equations, symbolic distillation was performed using PySR. In this study, symbolic distillation was performed at three levels.

First, a fragment-tier symbolic surrogate was fitted to the output of the fragment branch. In this surrogate, pooled fragment descriptors were used as the primary candidate variables, and some auxiliary global/gate descriptors recorded for symbolic distillation were also included as candidate variables.

$$f_{\mathrm{frag}}\approx {\varphi }_{\mathrm{frag}}\left(\mathrm{FragMean} ,\mathrm{FragMax} ,{\mathrm{FragSum}}_{*}\right).$$

Second, a molecular-tier symbolic surrogate was fitted to the output of the molecular branch. In this surrogate, molecular descriptors were used as the primary candidate variables, and auxiliary descriptors recorded for analysis were also included as candidate variables.

$$f_{\mathrm{mol}}\approx {\varphi }_{\mathrm{mol}}\left(\mathrm{MolDescriptors}\right)$$

Third, the final integration law of linear_stack was confirmed.

$$y=a{f}_{\mathrm{frag}}+b{f}_{\mathrm{mol}}+c$$

Thus, the fragment- and molecular-tier symbolic surrogates are diagnostic surrogates for approximating the corresponding branch output logits at low complexity; they are not strict mechanistic reconstructions constrained only to the input features of each branch. Therefore, the variables in branch-level symbolic surrogates do not represent the pure input features of the corresponding branch and may include auxiliary variables selected from the diagnostic feature pool recorded for symbolic distillation to approximate branch output logits at low complexity.

Because the stacking law is the model structure, it can be represented with a fidelity of R² = 1.0. In contrast, the fragment-tier and molecular-tier symbolic surrogates approximate internal neural computations within the branches using low-complexity equations, and their fidelity is target-dependent.

4.13. PySR Settings and Interpretability Constraints

During symbolic regression, priority was given to obtaining chemically readable low-complexity equations. Therefore, the PySR search space was deliberately conservative [27,34]. In the primary analysis, the binary operators were addition, subtraction, multiplication, and division, and the unary operators were square and square root. An upper limit ($\times /)$was imposed on equation complexity to prevent the selection of overly complex equations.

This setting was not intended to fully recover all nonlinearities of the internal neural model. Although allowing more flexible operators or a higher complexity budget can improve fidelity, it can also reduce chemical readability. Therefore, the symbolic surrogate equations obtained in this study were interpreted as low-complexity approximations of predictive evidence.

The PySR symbolic regression settings used for stacking law distillation and branch-level distillation are summarized in Supplementary Table S4.

4.14. Fragment Feature Ablation

To evaluate whether the fragment-level LogP/TPSA descriptors are redundant with global MolLogP/TPSA, a fragment feature ablation study was conducted. In the primary analysis, Frag_LogP_Sum and Frag_TPSA_Sum were structurally removed from the fragment input vector, and the model input dimensions and scaler were updated before retraining. This structural removal evaluated a condition under which the relevant features were completely excluded from the model structure. Additional conditions, including topology-only, topology plus electronic descriptors, and no local physicochemical descriptors, were compared to assess the contributions of the local physicochemical features.

4.15. Standardized Evidence Contribution Analysis

To compare the relative contributions of fragment- and molecular-level evidence for each endpoint, the standardized evidence contributions were calculated rather than directly comparing raw stacking coefficients. For each seed and endpoint, the fragment and molecular contributions were defined as C_frag = a * SD(f_frag) and C_mol = b * SD(f_mol), respectively, where a and b represent the learned linear stacking coefficients, and SD(f_frag) and SD(f_mol) denote the corresponding seed-specific standard deviations of the branch logits. SD(f_frag) and SD(f_mol) were calculated from the branch logits for the held-out scaffold test molecules for each seed and endpoint. The fragment evidence share was defined as S_frag = |C_frag| / (|C_frag| + |C_mol|), and the molecular evidence share was defined analogously. This analysis was used only to quantify the relative contributions of the two evidence logits and was not used for model training.

4.16. Statistical Analysis

Using the seed-wise results obtained from the production 10-seed evaluation, paired comparisons were performed between linear_stack and each comparator model, including current_gate, fragment_only, and molecular_only.

For each comparison, the seed-wise paired differences in ROC-AUC and PR-AUC were calculated, and the mean differences were reported as Delta ROC-AUC and Delta PR-AUC. In statistical testing, a two-sided paired t-test was used as the primary test, whereas a Wilcoxon signed-rank test was used as an auxiliary test[35,36,37,38,39]. The number of seeds in which linear_stack outperformed each comparator was also recorded as wins. A comparison of linear_stack and current_gate demonstrated that no significant difference was observed in ROC-AUC, whereas linear_stack exhibited numerically higher PR-AUC even though the difference was not statistically significant. Therefore, the superiority of linear_stack over current_gate is limited to its simpler and more explicit final evidence integration while maintaining comparable discrimination performance.

Considering the class imbalance in Tox21, PR-AUC and ROC-AUC were treated as the primary metrics. Threshold-dependent metrics, such as MCC and balanced accuracy, were recorded as auxiliary metrics, and the selection of the primary model was based on ROC-AUC and PR-AUC. Because multiple paired comparisons were performed, the Benjamini–Hochberg false discovery rate correction was applied as auxiliary analysis. Cohen’s dz was calculated as the effect size for paired differences. The Benjamini–Hochberg correction was applied to the six primary paired t-test p-values for ROC-AUC and PR-AUC. These adjusted p-values and Cohen’s dz are shown in Supplementary Table S8. These corrections and effect sizes guided the significance interpretations, ensuring that the primary conclusions were not solely dependent on unadjusted p-values.

5. Conclusions

This study proposes a fragment–molecular evidence stacking framework that separately learns a fragment-level graph model and a molecular-level descriptor model and then integrates their evidence logits via two-logit linear stacking. In the proposed framework, the fragment branch outputs an evidence logit derived from local structural information; the molecular branch outputs an evidence logit derived from whole-molecule descriptors.

The final prediction is explicitly expressed as y = a f_frag + b f_mol + c, where f_frag and f_mol denote the fragment- and molecular-branch evidence logits, respectively.

In the production 10-seed evaluation under Tox21 scaffold splits, linear_stack exhibited no statistically significant difference in ROC-AUC compared with current_gate; however, it exhibited a numerically higher mean PR-AUC even though the difference was not statistically significant. In contrast, linear_stack significantly improved both ROC-AUC and PR-AUC compared with fragment_only and molecular_only. These results indicate that although linear_stack does not consistently outperform complex feature-based gating, it maintains superior predictive performance over single-branch baselines while explicitly separating and integrating fragment- and molecular-level evidence.

Furthermore, standardized evidence contribution analysis demonstrated that under the coefficient-weighted branch-logit standard deviation standardization defined in this study, the mean fragment evidence share exceeded 0.5 for all endpoints. However, the magnitude differed by endpoint, and the molecular evidence contribution was significant at endpoints with shares close to 0.5. This result supports the interpretation of linear_stack as a framework that integrates fragment- and molecular-level evidence in different proportions across endpoints rather than as a simple fragment-only classifier.

The performance decreased after the structural removal of fragment-level LogP/TPSA descriptors, indicating that local hydrophobicity and polarity belonging to the same chemical axes as MolLogP/TPSA provide auxiliary predictive data not fully captured by global descriptors alone when handled as local distributions on the fragment graph. This result supports the importance of modeling not only the presence of fragments but also the local physicochemical properties carried by those fragments.

In symbolic distillation, the final stacking law was explicitly represented as the model structure. The obtained fidelity of R² = 1.0 reflects the fact that the final integration layer is linear by design. In contrast, the symbolic compression of the fragment and molecular branches was target-dependent, and not all internal mappings could be explained by low-complexity equations. Therefore, the interpretability of the proposed framework lies in hierarchically separating and analyzing fragment evidence, molecular evidence, and final evidence integration rather than in converting the entire model into a complete mechanistic equation.

Although the KAN-style head exhibited potential usefulness as a branch-level nonlinear approximator, it was not required for the final evidence integration. Thus, the primary contribution of this study is the proposed framework that separates fragment- and molecular-level evidence and integrates them via explicit two-logit linear stacking rather than the KAN.

Overall, these results demonstrate that fragment–molecular evidence stacking is a useful framework for toxicity endpoint prediction in Tox21 because it maintains competitive predictive performance while enabling explicit analysis of the final integration of fragment- and molecular-level evidence. Future work should examine the generalizability of the proposed framework through validation on external toxicity datasets, stability assessment of endpoint-wise evidence contributions, and reproducibility analysis of symbolic surrogate equations.