We assessed our method using two publicly available datasets: O'Neil (O'Neil et al., 2016) and NCI-ALMANAC (Holbeck et al., 2017). The O'Neil dataset comprised 23,062 drug combination samples involving 38 drugs and 39 human cancer cell lines. The NCI-ALMANAC dataset was relatively larger, containing 304,549 data points across 104 drugs and 60 cell lines. The synergy value for each sample is represented by the Loewe and combination scores for O'Neil and NCI-ALMANAC, respectively. The characteristics of the cell lines were represented by 651 gene expression values obtained from the COSMIC database (Forbes et al., 2015). Following established preprocessing steps (Liu et al., 2022), the final datasets included 18,950 and 74,139 drug-drug-cell line combinations for O'Neil and NCI-ALMANAC, respectively. Figure 1 depicts the distribution of synergy scores for both datasets. Notably, the left side of the distribution, centered around 30, constitutes more than half of the dataset. These values correspond to the negative pairs that exhibit either additive or antagonistic effects, indicating that a significant portion of the drug combinations do not show a synergistic benefit over the individual effects of the drugs. This observation underscores the complexity of identifying truly synergistic drug pairs and highlights the importance of systematic screening and computational approaches to optimize drug combination therapies.

To train and evaluate the model, we began by randomly selecting 90% of drug pairs and cell lines from each dataset to conduct three different experimental settings: random setting, cold cell line setting, and cold drug pair setting. The remaining 10% of the samples were set aside as an independent test set to evaluate generalization performance. For the random splitting setting, we divided the samples into five equal subsets. One subset served as the test set, while the remaining four were further split into training and validation sets in a 9:1 ratio. In the cold cell line setting, all the unique cell lines were divided into five equal groups randomly. The related samples which contain the cell line from one of these groups were used for testing, while the remaining samples were split into a 9:1 ratio as the training set and validation set. This ensured that the test set included only cell lines not present in the training set. For the cold drug pair setting, drug pairs were similarly partitioned into five equal groups. Four groups were used for training, with the test set containing only those drug pairs not seen during training. This ensured that the model was tested on an entirely new pair of drugs.

In this study, we formulate the synergy prediction problem as a regression task. Let X train = { d 1 i , d 2 i , c i } i ‐ 1 N denote the set of the training samples where d 1 i , d 2 i denote the drug pair and c i is the cell lines, and N denotes the number of training samples. Also, the corresponding synergy effect is represented by the label Y t r a i n = { y i } i = 1 N . The paper aims at learning a drug combination function f · , given a drug pair and a cell line, f · can generate the target value y ^ .

A variety of molecular representations have been employed for drug combination prediction tasks. Fingerprints, such as ECFP and MHFP, are commonly used to encode compound structures. In this study, we adopted the MinHashed Atom-Pair fingerprint extended to four bonds (MAP4) as our molecular representation. MAP4 offers a versatile approach to representing diverse chemical structures.

Gene expression profiles have been commonly employed to represent cell lines in drug combination prediction tasks. In this study, we utilized gene expression data extracted from COSMIC, represented as 651-dimensional vectors ( g ), where each element corresponds to the expression level of a specific gene. In the most of the deep learning-based models treat the gene expression g as a vector which does not satisfy the biomedical meaning which each gene expression should be treated separately. In the bio-mechanism of drug synergy, only a part of genes contributes to the synergy effect. So, we treat the 651-dimensional vectors ( g ) as a matrix H = { h j } j = 1 651 ∈ R 1 × 651 . In the following work, we use the CNN to extract the important genes to simulate the bio-mechanism.

The weighted gene expression representation of a cell line is fed into a cell line feature encoder to learn abstract cell line representations. This encoder comprises three convolutional layers interleaved with pooling layers. The initial convolutional layer transforms the input into feature maps, which are subsequently downsampled using max-pooling. This process is repeated three times.

The MAP4 vector representing a drug is input into a drug feature encoder to extract high-level abstract features. The encoder consists of two fully connected (FC) layers followed by Gaussian Error Linear Units (GELU) (Hendrycks and Gimpel, 2016) and batch normalization. The resulting features serve as essential inputs for subsequent fusion operations. The formulation of the drug feature encoder can be summarized as follows (Equation 1): F o i = B N G E L U F C 1024 B N G E L U F C 2048 x o i (1) where x o i is one of the input features and F o i denotes the corresponding generated feature. BN represents 1day batch normalization. F C n · represents an FC layer with n neurons. During the feature extraction stage, we project drug features and cell line feature into the same dimension to obtain higher-quality information for use in the subsequent modules.

