MultiChem: predicting chemical properties using multi-view graph attention network

Abstract

Background

Understanding the molecular properties of chemical compounds is essential for identifying potential candidates or ensuring safety in drug discovery. However, exploring the vast chemical space is time-consuming and costly, necessitating the development of time-efficient and cost-effective computational methods. Recent advances in deep learning approaches have offered deeper insights into molecular structures. Leveraging this progress, we developed a novel multi-view learning model.

Results

We introduce a graph-integrated model that captures both local and global structural features of chemical compounds. In our model, graph attention layers are employed to effectively capture essential local structures by jointly considering atom and bond features, while multi-head attention layers extract important global features. We evaluated our model on nine MoleculeNet datasets, encompassing both classification and regression tasks, and compared its performance with state-of-the-art methods. Our model achieved an average area under the receiver operating characteristic (AUROC) of 0.822 and a root mean squared error (RMSE) of 1.133, representing a 3% improvement in AUROC and a 7% improvement in RMSE over state-of-the-art models in extensive seed testing.

Conclusion

MultiChem highlights the importance of integrating both local and global structural information in predicting molecular properties, while also assessing the stability of the models across multiple datasets using various random seed values.

Implementation

The codes are available at https://github.com/DMnBI/MultiChem.

Introduction

Accurate characterization of the molecular properties of chemical compounds is crucial for identifying lead candidates and evaluating cellular and tissue-level interactions in drug development. With over 100 million molecules available, experimentally screening desirable compounds within this vast chemical space is challenging, time-consuming, and costly. In this context, in silico evaluation provides a more efficient alternative, contingent on the development of accurate predictive models.

In recent years, machine learning methods have become widely adopted in drug discovery and molecular property prediction [43, 66]. These methods, grounded in quantitative structure–activity relationship (QSAR) principles, predict molecular activity based on structural similarities by carefully selecting relevant variables and inference functions [21]. Traditional machine learning models, such as support vector machine [7] and random forests [2], have demonstrated strong performance in estimating relationships between molecular structures and their properties.

In QSAR modeling, molecular fingerprints have become a standard for representing molecular structures [43, 66]. These fingerprints capture features such as electronic, geometric, and steric properties. Commonly used fingerprints, like PubChem fingerprint [30], MACCS Key [8], and extended-connectivity fingerprints (ECFP) [49], can be generated efficiently using tools like PaDel [65], RDKit [32], and CDK [53]. However, these approaches depend on predefined substructures, which may overlook other relevant molecular features [26], highlighting the need for more comprehensive molecular representations.

Deep learning has recently gained prominence across various fields due to its ability to automatically extract features without manual intervention [13]. In drug discovery, graph-based and sequence-based neural networks have shown particular promise by capturing a broader range of substructures than traditional fingerprints [27, 28, 37, 40, 44, 45, 56, 59, 67, 70]. Graph-based models are adept at distinguishing locally positioned substructures, while sequence-based models capture global patterns from distant substructures.

Several graph neural networks (GNNs) have been introduced for molecular property prediction, including graph convolutional network (GCN) [31], message-passing neural network (MPNN) [19], and graph attention network (GAT) [3, 57]. These models focus on individual nodes within molecular graphs to accurately capture local structural features. Studies have shown that GCNs and MPNNs perform comparably to traditional machine learning approaches [6, 9, 19, 29, 31, 39, 51, 55]. GATs and their variants, incorporating attention mechanisms, have shown better performance by assigning varying influence probabilities within the graph structure [34, 62, 68].

Building on the success of atom-centered GNNs, several variant models have been developed to enhance prediction accuracy [15, 36, 38, 54, 64, 69]. Among these, the directional message-passing neural network (DMPNN) [63] is notable for its focus on bond information, treating atoms and bonds are equally important to better differentiate non-isomorphic structures. Unlike atom-centered GNNs, the bond-centered DMPNN utilizes a line graph where the atoms and bonds are assigned to edges and nodes, respectively [18, 60, 63]. Inspired by the success of both atom- and bond-centered GNNs, several studies have attempted to combine atom and bond representations in a single framework [5, 15, 40, 52, 59]. In these studies, either a bond-centered GNN was added to an atom-centered GNN, or both GNNs were used together to exchange information. However, these approaches perform simple summation or averaging to represent atoms and bonds rather than employing more sophisticated attention, potentially limiting their effectiveness.

Sequence-based models have also been employed in molecular property prediction [20, 25, 50, 58]. These models primarily use SMILES sequences to represent molecular structures, enabling them to extract global structural features from entire sequences. Notably, BERT has been adapted for molecular property prediction using fully connected graph approaches [58]. Similar to GNNs, sequence-based models often outperform traditional methods.

