Assigning Confidence to Molecular Property Prediction

Introduction: Computational modeling has rapidly advanced over the last decades, especially to predict molecular properties for chemistry, material science and drug design. Recently, machine learning techniques have emerged as a powerful and cost-effective strategy to learn from existing datasets and perform predictions on unseen molecules. Accordingly, the explosive rise of data-driven techniques raises an important question: What confidence can be assigned to molecular property predictions and what techniques can be used for that purpose? Areas covered: In this work, we discuss popular strategies for predicting molecular properties relevant to drug design, their corresponding uncertainty sources and methods to quantify uncertainty and confidence. First, our considerations for assessing confidence begin with dataset bias and size, data-driven property prediction and feature design. Next, we discuss property simulation via molecular docking, and free-energy simulations of binding affinity in detail. Lastly, we investigate how these uncertainties propagate to generative models, as they are usually coupled with property predictors. Expert opinion: Computational techniques are paramount to reduce the prohibitive cost and timing of brute-force experimentation when exploring the enormous chemical space. We believe that assessing uncertainty in property prediction models is essential whenever closed-loop drug design campaigns relying on high-throughput virtual screening are deployed. Accordingly, considering sources of uncertainty leads to better-informed experimental validations, more reliable predictions and to more realistic expectations of the entire workflow. Overall, this increases confidence in the predictions and designs and, ultimately, accelerates drug design.


I. INTRODUCTION
While data-driven modeling has made large advances in image processing and speech recognition, groundbreaking contributions in drug discovery are harder to find. Often, only limited amounts of reliable in vitro and in vivo data suitable for supervised learning tasks are available. Additionally, uncertainties of experimental data, for instance regarding effectiveness and toxicity, are common and significant as the experimental design is complex, prone to noise and the outcome is dependent on many parameters, in particular compound dosage. Hence, these uncertainties are propagated to property prediction workflows when data-driven models are trained on that data.
While it is evident that machine learning (ML) prediction accuracy and data quality are intrinsically related, there are more subtle aspects that are sometimes neglected.
Dataset size, composition, coverage and training-test split have considerable influence on the final model performance but are insufficient to assess uncertainty alone without suitable models and methodologies for that purpose. For instance, using random or scaffold-based train/test splits does not lead to reliable measures of model performance and generalizability. Other common problems include high dataset bias, small sample size or low chemical diversity. However, some of these issues currently do not have a straightforward remedy as there is no commonly accepted way to describe chemical subspaces rendering it difficult to uncover all inherent biases.
Nevertheless, when the underlying data is not investigated appropriately model performance and generalizability tend to be overestimated leading to a significant number of false predictions and, ultimately, low confidence in the property prediction workflows.
In this article, we discuss methodologies and datasets for chemical properties important in drug design. In particular, we focus on the main sources of uncertainties for both simulation-based and data-driven property prediction. In doing so, we discuss uncertainties inherent in datasets, outputs of data-driven models, input features, and simulation of binding affinities. In addition, we dis-cuss the importance of uncertainty in generative models, especially for property-based molecule design. Finally, we close with our personal opinion on the most important aspects of uncertainty and confidence focusing on important problems and future avenues that will lead to higher predictive ability and models that naturally take uncertainty into account.