We refer to these generated drug pair features as F d 1 , F d 2 ∈ R D and cell line feature as F c ∈ R L × D .

Given the drug pair features F d 1 , F d 2 ∈ R D and extracted feature of weight gene expression of cell line F c ∈ R L × D , we first use a global max pooling operation to obtain the global cell line feature as G c ∈ R D . We treat the drug pair features and global cell line feature as whole global features G = F d 1 , F d 2 , G c ∈ R 3 × D and the F c as the local cell line feature. The deep feature aggregation module can be decomposed into two parts: 1) global feature fusion, and 2) global to local feature fusion. Details are discussed as follows:

Global feature fusion: This process aims to integrate drug and early cell line features, followed by reinforcing the fused global features back into the local cell features. We employ a transformer encoder for global feature fusion. The core idea of the transformer encoder is the attention mechanism. An attention function maps queries ( Q ), keys ( K ), and values ( V ) to an output o as follows (Equation 2): A t t e n t i o n Q , K , V = s o f t m a x Q K T d V (2) where d is the dimensionality of the query vector.

The multi-head attention mechanism consists of multiple attention heads, with each head conducting a linear transformation on the input vectors before performing the attention operation. Each attention head has its own set of trainable parameters, allowing it to potentially model an independent relationship between the input vectors. This is achieved by utilizing different parameters in the linear transformation step.

Then, for the h head, three weight matrices W Q h , W K h , W V h ∈ R d f × d p are used to project Q, K, and V, respectively, to a lower dimension d p ; then, an attention function is performed (Equation 3). A h = A t t e n t i o n Q h , K h , V h (3) wherein Q h = Q W Q h , K h = K W K h ,   V h = V W V h .

Then, the output of the multi-head attention mechanism is the linear transformation of the concatenation of the output vectors acquired from the attention heads (Equation 4): M u l t i H e a d Q , K , V = C o n c a t A 1 , A 2 , ⋯ , A H W O (4) where H is the number of heads and W O is a trainable weight matrix.

Besides the attention mechanism, the transformer encoder also contains the residual and feed-forward neural network. Formally, the global feature fusion can be defined as follows (Equation 5): A G = LN G + M u l t i H e a d G , G , G F G = LN A G + F F N A G (5) where LN · and F F N · represent layer normalization and feed-forward neural network, respectively.

Global to local cell line feature fusion: Inspired by recent findings that drugs can influence the synergistic or antagonistic effects of drug combinations through modulating key gene expression (Wu et al., 2023), we incorporate a global-to-local cell line feature fusion network to simulate drug-induced gene regulation effects. The local cell line feature is enhanced through multi-head attention where global features ( F G )are incorporated using a Transformer decoder. This enables adaptive re-weighting of gene expressions based on cross-tissue biological patterns, with layer normalization and residual connections stabilizing feature refinement. This process can be mathematically expressed as follows (Equation 6): F G = LN F c + MultiHead F c , F c , F c C G = LN F c + MultiHead F c , F G , F G F c = LN C c + F F N C c (6)

The final synergy value of a drug combination is predicted using the output of the global feature fusion network ( F G ) and the global-to-local cell line feature fusion network ( F c ). Specifically, F G is flattened into a 1D vector, and global max pooling is applied to F c to obtain another 1D vector. These vectors are then fed into separate multi-layer fully connected layers to refine their abstract features. Finally, the refined features are concatenated and passed through a final FC layer to predict the synergy value y ^ .

Given a training dataset, that contains N samples with ground-truth synergy scores y and the corresponding values y ^ predicted by our method, we can train the deep learning model in an end-to-end fashion using the mean squared error (MSE) loss as the loss function.