Data scarcity is a common challenge in molecular property prediction, leading to adoption of multi-task learning and transfer learning to address this issue. Multi-task learning improves performance by sharing information across related tasks [43, 61]. Transfer learning, widely used in both graph- and sequence-based models, enables the leveraging of knowledge from larger datasets to enhance predictive accuracy on smaller datasets [16, 35, 50, 58].

In this work, we developed a model that generates a comprehensive molecular representation by simultaneously capturing local and global structural information. Our approach incorporates two graph encoders to effectively capture local information from atoms and bonds. These local features are jointly aggregated using a graph attention mechanism, which represents a key contribution of our study. Unlike traditional aggregation methods [34, 57, 62], our mechanism more effectively captures structural information from atoms and bonds, making it particularly well-suited for molecular property prediction. To complement the local features, we applied a multi-head attention mechanism to capture global relationships between distant substructures across diverse molecular topologies. By directly applying multi-head attention to the local features from our graph encoders, our model emphasizes both substructural relationships while simultaneously extracting global relationships. This multi-view approach integrates local and global information, enabling a more refined representation of molecular structures. When we assessed the impact of individual components within our model, we observed AUROC improvements ranging from 1.0% to 5.1% when comparing the model to the other models lacking specific components. In addition, the final model with ensemble method showed a 2.3% increase in AUROC score.

In our experiments, several models showed substantial performance differences between the different datasets. These inconsistencies may be due to structural differences in the scaffold distributions of the datasets, suggesting that each model focused on different molecular regions during training. To address the challenges of data scarcity and imbalance, we also incorporated an ensemble method, which has also shown promising results in the biological fields [12, 14]. Ensemble methods reduce variance and bias, mitigating risks of overfitting. In particular, bagging in our method promotes model diversity by training on data subsets with varied distributions, improving generalization and stability by reducing variance across individual models.

We evaluated our method on nine benchmark datasets from the MoleculeNet, including both classification and regression tasks. Compared to state-of-the-art methods [16, 34, 35, 40, 50] on golden standard datasets, our model achieved better or comparable performances, with an average area under the receiver operating characteristic (AUROC) of 0.822 and a root mean squared error (RMSE) of 1.133, reflecting improvements of 3% in AUROC and 7% in RMSE over current leading methods.

Materials and methods

Datasets

We used classification and regression datasets from MoleculeNet [61] using the DeepChem [47] library. The classification datasets include six tasks of BBBP [42], Tox21 [23], ToxCast [48], SIDER [47], ClinTox [61], and BACE. The regression datasets include three tasks of ESOL, FreeSolv, and Lipophilicity regarding physical chemistry.

The BBBP dataset provides information on molecular penetration of the blood–brain barrier [42]. Tox21 includes qualitative screening data for nuclear receptor (NR) and stress response (SR) assays [23]. ToxCast comprises extensive in vitro high-throughput screening data [48]. SIDER contains adverse drug reaction data classified by MedDRA [47], and ClinTox includes FDA-approved drugs alongside those failed due to toxicity issues [61]. The BACE dataset provides biophysical properties related to binding for inhibitors of human beta-secretase 1 (BACE-1). The ESOL dataset contains 1,128 compounds with water solubility. FreeSolv provides compounds with experimental hydration-free energy in water. The Lipophilicity dataset contains experimental octanol/water distribution coefficient values. In classification tasks, each task classifies compounds as active, inactive, and inconclusive. For this study, we focus exclusively on the active and inactive outcomes to maintain accuracy and reliability. Detailed configurations of each dataset are provided in Table 1.

Table 1 Details of the datasets

Feature initialization

Our model utilizes two graph forms: atom graph and bond graph (Fig. 1A), designed for the atom-centered and bond-centered sub-models, respectively. The atom graph represents atoms as nodes and bonds as edges, with the adjacency matrix indicating the bond presence. The bond graph assigns bonds as nodes and atoms as edges, with the adjacency matrix indicating direct connections between two bonds via an atom. This dual-input approach ensures our model equally considers both atoms and bonds.

Fig. 1

MultiChem: a practical tool for molecular property prediction. A Inputs. We employed two molecular graphs for our model. One is the Atom graph, which assigns atoms and bonds as nodes and edges. The other is the Bond graph, which allocates bonds and atoms as nodes and edges. B GNN for atom. We modified the graph attention network to incorporate the bond information from the other GNN and used it for the Atom graph. C GNN for bond. Similarly, we adapted the graph attention network to reflect the atom information from the GNN for the atom and exploited it to the Bond graph. D Attention. We adopted a self-attention mechanism between the nodes in the molecule to capture the long-distance features