II. DATASET UNCERTAINTY
One of the most common challenges for data-driven molecular property prediction is dataset size and composition. For many important properties in drug discovery, especially pharmacologic properties of molecules such as absorption, distribution, metabolism, excretion and toxicity (ADMET), only a limited amount of high-quality data is available and it is usually only available for certain classes of molecules, which inherently introduces biases. To provide an overview of the amount and type of data available, Table 1 lists popular datasets for molecular property prediction. Accordingly, in this section, we discuss dataset characteristics that are important to consider for data-driven property prediction tasks along with methods and procedures to minimize their impact on property prediction performance.
In the framework of supervised learning, models learn relationships between input and output from the training data and, subsequently, predict properties for unlabelled input data, and the corresponding prediction ability is termed generalization. When predictors perform well on the training set, but performance on the unlabelled data is poor, the model is overfitted. This is a very common problem in ML, especially when the training set size is small (Figure 1-A). Learning theory states that collecting more data improves generalization if the data is independent and identically distributed [20]. However, real data is not a uniform sample of chemical space. Typically, molecules in a dataset are collected under specific criteria such as the number of atoms, the constituent elements, the similarity to known molecules or the availability of synthetic procedures introducing bias (Figure 1-B). In addition, there is an inherent bias in both industry and academia to publish only successful experiments which can be problematic for assembling training data as negative results are as important as positive ones for data-driven modeling.
First, the uncertainty of model predictions for unlabelled molecules is largely dependent on the difference of the structure and property distributions to the training data. To quantify this difference, in the field of quantitative structure-activity relationships (QSARs) [21], the concept of applicability domain (AD) is widely applied. The AD of a QSAR model is defined as "the response and chemical structure space in which the model makes predictions with a given reliability" [22] and makes the inherent bias of a model comprehensible (Figure 1-C). There are several methods to estimate ADs [23]. A common approach is to consider molecules with descriptors within a certain distance to the mean of the training data to be inside the AD. Notably, predictions for molecules inside the AD are considered reliable, outside the AD reliability cannot be guaranteed.
Secondly, a significant bias in the training data can cause models only to learn the inherent bias rather than physically meaningful relationships (Figure 1-B). When users are unaware of this bias, it inadvertently leads to overly optimistic conclusions concerning model performance. Recently [10], it was pointed out that there is significant hidden bias in the widely used dataset for structure-based virtual screening, the Directory of Useful Decoys: Enhanced (DUD-E) [24]. DUD-E includes two types of molecules: ligands, small molecules that bind to a receptor, and decoys, molecules with similar physical properties but dissimilar chemical structures compared to the ligands. Using DUD-E as a dataset for training, the performance of receptor-ligand and ligand-only ML models were compared to determine whether the receptor-ligand models also learned from the actual interactions or whether they just learned the inherent ligand bias. Ideally, the performance of these two types of models should be different as the binding energies are a function of both ligand and receptor structures, and the receptor structure can generally not be inferred from the ligand alone. However, the performance was found to be equivalent, the receptor-ligand models did not learn from the receptor structure, and it was concluded that "models may only learn the inherent bias in the dataset rather than physically meaningful features" [10]. This calls for general methods to test what property predictors actually learn as it would naturally lead to more powerful, generalizable models and to more confidence in the models. To minimize hidden biases, the importance of dataset diversity for the generalizability of ML models has been pointed out before [25].
Finally, small datasets pose a particular challenge as they are inherently biased due to their limited size and, additionally, overfitting is to be expected. An important strategy to reduce both bias and overfitting in ML models when data is sparse is transfer learning (Figure 1-D) [23,26]. When using that approach, models are first pre-trained on large, sometimes unlabelled, datasets, and subsequently trained on the property of interest. Importantly, transfer learning has been shown to improve the prediction performance for models trained on small datasets significantly [27,28]. For instance, recently [27], a language-based ML property prediction model was pretrained first via self-supervised learning on one million unlabeled molecules curated from ChEMBL. This pretraining improved the prediction performance after subsequent training on small datasets of molecular activity consistently, regardless of the specific property prediction FIG. 1. Common problems with datasets for data-driven molecular property prediction and strategies to improve confidence in the models. A) Small training datasets hamper generalization. B) Biased training datasets lead to models only learning the inherent bias rather than meaningful relationships. C) The applicability domain of a model defines the molecules for which the property predictions are expected to be reliable. D) Transfer learning via pre-training on large (unlabelled) datasets improves prediction accuracy for property prediction based on small training datasets.
task and the training dataset size. Notably, the smallest dataset only contained 642 samples [27].

III. UNCERTAINTY IN THE OUTPUTS
Supervised ML largely relies on an inductive inference approach to derive general rules from a finite number of training examples. Therefore, predictive models obtained using inductive inference are never formally correct and inherently uncertain.
There are often additional sources of uncertainty that stem from noise or imprecision on target measurements. For instance, molecular properties like binding affinity are often measured with an uncertainty, e.g. 1 µM ± 10%. The latter type is referred to as aleatoric uncertainty, while the former is referred to as epistemic or systematic uncertainty and can in theory be reduced in light of more training observations or increased knowledge about the correct model. Both types of uncertainty are important, although they are often not treated separately and the formal distinction between the two remains an active area of research [29]. Many ML methods are available for representing uncertainty. They can be categorized with respect to the way uncertainty is represented and whether or not they allow for differentiation between the two types (aleatoric and epistemic). Here we divide the discussion of methods derived from frequentist statistics and ones based on Bayesian inference.
A widely used class of frequentist methods are ensemble strategies which typically use a large set of distinct predictors as opposed to a single model. Ensembles of tree-based models such as random forests [30] are commonly used to produce frequentist probability estimates based on the relative output frequency of the ensemble members [31,32]. Ensemble learning can introduce variance through multiple ensemble models with different parameter initializations (ensembling) or through the use of randomly sampled training sets (bagging or boosting, Figure 2-a). Ensembles of neural networks can also provide computationally efficient, Database of commercially available compounds together with very simple estimated molecular properties for virtual screening.

