Bayesian semi-supervised learning for uncertainty-calibrated prediction of molecular properties and active learning

TL;DR

This paper introduces Bayesian semi-supervised graph convolutional neural networks for predicting molecular properties, emphasizing uncertainty quantification and active learning to improve reliability and data efficiency in drug discovery.

Contribution

It presents a novel Bayesian semi-supervised approach that enhances uncertainty estimation and active learning capabilities in molecular property prediction.

Findings

Uncertainty estimates are statistically principled and reliable.

Semi-supervised learning improves performance with limited data.

Active learning accelerates discovery by selecting informative samples.

Abstract

Predicting bioactivity and physical properties of small molecules is a central challenge in drug discovery. Deep learning is becoming the method of choice but studies to date focus on mean accuracy as the main metric. However, to replace costly and mission-critical experiments by models, a high mean accuracy is not enough: Outliers can derail a discovery campaign, thus models need reliably predict when it will fail, even when the training data is biased; experiments are expensive, thus models need to be data-efficient and suggest informative training sets using active learning. We show that uncertainty quantification and active learning can be achieved by Bayesian semi-supervised graph convolutional neural networks. The Bayesian approach estimates uncertainty in a statistically principled way through sampling from the posterior distribution. Semi-supervised learning disentangles representation learning and regression, keeping uncertainty estimates accurate in the low data limit and allowing the model to start active learning from a small initial pool of training data. Our study highlights the promise of Bayesian deep learning for chemistry.

Tables

Table 1. Table 1: The performance of the baseline XGBoost model on the datasets considered in this paper.

Dataset	Number of data points	${𝐑}^{\mathrm{𝟐}}$	RMSE
$\mathrm{𝙵𝚛𝚎𝚎𝚂𝚘𝚕𝚟}$	643	0.765	1.90
$\mathrm{𝙴𝚂𝙾𝙻}$	1128	0.704	1.14
$\mathrm{𝙲𝚊𝚝𝚂}$	595	0.654	0.367
$\mathrm{𝙼𝚎𝚕𝚝𝚒𝚗𝚐𝙿𝚘𝚒𝚗𝚝}$	3025	0.725	50.2
$\mathrm{𝚙𝟺𝟻𝟶}$	8817	0.291	0.996
$\mathrm{𝙼𝚊𝚕𝚊𝚛𝚒𝚊}$	13417	0.499	0.655

Table 2. Table 2: Combining semi-supervised learning with with SVGD yields uncertainty estimates that are strongly correlated with actual model error. The table shows the Spearman correlation coefficient between the variance predicted by the model and the error on test data points.

Dataset	Graph convolution with Dropout	Semi-supervised with Dropout	Semi-supervised with SVGD
$\mathrm{𝙵𝚛𝚎𝚎𝚂𝚘𝚕𝚟}$	$0.531\pm 0.061$	$0.439\pm 0.093$	$0.688\pm 0.053$
$\mathrm{𝙴𝚂𝙾𝙻}$	$0.112\pm 0.035$	$0.306\pm 0.079$	$0.553\pm 0.026$
$\mathrm{𝙲𝚊𝚝𝚂}$	$0.049\pm 0.036$	$0.066\pm 0.044$	$0.310\pm 0.019$
$\mathrm{𝙼𝚎𝚕𝚝𝚒𝚗𝚐𝙿𝚘𝚒𝚗𝚝}$	$0.192\pm 0.016$	$0.284\pm 0.035$	$0.337\pm 0.013$
$\mathrm{𝚙𝟺𝟻𝟶}$	$0.167\pm 0.015$	$0.185\pm 0.049$	$0.213\pm 0.010$
$\mathrm{𝙼𝚊𝚕𝚊𝚛𝚒𝚊}$	$0.315\pm 0.028$	$0.317\pm 0.031$	$0.378\pm 0.019$

Equations

$${\mathbf{m}}_v^t=\sum _{w\in \mathcal{N}(\mathbf{v})}({\mathbf{h}}_w^{t-1},{\mathbf{e}}_{vw}),$$

$${\mathbf{h}}_v^t=\sigma ({\mathbf{H}}_{t-1}^{\mathrm{deg}(v)}{\mathbf{m}}_v^t),$$

$$y=f\left(\sum _{v,t}\mathrm{softmax}({\mathbf{W}}_t{\mathbf{h}}_v^t)\right)$$

$$L=\sum _{m\in \mathcal{M}}\sum _{v\in m}\sum _{t=1}^T\log P({\mathbf{h}}_v^t|{\mathbf{u}}_m),$$

$$P({\mathbf{h}}_v^t|{\mathbf{u}}_m)=\frac{\exp \left(({\mathbf{h}}_v^t{)}^T{\mathbf{u}}_m\right)}{{\sum }_{n\in \mathcal{M}}\exp \left(({\mathbf{h}}_v^t{)}^T{\mathbf{u}}_n\right)},$$