To initialize features, we represented the molecule with three matrices (node vectors, edge vectors, and adjacency matrix), a common method in deep learning for drug discovery [63]. Node vectors are based on atomic properties such as atom type, the number of bonds, and mass. Edge vectors include bond properties like bond type, ring presence, and conjugation. Properties are transformed into feature vectors using one-hot encoding or scaling, resulting in a node vector of length 127 and an edge vector of length 12. For example, atom type and mass are converted to a vector of length 100 by one-hot encoding and a vector of length one by scaling, respectively. Bond type is transformed into a vector of length four by one-hot encoding. Feature configurations are detailed in Table 2, with properties extracted from the SMILES string using RDKit [32].

Table 2 Feature initialization

Model construction

In the MultiChem model, we incorporate three key components: a graph attention mechanism to integrate atom and bond information, a multi-head attention mechanism to capture both local and global features, and an ensemble method to enhance model stability. Since molecules are composed of atoms connected by covalent bonds, accurately predicting molecular properties requires models that account for both atomic and bonding information. Recent studies [5, 15, 40, 52, 59] have shown that models integrating both atomic and bond information provide more accurate predictions and richer representations compared to those focusing primarily on atoms and treating bonds merely as connections.

In our approach, we incorporate distinct graphs for atoms and bonds and iteratively updating both through a process of mutual reinforcement, where nodes are updated based on edge information and edges are updated based on node information. Further, we advance the integration of atom and bond information by employing a graph attention mechanism, which is known for its flexibility in learning structural information by assigning different weights to neighboring nodes, thus identifying crucial substructures. Our method uses two graph attention encoders—one for atoms and one for bonds–that are interconnected by sharing certain hidden states of atoms and bonds, updated using attention weights derived from neighboring atom and bond hidden states.

After extracting local features with the graph attention encoders, we apply a multi-head attention mechanism to the combined features, aiming to improve predictive power. Multi-head attention is recognized for its reliability and stability compared to single attention mechanisms. While GNNs are typically limited to a predefined number of hops, multi-head attention allows for capturing useful information between distant nodes. By applying multi-head attention to the outputs from our local graph-based module, the model can extract global features from distant parts of the graph, complementing the local features provided by the graph encoders. Further, to integrate local and global features effectively, we used a residual connection between the local and global modules, which helps prevent overfitting in the multi-head attention layers and preserves earlier local information.

To further improve robustness and performance, we adopted bagging, an ensemble method, during model training and inference. Ensemble methods like bagging reduce variance and bias, mitigating risks of overfitting and underfitting. In particular, bagging in our method promotes model diversity by training on data subsets with varied distributions, improving generalization and stability by reducing variance across individual models. Bagging is well-suited for our method because we can generate diverse training data by employing balanced scaffold splitting method, as described in later section. For final inference in both classification and regression tasks, we used soft voting, which generally produces more reliable results than hard voting by averaging the predictions from multiple models.

This combined approach—graph attention for interactive atom-bond representation, multi-head attention for global feature extraction, and ensemble bagging for robustness–provides a comprehensive framework for accurate molecular property prediction.

Sub-model for the node

The node sub-model (Fig. 1B) is based on the graph attention network [2, 57], applying attention to connected nodes in the message-passing layers of the message-passing neural network. Message-passing involves exchanging information (i.e., message) between two connected nodes, and we incorporated the attention method into this process. The initial hidden vector of the ith node as ${h}_{i}^{0}$, is derived from the initial node vector ${x}_{i}$ through the function ${f}_{init}$ (Eq. 1), which includes learning parameters, a bias, and the exponential linear unit (ELU) activation function.

The message vector ${m}_{j,i}^{l}$ in the message-passing layer is derived from the node vectors ${h}_{j}^{l-1}$, ${h}_{i}^{l-1}$, and the edge vector ${h}_{ji}^{l-1}$ from the previous layer using the function ${f}_{message}$ (Eq. 2). The hidden vector ${h}_{ji}^{l}$ is derived from the directional edge from the jth node to the ith node in the lth hidden layer in the edge sub-model, described in the next section. The attention coefficient ${\alpha }_{j,i}^{l}$ (Eq. 3) is calculated by applying the softmax function to the message vectors ${m}_{j,i}^{l}$, with j in the range of the neighbors of the ith node (N(i)). The node update function (Eq. 4) is formulated to obtain the node vector ${h}_{i}^{l}$, considering the edge vector ${h}_{ji}^{l-1}$, the attention score ${\alpha }_{j,i}^{l}$, and the node vector ${h}_{j}^{l-1}$. Finally, multi-head attention with k heads, a more stable method than single attention, is applied to the node update function (Eq. 5).

$$h_i^0=f_{init}\left(x_i\right)$$

(1)