billion
Inherently biased by currently synthesizable chemical space. Consequently, the molecular shapes have been shown to be highly biased against sphere-like molecules.

thousand
Biased towards small molecules only containing the elements C, H, N, O and F.
PubChemQC [3,4] Geometries and electronic properties of molecules with short string representations taken from PubChem.

million
Biased towards small molecules that have been reported in the literature before.
Tox21 [5] Toxicologic properties of molecules with respect to 12 different assays 13 thousand Biased towards environmental compounds and approved drugs.

ToxCast [6]
High-throughput screening and computational data for the toxicology of molecules from industry, consumer products and the food industry based on cell assays.

thousand
Biased towards molecules used in industry, consumer products and the food industry.
ClinTox [7] Drugs and drug candidates that made it to clinical trials and were either approved or failed.

thousand
Biased towards drugs that made it to clinical trials.
SIDER [8] Recorded adverse drug reactions of marketed drugs. 1.4 thousand Biased towards marketed drugs.
ChEMBL [9] Bioactive small molecules and their activities extracted from the literature, from clinical trials and from other databases.

million
Biased towards compounds for which bioactivity was published in the scientific literature.

DUD-E [10]
Ligand binding affinities against 102 distinct target proteins with both strong and weak binders.

thousand
Biased towards molecules that have been synthesized and evaluated for binding affinity.
AqSolDB [11] Aqueous solubility data of organic molecules taken from 9 different datasets. 10 thousand Biased towards organic molecules with relatively high aqueous solubility.
Olfaction Prediction Challenge [11] Olfactory perception of organic molecules at different concentrations.

thousand
Biased towards small and volatile organic molecules. Results biased by familiarity of smells.
FreeSolv [12] Experimental and computed hydration free energies of small and neutral molecules. 0.6 thousand Biased towards small and neutral molecules that have been studied in the literature both computationally and experimentally for hydration free energies.
ESOL [13] Experimental aqueous solubility combining datasets for small molecules from the literature, for medium-sized molecules used as pesticides and larger proprietary compounds from the pharmaceutical industry.

thousand
The sub-groups each have a different bias as they each have different application domains.

thousand
Biased towards molecules with distribution coefficients between -10 and 10.

million
Biased towards molecules of interest and molecules that are synthesizable.

thousand
Biased towards complexes with available crystal structures.
BBBP [19] The blood-brain penetration partition coefficient for molecules collected from the literature.