$$y_i=F({\mathbf{x}}_i,\mathbf{\theta })+{ϵ}_i.$$

$$P(\mathbf{\theta }|\{{\mathbf{x}}_i\},\{y_i\})=\frac{1}{Z}P(\{y_i\}|\mathbf{\theta },\{{\mathbf{x}}_i\})P(\mathbf{\theta }),$$

$$⟨\hat{y}⟩=\int P(\mathbf{\theta }|\{{\mathbf{x}}_i\},\{y_i\})F(\hat{\mathbf{x}},\mathbf{\theta })\mathrm{d}\mathbf{\theta }.$$

$$\mathrm{var}(\hat{y})=\int P(\mathbf{\theta }|\{{\mathbf{x}}_i\},\{y_i\})(⟨\hat{y}⟩-F(\hat{\mathbf{x}},\mathbf{\theta }){)}^2\mathrm{d}\mathbf{\theta }.$$

$$-\log P(\mathbf{\theta }|\{{\mathbf{x}}_i\},\{y_i\})=\sum _i\left(\frac{1}{2{\sigma }_i^2}(y_i-F(\hat{\mathbf{x}},\mathbf{\theta }){)}^2+\frac{1}{2}\log {\sigma }_i^2\right)-\log P(\mathbf{\theta }),$$

$$q(\mathbf{\theta })=\prod _{m=1}^M{\mathbf{\Theta }}_{\mathbf{m}}{z}_m,z_m\sim \mathrm{Bernoulli}(p),$$

$$⟨y⟩=\frac{1}{N}\sum _{m=1}^N{y}^m,\mathrm{var}(y)=\frac{1}{N}\sum _{m=1}^N(y^m-⟨y⟩{)}^2+\frac{1}{N}\sum _{m=1}^N({\sigma }^m{)}^2.$$

$${\mathbf{\theta }}_i^{t+1}={\mathbf{\theta }}_i^t+\eta \varphi ({\mathbf{\theta }}_i^t),$$

$$\varphi (\mathbf{\theta })=\frac{1}{N}\sum _{j=1}^N\left[k({\mathbf{\theta }}_j^t,\mathbf{\theta }){\nabla }_{{\mathbf{\theta }}_j^t}\log P({\mathbf{\theta }}_j^t|\{{\mathbf{x}}_i\},\{y_i\})+{\nabla }_{{\mathbf{\theta }}_j^t}k({\mathbf{\theta }}_j^t,\mathbf{\theta })\right].$$

Full text

I Introduction

Predicting physiological properties and bioactivity from molecular structure – quantitative structure-property relationships (QSPR) – underpins a large class of problems in drug discovery. Classical QSPR workflows Cherkasov et al. (2014) separate descriptor generation – mapping a 2D Randić (1991); Ivanciuc (2000); Durant et al. (2002); Rogers et al. (2005) or 3D molecular structure Cramer et al. (1988); Verma et al. (2010) into a vector of real numbers using some handcrafted rules – and machine learning method that connects descriptors to property. Pioneering advances in machine learning such as graph neural networks directly take a molecular graph as input and infer the optimal structure-to-descriptor map from data Scarselli et al. (2009); Duvenaud et al. (2015), outperforming classical machine learning methodologies with handcrafted descriptors Wu et al. (2018).

Nonetheless, graph neural networks are usually developed using frequentist maximum likelihood inference, with the benchmark being the mean error on a test set. However, if the goal of QSPR is to replace mission-critical but expensive experiments, a low mean error is insufficient: the user needs to have an estimate of uncertainty and know when the model is expected to fail. This is because typically only a small number of top-ranked predictions are selected to test experimentally, thus outliers can ruin a discovery campaign. Moreover, cost limits the number of experiments that can be run, thus an approach that judiciously designs the training set to maximise information gained is needed.

Uncertainty quantification, or domain applicability, has been extensively considered in the QSPR literature but not in the context of graph neural networks and not in a statistically complete way. Previous works estimate uncertainty of prediction as the distance in descriptor space between the input molecule and the training set, or training an ensemble of models and evaluating the variance Sheridan et al. (2004); Sushko et al. (2010); Sheridan (2012); Toplak et al. (2014). More recent works consider conformal regression Norinder et al. (2014); Svensson et al. (2017), which trains two models, one for the molecular property and one for the error. However, there are two sources of uncertainty: Epistemic uncertainty arises due to insufficient data in the region of chemical space that the model is asked to make predictions on. Aleatoric uncertainty arises due to noise in the measurements themselves (e.g. noisy biochemical assays) Kendall and Gal (2017). Distance to the training set and variance within a model ensemble approximately capture epistemic uncertainty, whilst employing an ancillary model for prediction error approximately captures aleatoric uncertainty. We will show that the Bayesian statistical framework captures both sources of uncertainty in a unified and statistically principled manner.

