Analyzing Learned Molecular Representations for Property Prediction
Advancements in neural machinery have led to a wide range of algorithmic solutions for molecular property prediction. Two classes of models in particular have yielded promising results: neural networks applied to computed molecular fingerprints or expert-crafted descriptors, and graph convolutional neural networks that construct a learned molecular representation by operating on the graph structure of the molecule. However, recent literature has yet to clearly determine which of these two methods is superior when generalizing to new chemical space. Furthermore, prior research has rarely examined these new models in industry research settings in comparison to existing employed models. In this paper, we benchmark models extensively on 19 public and 16 proprietary industrial datasets spanning a wide variety of chemical endpoints. In addition, we introduce a graph convolutional model that consistently matches or outperforms models using fixed molecular descriptors as well as previous graph neural architectures on both public and proprietary datasets. Our empirical findings indicate that while approaches based on these representations have yet to reach the level of experimental reproducibility, our proposed model nevertheless offers significant improvements over models currently used in industrial workflows.
MIT]Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA 02139, United States MIT]Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA 02139, United States MIT]Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA 02139, United States MIT2]Department of Chemical Engineering, MIT, Cambridge, MA 02139, United States BASF]BASF SE, Ludwigshafen 67063, Germany Amgen]Amgen Inc., Cambridge, MA 02141, United States Amgen]Amgen Inc., Cambridge, MA 02141, United States Amgen]Amgen Inc., Cambridge, MA 02141, United States Novartis]Novartis Institutes for BioMedical Research, Cambridge, MA 02139, United States BASF]BASF SE, Ludwigshafen 67063, Germany BASF]BASF SE, Ludwigshafen 67063, Germany BASF]BASF SE, Ludwigshafen 67063, Germany MIT]Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA 02139, United States MIT2]Department of Chemical Engineering, MIT, Cambridge, MA 02139, United States MIT]Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA 02139, United States
Molecular property prediction, one of the oldest cheminformatics tasks, has received new attention in light of recent advancements in deep neural networks. These architectures either operate over fixed molecular fingerprints common in traditional QSAR models, or they learn their own task-specific representations using graph convolutions duvenaud2015convolutional, Wu_2018, kearnes2016molecular, gilmer2017neural, li2015gated, kipf2016semi, defferrard2016convolutional, bruna2013spectral, coley_convolutional_2017, schutt2017quantum, battaglia2016interaction. Both approaches are reported to yield substantial performance gains, raising state-of-the-art accuracy in property prediction.
Despite these successes, many questions remain unanswered. The first question concerns the comparison between learned molecular representations and fingerprints or descriptors. Unfortunately, current published results on this topic do not provide a clear answer. Wu_2018Wu_2018 demonstrate that convolution-based models typically outperform fingerprint-based models, while experiments reported in mayr2018chemblmayr2018chembl report the opposite. Part of these discrepancies can be attributed to differences in evaluation setup, including the way datasets are constructed. This leads us to a broader question concerning current evaluation protocols and their capacity to measure the generalization power of a method when applied to a new chemical space, as is common in drug discovery. Unless special care is taken to replicate this distributional shift in evaluation, neural models may overfit the training data but still score highly on the test data. This is particularly true for convolutional models that can learn a poor molecular representation by memorizing the molecular scaffolds in the training data and thereby fail to generalize to new ones. Therefore, a meaningful evaluation of property prediction models needs to account explicitly for scaffold overlap between train and test data in light of generalization requirements.
In this paper, we aim to answer both of these questions by designing a comprehensive evaluation setup for assessing neural architectures. We also introduce an algorithm for property prediction that outperforms existing strong baselines across a range of datasets. The model has two distinctive features: (1) It operates over a hybrid representation that combines convolutions and descriptors. This design gives it flexibility in learning a task specific encoding, while providing a strong prior with fixed descriptors. (2) It learns to construct molecular encodings by using convolutions centered on bonds instead of atoms, thereby avoiding unnecessary loops during the message passing phase of the algorithm.
We extensively evaluate our model and other recently published neural architectures with over 850 experiments on 19 publicly available benchmarks from Wu_2018Wu_2018 and mayr2018chemblmayr2018chembl and on 16 proprietary datasets from Amgen, Novartis, and BASF (Badische Anilin und Soda Fabrik). Our goal is to assess whether the models’ performance on the public datasets and their relative ranking are representative of their ranking on the proprietary datasets. We demonstrate that under a scaffold split of training and testing data, the relative ranking of the models is consistent across the two classes of datasets. We also show that a scaffold-based split of the training and testing data is a good approximation of the temporal split commonly used in industry in terms of the relevant metrics. By contrast, a purely random split is a poor approximation to a temporal split, confirming the findings of sheridan2015relativesheridan2015relative. To put the performance of current models in perspective, we report bounds on experimental error and show that there is still room for improving deep learning models to match the accuracy and reproducibility of screening results.
Building on the diversity of our benchmark datasets, we explore the impact of molecular representation with respect to the dataset characteristics. We find that a hybrid representation yields higher performance and generalizes better than either convolution-based or fingerprint-based models. We also note that on small datasets (up to 1000 training molecules) fingerprint models can outperform learned representations, which are negatively impacted by data sparsity. Beyond molecular representation issues, we observe that hyperparameter selection plays a crucial role in model performance, consistent with prior workshahriari2016bayes. We show that Bayesian optimization yields a robust, automatic solution to this issue. The addition of ensembling further improves accuracy, again consistent with the literaturedietterich2000ensemble.
Our experiments show that our model achieves consistently strong out-of-the-box performance and even stronger optimized performance across a wide variety of public and proprietary datasets. Our model achieves comparable or better performance on 11 out of 19 public datasets and on 15 out of 16 proprietary datasets compared to all baseline models. Furthermore, no single baseline model is clearly superior across the remaining 8 public datasets, and the relative performance of the baseline models often varies from dataset to dataset, whereas our model is consistently strong across datasets. These results indicate that our model, and learned molecular fingerprints in general, are applicable and ready to be used as a powerful tool for chemists actively working on drug discovery.
Since the core element of our model is a graph encoder architecture, our work is closely related to previous work on graph encoders, such as those for social networkskipf2016semi, hamilton2017inductive or for chemistry applicationsscarselli2009graph, bruna2013spectral, duvenaud2015convolutional, henaff2015deep, dai2016discriminative, defferrard2016convolutional, lei2017deriving, coley_convolutional_2017, kusner2017grammar, gomez2018automatic, jin2018junction, jin2018learning.
Common approaches to molecular property prediction today involve the application of well-known models like support vector machinescortes1995support or random forestsbreiman2001random to expert-engineered descriptors or molecular fingerprints, such as the Dragon descriptorsmauri2006dragon or Morgan (ECFP) fingerprintsrogers2010extended. One direction of advancement is the use of domain expertise to improve the base feature representation of molecular descriptorsswamidass2005kernels, cao2013chemopy, durant2002reoptimization, moriwaki2018mordred, mauri2006dragon to drive better performancemayr2018chembl. Additionally, many studies have leveraged explicit 3D atomic coordinates to improve performance furtherWu_2018, schutt2017schnet, kondor2018covariant, faber2017machine, feinberg2018potentialnet.
The other main line of research is the optimization of the model architecture, whether the model is applied to descriptors or fingerprintsmayr2018chembl, Lee201810847 or is directly applied to SMILESweininger1988smiles stringsmayr2018chembl or the underlying graph of the moleculeduvenaud2015convolutional, Wu_2018, kearnes2016molecular, gilmer2017neural, li2015gated, kipf2016semi, defferrard2016convolutional, coley_convolutional_2017, bruna2013spectral, schutt2017quantum, battaglia2016interaction. Our model belongs to the last category of models, known as graph convolutional neural networks. In essence, such models learn their own expert feature representations directly from the data, and they have been shown to be very flexible and capable of capturing complex relationships given sufficient datagilmer2017neural, Wu_2018.
In a direction orthogonal to our own improvements, ishiguro2019graphishiguro2019graph also make a strong improvement to graph neural networks. liu2018chemiliu2018chemi also evaluate their model against private industry datasets, but we cannot compare against their method directly owing to dataset differencesliu2018chemi.
The property prediction models most similar to our own are encapsulated in the Message Passing Neural Network (MPNN) framework presented in gilmer2017neuralgilmer2017neural. We build upon this basic framework by adopting a message-passing paradigm based on updating representations of directed bonds rather than atoms. Additionally, we further improve the model by combining computed molecule-level features with the molecular representation learned by the MPNN.
We first summarize MPNNs in general using the terminology of gilmer2017neuralgilmer2017neural, and then we expand on the characteristics of Directed MPNN (D-MPNN)dai2016discriminative used in this paper. (D-MPNN is originally called structure2vec in dai2016discriminativedai2016discriminative. In this paper, we refer to it as Directed MPNN to show it is a variant of the generic MPNN architecture.)
3.1 Message Passing Neural Networks
An MPNN is a model which operates on an undirected graph with node (atom) features and edge (bond) features . MPNNs operate in two phases: a message passing phase, which transmits information across the molecule to build a neural representation of the molecule, and a readout phase, which uses the final representation of the molecule to make predictions about the properties of interest.
More specifically, the message passing phase consists of steps. On each step , hidden states and messages associated with each vertex are updated using message function and vertex update function according to:
where is the set of neighbors of in graph , and is some function of the initial atom features . The readout phase then uses a readout function to make a property prediction based on the final hidden states according to
The output may be either a scalar or a vector, depending on whether the MPNN is designed to predict a single property or multiple properties (in a multitask setting).
During training, the network takes molecular graphs as input and makes an output prediction for each molecule. A loss function is computed based on the predicted outputs and the ground truth values, and the gradient of the loss is backpropagated through the readout phase and the message passing phase. The entire model is trained end-to-end.
3.2 Directed MPNN
The main difference between the Directed MPNN (D-MPNN)dai2016discriminative and the generic MPNN described above is the nature of the messages sent during the message passing phase. Rather than using messages associated with vertices (atoms), D-MPNN uses messages associated with directed edges (bonds). The motivation of this design is to prevent totters mahe2004extensions, that is, to avoid messages being passed along any path of the form where for some . Such excursions are likely to introduce noise into the graph representation. Using Figure 1 as an illustration, in D-MPNN, the message will only be propagated to nodes 3 and 4 in the next iteration, whereas in the original MPNN it will be sent to node 1 as well, creating an unnecessary loop in the message passing trajectory. Compared to the atom based message passing approach, this message passing procedure is more similar to belief propagation in probabilistic graphical modelskoller2009probabilistic. We refer to dai2016discriminativedai2016discriminative for futher discussion about the connection between D-MPNN and belief propagation.
The D-MPNN works as follows. The D-MPNN operates on hidden states and messages instead of on node based hidden states and messages . Note that the direction of messages matters (i.e., and are distinct from and ). The corresponding message passing update equations are thus
Observe that message does not depend on its reverse message from the previous iteration. Prior to the first step of message passing, we initialize edge hidden states with
where is a learned matrix, is the concatenation of the atom features for atom and the bond features for bond , and is the ReLU activation functionnair2010relu.
We choose to use relatively simple message passing functions and edge update functions . Specifically, we define and we implement with the same neural network on every step,
where is a learned matrix with hidden size . Note that the addition of on every step provides a skip connection to the original feature vector for that edge.
Finally, we return to an atom representation of the molecule by summing the incoming bond features according to
where is a learned matrix.
Altogether, the D-MPNN message passing phase operates according to
for , followed by
The readout phase of the D-MPNN is the same as the readout phase of a generic MPNN. In our implementation of the readout function , we first sum the atom hidden states to obtain a feature vector for the molecule
Finally, we generate property predictions where is a feed-forward neural network.
3.3 Initial Featurization
Our model’s initial atom and bond features are listed in Tables 1 and 2, respectively. The D-MPNN’s initial node features are simply the atom features for that node, while the D-MPNN’s initial edge features are the bond features for bond . All features are computed using the open-source package RDKitlandrum2006rdkit.
|Atom type||Type of atom (ex. C, N, O), by atomic number.||100|
|# Bonds||Number of bonds the atom is involved in.||6|
|Formal charge||Integer electronic charge assigned to atom.||5|
|Chirality||Unspecified, tetrahedral CW/CCW, or other.||4|
|# Hs||Number of bonded Hydrogen atom.||5|
|Hybridization||sp, sp2, sp3, sp3d, or sp3d2.||5|
|Aromaticity||Whether this atom is part of an aromatic system.||1|
|Atomic mass||Mass of the atom, divided by 100.||1|
|Bond type||Single, double, triple, or aromatic.||4|
|Conjugated||Whether the bond is conjugated.||1|
|In ring||Whether the bond is part of a ring.||1|
|Stereo||None, any, E/Z or cis/trans.||6|
3.4 D-MPNN with Features
Next, we discuss further extensions and optimizations to improve performance. Although an MPNN should ideally be able to extract any information about a molecule that might be relevant to predicting a given property, two limitations may prevent this in practice. First, many property prediction datasets are very small, i.e., on the order of only hundreds or thousands of molecules. With so little data, MPNNs are unable to learn to identify and extract all features of a molecule that might be relevant to property prediction, and they are susceptible to overfitting to artifacts in the data. Second, most MPNNs use fewer message passing steps than the diameter of the molecular graph, i.e. , meaning atoms that are a distance of greater than bonds apart will never receive messages about each other. This results in a molecular representation that is fundamentally local rather than global in nature, meaning the MPNN may struggle to predict properties that depend heavily on global features.
In order to counter these limitations, we introduce a variant of the D-MPNN that incorporates 200 global molecular features that can be computed rapidly in silico using RDKit. The neural network architecture requires that the features are appropriately scaled to prevent features with large ranges dominating smaller ranged features, as well as preventing issues where features in the training set are not drawn from the same sample distribution as features in the testing set. To prevent these issues, a large sample of molecules was used to fit cumulative density functions (CDFs) to all features. CDFs were used as opposed to simpler scaling algorithms mainly because CDFs have the useful property that each value has the same meaning: the percentage of the population observed below the raw feature value. Min-max scaling can be easily biased with outliers and Z-score scaling assumes a normal distribution which is most often not the case for chemical features, especially if they are based on counts.
The CDFs were fit to a sample of 100k compounds from the Novartis internal catalog using the distributions available in the scikit-learn packagepedregosa2011scikit, a sample of which can be seen in Figure 2. One could do a similar normalization using publicly available databases such as ZINCirwin2005zinc and PubChemkim2015pubchem. scikit-learn was used primarily due to the simplicity of fitting and the final application. However, more complicated techniques could be used in the future to fit to empirical CDFs, such as finding the best fit general logistic function, which has been shown to be successful for other biological datasetssartor2008lrpath. No review was taken to remove odd distributions. For example, azides are hazardous and rarely used outside of a few specific reactions, as reflected in the fr_azide distribution in Figure 2. As such, since the sample data was primarily used for chemical screening against biological targets, the distribution used here may not accurately reflect the distribution of reagents used for chemical synthesis. For the full list of calculated features, please refer to the Supporting Information.
To incorporate these features, we modify the readout phase of the D-MPNN to apply the feed-forward neural network to the concatenation of the learned molecule feature vector and the computed global features ,
This is a very general method of incorporating external information and can be used with any MPNN and any computed features or descriptors.
3.5 Hyperparameter Optimization
The performance of MPNNs, like most neural networks, can depend greatly on the settings of the various model hyperparameters, such as the hidden size of the neural network layers. Thus to maximize performance, we perform hyperparameter optimization via Bayesian Optimizationshahriari2016bayes using the Hyperopthyperopt Python package. We specifically optimize our model’s depth (number of message-passing steps), hidden size (size of bond message vectors), number of feed-forward network layers, and dropout probability.
A common technique in machine learning for improving model performance is ensembling, where the predictions of multiple independently trained models are combined to produce a more accurate predictiondietterich2000ensemble. We apply this technique by training several copies of our model, each initialized with different random weights, and then averaging the predictions of these models (each with equal weight) to generate an ensemble prediction.
Since prior work did not report performance using ensembling, all direct comparisons we make to prior work use a single D-MPNN model for a fair comparison. However, we also report results using an ensemble to illustrate the maximum possible performance of our model architecture.
We implement our model using the PyTorchpaszke2017pytorch deep learning framework. All code for the D-MPNN and its variants is available in our GitHub repository chemprop_code. Code for computing and using the RDKit feature CDFs is available in the Descriptastorus packagedescriptastorus. Additionally, a web demonstration of our model’s predictive capability on public datasets is available onlinechemprop_website.
We test our model on 19 publicly available datasets from Wu_2018Wu_2018 and mayr2018chemblmayr2018chembl. These datasets range in size from less than 200 molecules to over 450,000 molecules. They include a wide range of regression and classification targets spanning quantum mechanics, physical chemistry, biophysics, and physiology. Detailed descriptions are provided in Table 3.
|QM7, QM8, QM9||Quantum Mechanics||Computer-generated quantum mechanical properties|
|ESOL||Physical Chemistry||Water solubility|
|FreeSolv||Physical Chemistry||Hydration free energy in water|
|Lipophilicity||Physical Chemistry||Octanol/water distribution coefficients|
|PDBbind||Biophysics||Protein binding affinity|
|PCBA||Biophysics||Assorted biological assays|
|MUV||Biophysics||Assorted biological assays|
|HIV||Biophysics||Inhibition of HIV replication|
|BACE||Biophysics||Inhibition of human -secretase|
|BBBP||Physiology||Ability to penetrate the blood-brain-barrier|
|SIDER||Physiology||Side effects of drugs|
Summary statistics for all the datasets are provided in Table 4, and further details on the datasets are available in Wu_2018Wu_2018, with the exception of the ChEMBL dataset which is described in mayr2018chemblmayr2018chembl. Additional information on the class balance of the classification datasets is provided in the Supporting Information. Although most classification datasets are reasonably balanced, the MUV dataset is particularly unbalanced, with only 0.2% of molecules classified as positive. This makes our model unstable, leading to the wide variation in performance on this dataset in the subsequent sections.
It is worth noting that for some datasets, the number of compounds in Table 4 does not precisely match the numbers from Wu_2018Wu_2018. This is because Wu_2018Wu_2018 included duplicate molecules in that count while we count the unique number of molecules. Additionally, we left out one or two molecules which could not be processed by RDKitlandrum2006rdkit. However, the impact of removing these molecules is negligible on overall model performance. Furthermore, we have fewer molecules in QM7 because we used SMILES strings generated by Wu_2018Wu_2018 from the original 3D coordinates in the dataset, but the SMILES conversion process failed for molecules. For this reason, we do not directly compare our model’s performance on QM7 to the QM7 performance numbers reported by Wu_2018Wu_2018.
|Dataset||# Tasks||Task Type||# Compounds||Metric|
4.2 Experimental Procedure
Cross-Validation and Hyperparameter Optimization.
Since many of the datasets are very small (two thousand molecules or fewer), we use a cross-validation approach to decrease noise in the results both while optimizing the hyperparameters and while determining final performance numbers. For consistency, we maintain the same approach for all of our datasets. Specifically, for each dataset, we use 20 iterations of Bayesian optimization on 10 randomly-seeded 80:10:10 data splits to determine the best hyperparameters, selecting hyperparameters based on validation set performance. We then evaluate the model by retraining using the optimal hyperparameters and checking performance on the test set. Due to computational cost, we only use 3 splits for HIV, QM9, MUV, PCBA, and ChEMBL. When we run the best model from mayr2018chemblmayr2018chembl for comparative purposes, we optimize their model’s hyperparameters with the same splits, using their original hyperparameter optimization script.
We evaluate all models on random and scaffold-based splits as well as on the original splits from Wu_2018Wu_2018 and mayr2018chemblmayr2018chembl. The one exception is the model of mayr2018chemblmayr2018chembl, which we only ran on scaffold-based splits, due to the large computational cost of optimizing their model. Results on scaffold-based splits are reported below while results on random splits are presented in the Supporting Information.
Our scaffold split is similar to that of Wu_2018Wu_2018. Molecules are partitioned into bins based on their Murcko scaffold calculated by RDKitlandrum2006rdkit. Any bins larger than half of the desired test set size are placed into the training set, in order to guarantee the scaffold diversity of the validation and test sets. All remaining bins are placed randomly into the training, validation, and test sets until each set has reached its desired size. As this latter process involves randomly placing scaffolds into bins, we are able to generate several different scaffold splits for evaluation.
None of our splits on classification datasets are stratified; we do not enforce class balance. Compared to random splits, the scaffold splits have more class imbalance on average, but are not excessively imbalanced; we analyze this class balance quantitatively in the Additional Dataset Statistics section of the Supporting Information.
Compared to a random split, a scaffold split is a more challenging and realistic evaluation setting as shown in Figures 11 and 13. This allows us to use a scaffold split as a proxy for the chronological split present in real-world property prediction data, where one trains a model on past data to make predictions on future data, although chronological splits are still preferred when available. However, as chronological information is not available for most public datasets, we use a scaffold-based split for all evaluations except for our direct comparison with the MoleculeNet models from Wu_2018Wu_2018, for which we use their original data splits.
We compare our model to the following baselines:
The best model for each dataset from MoleculeNet by Wu_2018Wu_2018
The best model from mayr2018chemblmayr2018chembl, a feed-forward neural network on a concatenation of assorted expert-designed molecular fingerprints.
Random forest on binary Morgan fingerprints.
Feed-forward network (FFN) on binary Morgan fingerprints using the same FFN architecture that our D-MPNN uses during its readout phase.
FFN on count-based Morgan fingerprints.
FFN on RDKit-calculated descriptors.
The models in MoleculeNet by Wu_2018Wu_2018 include MPNNgilmer2017neural, Weavekearnes2016molecular, GraphConv, kernel ridge regression, gradient boostingfriedman2001greedy, random forestbreiman2001random, logistic regressioncramer2010logreg, directed acyclic graph modelslusci2013deep, support vector machinesvapnik1995svm, Deep Tensor Neural Networksschutt2017quantum, multitask networksma2015deep, bypass networksramsundar2017multitask, influence relevance votingswamidass2009influence, and/or ANI-1smith2017ani, depending on the dataset. Full details can be found in Wu_2018Wu_2018. For the feed-forward network model from mayr2018chemblmayr2018chembl, we modified the authors’ original code with their guidance in order to run their code on all of the datasets, not just on the ChEMBL dataset they experimented with. We tuned learning rates and hidden dimensions in addition to the extensive hyperparameter search already present in their code.
5 Results and Discussion
In the following sections, we analyze the performance of our model on both public and proprietary datasets. Specifically, we aim to answer the following questions:
How does our model perform on both public and proprietary datasets compared to public benchmarks, and how close are we to the upper bound on performance represented by experimental reproducibility?
How should we be splitting our data, and how does the method of splitting affect our evaluation of the model’s generalization performance?
What are the key elements of our model, and how can we maximize its performance?
In the following sections, all results using root-mean-square error (RMSE) or mean absolute error (MAE) are displayed as plots showing change relative to a baseline model rather than showing absolute performance numbers. This is because the scale of the errors can differ drastically between datasets. All results using , area under the receiver operating characteristic curve (ROC-AUC), or area under the precision recall curve (PRC-AUC) are displayed as plots showing the actual values. For RMSE and MAE, lower is better, while for , ROC-AUC, and PRC-AUC, higher is better. Table 4 indicates the metric used for each dataset. Tables showing the exact performance numbers for all experiments can be found in the Supporting Information. Note that the error bars on all plots show the standard error of the mean across multiple runs, where standard error is defined as the standard deviation divided by the square root of the number of runs.
We evaluate statistical significance using two statistical tests: a one-sided Wilcoxon signed-rank test and a one-sided Welch’s t-test. While the Wilcoxon test is stronger, as it is a paired test comparing performance molecule-by-molecule, it requires knowing per-molecule predictions, which we do not have easy access to for the models from MoleculeNetWu_2018 and mayr2018chemblmayr2018chembl. Furthermore, comparisons between data split types inherently involves comparing performance on different test molecules, meaning a per-molecule test is not possible. Therefore, for these comparison we use the weaker Welch’s t-test and for all other comparisons we use the Wilcoxon test. When using the Wilcoxon test for regression datasets, we directly compare test errors molecule-by-molecule. For the classification datasets, we divide all the test molecules into 30 equal parts, compute AUC on each part, and then use the Wilcoxon test on these AUC values. This subdivision of the test molecules into 30 parts gives the Wilcoxon test more strength than evaluating directly on the original 3 or 10 test cross-validation folds while still keeping each part large enough to result in a meaningful AUC computation. We define statistical significance as p-value less than 0.05.
Additionally, note that in all figures and tables, “D-MPNN” refers to the base D-MPNN model, “D-MPNN Features” refers to the D-MPNN with RDKit features, “D-MPNN Optimized” refers to the D-MPNN with RDKit features and optimized hyperparameters, and “D-MPNN Ensemble” refers to an ensemble of five D-MPNNs with RDKit features and optimized hyperparameters.
5.1 Comparison to Baselines
After optimizing our model, we compare our best single (non-ensembled) model on each dataset against models from prior work.
5.1.1 Comparison to MoleculeNet
We first compare our D-MPNN to the best model from MoleculeNetWu_2018, Ramsundar-et-al-2019 on the same datasets and splits on which Wu_2018Wu_2018 evaluate their models. We were unable to reproduce their original data splits on BACE, Toxcast, and QM7, but we have evaluated our model against their original splits on all of the other datasets. The splits are a mix of random, scaffold, and time splits, as indicated in Figure 3.
Overall on the 10 datasets where the MoleculeNet models only use 2D information, i.e. all datasets except the QM and PDBbind datasets, our D-MPNN is significantly better than the best MoleculeNet models on 5 datasets, is not significantly different on 3 datasets, and is significantly worse on 2 datasets. This indicates that D-MPNN tends to outperform even the best MoleculeNet models, with the added benefit that the D-MPNN model architecture is the same for every dataset while the best MoleculeNet model architecture differs between datasets.
Furthermore, we note that there are two cases in which our D-MPNN may underperform. The first is the MUV dataset, which is large but extremely imbalanced; only of samples are labeled as positives. Wu_2018Wu_2018 also encountered great difficulty with this extreme class imbalance when experimenting with the MUV dataset; all other datasets we experiment on contain at least % positives (see the Supporting Information for full class balance information). The second exception is when there is auxiliary 3D information available, as in the three variants of the PDBbind dataset and in QM9. The current iteration of our D-MPNN does not use 3D coordinate information, and we leave this extension to future work. Thus it is unsurprising that our D-MPNN model underperforms models using 3D information on a protein binding affinity prediction task such as PDBbind, where 3D structure is key. Nevertheless, our D-MPNN model outperforms the best graph-based method in MoleculeNet on PDBbind and QM9. Moreover, we note that on another dataset that provides 3D coordinate information, QM8, our model outperforms the best model in MoleculeNet with or without 3D coordinates.
5.1.2 Comparison to mayr2018chemblmayr2018chembl
In addition, we compare D-MPNN to the baseline from mayr2018chemblmayr2018chembl in Figure 4. We reproduced the features from their best model on each dataset using their scripts or equivalent packageslsc_experiments. We then ran their code and hyperparameter optimization directly on the classification datasets, and we modified their code to run on regression datasets with the authors’ guidancelsc_experiments. On most classification datasets, we obtain similar performance to mayr2018chemblmayr2018chembl. On regression datasets, the baseline from mayr2018chemblmayr2018chembl performs poorly in comparison, despite extensive tuning. We hypothesize that this poor performance on regression in comparison to classification is the result of a large number of binary input features to the output feed-forward network; this hypothesis is supported by the similarly poor performance of our Morgan fingerprint FFN baseline. In addition, their method does not employ early stopping based on validation set performance and therefore may overfit to the training data in some cases; this may be the source of some numerical instability.
Overall, our D-MPNN is significantly better than the mayr2018chemblmayr2018chembl model on 8 datasets, is not significantly different on 10 datasets, and is significantly worse on 1 dataset. This indicates that D-MPNN generally outperforms the mayr2018chemblmayr2018chembl model, especially on regression datasets.
5.1.3 Out-of-the-Box Comparison of D-MPNN to Other Baselines
For our final baseline comparison, we evaluate our model’s performance “out-of-the-box,” i.e. using all the default settings (hidden size = 300, depth = 3, number of feed-forward layers = 2, dropout = 0) without any hyperparameter optimization and without any additional features. For this comparison, we compare to a number of simple baseline models that use computed fingerprints or descriptors:
Random forest (RF) with 500 trees run on Morgan (ECFP) fingerprints using radius 2 and hashing to a bit vector of size 2048.
Feed-forward network (FFN) on Morgan fingerprints.
FFN on Morgan fingerprints which use substructure counts instead of bits.
FFN on RDKit descriptors.
The parameters of the simple baseline models are also out-of-the-box defaults. We make this comparison in order to demonstrate the strong out-of-the-box performance of our model across a wide variety of datasets. Finally, we include the performance of the automatically optimized version of our model as a reference.
Figure 5 shows that even without optimization, our D-MPNN provides an excellent starting point on a wide variety of datasets and targets, though it can be improved further with proper optimization.
5.2 Proprietary Datasets
We also ran our model on several private industry datasets, verifying that our model’s strong performance on public datasets translates to real-world industrial datasets.
We ran our model along with mayr2018chembl’smayr2018chembl model and our simple baselines on four internal Amgen regression datasets. The datasets are as follows.
Rat plasma protein binding free fraction (rPPB).
Solubility in 0.01 M hydrochloric acid solution, pH 7.4 phosphate buffer solution, and simulated intestinal fluid (Sol HCL, Sol PBS, and Sol SIF respectively).
Rat liver microsomes intrinsic clearance (RLM).
Human pregnane X receptor % activation at 2uM and 10uM (hPXR).
In addition, we binarized the hPXR dataset according to Amgen’s recommendations in order to evaluate on a classification dataset. Details of the datasets are shown in Table 5. Throughout the following, note that rPPB is in logit while Sol and RLM are in .
|Category||Dataset||# Tasks||Task Type||# Compounds||Metric|
For each dataset, we evaluate on a chronological split. Our model outperforms the baselines on 4 out of the 5 datasets, as shown in Figure 6. Thus our D-MPNN’s strong performance on scaffold splits of public datasets can translate well to chronological splits of private industry datasets.
We ran our model on 10 highly related quantum mechanical datasets from BASF. Each dataset contains 13 properties calculated on the same 30,733 molecules, varying the solvent in each dataset. Dataset details are in Table 6.
|Category||Dataset||Tasks||Task Type||# Compounds||Metric|
|Quantum Mechanics||Ethyl acetate||13||regression||30,733||R2|
For these datasets, we used a scaffold-based split because a chronological split was unavailable. We found that the model of mayr2018chemblmayr2018chembl is numerically unstable on these datasets, and we therefore omit it from the comparison below. Once again we find that our model, originally designed to succeed on a wide range of public datasets, is robust enough to transfer to proprietary datasets as shown in Figure 7.
|Category||Dataset||Tasks||Task Type||# Compounds||Metric|
5.3 Experimental Error
As a final “oracle” baseline, we compare our model’s performance with to an experimental upper upper bound: the agreement between multiple runs of the same assay, which we refer to as the experimental error. Figure 9 shows the R2 of our model on the private Amgen regression datasets together with the experimental error; in addition, this graph shows the performance of Amgen’s internal model using expert-crafted descriptors. Both models remain far less accurate than the corresponding ground truth assays. Thus there remains significant space for further performance improvement in the future.
5.4 Analysis of Split Type
We now justify our use of scaffold splits for performance evaluation. The ultimate goal of building a property prediction model is to predict properties on new chemistry in order to aid the search for drugs from new classes of molecules. On proprietary company datasets, performance on new chemistry is evaluated using a chronological split of the data, i.e., everything before a certain date serves as the training set while everything after that date serves as the test set. This approximates model performance on molecules that chemists are likely to investigate in the future. Since chronological data is typically unavailable for public datasets, we investigate whether we can use our scaffold split as a reasonable proxy for a chronological split, following the work of sheridan2015relativesheridan2015relative.
Figure 10 provides motivation for a scaffold split approach. As illustrated in the figure, train and test sets according to a chronological split share fewer molecular scaffolds than train and test sets split randomly. Since our scaffold split enforces zero molecular scaffold overlap between the train and test set, it should ideally provide a split that is at least as difficult as a chronological split.
As illustrated in Figures 11, 12, and 13, performance on our scaffold split is on average closer to performance on a chronological split on proprietary datasets from Amgen and Novartis and on the public PDBbind datasets. However, the results are noisy due to the nature of chronological splitting, where we only have a single data split, as opposed to random and scaffold splitting, which both have a random component and can generate different splits depending on the random seed. We can alleviate the problem with noise in chronological datasets by using a sliding time window to get different equally-sized splits, at the cost of significantly decreasing the dataset size. We report results on such sliding window splits in the Supporting Information, as the conclusions from these splits are qualitatively similar to those in the main paper.
Figure 14 shows the difference between a random split and a scaffold split on the publicly available datasets, further demonstrating that a scaffold split generally results in a more difficult, and ideally more useful, measure of performance. Therefore, all of our results are reported on a scaffold split rather than a random split in order to better reflect the generalization ability of our model on new chemistry. Nevertheless, it should be emphasized that the chronological split is still the ideal split on which to evaluate when it is available. Thus we additionally report the results on chronological splits on all datasets where they are available.
Overall, our results confirm the findings of sheridan2015relativesheridan2015relative that scaffold and chronological splits are more difficult than random splits, and hence scaffold splits should be preferred over random splits during evaluation. Our findings differ somewhat from those in sheridan2015relativesheridan2015relative in that we find some evidence that chronological splits may actually be harder than scaffold splits. However, owing to the small number of datasets where chronological splits are available, further investigation is necessary on this point, ideally on a larger range of datasets.
Finally, we analyze and justify our modeling choices and optimizations.
5.5.1 Message Type
The most important distinction between our D-MPNN and related work is the nature of the messages being passed across the molecule. Most prior work uses messages centered on atoms whereas our D-MPNN uses messages centered on directed bonds. To isolate the effect of the message passing paradigm on property prediction performance, we implemented message passing on undirected bonds and on atoms as well, as detailed in the Supporting Information and in our code. Figure 15 illustrates the differences in performance between these three types of message passing. While on average the method using directed bonds outperforms the alternatives, the results are largely not statistically significant, so more investigation is warranted on this point.
5.5.2 RDKit Features
Next, we examined the impact of adding additional molecule-level features from RDKit to our model. Figure 16 shows the effect on model performance. The results appear to be highly dataset-dependent. Some datasets, such as QM9 and ESOL, show marked improvement with the addition of features, while other datasets, such as PCBA and HIV, actually show worse performance with the features. We hypothesize that this is because the features are particularly relevant to certain tasks while possibly confusing and distracting the model on other tasks. This implies that our model’s performance on a given dataset may be further optimized by selecting different features more relevant to the task of interest.
Another interesting trend is the effect of adding features to the three PDBbind datasets. The features appear to help on all three datasets, but the benefit is much more pronounced on the extremely small PDBbind-C (core) dataset than it is on the larger PDBbind-R (refined) and PDBbind-F (full) datasets. This indicates that the features may help compensate for the lack of training data and thus may be particularly relevant in low-data regimes. In particular, we hypothesize that the features may help to regularize a representation derived from a small dataset: because the features are derived from more general chemical knowledge, they implicitly provide the model some understanding of a larger chemical domain. Thus, it is worthwhile to consider the addition of features both when they are particularly relevant to the task of interest and when the dataset is especially small.
5.5.3 Hyperparameter Optimization
To improve model performance, we performed Bayesian Optimization to select the best model hyperparameters for each dataset. Figure 17 illustrates the benefit of performing this optimization, as model performance improves on virtually every dataset. Interestingly, some datasets are particularly sensitive to hyperparameters. While most datasets experience a moderate 2-5% improvement in performance following hyperparameter optimization, the quantum mechanics datasets (QM7, QM8, and QM9) and PCBA see dramatic improvements in performance, with our D-MPNN model performing 37% better on QM9 after optimization.
To maximize performance, we trained an ensemble of models. For each dataset, we selected the best single model—i.e. the best hyperparameters along with the RDKit features if the features improved performance—and we trained five models instead of one. The results appear in Figure 19. On most datasets, ensembling only provides a small 1-5% benefit, but as with hyperparameter optimization, there are certain datasets, particularly the quantum mechanics datasets, which especially benefit from the effect of ensembling.
While each of the latter three optimizations (RDKit descriptors, hyperparameter optimization, and ensembling) on its own has limited benefits, altogether they significantly improve the model’s performance on every dataset except MUV.
5.5.5 Effect of Data Size
Finally, we analyze the effect of data size on the performance of our model, using the ChEMBL dataset. ChEMBL is a large dataset of 456,331 molecules on 1,310 targets, but is extremely sparse: only half of the 1,310 targets have at least labels. For this analysis, we use the original scaffold-based split of mayr2018chemblmayr2018chembl, containing 3 cross-validation folds. From Figure 20, we hypothesize that our D-MPNN struggles on low-label targets in comparison to this baseline. As our D-MPNN model does not use any human-engineered fingerprints or descriptors and must therefore learn its features completely from scratch based on the input data, it would be unsurprising if the average ROC-AUC score of D-MPNN is worse than that of the feed-forward network running on human-engineered descriptors in mayr2018chemblmayr2018chembl. When we filter the ChEMBL dataset by pruning low-data targets at different thresholds, we find that our D-MPNN indeed may outperform the best model of mayr2018chemblmayr2018chembl at larger data thresholds (Figure 20), though our results are not fully conclusive.
6 Conclusion and Future Work
|Baseline||D-MPNN is better||D-MPNN is same||D-MPNN is worse||# datasets|
|RF on Morgan||9||1||4||15|
|FFN on Morgan||14||5||0||19|
|FFN on Morgan Counts||15||4||0||19|
|FFN on RDKit||8||5||4||19|
In this paper, we performed an extensive comparison of molecular property prediction models based on either fixed descriptors or learned molecular representations by performing over 850 experiments on 19 public and 16 proprietary datasets. Table 8 shows a summary of how our D-MPNN compares to each of the baseline models. Our model consistently matches or outperforms each baseline individually, and across all baselines, our model achieves comparable or better performance on 11 of the 19 public datasets: QM7, QM8, QM9, ESOL, FreeSolv, Lipophilicity, BBBP, PDBbind-F, PCBA, Tox21, and ClinTox. On the remaining 8 datasets, no single baseline model is consistently superior. Furthermore, our model’s strong results transfer to proprietary datasets, where our model outperforms the random forest, feed-forward neural network, and mayr2018chemblmayr2018chembl models on 15 out of the 16 datasets. The strong performance of our model over these baselines, many of which use computed fingerprints or descriptors, demonstrate that learned molecular representations are indeed ready for “prime time” use in industrial property prediction settings.
Nevertheless, several avenues for future research remain. When analyzing the performance of our D-MPNN, we found that it typically underperforms when either 1) the other models incorporate 3D information, as in MoleculeNet’s best QM and PDBbind models, 2) the dataset is especially small, as in the PDBbind-C dataset, or 3) the classes are particularly imbalanced, as in the MUV dataset. One avenue of improvement is the incorporation of additional 3D information into our model, which currently includes only a very restricted and naive representation of such features. Another potential improvement is a principled pretraining approach, which some authors have already begun to explorenavarin2018pre, goh2018using. Such an approach could enable models to transfer learning from large chemical datasets to much smaller datasets, thereby improving performance in limited data settings. Another direction for future research is to determine how to adapt models and training algorithms to classification datasets with extreme class imbalance. Finally, in addition to these potential improvements, our analysis of how estimation of model generalizability is affected by split type opens the door to future work in uncertainty quantification and domain of applicability assessment.
We thank the Machine Learning for Pharmaceutical Discovery and Synthesis (MLPDS) consortium, Amgen Inc., BASF, and Novartis for funding this research. The MIT authors are funded by MIT, in particular the MLPDS consortium; the authors at companies are funded by their respective organizations. This work was supported by the DARPA Make-It program under contract ARO W911NF-16-2-0023. The authors declare no competing financial interest.
We would like to thank the other members of the computer science and chemical engineering groups in the Machine Learning for Pharmaceutical Discovery and Synthesis consortium for their helpful feedback throughout the research process. We would also like to thank the other industry members of the consortium for useful discussions regarding how to use our model in a real-world setting. We thank Nadine Schneider and Niko Fechner at Novartis for helping to analyze our model on internal Novartis data and for feedback on the manuscript. We thank Ryan White, Stephanie Geuns-Meyer, and Florian Boulnois at Amgen Inc. for their help enabling us to run experiments on Amgen datasets. In addition, we thank Lior Hirschfeld for his work on our web-based user interfacechemprop_website. Finally, we thank Zhenqin Wu for his aid in recreating the original data splits of Wu_2018Wu_2018, and we thank Andreas Mayr for helpful suggestions regarding adapting the model from mayr2018chemblmayr2018chembl to new classification and regression datasets.
In our Supporting Information, we provide links to our code and to a demonstration of our web-based user interface. We also provide further comparisons of model performance on both scaffold-based and random splits of the data, and we provide tables with all raw performance numbers (including p-values) which appear in the charts in this paper. In addition, we analyze the class balance of classification datasets, and we provide a list of the RDKit calculated features used by our model.