thousand
Biased towards molecules studied in the literature for blood-brain penetration.
readily parallelizable uncertainty estimates [33,34]. However, care must be taken using these methods as the uncertainties are often biased, poorly calibrated, and are only concerned with aleatoric uncertainty. Heteroscedastic aleatoric uncertainty can also be modeled using a single neural network trained with frequentist maximum likelihood inference that outputs a probability distribution ( Figure 2-b) [35]. Conformal prediction (CP) has recently received significant attention, in which one ML model is trained to predict a property and another to predict the uncertainty (Figure 2 d) [36][37][38][39]. CP methods are based on a rigorous mathematical framework and allow users to select intuitive confidence levels for their predictions, e.g. selecting a confidence level of 0.9 means that at most 10% of predictions will be outside the predicted range.
Methods inspired by Bayesian inference seek to update a prior distribution p(θ) on the space of all possible models in light of training data, D, yielding a posterior distribution p(θ|D). Given a molecular representation for prediction, x * , the predictive posterior distribution p(y * |x * , D) is in theory obtained by averaging all possible models weighted by their respective posteriors. Bayesian methods have experienced a recent resurgence in drug design applications due to increasing computational Gaussian processes (GPs) [40] allow for Bayesian inference in a non-parametric way, where the prior is determined by a mean and kernel function. The prior may be updated with observed data and for regression with Gaussian noise, the predictive posterior distribution is Gaussian with mean µ and variance σ 2 . The predictive variance σ 2 represents the total uncertainty but can be decomposed into the variance of the error term σ 2 which represents aleatoric uncertainty. Epistemic uncertainty is the difference between σ 2 − σ 2 and is determined by the hyperparameters of the kernel function such as the characteristic length scale. GPs are commonly used tools for property prediction in the low data regime (i.e. below 1000 data points) due to their known robustness to overfitting [40][41][42]. Hie et al. [43] used a GP trained on less than 100 compounds to screen a large library for candidates with nanomolar affinity for diverse kinases and whole-cell growth inhibition of Mycobacterium tuberculosis [43]. The authors show how uncertainty estimation can enable machine-guided discovery. However, the traditional implementation of GPs incurs cubic scaling with training set size, rendering them less desirable when training data is abundant.
Bayesian neural networks (BNNs) place probability distributions on neural network parameters w and the analytically intractable posterior distribution p(w|D) is typically estimated using variational Bayesian inference [44]. The predictive posterior distribution is approximated using Monte Carlo sampling, which intuitively involves averaging the predictions of multiple parameter configurations sampled from the posterior (Figure 2-c).
Replacing point-estimated network parameters with distributions allows for the quantification of epistemic uncertainty as well as heteroscedastic aleatoric uncertainty. Recently, such networks have been combined with state-of-the-art representation learning to give uncertain predictions of the physiological properties of small molecules [45].
For classification tasks, confidence calibration seeks to produce probability estimates which correctly reproduce the true correctness likelihood. In other words, given N predictions each with confidence 0.9 about a label, we would expect 90% of the predictions to give the correct label. Scalar valued confidence statistics such as expected calibration error and maximum calibration error have also been proposed. Many calibration methods such as isotonic regression or temperature scaling are based on a learned post-processing step in which a calibrated probability q i is derived from the model's output probability p i * [46]. Confidence calibration schemes can be applied as a scalable post-processing step to any predictive learning algorithm, and their ideas can be extended to regression tasks.
Confidence or uncertainty estimation in molecular property prediction is crucial for tasks that involve data-driven decision making where valuable resources or human risk is at stake. For example, BNNs have been employed to give intuitive probabilistic statements about the severity of the risk of drug-induced liver injury [47,48]. Predicting clinical outcomes from pre-clinical data and using uncertainty estimates to quantify and communicate risks can preserve time, resources, and patient well-being [49]. Predictive uncertainty is central to the field of active learning [50], which has been known in the context of drug discovery for almost two decades. Active learning strategies balance exploitation of the knowledge of a predictive model to retrieve active candidate compounds with the exploration of regions in compound space where the model is highly speculative [45,51,52]. Explorative behaviour of an active learner is paramount when structural novelty is prioritized. Elucidation of error sources in training data remains a critical application for uncertain predictive models. Ryu et al. modeled aleatoric and epistemic uncertainty separately using Bayesian graph convolutional neural networks for a variety of molecular property prediction tasks including classification of bioactivity and toxicity. Analysis of predictive aleatoric uncertainty allowed for the identification of poor quality training data for some molecular candidates [53].
In the near future, we believe predictive uncertainty will become an indispensable component of molecular property prediction.
Researchers will continue to systematically benchmark uncertainty quantification methods on various datasets [54]. We anticipate software packages that facilitate uncertain predictions will continue to appear [55].

IV. UNCERTAINTY IN THE INPUT FEATURES
Contrary to uncertainty in the properties being predicted, imprecision in the input features is a source of uncertainty that is rarely considered. Property prediction in drug design often relies on representations encoding a specific molecule, such as SMILES strings [58], or on three-dimensional structural data. However, these static descriptions of molecular systems are not representative of their dynamic nature in solution. In this section, we discuss how input-related uncertainty might arise when predicting chemical properties from molecular structures, and how some of the latest ML models might be used to capture this uncertainty in the confidence of the predictions.
One form of input uncertainty is due to the presence of multiple protonation and tautomeric states in solution. String and graph representations of molecules are static and unable to capture these multiple states and their probabilities in solution. SMILES strings representing two different tautomers are likely to return different property predictions when they in fact should be the same given the behaviour of these molecules in solution. This represents a challenge to many compound property predictions and cheminformatics tools and may be seen as a form of input uncertainty, given that the chemical species in solution can be modeled only approximately. Similarly, protonation states depend on the chemical environment of the molecule, there may not be one dominant state in solution and, finally, these may change upon interaction with a binding partner. All these aspects contribute to creating uncertainty about the input representation used in ML prediction tasks.
Three-dimensional structural data is often used to predict properties like drug binding affinity, selectivity, and mutational effects on protein stability and drug resistance. This data is typically derived by X-ray crystallography and, increasingly, cryo-electron microscopy (cryo-EM). Despite providing abundant information on protein-ligand interactions, which can be used to drive drug design, these structures are only a single snapshot of all possible protein-ligand conformations in solution. The protonated form of histidine has a pKa of approximately 6.0, which means that it is generally expected to be found in its neutral form at physiological pH. Yet, because of the local environment in protein binding pockets, it is often unclear which of the two is the major species, and they may both co-exist in non-negligible proportions. On the right-hand side, some of the tautomeric states of warfarin [57] are shown.
As such, this three-dimensional data is a single datapoint of a much broader input probability distribution. Some structure-based modeling approaches, like ensemble docking, try to take this uncertainty into account by sampling multiple conformations with molecular dynamics simulations. In such a way, ensemble docking tries to capture the distribution/uncertainty of protein-ligand structures. By generating an "ensemble" of drug target confirmations, the docking results are less likely to be overfitted to a specific protein conformation. However, fast scoring functions for ligand binding affinities are typically trained on databases of static structures, such as the PDBbind database, and do not capture this type of uncertainty in their predictions. In addition to the statistical uncertainty caused by conformational variability, differences between protein structures in crystallo and in vivo due to packing effects and cryo temperatures may result in systematic biases.
Ideally, one would like uncertainty about multiple possible chemical and conformational states in solution to be captured by the ML model and reflected in the confidence of the property predictions. Currently, only a handful of specialized ML approaches are able to capture uncertainty in the input features. Gaussian process latent variable models, together with other similar kernelbased approaches [59,60], are examples of such models [61,62]. Recent work has also started to investigate the effect of molecular ensembles in property predictions [63] with ML models. However, we are unaware of widely adopted supervised ML models that can either take into account known input uncertainty, or that can infer it directly from data. We believe this to be an underexplored research direction with much potential for practical applications. Such models could become an indispensable component of molecular property prediction, being able to capture or infer feature uncertainty in structural and chemical data.