Active learning strategies have been considered in the drug discovery literature Reker and Schneider (2015); Reker et al. (2017). However, those pioneering works considered a priori defined molecular descriptors, and estimate uncertainty via variance within an ensemble of models. Notwithstanding the shortcomings with incomplete modelling of uncertainty discussed above, employing graph neural networks in active learning presents unique opportunities and challenges: High model accuracy in the big data limit comes at the cost of being data-hungry. As the descriptor is fully data-driven, the model cannot estimate how “far” a compound is from the test set in the low-data limit, leading to poor uncertainty estimate and breaking down the active learning cycle. Low-data drug discovery has been considered in the context of one-shot learning Altae-Tran et al. (2017) which estimates distance in chemical space by pulling data from related tasks. Nonetheless, this approach requires a priori knowledge on which tasks are related. Works on generative molecular design overcome this problem Gómez-Bombarelli et al. (2018) by starting the active learning cycle with $<1000$ quantitative measurements, which impose a significant upfront experimental cost.

In this paper, we combine Bayesian statistics – a principled framework for uncertainty estimation – with semi-supervised learning which learns the representation from unlabelled data. We show that Bayesian semi-supervised graph convolutional neural networks can robustly estimate uncertainty even in the low data limit and drive an active learning cycle, and overcome dataset bias in the training set. Further, we demonstrate that the quality of posterior sampling is directly related to accuracy of the uncertainty estimates. As different Bayesian inference methods can be mixed and matched with different models, our study opens up a new dimension in the design space of uncertainty-calibrated QSPR models.

II Methods and Data

A machine learning method has two independent components: model and inference. The model is a function with parameters that relate the input to the output. Inference pertains to the methodology by which the model parameters are inferred from data. In terms of model, we focus on graph convolutional neural network models that take molecular graphs as input. In terms of inference, we focus on the Bayesian methodology.

II.1 Supervised graph convolutional neural network

Our baseline model is the graph convolutional fingerprint model Duvenaud et al. (2015). The salient idea is the message passing operation Gilmer et al. (2017), which creates a vector that summarises the local atomic environment around each atom while respecting invariance with respect to atom relabelling. A molecule is described by a graph, where the nodes are atoms and the edges are bonds. Atom $i$ is described by a vector of atomic properties $\mathbf{x}_{v}$, and a bond connecting $i$ and $j$ is described by bond properties $\mathbf{e}_{vw}$. The algorithm is iterative: at step $t$, each atom has a hidden state $\mathbf{h}_{v}^{t}$, which depends on “messages” $\mathbf{m}_{v}^{t}$ received from surrounding atoms as well as $\mathbf{h}_{v}^{t-1}$. The hidden states can be interpreted as descriptors of local atomic environment, and the messages allow adjacent atoms to comprehend the environment of its neighbours. Each atom is initialised to its atomic features, $\mathbf{h}_{v}^{0}=\mathbf{x}_{v}$, and

[TABLE]

and

[TABLE]

where $\mathcal{N}(\mathbf{v})$ denotes the set of atoms bonded to atom $v$, $\sigma(\cdot)$ is the sigmoid function, $\mathbf{H}_{t}^{N}$ is a learned matrix for each step $t$ and vertex degree $N$. The algorithm is run $T$ times, with $T$ being a hyperparameter. In the final step, the output is given by a multilayer neural network $f(\cdot)$ that takes a weighted average of the hidden states at each step as input return a prediction,

[TABLE]

where $\mathbf{W}_{t}$ are learned readout matrices, one for each step $t$.