We use an ensemble method to further enhance generalization and quantify the uncertainty of the CDFA. Specifically, we trained M distinct model replicas. Each replica shares the same neural network architecture and settings but uses a different initial random seed for parameter initialization. This ensures that while the models are structurally identical, they develop unique parameter values during training, leading to diverse predictions and a more robust uncertainty estimation. For every input drug combination, each model generates a predicted synergy value, denoted as y ^ k d 1 i , d 2 i , c i . The final synergy prediction, μ d 1 i , d 2 i , c i , is determined by averaging these individual predictions. Meanwhile, the uncertainty associated with this prediction, σ d 1 i , d 2 i , c i , is quantified by calculating the standard deviation of the individual predictions from the ensemble.

Calibration errors (Mervin et al., 2021) in probability estimates compromise reliability by creating discrepancies between predicted and true probabilities. Specifically, they refer to the discrepancy between the model’s predicted confidence and the actual observed frequency of correctness at that confidence level. For example, if a model assigns 80% confidence to a set of predictions, but only 70% of them are correct, this indicates a calibration error in that confidence range. Such miscalibration reduces the effectiveness of uncertainty estimates as indicators of trustworthiness in predictions.

To address this issue, a common strategy is to learn a recalibration function that adjusts the predicted uncertainties to better align with the true underlying probabilities. The recalibration function is often a non-linear uncertainty scaling function, learned using a hold-out validation dataset to create a calibration map, and is often assessed using metrics like Expected Calibration Error (ECE). In our method, we adopt a simple yet effective single-parameter scaling approach that adjusts only the uncertainty component σ · . We achieve this by multiplying σ · with a scaling factor r , while keeping the predicted synergy value μ · unchanged. This choice is motivated by the fact that μ · , as the model’s point estimate, already captures the optimal synergy prediction and should not be altered during post-hoc calibration. Instead, we rescale σ · by a positive scalar factor r , resulting in the recalibrated output μ · , r σ · . The scaling factor r is optimized using Brent’s method (Brent, 1971) to ensure that the recalibrated uncertainties accurately reflect the true probability of correctness. The objective is to minimize the miscalibration, quantified by ECE, on a separate validation set. This optimization ensures that the adjusted uncertainties more accurately reflect the true likelihood of correct predictions across confidence levels. The result is an uncertainty estimate that is better aligned with the model’s empirical behavior and more trustworthy for downstream decision-making.

CDFA is an ensemble deep learning framework for predicting the potential synergy effects of drug combinations based on the drugs’ molecular information and the cells’ gene expression. The overall architecture of CDFA is shown in Figure 2. It consists of three main components: the feature encoders for the drug pair and cell line, the feature aggregation module, and the synergy prediction module. First, MAP4 is used to represent diverse chemical structures of the paired drugs. Gene expression profiles are employed to represent cell lines in drug combination prediction tasks. Then, feature encoders are used to extract these three types of features separately. A novel feature aggregation network is involved based on the Transformer which tries to capture the intricate interactions between drug pairs and cell lines. Finally, the aggregated features are connected to another synergy prediction module. The subsequent sections of this section provide detailed evidence of the superiority of this computational framework.

To evaluate CDFA’s performance, we compared it with nine existing drug combination synergy prediction models: HypergraphSynergy, DeepSynergy, DTF, CombFM, Celebi’s method, PermuteDDS, MatchMaker, GTextSyn and MMGCSyn. We employed three common regression evaluation metrics to assess the performance of these methods: root mean squared error (RMSE), coefficient of determination ( R 2 ), and Pearson’s Correlation Coefficient (PCC).

As Table 1 shows, we compared CDFA’s performance with several models using the O'Neil dataset across three different experimental setups. In the random split scenario, where data is divided without specific constraints, the CDFA model outshone others with the lowest RMSE at 13.522, alongside the highest R 2 at 0.651 and PCC at 0.808. When tested on unseen cell lines (cold cell line setting), HypergraphSynergy led with the highest R 2 of 0.252 and the lowest RMSE of 19.537. However, CDFA maintained a competitive edge despite not leading in every metric. For the cold drug pair setting, where models predict outcomes for drug combinations not encountered during training, CDFA performed exceptionally well, achieving the lowest RMSE (15.976), highest R 2 (0.511), and a PCC of 0.717, demonstrating its strength in handling unseen drug pairs.