V. UNCERTAINTY IN THE BINDING AFFINITIES
A key property of interest in drug discovery is the binding affinity of a ligand to a target receptor. Of the many physics-based models that have been developed and refined to predict binding affinity, here we focus on two widely-used techniques that utilize the three-dimensional structures of ligand and receptor, and a molecular-mechanics energy function to compute ligand affinities: molecular docking and free energy perturbation (FEP) methods.
Large uncertainties in docking estimates can be attributed to conceptual deficiencies in the scoring function, for example, a lack of proper accounting of solvation thermodynamics, which may require incorporation of explicit water molecules against an ensemble of receptor poses to better represent the physical system [73]. Docking scores are also made uncertain by implementation choices made for algorithmic efficiency. For example, flexible ligand or receptor degrees of freedom are usually restricted to enable efficient conformational searching. Recent studies using deep learning neural-network-based scoring functions show promise in reducing prediction uncertainties [72,74]. Regardless of these inaccuracies, docking approaches are likely to remain an essential component of virtual screening, as they can reduce and enrich the chemical space of potential binders, and provide initial bound-pose conformations for seeding molecular dynamics (MD) simulations.
In contrast to docking, FEP methods are much more expensive, requiring all-atom MD simulations to perform statistical sampling of ligand-receptor configurations [75][76][77][78]. The benefit is a best-in-class prediction of ligand binding affinity, achieving accuracies within chemical accuracy, i.e. 1.0 kcal/mol [79,80]. Assigning uncertainties to such predictions, however, requires consideration of the many error sources these methods introduce. These can be categorized as arising from both sampling and scoring inaccuracies.

A. Uncertainty in FEP Estimates due to Sampling Error
To understand uncertainty from statistical sampling error, consider that using MD simulation naïvely to predict ligand affinity would require equilibrium sampling of multiple reversible binding transitions which typically exceed the timescales of state-of-the-art MD, even on special-purpose high-performance computers [81]. FEP and related methods circumvent this problem by performing simulations in multiple thermodynamic ensembles for a series of alchemical intermediates and analyzing the resulting samples with a statistical free energy estimator. Absolute binding free energy (ABFE) is computed along alchemical paths that decouple all non-bonded interactions in solution of the ligand, including the binding to a receptor. Relative binding free energy (RBFE) is computed along with an alchemical transformation from one molecular topology to another.
Statistical free energy estimators that utilize multiensemble sampling include exponential averaging (EXP) [82], thermodynamic integration (TI) [83], the weighted-histogram method (WHAM) [84], the Bennett acceptance ratio (BAR) method [85], the multistate Bennett acceptance ratio (MBAR) method [86], and newer dynamical estimators like TRAM [87] and DHAM [88]. Uncertainties from statistical sampling can influence the confidence of calculated Gibbs free energies in complicated ways. Hence, regardless of which free estimator is used, a non-parametric bootstrap is preferred to estimate uncertainty in Gibbs free energies [76]. Similarly, uncertainties can be calculated over block-averages of a long and equilibrated simulation trajectory, or can alternatively be estimated across multiple independent simulations, in some cases initiated with different ligand starting poses. For FEP, which involves computing free energy differences between the ligand-bound state and the free ligand in solution, overall uncertainties are often obtained by summing individual uncertainties for each ensemble in quadrature [71,76].
Different free energy estimators have unique bias and variance due to finite sampling, but they all depend on having sufficient thermodynamic overlap between neighbouring ensembles along the thermodynamic path. For this reason, relative free energy estimation is often more accurate for pairs of similar molecules, since the alchemical transformation involves a smaller number of atoms, and hence inherently has a better overlap. Achieving good thermodynamic overlap is the main consideration guiding choices about the number and type of alchemical intermediates, which may require position or angular restraints on the ligand to focus sampling on particular ligand poses [71].
Most free-energy estimators assume that samples are statistically independent and drawn from the equilibrium distribution. Sufficient equilibration can be notoriously slow to achieve for some systems, e.g. ligands with "hidden barriers" between conformational states, a problem which Hamiltonian replica exchange and other enhanced sampling methods have been designed to address [16].
To avoid statistically correlated input, which results in artificially low uncertainty estimates, recommended strategies discard initial unequilibrated samples, and subsample snapshots at a frequency (2τ + 1), where τ is a computed correlation time [89]. An algorithm to automatically detect simulation equilibration in this way has been developed [90].
Related to FEP is a class of non-equilibrium work (NEW) free energy methods, based on the Crooks fluctuation theorem [91], which drives a system between alchemical endpoints and analyzes the distribution of non-equilibrium work values [77]. These methods can achieve comparable accuracies to FEP with fewer intermediates [92][93][94], yielding uncertainties that depend on the variance of the work distributions, in turn depending on the timescale of the alchemical schedule. Sampling uncertainty can also be introduced by artifacts from the implementation of MD. For example, uncertainties on thermodynamic and kinetic properties can depend greatly on solvation artifacts that depend on simulation box size [95]. Furthermore, surprisingly, variation in computed energies across different simulation packages [96], integrators, and thermostats may be non-negligible [97].