$${m}_{j,i}^{l} = {f}_{message}\left({h}_{j}^{l-1},{h}_{i}^{l-1},{h}_{ji}^{l-1}\right)$$

(2)

$${\alpha }_{j,i}^{l} = {Softmax}_{j \in N\left(i\right)}\left({m}_{j,i}^{l}\right)$$

(3)

$$h_i^l=\sum\nolimits_{j\in N\left(i\right)}\alpha_{j,i}^l{\ast f_{update}(h}_j^{l-1},h_{ji}^{l-1})$$

(4)

$${h}_{i}^{l} = \frac{1}{K}\sum_{k = 1}^{K}{\sum }_{j \in N\left(i\right)}{\alpha }_{j,i}^{l,k}{*{f}_{update}^{k}(h}_{j}^{l-1,k},{h}_{ji}^{l-1,k})$$

(5)

Sub-model for the edge

The edge sub-model (Fig. 1C), also based on the graph attention network, was constructed similarly to the node sub-model. We applied the attention mechanism to the directional message-passing neural network [63] to incorporate hidden states and messages for the edges, creating a novel model. By focusing on bonds, our model captures features distinct from those obtained by the node model, making it valuable for bond-centric characteristics. This initial hidden vector of the bond between the ith and jth nodes, ${h}_{ij}^{0}$, is derived from the initial edge vector ${e}_{ij}$, node vector ${x}_{i}$, and node vector ${x}_{j}$ through the function ${f}_{init}$ (Eq. 6). The hidden vector of the ith node in the lth hidden layer, ${h}_{i}^{l}$, is inherited from the node sub-model.

With the initial states defined, we constructed the novel message function ${f}_{message}$. The message ${m}_{ki,ij}^{l}$ is calculated using the edge ${h}_{ki}^{l-1}$, the edge ${h}_{ij}^{l-1}$, and the node ${h}_{i}^{l-1}$ (Eq. 7). Here, ki represents the edge from the kth node to the ith node, ij represents the edge from the ith node to the jth node, and l denotes the layer depth. These messages derive the attention score ${\alpha }_{ki,ij}^{l}$ through a softmax function (Eq. 8). This function operates among the edges (E(ij)), directly connected to the edge ij. The aggregating function (Eq. 9) leverages attention scores, edges, and nodes. Finally, we introduced the multi-head attention with n heads to this sub-model (Eq. 10).

$${h}_{ij}^{0} = {f}_{init}\left({x}_{i},{x}_{j},{e}_{ij}\right)$$

(6)

$${m}_{ki,ij}^{l} = {f}_{message}\left({h}_{ki}^{l-1},{h}_{ij}^{l-1},{h}_{i}^{l-1}\right)$$

(7)

$${\alpha }_{ki,ij}^{l} = {Softmax}_{ki \in E(ij)}\left({m}_{ki,ij}^{l}\right)$$

(8)

$${h}_{ij}^{l} = {\sum }_{ki \in E(ij)}{\alpha }_{ki,ij}^{l}{*{f}_{update}(h}_{ki}^{l-1},{h}_{i}^{l-1})$$

(9)

$$h_{ij}^l=\frac1N\sum_{n=1}^N\sum\nolimits_{ki\in E(ij)}\alpha_{ki,ij}^{l,n}{\ast f_{update}^n(h}_{ki}^{l-1,n},h_i^{l-1,n})$$

(10)

Multi-head attention network

We built our model on the graph neural network that extracts the features focusing on the nodes and edges, typically capturing local structural features. To incorporate local and global crucial features, we adopted a multi-head attention layer known for capturing distant nodes. Before applying this method, we merged the hidden states from the node and edge sub-models. As shown in Eq. 11, we aggregated the edges connected to the ith node and the node ${h}_{i}^{L}$ in the last (Lth) layer of the sub-models, obtaining a new initial hidden state ${h}_{i}^{t0}$ using the function ${f}_{update}$ (Fig. 1D).

Once the new hidden states containing the local structural features were initialized, we applied the multi-head attention method to our model. The technique allowed us to determine the contributions of these local features by using the softmax function to calculate the importance of each sub-structure. In Eq. 12, the attention score ${\alpha }_{i,j}^{t}$ representing the importance of the jth node relative to the ith node, is derived using the softmax function on the hidden states of nodes ${h}_{i}^{t-1}$ and ${h}_{j}^{t-1}$, where t indicates the layer depth. The hidden state ${h}_{i}^{t}$ is then systematically calculated by summing all the nodes in the graph (N(G)) with their attention scores (Eq. 13). Finally, we applied the multi-head attention mechanism with k heads to our model (Eq. 14).

$${h}_{i}^{t0} = {f}_{update}({h}_{i}^{L},\sum_{j \in N(i)}{h}_{ji}^{L})$$