We use the implementation reported in the repository 111$\mathtt{https://github.com/debbiemarkslab/neural-fingerprint-theano}$. In all experiments, we consider $T=3$, hidden layer at each level has 128 units, fingerprint length 256 (i.e. $\mathbf{W}_{t}\in\mathbb{R}^{128\times 256}$), and $f(\cdot)$ is a two layer neural network with 128 units each and $relu$ units.

II.2 Semi-supervised graph convolutional neural network

The fully supervised approach learns molecular descriptors directly from data. This is an advantage if one has a lot of data but a disadvantage in the data-limited settings such as active learning applications, where the objective is to design informative experiments starting from a small pool of initial training data.

The insight behind the semi-supervised approach is that significant amount of chemical knowledge is contained within the molecular structures themselves, without any associated molecular properties (i.e. unlabelled data). Thermodynamic stability puts constrains on what bonds are possible, and tends to put certain bonds near each other, forming persistent chemical motifs. For example, just by looking at drug molecules, one would immediately spot ubiquitous motifs such as amide group, benzene rings etc., and some motifs often occur together as scaffolds Welsch et al. (2010). The key assumption is that those persistent chemical motifs contribute to the molecular property that we want to predict. We can make mathematical progress by constructing a descriptor akin to Equation (1) - (3). However, the objective is no longer trying to fit a particular property. Rather, the hidden states $\mathbf{h}_{v}^{t}$, which summarises the atomic environment around atom $v$ within radius $t$, are constructed such that they are predictable from the hidden states of the surrounding atoms. Therefore, the model learns a descriptor that clusters similar environments.

Specifically, we use the semi-supervised approach developed by Ngyen et al. Nguyen et al. (2017), which builds on the paragraph vector approach in natural language processing Le and Mikolov (2014). Given a set of molecular structures $\mathcal{M}$, the hidden states $\mathbf{h}_{v}^{t}$ maximise the log-likelihood

[TABLE]

[TABLE]

where $\mathbf{u}_{n}$ is the molecular identifier, obtained by maximising Equation (4), with $\mathbf{h}_{v}^{t}$ defined by Equation (1)-(2). We can interpret $\mathbf{u}_{n}$ as a vector that describes the “type” of molecule, and the objective encourages the hidden states $\mathbf{h}_{v}^{t}$ to take values such that similar molecules have similar atomic environments.

After finding parameters that maximise the objective (4), $\{\mathbf{h}_{v}^{t}\}$ are then passed to a neural network, Equation (3). The parameters of the neural network as well as the readout matrices $\mathbf{W}_{t}$ are learned in a supervised manner. Note that this formalism infers descriptors using unsupervised learning and uses supervised learning to relate descriptors to molecular properties.

We use the implementation reported in the Github repository222$\mathtt{https://github.com/pfnet-research/hierarchical-molecular-learning}$ accompanying ref Nguyen et al. (2017). In all experiments, we consider $T=3$, hidden layer at each level has 128 units, fingerprint length 256 (i.e. $\mathbf{W}_{t}\in\mathbb{R}^{128\times 256}$), and $f(\cdot)$ is a two layer neural network with 128 units each and $relu$ units.

II.3 Bayesian deep learning

In Bayesian inference, the aim is to determine the distribution of model parameters that conforms to the data, the so-called posterior distribution. Let $\boldsymbol{\theta}$ be model parameters, $\mathbf{x}_{i}$ the dependent variables and $y_{i}$ the independent variable, such that

[TABLE]

where $\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}_{i})$ is the measurement noise. Bayes theorem states the posterior distribution, $P(\boldsymbol{\theta}|\{\mathbf{x}_{i}\},\{y_{i}\})$, is related to the likelihood, $P(\{y_{i}\}|\boldsymbol{\theta},\{\mathbf{x}_{i}\})$ and the prior $P(\boldsymbol{\theta})$ via

[TABLE]

where $Z$ is a normalising constant. The prediction for an unknown input $\mathbf{\hat{x}}$ is obtained by averaging over the posterior

[TABLE]

The uncertainty of model predictions can be readily derived from this Bayesian formalism. There are two types of uncertainties Kendall and Gal (2017). First, the epistemic uncertainty, is given by the variance of the prediction with respect to the posterior

[TABLE]

Second, the aleatoric uncertainty, is the intrinsic noise of the measurement $\sigma_{i}^{2}$. This aleatoric noise can depend on the input, $\sigma_{i}^{2}=\sigma(\mathbf{x})^{2}$, as certain areas of the chemical space can be intrinsically more variable.

We note that the log posterior is, up to a constant,

[TABLE]

which is exactly the mean-squared loss if $\sigma_{i}$ is constant, with $\log P(\boldsymbol{\theta})$ being the regulariser. Therefore, maximum likelihood inference is a special case of Bayesian inference.

The Bayesian formalism is easy to state but computationally expensive. The numerical bottleneck is the numerical evaluation of the high dimensional integrals (8)-(9). A plethora of approximate numerical methods have been developed in the literature to overcome this bottleneck. However, there is no free lunch, and methods which approximate the posterior well are usually computationally expensive. In this paper, we will consider two approximate methods spanning the cost-accuracy spectrum.

II.3.1 Dropout variational inference

Variational inference seeks to approximate the posterior distribution by a distribution that is much easier to sample from. Refs Gal and Ghahramani (2015); Kendall and Gal (2017) show that a popular way to regularise neural networks – dropout – is equivalent to approximate Bayesian inference. The algorithm is simple: The neural network is forked at the last layer to have two outputs, the predicted aleatoric uncertainty $\sigma_{i}^{2}$ and dependent variable $y_{i}$, and trained to minimise the loss (10). However, during training, each unit has a probability $p$ of being set to 0.