B. Uncertainty in FEP Estimates due to Scoring Error
The workhorse of FEP is a molecular mechanics energy function, called a force field, which has usually been designed to approximate a quantum mechanical potential energy surface, and tuned to reproduce experimental observables. Iterative cycles of force field development, parameterization and assessment are thus vital for reducing the uncertainty introduced by force field inaccuracy [94,98,99]. With new machine learning advances have come new ways to improve force fields for the prediction of ligand binding free energies. One product of the recently-launched Open Force Field Initiative is a new generation of force fields that eschews atom-types in favour of direct chemical perception, utilizing the SMIRKS-native Open Force Field (SMIRNOFF) scheme for parameter assignment [100]. Parameterization of Open Force Field v1.0.0, code-named Parsley [38], was made possible by the ForceBalance algorithm, a Bayesian algorithm for automatic and reproducible parameterization [101][102][103]. Polarizable force fields such as AMOEBA continue to make gains in accuracy for predictions of Gibbs free energy of binding [104,105], and there are examples of neural network-based potentials, which aim to reproduce molecular properties with DFT accuracy at force field cost [106,107].
Finally, uncertainty can also result from the assumptions of a selected model. For example, most fixed-charge force fields require the pre-selection of particular protonation states for both ligand and ionizable side chains of a receptor. Depending on the predicted pK a of acidic or basic groups, which may differ significantly between bound and unbound states, multiple simulations may be required to explore these possibilities. Overall, there are several factors to consider when designing a study to compute binding affinities for drug-like molecules. New methodologies and best practices can often be found through blind assessment challenges such as SAMPL [108] (Statistical Assessment of Modeling of Proteins and Ligands), which are created annually to assess predictions of binding affinities among other properties important to drug discovery.