(11)

$${\alpha }_{i,j}^{t} = {Softmax}_{j \in N(G)}\left({f}_{query}({h}_{i}^{t-1}), {{f}_{key}({h}_{j}^{t-1})}^{T}\right)$$

(12)

$${h}_{i}^{t} = {\sum }_{j \in N(G)}{\alpha }_{i,j}^{t}*{f}_{value}({h}_{j}^{t-1})$$

(13)

$${h}_{i}^{t} = \frac{1}{K}\sum_{k = 1}^{K}{\sum }_{j \in N(G)}{\alpha }_{i,j}^{t,k}{*{f}_{update}^{k}(h}_{j}^{t-1,k})$$

(14)

Readout phase

Following the multi-head attention layers, we adopted the readout method in our model to extract the hidden feature $H$ of the molecule. For all nodes in the graph (N(G)), ${h}_{i}^{T}$ from the last (= Tth) multi-head attention layer was averaged (Eq. 15). Finally, the feed-forward neural network was applied to the output of the readout phase, yielding the predicted values.

$$H = \frac{1}{|N\left(G\right)|}{\sum }_{i \in N(G)}{h}_{i}^{T}$$

(15)

Multi-task learning and ensemble method

We employed binary cross-entropy (BCE) and mean squared error (MSE) loss functions for classification and regression tasks, respectively. For classification tasks, we adopted a multi-task learning approach, training a model on multiple tasks simultaneously [4]. This approach leverages larger datasets combined from multiple tasks, which is advantageous when individual tasks have limited data. Multi-task learning also acts as a regularization effect, reducing the risk of overfitting to any specific task. For instance, on the Tox21 dataset, consisting of 12 tasks in NR and SR, multi-task learning yielded better performance than single-task learning. We further enhanced our model through bagging, an ensemble method that improves generalization and mitigates overfitting. In bagging, multiple training datasets are created by sampling subsets with varying distributions, addressing data limitations and providing a diverse training base. Independent models are then trained on these samples, and the final prediction is made by combining the outputs from all trained models. For this ensemble inference, we used soft voting to calculate the final predictions. Unlike hard voting, which selects the prediction favored by the majority of models, soft voting averages the outputs from each model to produce the final prediction. As shown in Eq. 16, the result of sample i is derived by averaging the outputs across multiple models ${M}_{k}$.

$${\widehat{y}}_{i}^{t} = \frac{1}{K}{\sum }_{k = 1}^{K}{M}_{k}({x}_{i}^{t})$$

(16)

Training and evaluation data

Similar to previous studies, we employed the balanced scaffold splitting method for evaluation, which introduces a more rigorous and scientifically meaningful challenge compared to random or scaffold splitting [35, 50]. The balanced scaffold splitting method divides molecules into disjoint scaffold sets based on their structural scaffolds. These scaffold sets are then categorized into large and small sets according to a specified size threshold (e.g., half of the size of a test dataset). To ensure the representativeness of the training set and reduce bias in the validation and test datasets, large scaffold sets are consistently allocated to the training dataset.

To train and evaluate our ensemble model, we initially split a dataset using balanced scaffold splitting with a random seed, allocating 10% for testing and 90% for subsequent processing. We then applied a new random seed, referred to as the e-seed, distinct from the existing seed, to split the remaining data. Once the e-seed was determined, we further divided this data into 80% for training and 10% for validation. This approach enabled us to obtain multiple training-validation dataset pairs using different e-seed values, each pair containing slightly different scaffolds in the training dataset. This variability enabled us to train multiple models independently, supporting our ensemble model.

Model training

We trained our ensemble model using datasets created with balanced scaffold splitting. Before training, graphs and line graphs were extracted from the molecules in the datasets using RDKit [32]. Our model was implemented using Pytorch [1], PyTorch Lightning [11], PyTorch Geometric [17], and TorchVision [41]. PyTorch Geometric and TorchVision were specifically used to implement the graph and multi-head attention layers.

We trained our model by minimizing BCE and MSE loss functions for classification and regression tasks, respectively, on the preprocessed datasets. To prevent overfitting and improve model efficiency, we applied an early stopping based on the minimum validation loss [33, 46]. We trained our model on a single NVIDIA Tesla V100-DGXS-32 GB GPU and a CUDA 11.3 environment.

After training, we obtained multiple models from different e-seed splits, allowing us to make predictions on the test dataset using an ensemble approach with soft voting. In this method, the final prediction for the test set was obtained by averaging the outputs from each model. Final evaluation metrics, including AUROC for classification tasks and RMSE for regression tasks, were computed using Scikit-Learn [10].