For a neural network with $M$ units, refs Gal and Ghahramani (2015); Kendall and Gal (2017) show that the above algorithm is approximately equal to finding parameters $\boldsymbol{\theta}=(\boldsymbol{\Theta_{1}},\boldsymbol{\Theta}_{2}\cdots\boldsymbol{\Theta}_{M})$ that fit the distribution

[TABLE]

to the posterior distribution $P(\boldsymbol{\theta}|\{\mathbf{x}_{i}\},\{y_{i}\})$, where $\Theta_{m}$ is the parameter vector associated with the $m^{th}$ unit.

Distribution (11), although not the same as the true posterior distribution, is significantly easier to sample: In the prediction phase, the model is run $N$ times, and akin to the training phase each unit has probability $p$ of being set to 0. The final prediction and total uncertainty is taken to be the mean over $N$ different predictions of depending variable and variance, $\{y^{i},(\sigma^{i})^{2}\}_{i=1}^{N}$,

[TABLE]

The first term in Equation (12) is the epistemic uncertainty and the second term is the aleatoric uncertainty.

In our numerical experiments, dropout is applied to every unit that is trained using supervised learning, i.e. every unit in the supervised graph convolutional neural network is subjected to dropout, whereas for the semisupervised case the layers on top of the hidden states are trained with dropout.

II.3.2 Stein Variational Gradient Descent

Rather than fitting a distribution to the posterior, Stein Variational Gradient Descent (SVGD) Liu and Wang (2016) directly draws samples from the posterior via gradient descent. Specifically, let $\{\boldsymbol{\theta}^{0}_{i}\}_{i=1}^{N}$ be parameters randomly and independently initialised in parameter space. We want to evolve parameters such that, after $T$ steps, $\{\boldsymbol{\theta}^{T}_{i}\}_{i=1}^{N}$ are $N$ independent samples drawn from $P(\boldsymbol{\theta}^{t}_{j}|\{\mathbf{x}_{i}\},\{y_{i}\})$. Ref Liu and Wang (2016) shows that the following dynamical system does the trick:

[TABLE]

where

[TABLE]

and $k(\cdot,\cdot)$ is a generic kernel function and $\eta$ is the learning rate. Equations (13)-(14) can be interpreted as free energy minimisation of an interacting particle system: a “particle” (parameter vector) is subjected to a “force” $\phi(\boldsymbol{\theta})$, which drives particles to regions of low energy (low loss), whilst forcing the particles apart to maximise entropy. The total uncertainty is evaluated also with Equation (12), except $\{y^{m}\}_{m=1}^{N}$ are predictions from different model parameters $\{\boldsymbol{\theta}_{i}\}_{i=1}^{N}$.

The key advantage of Equation (13) - (14) is that it is a well-defined approximation: Frequentist inference (c.f. Equation (10)) is recovered if $N=1$, whereas when $N\rightarrow\infty$ the system exactly samples from the posterior. Therefore, for finite $N$, the algorithm interpolates between frequentist and full Bayesian inference. The computational cost and memory demands increase with $N$, and in this paper we use $N=50$.

To illustrate the computational demands of Stein Variational Gradient Descent, Figure 1 shows the wall clock time, on a Nvidia P100 GPU, as a function of the number of gradient updates steps for graph convolution with dropout, semi-supervised with dropout, and semi-supervised with SVGD on the melting point dataset discussed below. Both SVGD and Stochastic Gradient Descent use back-propagation to optimize the neural network parameters. For models trained using Stochastic Gradient Descent, the computational complexity of back-propagation at each iteration is $O(BM)$, where $B$ is the number of training samples at each iteration and $M$ is the number of parameters in the model. The semi-supervised model has less parameters than the fully supervised model, thus the wall-clock time is less per iteration. In SVGD, we need to update $N$ Stein particles per iteration, thus wall clock time per iteration scales as $O(BMN)$.

II.3.3 Architectures and hyperparameters

As the objective of this paper is to demonstrate the types of chemical problems that Bayesian deep learning can tackle, we adopt common parameters for graph convolutional neural networks taken from the literature rather than performing extensive hyperparameter optimisation and neural architecture search. For both supervised and semi-supervised graph convolutional neural networks, we keep the number of hidden layers the same as the original implementation in the GitHub repositories cited above. Following the implementations in the repositories we keep the dimension of the fingerprint twice of $n_{h}$, the number of neurons in the hidden layers, and only optimise the $n_{h}$ by a grid search over $\{32,64,128,256\}$, choosing the value of $n_{h}$ with the best averaged 5-fold cross-validation root mean squared error over all the datasets. This leads to the architecture of $T=3$, $N=50$, hidden units = 128 and two layers.

II.4 Datasets