As shown in Table 2, the consistent superiority of CDFA has also been demonstrated on the NCI-ALMANAC dataset. In the random split setup, CDFA exhibited the best performance with the lowest RMSE of 41.893, highest R 2 of 0.552, and highest PCC of 0.746. Under the cold cell line condition, HypergraphSynergy performed best with an RMSE of 53.398, R 2 of 0.273, and PCC of 0.538. In the cold drug pair scenario, CDFA once again stood out, achieving the lowest RMSE (50.522), sub-optimal R 2 (0.346), and highest PCC (0.593), underscoring its effectiveness in predicting responses for novel drug combinations.

The 10% of the samples of the O'Neil and NCI-ALMANAC datasets were set aside as an independent test set to evaluate these models’ generalization performance. In the independent test data section of the O'Neil and NCI-ALMANAC datasets, the superior performance of CDFA has been once again proven. As Table 3 shows, it illustrates the performance of various methods when applied to the independent test datasets. For the O'Neil dataset, CDFA demonstrates superior accuracy with the lowest RMSE of 15.111 and the highest R 2 of 0.660. PermuteDDS trails closely behind with similarly strong results, showing almost no difference from CDFA. On the NCI-ALMANAC dataset, CDFA retains its leadership by achieving the best RMSE at 42.307 and the highest R 2 value at 0.508, confirming its robustness in both precision and explanatory capability. Although PermuteDDS performs well, it still lags slightly behind CDFA across all metrics. The remaining methods exhibit higher RMSE figures and lower R 2 values, suggesting they are less precise and less effective compared to our method.

Overall, CDFA consistently demonstrated strong performance, particularly excelling in the random split and cold drug pair settings. However, the poor performance of all methods in the cold cell line setting suggests that future research should focus on improving models' ability to generalize to new cell lines.

Both previous studies and our own experiments have consistently demonstrated that model performance deteriorates significantly under the cold cell-line scenario, where test cell lines are entirely disjoint from those seen during training. This setting introduces substantial biological variability, making it difficult to disentangle whether performance degradation arises from tissue-specific effects or from the challenge of generalizing to unseen cell-line profiles.

To avoid this confounding factor, we also conducted a tissue-specific analysis on the O’Neil and NCI-ALMANAC datasets. The O’Neil dataset is built on testing 38 drugs on 39 cell lines representing multiple cancer types from six tissue origins. The NCI-ALMANAC dataset covers 104 drugs in 60 cell lines from nine tissue origins. As illustrated in Figures 3, 4, our analysis employs raincloud plots to visualize the distribution of MSE for the two independent test datasets. These plots combine box plots with kernel density estimates ('clouds') to visualize both the shape and central tendency of the error distributions, with outliers indicated by diamond markers.

Our analysis reveals that although the MSE values of the median, second quartile, and third quartile are low, almost all tissues included by the two datasets have MSE values exceeding 500 and 2000, respectively. This suggests that while there is a small number of higher error values across most tissues, the central tendency of the error distribution may be relatively low. This pattern indicates that the model can achieve efficient prediction across different tissues. Our analysis confirms that the presence of high-error predictions—though limited in quantity—reveals significant variability in model performance. Such findings highlight the need for further investigation into the factors contributing to these higher errors and suggest that improvements in model accuracy and consistency are necessary for more reliable predictions across different tissue types.

Figure 5 displays the calibration curves of CDFA under various settings for both the O'Neil and NCI-ALMANAC datasets. The figures are organized from left to right, representing random splits, cold cell line settings, and cold drug pair settings, respectively. The first row showcases the O'Neil dataset, whereas the second row pertains to the NCI-ALMANAC dataset. The space between the calibration curves and the diagonal line represents the miscalibration area, which quantifies the extent of uncertainty calibration. As illustrated in Figure 5, CDFA’s recalibration algorithm successfully shifts the calibration curves closer to the diagonal line, thereby reducing the miscalibration area and improving the reliability of the predictions.

Figure 6 illustrates the relationship between prediction error and uncertainty, with uncertainty measured as the standard deviation (std). In this figure, red points indicate errors that do not fall within two standard deviations, while black and blue points represent errors that fall within one and two standard deviations, respectively. It is evident that the majority of the observed errors lie within two standard deviations, reflecting a reasonable alignment between the model’s predicted uncertainty and its actual prediction error.