During hyperparameter tuning, we applied early stopping method to prevent overfitting, with a maximum epoch of 1000 and patience of 30 epochs. Patience represents the number of epochs that training continues after a new minimum loss value is found. We optimized our model on validation datasets with various hyperparameters: batch size, layer size, layer depth, learning rate, weight decay, and dropout rate. Detailed hyperparameter configurations are provided in Supplementary Table 1.

We experimented with different batch sizes to assess their effects on overfitting, underfitting, and generalization. The optimal batch sizes were found to be 64 for smaller datasets (e.g., BBBP) and 256 for larger datasets (e.g., ToxCast). We also tested different sizes and depths for the graph and attention layers, which significantly affected model performance. The best performance was achieved with a graph layer depth of three, an attention layer depth of one, and a layer size of 128 for both layers.

For optimization, we adopted the Adam optimizer, known for its stability and speed. We adjusted two key parameters–learning rate and weight decay. Weight decay, similar to L2 regularization, was effective in preventing overfitting. The optimal model configuration was found with a learning rate of 1e-3 and a weight decay of zero. Additionally, dropout was applied to further mitigate overfitting, with the highest performance observed at a dropout rate of 0.3.

Methods for performance comparison

We conducted a comprehensive comparison of our model with various existing methods of different approaches, including graph-based [19, 29, 31, 34, 39, 40, 51, 62, 63], SMILES-based [58], and pre-training methods [16, 22, 25, 35, 50, 58], using the MoleculeNet dataset. This thorough evaluation ensured the fairness and reliability of our model.

Firstly, we compared our model with nine graph-based methods built on the message-passing neural network framework, including AttentiveFP [62], DMPNN [63], TrimNet [34], and CD-MVGNN [40]. AttentiveFP integrates a small graph into a large one using a graph attention mechanism to extract the critical substructures. DMPNN treats bonds as nodes to solve the tottering problem, while TrimNet exploits both bond and atom features to enhance the representation of molecular features. CD-MVGNN integrates atom- and bond-based GNNs with disagreement loss, mitigating the difference in predictions between these GNNs.

Secondly, we evaluated our model against the SMILES-BERT model [58], which uses the BERT architecture with SMILES string as input to represent chemical structures.

Lastly, we compared our model with six pre-training methods, such as GROVER [50], MPG [35], and KANO [16]. These methods were pre-trained on a large datasets like ZINC [24] and fine-tuned on small datasets like MoleculeNet. GROVER employs multi-level pre-tasks predicting masked information and motifs, MPG used self-supervised learning with the tasks of predicting the masked atom information and pairwise half-graph discrimination. KANO introduces contrastive learning between the augmented molecular graph for the knowledge graph of elements and the enhanced molecular graph for functional groups.

Results

We evaluated the accuracy of our model by comparing it to state-of-the-art models. To assess the stability and generalization of the models, we tested our model and several models across different data characteristics using various random seed values. Additionally, we analyzed the effectiveness of our ensemble method in handling data scarcity. An ablation study provided insight into the contribution of each model component, offering a comprehensive understanding of their roles and impact.

Performance evaluation on MoleculeNet dataset

We evaluated our model with eight state-of-the-art graph-based models and six pre-trained models [16, 19, 22, 25, 29, 31, 34, 35, 39, 50, 51, 58, 62, 63]. These methods were tested on two datasets: dataset-a, initially used to validate graph-based methods [50], and dataset-b, initially used to evaluate pre-trained models [35]. Both datasets were split using different random seeds for scaffold splitting, resulting in varying scaffold distributions across the training, validation, and test sets. All models, including ours, were evaluated on consistent test datasets to ensure fair comparison. We measured AUROC for classification tasks and RMSE for regression tasks, averaging results over three random seed values for comparison.

On dataset-a, initially used in GROVER [50], our model achieved an average AUROC of 0.821 across six classification tasks, placing second highest. GROVER led with an average AUROC of 0.834, followed by KANO with 0.814 as the third (Table 3). For regression tasks, GROVER achieved the lowest average RMSE scores of 1.544 and 0.560 on FreeSolv and Lipo, respectively, while our method reached the lowest RMSE of 0.746 on ESOL. On dataset-b, initially used in MPG [35], our model showed strong performance with an average AUROC of 0.851 on six classification tasks, followed by MPG and the KANO with AUROCs of 0.841 and 0.833, respectively. In regression tasks on FreeSolv, ESOL, and Lipo datasets, MPG, KANO, and our method achieved the lowest RMSE scores of 1.269, 0.405, and 0.545, respectively (Table 4).

Table 3 Performance comparison with current state-of-the-art methods on dataset-a

Table 4 Performance comparison with current state-of-the-art methods on dataset-b