We consider a set of common regression benchmarks for physical properties prediction and bioactivity prediction. The $\mathtt{MeltingPoint}$ dataset is a collection of 3,025 melting point measurements of drug-like molecules used in a benchmark study Coley et al. (2017). The $\mathtt{ESOL}$ dataset is a set of 1,128 measured aqueous solubilities of organic molecules Delaney (2004), and the $\mathtt{FreeSolv}$ dataset is a set of 643 hydration free energy measurements small molecules in water Mobley and Guthrie (2014). The $\mathtt{ESOL}$ and $\mathtt{FreeSolv}$ datasets are used in the MoleculeNet benchmark Wu et al. (2018). The $\mathtt{CatS}$ dataset comprises half-maximal inhibitory concentration $\log\mathrm{IC}_{50}$ measurements of 595 molecules against Cathepsin S, taken from the D3R Grand Challenge 3 and 4 Gaieb et al. (2019). The $\mathtt{Malaria}$ dataset is a set of in vitro half-maximal effective concentration ($\log\mathrm{EC}_{50}$) measurement of 13,417 molecules against a sulfide-resistant strain of P. falciparum, the parasite that causes malaria Gamo et al. (2010) used in benchmark Duvenaud et al. (2015). The $\mathtt{p450}$ dataset is a dataset of half-maximal effective concentration measurements of 8,817 molecules against Cytochrome P450 3A4, a key enzyme for metabolism and clearance of xenobiotics, taken from the PubChem assay AID 1851.

To give the reader a sense of how “hard” the different datasets are, we consider a simple baseline model of XGBoost Chen and Guestrin (2016) on ECFP6 fingerprints Rogers and Hahn (2010). We split the data into 80/10/10 (training/validation/testing). We take $\mathrm{MaxDepth}=5$, $\mathrm{LearningRate}=0.01$, and optimise the number of estimators ($\mathrm{nEstimators}=[50,100,150,200,250]$) using the validation set. Table 1 shows the coefficient of determination $R^{2}$ and root mean squared error (RMSE). Judging from the coefficient of determination, the dataset difficulty is (from easiest to hardest) $\mathtt{FreeSolv}<\mathtt{Melting}<\mathtt{ESOL}<\mathtt{CatS}<\mathtt{Malaria}<\mathtt{p450}$.

III Results

III.1 Uncertainty quantification

We first consider how well can the model estimate its own uncertainty given the full dataset, split into training (80%) and test (20 %) sets. The quality of the uncertainty estimate is operationalised by asking what is the model accuracy when the most uncertain predictions are removed, with uncertainty quantified by the variance computed from Equation (12). Our baseline method is graph convolution with dropout, which has recently been implemented in the DeepChem package as a feature Wu et al. (2018), although to our knowledge this is the first study that benchmarks Bayesian graph convolutional neural networks in terms of uncertainty quantification.

Figure 2 shows that semi-supervised learning with Stein Variational Gradient Descent accurately estimates uncertainty and significantly outperforms the baseline on every dataset. The plots show how the test set error varies as a function of confidence percentile – i.e. what is the error if we only consider the top-$n$ % of compounds in the test set ranked by confidence (note that confidence is inverse of the uncertainty, quantified by Equation (12)); the shaded region is one standard deviation, estimated by analysing 5 random partitions of the data into training and test sets. In every case, the error is a decreasing function of model confidence, thus the model successfully estimates which predictions are likely to be correct and which predictions are outliers.

Another metric that we can evaluate is the shape of the confidence-error curve. For $\mathtt{ESOL}$ and $\mathtt{FreeSolv}$, the error is a steeply decreasing at the low confidence limit before plateauing, suggesting that most predictions are accurate but for a few outliers, which the Bayesian method can identify. The situation is different for $\mathtt{MeltingPoint}$, $\mathtt{Malaria}$ and $\mathtt{p450}$ – the error is slowly decreasing at the low confidence limit before sharply decreasing when it approaches the 100% confidence percentile limit (see also insets of Figure 2). This suggests that a few predictions are very accurate, and Bayesian method can pick out those accurate predictions amid many less accurate ones. We note that our Bayesian model is well-suited for virtual screening applications, where the challenge is ensuring that the top-ranked actives picked out by the model are indeed actives, since only a very small proportion of the compounds ranked will actually be screened experimentally (the ”early recognition problem”) Bender and Glen (2005); Truchon and Bayly (2007).

A lingering question whether the quality of the uncertainty estimate is due to a set of good descriptors (obtained via semi-supervised learning) or accuracy of the Bayesian methodology. Figure 2 also shows that replacing Stein Variational Gradient Descent with dropout significantly reduces model performance. At the same confidence percentile, Stein Variational Gradient Descent consistently outperforms dropout. This suggests that the quality of posterior sampling drastically impacts the quality of uncertainty estimation.