VI. UNCERTAINTY IN GENERATIVE MODELING
Generative models can produce synthetic data samples, some of them by learning from reference data. For molecular design, generative models are used to tackle the inverse design problem [109], where novel structures with desired properties are generated. Molecular property estimators are often used in conjunction with these models [110]. In this section, we discuss the importance of property prediction uncertainty in the context of generative modeling. In particular, we consider variational autoencoders (VAEs) [111], generative adversarial networks (GANs) [112], reinforcement learning (RL) [113], and genetic algorithms (GAs) [114].
VAEs possess the ability to take discrete structures such as molecules and convert them into high-dimensional latent representations by way of neural networks (referred to as encoders). Simultaneously, a mapping from the latent representation back to the discrete structure is trained using a separate neural network (referred to as a decoder). During training, VAEs learn to encode and decode correct molecules and simultaneously shape the latent space into a gaussian distribution for sampling. For molecular design, Gómez-Bombarelli et al. [115] introduced an additional neural network that predicts molecular properties from the latent space representation. During the training of this network, gradient optimization steps are taken for both the parameters of the encoder network and the property prediction network. This shapes the latent space based on the predicted property.
The formulation of VAEs as generative models brings about three sources of uncertainty for the associated property prediction model: (1) Decoding points from a latent space does not always lead to valid molecules. However, for every point in the latent space one can predict molecular properties via the prediction model. Hence, these regions can be viewed as uncertain spaces where predictions have to be disregarded due to the lack of a physical manifestation. (2) The encoding of structures to the latent space is stochastic. Kingma et al. [111] introduced a reparameterization step that ensures differentiability through the latent space for training the overall architecture and imposition of a gaussian distribution. As a result, identical molecules are typically encoded at a range of points in the latent space. While these encoded molecules should have the exact same predicted property, the range of positions tends to give rise to different predictions. Hence, there are multiple property predictions associated with any given molecule representing the corresponding uncertainty. (3) Discrete objects within VAEs, in particular molecular string representations and adjacency matrices, are usually decoded from the latent space as a time-series. These discrete objects are generated one at a time, by sampling over the distribution of all user-defined actions, i.e. units such as the allowed characters in a string. Typically, the unit with the largest probability is selected. However, multiple units can have identical probabilities. Consequently, decoding the same point from the latent space multiple times can result in distinct molecules. Additionally, as the same latent space point represents all these different molecules simultaneously, they will have the exact same predicted property.
GANs [112,116,117] are generative models with joint training of two competing networks, a generator and a discriminator. Like VAEs, they rely on a gaussian latent space to draw samples from and decode them using a neural network (referred to as generator network). The generated samples are compared to references from a dataset by another neural network (referred to as discriminator). The generator and discriminator compete with each other, the generator tries to produce samples that appear more realistic to the discriminator while the discriminator tries to improve its reliability to distinguish generated from reference samples. Similar to VAEs, in the traditional formulation of GANs, property predictors can be associated with the latent space. The sources of uncertainty are essentially identical to VAEs.
RL is a field where agents learn to navigate environments to maximize rewards. In molecular design, RL agents build chemical compounds [118][119][120], and rewards are provided for achieving desired molecular properties. Actions constitute making changes to chemical structures, while rewards are calculated on-the-fly based on the resulting structures. In contrast, GAs [121][122][123] are generative models that generate structures without gradient-based optimization. The algorithm is initiated with a set of molecules, i.e. the initial population, that are ranked by calculating their fitness. For molecular design, the fitness is typically assessed based on a set of target properties. The fittest molecules are carried over to the next generation, the unfit molecules are replaced with mutations and crossovers of the best. For both RL and GAs, to offset the often prohibitive computational cost of obtaining accurate property predictions, cheap approximations such as ML models are utilized to estimate rewards and fitness, respectively. By using cheap property predictors, in high-dimensional spaces, the generation of structures can be skewed to uncertain regions of the chemical space. For instance, certain structures can be adversarial examples or outliers for property predictors, where the model can be overly confident or uncertain about its prediction. Hence, the resulting molecule modifications are uncertain due to the uncertainty in the rewards. Notably, Thiede et al. [124] demonstrated that the uncertainty from prediction networks can be used directly in the reward function to guide the generation of structures outside of the confidence region of the prediction network for curiositydriven molecular design. Interestingly, Shen et al. [125] propose PASITHEA, where a property prediction neural network is directly used for generating molecules with deep dreaming. An interesting direction would be the incorporation of uncertainty within property predictors for structure generation via the deep dreaming algorithm.
Generative models coupled with property predictors are yet to achieve wide adoption due to their complex setup, their high computational costs and the need to adapt them to specific problems. When used for properties that are inherently complicated and hard to predict, the role of uncertainty and prediction confidence will be very important. In addition, in models with latent spaces such as VAEs and GANs, the choice of molecular representation influences both the validity of the latent space and thereby the prediction confidence. For instance, the use of VAEs and GANs was demonstrated in conjunction with the SMILES [58] representation, where the probability to produce invalid strings is high. Consequently, both VAEs and GANs trained on SMILES can produce invalid molecules. Importantly, this problem has recently been solved by the SELFIES representation [126], which guarantees molecular validity by design. While SELF-IES increases model confidence by generating only valid molecules, it does not remove uncertainty from inherent model stochasticities such as probabilistic encoding and decoding. Furthermore, confidence in molecular property predictors can be increased by extending datasets used for training and validation of property predictors systematically. Recently, STONED [127], a simple generative model producing molecules at high speed via systematic navigation within SELFIES, has been introduced. Accordingly, based on an initial dataset, it can be used to diversify the molecule space systematically leading to better property prediction models.