Interestingly, several models showed substantial performance differences between the two datasets. Our method and KANO demonstrated stable performance across both datasets, whereas GROVER and MPG showed significant performance variations. These inconsistencies may be due to structural differences in the scaffold distributions of the two datasets, suggesting that each model focused on different molecular regions during training. To address this and ensure fair comparisons, we conducted an extensive evaluation using 30 random seed values, detailed in Sect. "Performance evaluation with multiple random seed values".

Performance evaluation with multiple random seed values

The performance variations observed between dataset-a and dataset-b for several methods (Sect. "Performance evaluation on MoleculeNet dataset") underscore the need to assess robustness across different dataset compositions. Thus, using 30 random seed values, we compared our model to prominent methods, including GROVER, MPG, CD-MVGNN, and KANO, on nine benchmark datasets. GROVER [50] and MPG [35] were selected due to their notable performance difference, KANO [16] as a recent advancement in the field, and CD-MVGNN [40] for its multi-view approach that considers both atoms and bonds for non-isomorphic graph differentiation.

As shown in Table 5, our method attained higher average AUROC scores of 0.957, 0.740, 0.628, and 0.919 on the BBBP, ToxCast, SIDER, and ClinTox, respectively. In the regression tasks, our method also obtained the lowest average RMSE scores of 0.783 and 0.597 on the ESOL and Lipo datasets. Meanwhile, KANO and CD-MVGNN outperformed other models on Tox21 and BACE datasets, with AUROC scores of 0.837 and 0.874, respectively. KANO also achieved the best score of 1.970 in the regression task on the FreeSolv dataset. These results demonstrate the robustness of our model in handling data variations, showing the effectiveness of the ensemble method, particularly for small datasets. Statistical significance was validated through Student’s t-test and Welch’s t-test (Supplementary Fig. 1 and 2).

Table 5 Performance comparison with sufficient random seed values

Among the five models evaluated on nine datasets, our model ranked the best on six datasets and the second-best on two others. These highlight the importance of developing more generalized and integrated models. Additionally, CD-MVGNN displayed the best performances on one dataset and the second-best performances on the two. These results suggest that multi-view approaches consistently improved the prediction of molecular properties. KANO presented the best scores on two datasets and the second-best scores on the three, demonstrating the effectiveness of the pretraining and knowledge- guided learning. In this context, the balanced performance of our model across multiple datasets highlights its potential for a wide range of molecular property prediction tasks, and our method could promise better predictions with pretraining in future work.

Contribution of the ensemble model

We applied an ensemble method to our model to address data scarcity on the nine MoleculeNet datasets, which contain only a few thousand samples. To validate the effectiveness of the ensemble method, we compared our ensemble model to single models across 30 random seed values. We used balanced scaffold splitting to create six single models with distinct e-seed values (from 0 to 5), training each model on a unique training and validation dataset (Figs. 2 and 3). The ensemble model was then constructed by aggregating these six models using soft voting.

Fig. 2

Performance comparison of the ensemble and individual models on classification tasks. We evaluated our ensemble and individual models across 30 random seeds on the classification datasets to validate the robustness and generalization of the ensemble method. Each e-seed n model was independently trained on different training and validation datasets splitting with different e-seed values. The Ensemble model used soft voting over the six e-seed n models. The test dataset was fixed, and all models were evaluated on the same test dataset on a specific seed value. Notable variations in AUROC values were observed among the e-seed n models. The Ensemble model consistently achieved the best performances across all datasets. Significant improvements of 0.009, 0.018, 0.020, 0.016, 0.022, and 0.022 AUROC over the second-highest scores were observed in BBBP, Tox21, ToxCast, SIDER, ClinTox, and BACE, respectively

Fig. 3

Performance comparison of the ensemble and individual models on regression tasks. We evaluated our ensemble and individual models across 30 random seeds on the regression datasets. Each e-seed n model was independently trained, and the Ensemble model used soft voting over the six e-seed n models. The Ensemble model consistently outperformed across all datasets. Significant improvements of 0.202, 0.060, and 0.047 RMSE over the second-lowest scores were observed in FreeSolv, ESOL, and Lipo, respectively

Our results show that the ensemble model consistently outperformed each individual model. For classification tasks, the ensemble model achieved higher AUROC scores than all single models (Fig. 2). The ensemble model improved AUROC scores by 1.4%, 2.4%, 3.0%, 3.3%, 2.9%, and 3.4% on average for the BBBP, Tox21, ToxCast, SIDER, ClinTox, and BACE, respectively. For regression, the ensemble model achieved RMSE improvements of 13.8%, 10.3%, and 8.3% on average for the FreeSolv, ESOL, and Lipo datasets, respectively (Fig. 3). These results indicate that the ensemble method significantly enhances performance, especially in regression tasks. Detailed classification and regression results are presented in Supplementary Table 3 and 4.