The quality of uncertainty estimates can also be gauged by the correlation between the predicted uncertainty on test data points and the error that the model incurs. Table 2 shows the Spearman correlation coefficient between the predicted variance and model error. As expected, combining semi-supervised learning with Stein Variational Gradient Descent leads to method with the highest rank correlation. This result is consistent with Figure 2, which shows that semi-supervised learning with Stein Variational Gradient Descent has the lowest confidence-error curve. Moreover, the rank correlation between predicted uncertainty and error broadly (although not exactly) follows the “difficulty” of the data (c.f. Table 1) – $\mathtt{FreeSolv}$ and $\mathtt{ESOL}$ have the highest rank correlation, and $\mathtt{p450}$ has the lowest.

Previous works on domain applicability have focused on either building an auxiliary model to predict the error Norinder et al. (2014); Svensson et al. (2017), or estimating the uncertainty of a prediction via the distance of the input to the training set Sheridan et al. (2004); Sushko et al. (2010); Toplak et al. (2014). The former models aleatoric uncertainty whereas the latter approximately captures epistemic uncertainty. Our Bayesian method captures both sources of uncertainty in a statistically principled manner. However, our model also provides independent estimates of epistemic and aleatoric uncertainties. As such, we can ask the question: is knowing epistemic or aleatoric uncertainty alone sufficient to estimate whether a prediction is accurate?

Figure 3 shows that the confidence-error curve for semi-supervised learning with Stein Variational Gradient Descent obtained by considering both epistemic and aleatoric uncertainty is below (or matches) that obtained by considering epistemic or aleatoric uncertainty alone. Considering both sources of uncertainty leads to much more accurate predictions at the high confidence limit for $\mathtt{ESOL}$ and $\mathtt{p450}$. Moreover, there is no consistent trend as to whether epistemic or aleatoric uncertainty is more important – for $\mathtt{ESOL}$, epistemic uncertainty is a better estimate of error than aleatoric uncertainty, whereas the opposite is true for $\mathtt{p450}$ and $\mathtt{CatS}$. As such, one cannot overlook epistemic or aleatoric uncertainty a priori, and our approach of combining both sources of uncertainty leads to an accurate uncertainty estimate.

III.2 Overcoming dataset bias

Our Bayesian methodology also overcomes dataset bias, which has be noted in the recent literature as the leading cause for overly optimistic results on benchmarks Wallach and Heifets (2018). Most ligand-based benchmarks are biased in the sense that the molecules reported are tightly clustered around a few important chemical scaffolds, such that when the dataset is randomly split into training set and test set, the molecules in the test set are structurally very similar to the training set. Therefore, a model that only memorises the training set will still achieve a high accuracy on the test set yet cannot generalise to other regions of chemical space. Methods such as scaffold splitting Wu et al. (2018) and attribution McCloskey et al. (2018) attempt to estimate what would be the true performance of the model if the dataset were not biased. However, bias is fundamental in chemical data – an ”uniform distribution” in chemical space does not exist because chemical space does not have a well-defined metric. Regardless of how one preprocesses the dataset or train the model, model predictions will always be awry for scaffolds that are not represented in the dataset. As such, rather than “unbiasing” the data, the practical question is whether the model can estimate whether it is likely to make a correct prediction for an unseen molecule given a biased training set.

To show that Bayesian uncertainty estimation overcomes dataset bias, we consider a toy problem where we know the ground truth and deliberately introduce bias: We consider the problem of predicting octanol-water partition coefficient ($\log P$) values, and use computed $\mathrm{ACD}\log P$ values as a surrogate. We construct a dataset comprising all molecules on ChEMBL with either a beta-lactam or a benzodiazepine scaffold. The dataset is obviously very biased as it only contains 2 scaffolds. We then train a model using Stein Variational Gradient Descent with semi-supervised learning, with the standard 8:2 split between training and test set on the biased dataset. Figure 4 (left) show that model is reasonably accurate on the biased test set. We now simulate how an user might unwittingly fall foul of dataset bias – suppose we use the model to predict $\log P$ of all molecules with a steroid scaffold on ChEMBL. Figure 4 (middle) show that the model performance, perhaps unsurprisingly, is poor. Steroids are not part of the training set, thus the model cannot predict its physiochemical properties. Bayesian uncertainty estimation provides a way out of this quandary – Figure 4 (right) shows that the estimated uncertainty of $\log P$ prediction on steroids is significantly greater than $\log P$ prediction on the test set of beta-lactams or a benzodiazepines. In other words, the model can inform the user when it is inaccurate, thus mitigating the impact of dataset bias.

III.3 Low data active learning