VII. CONCLUSION
In this review, we inspected common sources of uncertainties encountered in data-driven drug design workflows, some of which are well-studied in the literature, and some of which are typically neglected. Hidden biases in data can lead to overly confident property predictions that fail to manifest in the real world. Data-driven prediction models that explicitly incorporate uncertainty convey insights into the confidence of the model itself and provide a handle to assess the prediction quality. Accounting for uncertainty in the input features is underappreciated in ML research but is inherent to complex systems such as the chemistry of the human body. In contrast, error estimation in the prediction of binding free energies is one of the central methods in the field to assess the quality of simulation workflows that guides the development of efficient sampling techniques. Uncertainties in generative models are another underappreciated topic that needs to be addressed to improve the success of property-driven molecular design approaches. Overall, it is evident that assigning confidence in molecular prop-erty predictions is intimately tied to understanding the associated sources of error, and we believe that future progress in this field will require the development and application of even more methods accounting for uncertainties.

VIII. EXPERT OPINION
Assigning confidence to predictions is directly connected to uncertainty identification and quantification. When training ML models for predicting complex molecular properties, users need to be aware of potential biases of the training data when assessing generalizability. The recent use of DUD-E to train data-driven models serves as a cautionary tale for models learning to exploit biases in the data rather than learning the underlying physics of intermolecular interactions [24]. Most importantly, the smaller the training dataset size the larger the danger of hidden bias and the higher the importance of accounting for them. In that regard, transfer learning, especially in extremely sparse data regimes, requires more research to establish it as the method of choice for molecular property prediction.
Moreover, not only do small datasets lead to significant biases but they also lead to higher epistemic uncertainty because the model does not have enough information to learn the true relationship between input features and output labels [29]. However, expanding datasets, especially for properties relevant in drug design, can be prohibitively expensive and very time-consuming, in particular in the case of experimental data. Furthermore, not every new datum added to an existing dataset provides the same amount of additional information. Accordingly, developing systematic methods to identify data to be added to the existing dataset to maximize the improvement in prediction accuracy and the reduction in epistemic uncertainty is the subject of ongoing research.
Additionally, the aleatoric uncertainty [29] that is inherent in the data used to train ML models is at least as important as the epistemic one. One prominent example from the field of drug discovery discussed in this review is the use of estimated and simulated binding affinities to assess the interaction strength between ligands and receptors. Notably, binding affinities estimated via molecular docking are only of limited accuracy, but their computation efficiency allows assembling large datasets. The main question is whether the estimated property, regardless of how efficient it can be computed, conveys any information about the real relationship between input features and model output. When noise dominates over valuable information, the largest dataset will not provide useful information for ML models to be learned. In summary, it is important to build datasets and construct models that explicitly and consistently account for both aleatoric and epistemic uncertainties. While uncertainty estimates are relatively straightforward to obtain in general, they are only rarely benchmarked or evaluated during model training and future work needs to establish robust procedures for that purpose.
Another important and often underappreciated source of uncertainty, especially in complex systems, are the input features. This is particularly true for chemistry where molecular properties are generally conformer-dependent, and most molecules adopt much more than just one conformation under ambient conditions. Data-driven models, in contrast, are largely trained on 2D molecular representations that neglect conformations entirely. Accordingly, incorporating 3D structures and appropriate features in ML models is an active field of research [63]. The main challenge is to offset the added cost of generating representative 3D conformations with higher prediction accuracy in the final models. Currently, 2D molecular representations are still the state of the art in the field but we believe 3D features will become the standard in the near future.
One topic of paramount importance in drug design is simulating binding affinity between ligands and receptors. For that purpose, docking is a popular technique that provides crude estimates of interaction geometries and strengths using extremely efficient scoring approaches. However, this efficiency comes at the price of accuracy and, hence, cannot match the quality of FEP simulations, which are several orders of magnitude more computationally demanding. Importantly, even among FEP approaches there is an inevitable trade-off between cost and accuracy as the extent of sampling conducted directly controls the quality of the results. Accordingly, for extremely accurate results [71], FEP methods need to capture multiple binding transitions and this can increase the simulation cost by orders of magnitude making it less practical. Finally, a large portion of the uncertainties from both docking approaches and MD simulations stem from the scoring functions and the force fields used, respectively. While systematic optimization procedures for force field parameters are being developed and implemented, the development of force fields that not only provide energy but also uncertainty estimates has received less attention and we believe that future work in that direction is necessary to deliver more accurate models.
Finally, the generative models for molecular design introduced recently rely on efficient molecular property predictors to generate molecules with both desirable structure and function [110]. However, an explicit account and treatment of uncertainties, unfortunately, is largely neglected in the field. While accounting for property prediction uncertainties seems relatively straightforward with current methodologies, dealing with uncertainties in the molecular structure itself, and modifications thereof remains unexplored territory in comparison. Accord-ingly, we believe that future research in that direction is paramount and has the potential to lead to significant advances in computer-driven molecular design in general and drug design in particular.