Ablation study

We conducted a systematic ablation study to assess the impact of individual components within our model, using MoleculeNet physiology datasets from dataset-b. We first evaluated the complete model (M1 in Fig. 4), which integrates all components such as sub-model for node, sub-model for edge, multi-head attention, and ensemble method. We then tested a variant of the model without the ensemble method (M2 in Fig. 4). From the M2 model, we built and evaluated three additional models, each independently eliminating one component − either the sub-model for node, sub-model for edge, and multi-head attention.

Fig. 4

Ablation study evaluating the contribution of four components in our model. M1 (green) represents the full model, incorporating the atom graph, bond graph, attention mechanism, and ensemble method. M2 (light green) includes all components except the ensemble method. The next three models (light blue) each exclude one component-either the attention mechanism, bond graph, or atom graph. The last three models (pink) include only one component each (either the attention mechanism, bond graph, or atom graph). The full model, M1, achieved the highest performance, with an average AUROC score of 0.836 on MoleculeNet physiology datasets

The ablation study revealed that each component significantly contributed to the overall performance, with the full model achieving the best performance, underscoring the importance of incorporating all components (Fig. 4). The basic models containing only a single component (either the atom graph, bond graph, or attention) showed lower AUROC scores of 0.795, 0.801, and 0.777, respectively, which are approximately 5%, 4%, and 7% lower than our full model. Notably, the comparison between the M1 and M2 models showed a significant 2.3% increase in AUROC score attributed to the ensemble method. Additionally, we observed AUROC improvements ranging from 1.0% to 5.1% when comparing the M2 model to the other models lacking specific components. Detailed results on the physiology datasets are presented in Supplementary Table 5.

Discussion

As with many clinical prediction tasks, this study faced limitations due to the availability of experimentally labeled data, resulting in small sample sizes that can challenge the training of deep learning models and potentially compromise robustness and generalization. Despite these challenges, our multi-view model, which integrates atom, bond, and global features alongside an ensemble method, showed consistently stable performance across most experiments. These findings suggest that our approach to enhancing generalization could be valuable for clinical applications where data availability is limited. We believe this model holds promise for clinical studies, potentially facilitating reliable predictions in scenarios with constrained data resources.

In addition, the ensemble method has a limitation, increasing computational cost as the number of models increases. Based on our experiments, our results indicate that while our model uses higher computational resources than other architectures, it still performs efficiently within practical limits. For example, we measured training time for MPG [35] and our method; MPG is one of the pretrained methods, and our method exploits the ensemble method. To compare these two methods, we chose the BBBP dataset, which includes a small amount of data and one task, and the ToxCast dataset, which contains a relatively large amount of data and 617 tasks. On the BBBP dataset, MPG takes 27 min for finetuning, and our method takes 24 min for the ensemble. If a dataset is small, the ensemble method works efficiently through optimization, such as early stopping. For the ToxCast dataset, MPG takes 46 min, and our method takes one hour and 21 min. Following our anticipation, we observed that a large dataset could demand higher computational resources for the ensemble method, but it still has been within practical limits. Additionally, we can apply several potential optimizations to reduce computational costs, such as reducing the number of attention layers without significant loss of accuracy and leveraging sparse attention mechanisms that selectively focus on critical substructures when it was applied to larger datasets or more complex molecular structures.

To address these challenges associated with the ensemble method, future work should be focused on developing more efficient and practical ensemble approaches for deep learning models. Pretraining methods might be synergetic to the ensemble method by diminishing learning time with fine-tuning. Considering KANO achieved better performances in some of our experiments, adopting the ensemble method to a pre-trained model could promise to reduce computational costs and enhance performance. Additionally, replacing soft voting with attention-based networks could provide improved interpretability, which is essential for many clinical applications.

Conclusions

This study presents a graph-integrated model designed to enhance the prediction of molecular properties by effectively capturing both local and global substructural features. The model incorporates graph attention encoders and multi-head attention mechanisms to a deeper structural understanding. The graph attention networks enable our model to detect more specific structural properties by considering interactions between atoms and edges, while the multi-head attention network captures the relationships between both proximate and distant substructures. This constitution could lead our model to reflect locally and globally essential substructures. Furthermore, the multi-task learning method was adapted to increase performance by facilitating task cooperation and employed an ensemble method to enhance robustness and generalization.

We validated the effectiveness of our model through experiments across 30 random seed values on nine MoleculeNet datasets, comparing it with several state-of-the-art methods. Our model outperformed other models on four classification and two regression datasets, where the results demonstrate higher predictive accuracy. These results underscore the advantage of our ensemble method and highlight the importance of integrating local and global structural information in predicting molecular properties.