Having considered the quality of the uncertainty estimates in the data-abundant limit, our next question is whether we can estimate uncertainty in the low data limit and drive an active learning cycle. We consider the objective of obtaining a low model error with a small training set. The model is first trained from a small initial pool of data (25 % of the full training set, picked randomly), the model then selects a batch of molecules (2.5% of the full training set) that has the largest epistemic uncertainty to put into the training set, and then the model is retrained to suggest other additions, and the cycle continues. The test set is always 20% of the full dataset, held out at the beginning of the experiment. We note that other acquisition functions have been suggested in the literature Gal et al. (2017), and the objective function is problem-dependent Reker and Schneider (2015); Reker et al. (2017). Nonetheless, the goal of our experiment is to evaluate the quality of uncertainty estimate, thus we focus on a simple objective and acquisition functions. Further, as active learning requires constant retraining of the model, and Stein Variational Gradient Descent is significantly more computationally intensive than dropout, we will only consider dropout variational inference.

Figure 5 shows that semi-supervised learning significantly outperforms full supervised learning in the low data limit. The mean learning curves and error bars are obtained by analysing 20 active learning runs starting from random dataset splits. Moreover, in the case of full supervised learning, active learning is unable to deliver a better learning curve than random sampling, whereas for semi-supervised learning there is a sizeable gap between the learning curves of random sampling and active learning. This is because the full supervised method generates molecular descriptors directly from data. Therefore, in the low data limit, it is unable to learn descriptors that describe the structure of chemical space and chemical similarity between compounds, thus cannot generate meaningful uncertainty estimates to drive active learning.

The importance of choosing diverse compounds in the initial screen has been discussed extensively in the literature Huggins et al. (2011); Bakken et al. (2012); Paricharak et al. (2018), and the performance of our active learning method also depends on the chemical diversity in the initial screen. Figure 6 shows that active learning does not outperform random sampling when the initial training set biased and contain only a small number of scaffolds. We model scaffold bias by splitting the data using scaffold splitting implemented in $\mathtt{DeepChem}$ Wu et al. (2018), and consider the $\mathtt{Malaria}$ example where active learning most clearly outperforms random sampling in Figure 5. The underperformance of active learning is perhaps unsurprising – if the initial screen only consists of one scaffold, the knowledge that the model has on the other scaffolds would be minimal, i.e. the model is equally ignorant about all the other scaffolds. As such, randomly sampling the other scaffolds becomes a reasonable strategy.

IV Discussion and Conclusion

We propose a novel method to quantify uncertainty in molecular properties prediction. We show that our methodology significantly outperforms the baseline on a range of benchmark problems, both in terms of model accuracy and in terms of uncertainty estimates. Our method also overcomes dataset bias by returning a large uncertainty estimate when the test set is drawn from a different region of chemical space compared to the training set. Moreover, our methodology can drive an active learning cycle, maximising model performance while minimising the size of the training set. The key to the success of our method is the combination of semi-supervised learning and Bayesian deep learning. Semi-supervised learning allows us to learn informative molecular representation in the low data limit. Bayesian deep learning allows us to estimate aleatoric and epistemic uncertainty in a statistically principled manner. We exemplified our methodology on regression as it is generally more challenging than classification, although it can be readily extended to classification problems.

Our observation that the choice of Bayesian inference methodology significantly impacts the quality of the uncertainty estimate suggests an evident followup that probes the mathematical limit of Bayesian inference – i.e. benchmarking approximate inference techniques against importance sampling of the posterior using Markov Chain Monte Carlo till convergence, which is computationally expensive but mathematically exact. Moreover, we note that most approximate inference techniques in the literature have been benchmarked in terms of RMSE error or log-likelihood Mukhoti et al. (2018), rather than explicitly considering the quality of the uncertainty estimate in a manner relevant for chemoinformatics such as the confidence-error curve. An open question is the design of appropriate approximate inference techniques for graph convolutional neural networks that solves the trilemma between computational cost, model accuracy, and the quality of uncertainty estimate.

Another open question is whether the model has accurately disentangled aleatoric and epistemic uncertainty. Answering this question would require estimates of the ground truth aleatoric uncertainty, which is obtainable via repeating the experimental measurement and reporting the variance. Benchmark datasets which provide accurate experimental uncertainty estimates will be invaluable to advancing the Bayesian methodology.

Finally, our active learning methodology performs well when the initial screen covers diverse compounds. To successfully perform active learning on a scaffold-biased initial set, the model needs information on the bioactivity of those unseen scaffolds. We speculate that strategies such as multitask learning Ramsunda et al. (2017); Wenzel et al. (2019), which pools information from other cognate assays which have explored the unseen scaffolds, will be a fruitful avenue.

Acknowledgements.

AAL acknowledges the support of the Winton Programme for the Physics of Sustainability.