AMPL: A Data-Driven Modeling Pipeline for Drug Discovery

AMPL: A Data-Driven Modeling Pipeline for Drug Discovery

Amanda J. Minnich minnich2@llnl.gov    Kevin McLoughlin Lawrence Livermore National Laboratory    Margaret Tse    Jason Deng    Andrew Weber    Neha Murad GlaxoSmithKline    Benjamin D. Madej Frederick National Laboratory    Bharath Ramsundar Computable    Tom Rush    Stacie Calad-Thomson GlaxoSmithKline    Jim Brase    Jonathan E. Allen Lawrence Livermore National Laboratory
Abstract

One of the key requirements for incorporating machine learning into the drug discovery process is complete reproducibility and traceability of the model building and evaluation process. With this in mind, we have developed an end-to-end modular and extensible software pipeline for building and sharing machine learning models that predict key pharma-relevant parameters. The ATOM Modeling PipeLine, or AMPL, extends the functionality of the open source library DeepChem and supports an array of machine learning and molecular featurization tools. We have benchmarked AMPL on a large collection of pharmaceutical datasets covering a wide range of parameters. Our key findings include:

  • Physicochemical descriptors and deep learning-based graph representations are significantly better than traditional fingerprints to characterize molecular features

  • Dataset size is directly correlated to performance of prediction: single-task deep learning models only outperform shallow learners if there is enough data. Likewise, data set size has a direct impact of model predictivity independently of comprehensive hyperparameter model tuning. Our findings point to the need for public dataset integration or multi- task/ transfer learning approaches.

  • DeepChem uncertainty quantification (UQ) analysis may help identify model error; however, efficacy of UQ to filter predictions varies considerably between datasets and model types.

This software is open source and available for download at http://github.com/ATOMconsortium/AMPL.

Machine learning, deep learning, pipeline, data-driven modeling, QSAR, drug discovery, pharmacokinetics, liability assays
\abbreviations

IR,NMR,UV

1 Introduction

Discovery of new compounds to treat human disease is a multifaceted process involving the selection of chemicals with favorable pharmacological properties: a high potency to the desired target, elimination or minimization of safety liabilities, and a favorable pharmacokinetic (PK) profile. To address this challenge, the drug discoverer has a wealth of choices, with total “drug-like” chemical matter estimated between - unique molecules. However, evaluating the desirability of these molecules with respect to potency, pharmacokinetics, and safety liabilities is a time-consuming and expensive process. Many of these molecules require de novo synthesis, which is a rate-limiting step. Furthermore, evaluation of pharmacological properties both in vitro and especially in vivo is prohibitively expensive given the universe of possible choices. To aid in this design challenge, the field of computer-aided drug design has evolved to rapidly predict the properties of pharmacological matter in silico, allowing for rational selection of a feasible set compounds for synthesis and evaluation. These techniques generally fall into two categories: (1) structure-based drug design, which relies on knowledge of the target structure (i.e. docking, molecular dynamics, free energy perturbation) and (2) ligand-based drug design, which uses known properties of molecules to develop models of quantitative structure-activity relationships (QSAR). Ligand-based drug design generally relies on machine learning-based techniques to identify the link between structure and the property of interest. Recently, a proliferation of advanced machine learning techniques have shown great promise in increasing the predictability of QSAR models. A deep learning model first won the Merck Kaggle multi-activity challenge in 2014dahl_multi-task_2014, and since then these models have continued to show increased predictive accuracy over QSAR models based on classical machine learning techniques in many studiesgilmer_neural_2017. A recent example of success with deep learning is the paper by Feinberg et al. that compared the PotentialNet deep learning method with existing shallow learners on a wide array of pharmaceutically-relevant datasets feinberg_step_2019. These results showed dramatic improvements for deep learning based on temporal splits using data collected from a pharmaceutical company. Another evaluation showed that a directed message-passing neural network model can provide robust performance over a range of experimental datasets characterizing molecular propertiesyang_analyzing_2019. The authors provide an open-source deep learning software to go with this paper that has been tested on a wide range of parameters. However, this software does not include any type of modular pipeline that would allow for the incorporation of different models and chemical representations. Overall there has been a lack of publicly available suites of software tools that support a transparent and reproducible generation of a diverse array of deep and classical machine learning models, especially ones that can scale to model the large set of pharmaceutically-relevant parameters. A major advance towards this goal was made with the introduction of DeepChem deepchem, which supports the building of a variety of machine learning models for small molecule property prediction. DeepChem contains a variety of very helpful modules and tools, but has limitations in its ability to robustly train models from a wide selection of hyperparameters, and published performance evaluation is limited to a small number of public datasets with less diverse pharmaceutical relevance wu_moleculenet:_2018. In this paper we introduce a new small molecule property prediction pipeline, AMPL. This software was developed through the Accelerating Therapeutics Opportunities in Medicine (ATOM) Consortium as the ATOM Modeling PipeLine. The key contributions of this work are to automate deep learning training, particularly in hyperparameter search; to enable extensive performance benchmarking; and to apply AMPL to a large collection of pharmaceutically-relevant property-prediction datasets. Most notably, AMPL is available as open source to benefit the drug discovery community. The closest existing pipeline tools are BIOVIA Pipeline Pilot pipeline_pilot and KNIME knime. Pipeline Pilot is a license-based graphical tool for machine learning pipelining. It has capabilities for data cleaning, splitting, training, and model deployment, but are all mainly GUI-based, limiting the customizability of the software. Furthermore, it is only available for a licensing fee, so it does not target the open source community. In terms of free and open source software suites for data analytics, the main alternative is KNIME. This software provides an environment for creating general data flows to process data, use predictive models, and analyze complex datasets. An ecosystem of open source and commercial KNIME node extensions has developed which enable workflows for library analyses, virtual screening, model fitting and prediction. In contrast, AMPL is tightly focused on integrating modern machine learning methods with best practices for chemical activity and property prediction. Important issues with machine learning models, such as dataset characterization, model validation, and uncertainty quantification are addressed by AMPL in automated and reproducible ways. The code suite also provides high performance computing modules for model fitting, hyperparameter optimization, and predictions. AMPL currently targets job submission-based clusters to scale training runs; however, AMPL could be adapted to operate on other scalable platforms such as Spark in the future. Furthermore, AMPL is implemented as a modular and reusable Python library to allow for easy integration with other data science software platforms. An extensive set of experiments were conducted with AMPL, and key observations include:

  • Physicochemical descriptors and deep learning-based graph representations are significantly better than traditional fingerprints to characterize molecular features

  • Dataset size is directly correlated to performance of prediction: single-task deep learming models only outperform shallow learners if there is enough data. Likewise, data set size has a direct impact of model predictivity independently of comprehensive hyperparameter model tuning. Our findings point to the need for public dataset integration or multi- task/ transfer learning approaches.

  • DeepChem uncertainty quantification (UQ) analysis may help identify model error; however, efficacy of UQ to filter predictions varies considerably between datasets and model types.

The aim of this paper is to present the rigorous and transparent open source software pipeline AMPL to build global and local ‘baseline’ models for a wide array of molecular properties that are needed for in silico drug discovery. This new software will support reproducible training and testing protocols that enable the broader modeling community to evaluate and improve modeling approaches over time.

2 Methods

Figure 1 shows the overall architecture of AMPL. This end-to-end pipeline supports all functions needed to generate, evaluate, and save machine learning models: data ingestion and curation, featurization of chemical structures into feature vectors, training and tuning of models, storage of serialized models and metadata, and visualization and analysis of results. It also contains modules for parallelized hyperparameter search on high-performance computing (HPC) clusters.

Figure 1: Overview of AMPL

2.1 Data curation

AMPL includes several modules to curate data into machine learning-ready datasets. Functions are provided to represent small molecules with canonicalized SMILES strings using RDKit landrum_fingerprints_nodate and the MolVS package matt_swain_2017_260237, by default stripping salts and preserving isomeric forms. Data curation procedures are provided with AMPL as Jupyter notebooks noauthor_jupyterlab_nodate, which can be used as examples for curating new datasets. Procedures allow for averaging response values for compounds with replicate measurements and filtering compounds with high variability in their measured response values. AMPL also provides functions to assess the structural diversity of the dataset, using either Tanimoto distances between fingerprints, or Euclidean distances between descriptor feature vectors. Data ingestion and curation-related parameters include:

  • Unique human readable name for training file

  • Data privilege access group

  • Parameter for overriding the output files/dataset object names

  • ID for the metadata + dataset

  • Boolean flag for using an input file from the file system

  • Name of column containing compound IDs

  • Name of column containing SMILES strings

  • List of prediction task names

  • Number of classes for classification

  • User specified list of names of each class

  • Boolean switch for using transformation on regression output. Default is True

  • Response column normalization type

  • Minimum number of dataset compounds considered adequate for model training. A warning message will be issued if the dataset size is less than this.

2.2 Featurization

AMPL provides an extensible featurization module which can generate a variety of molecular feature types, given SMILES strings as input. They include:

  • Extended connectivity fingerprints (ECFP) with arbitrary radius and bit vector length rogers_extended-connectivity_2010

  • Graph convolution latent vectors, as implemented in DeepChem duvenaud_convolutional_2015

  • Chemical descriptors generated by the Mordred open source package moriwaki_mordred:_2018

  • Descriptors generated by the commercial software package Molecular Operating Environment (MOE) noauthor_chemical_nodate

  • User-defined custom feature classes

Because some types of features are expensive to compute, AMPL supports two kinds of interaction with external featurizers: a dynamic mode in which features are computed on-the-fly and a persistent mode whereby features are read from precomputed tables and matched by compound ID or SMILES string. In the persistent mode, when SMILES strings are available as inputs, the featurization module matches them against the precomputed features where possible, and computes features dynamically for the remainder. Because precomputed feature tables may span hundreds or thousands of feature columns for millions of compounds, the module uses the feather format noauthor_feather_nodate to speed up access. Featurized datasets for feature types that support persistent mode (currently, all except ECFP fingerprints and graph convolution format) are saved in the filesystem or remote datastore, so that multiple models can be trained on the same dataset. This also facilitates reproducibility of model results. Chemical descriptor sets such as those generated by MOE often contain descriptors that are exact duplicates or simple functions of other descriptors. In addition, large blocks of descriptors may be strongly correlated with one another, often because they scale with the size of the molecule. The featurization module deals with this redundancy by providing an option to remove duplicate descriptors and to scale a subset of descriptors by the number of atoms in the molecule (while preserving the atom count as a distinct feature). Factoring out the size dependency often leads to better predictivity of models. The featurization module can be easily extended to handle descriptors generated by other software packages, latent vectors generated by autoencoders, and other types of chemical fingerprints. In most cases, this can be accomplished by writing a small function to invoke the external feature generation software, and by adding an entry to a table of descriptor types, listing the generated feature columns to be used. In more complicated cases, one may need to write a custom subclass of one of the base featurization classes. Featurization-relevant input parameters include:

  • Type of molecule featurizer

  • Feature matrix normalization type

  • Boolean flag for loading in previously featurized data files

  • Type of transformation for the features

  • Radius used for ECFP generation

  • Size of ECFP bit vectors

  • Type of autoencoder, e.g. molvae, jt

  • Trained model HDF5 file path, only needed for MolVAE featurizer

  • Type of descriptors, e.g. MOE, Mordred

  • Max number of CPUs to use for Mordred descriptor computations. None means use all available

  • Base of key for descriptor table file

2.3 Dataset partitioning

AMPL supports several options for partitioning datasets for model training and evaluation, By default, datasets are split into 3 parts: a training set, a validation set (for parameter selection), and a holdout test set (for evaluation). Alternatively, AMPL offers a k-fold cross-validation option, to assess the performance impact of sampling from the modeled dataset. Under k-fold cross-validation, the holdout test set is selected first, and the remainder is divided into k-fold sets for training and validation. AMPL offers a number of dataset splitting algorithms, which offer different approaches to the problem of building models that generalize from training data to novel chemical space. It supports several of the methods included in DeepChem, including random splits, Butina clustering, Bemis-Murcko scaffold splitting, and a simple algorithm based on fingerprint dissimilarity wu_moleculenet:_2018. In addition, we implemented temporal splitting and a modified version of the asymmetric validation embedding (AVE) debiasing algorithm wallach_most_2018. We compared random splitting with Bemis-Murcko scaffold splitting for our benchmarking experiments. Input parameters related to data splitting include:

  • Type of splitter to use: index, random, scaffold, Butina, ave_min, temporal, fingerprint, or stratified

  • Boolean flag for loading in previously-split train, validation, and test CSV files

  • UUID for CSV file containing train, validation, and test split information

  • Choice of splitting type between k-fold cross validation and a normal train/valid/test split

  • Number of k-folds to use in k-fold cross validation

  • Type of splitter to use for train/validation split if temporal split used for test set (random, scaffold, or ave_min)

  • Cutoff Tanimoto similarity for clustering in Butina splitter

  • Cutoff date for test set compounds in temporal splitter

  • Column in dataset containing dates for temporal splitter

  • Fraction of data to put in validation set for train/valid/test split strategy

  • Fraction of data to put in held-out test set for train/valid/test split strategy

2.4 Model training and tuning

AMPL includes a train/tune/predict framework to create high-quality models. This framework supports a variety of model types from two main libraries: scikit-learn sklearn and DeepChem deepchem. Currently, specific input parameters are supported for:

  • Random forest models from scikit-learn

  • XGBoost models xgboost

  • Fully connected neural network models

  • Graph convolution neural network models graphconv

As with the featurization module, AMPL supports integration of custom model sub-classes. Parameters for additional models can be easily added to the parameter parser module. Model-relevant input parameters include:

  • Type of model to fit (neural network, random forest, or xgboost)

  • Prediction type (regression or classification)

  • Singletask or multitask model

  • Number of decision trees in the forest for random forest models

  • Max number of features to split random forest nodes

  • Number of estimators to use in random forest models

  • Batch size for neural network model

  • Optimizer type for neural network model

  • Optimizer specific for graph convolutional models, defaults to “adam”

  • Model batch size for neural network model

  • List of hidden layer sizes for neural network model

  • List of dropout rates per layer neural network model

  • List of standard deviations per layer for initializing weights for neural network model

  • The type of penalty to use for weight decay, either “l1” or “l2”

  • The magnitude of the weight decay penalty to use

  • List of initial bias parameters per layer for neural network model

  • Learning rate for dense neural network models

  • Epoch for evaluating baseline neural network model performance, if desired

  • Maximum number of training epochs for neural network model

  • Type of score function used to choose best epoch and/or hyperparameters

  • Boolean flag for computing uncertainty estimates for regression model predictions

2.4.1 Epoch selection for neural network models

Early stopping is an essential strategy to avoid overfitting of neural networks, thus the number of training epochs is one of the key hyperparameters that must be optimized. To implement early stopping, AMPL trains neural network models for a user-specified maximum number of epochs, evaluating the model on the validation set after each epoch, and identifies the epoch at which a specified performance metric is maximized. By default this metric is the coefficient of determination for regression models, and the area under the receiver operating characteristic curve (ROC AUC) for classification models.

2.4.2 Model persistence

Serialized models are saved after training and prediction generation are complete, along with detailed metadata to describe the model. This supports traceability and reproducibility, as well as model sharing. AMPL supports saving models and results either using the file system or optionally through a collection of database services. The metadata can be stored in a mongoDB database mongo or as JSON files. AMPL has functions for saving models and loading pre-built models for prediction generation.

2.5 Model performance metrics

AMPL calculates a variety of performance metrics for predictions on the training, validation and test sets. Metrics may be saved in a mongoDB database or in JSON files. For regression models, we calculate:

  • Coefficient of determination (). This is calculated using sklearn’s metrics function. Note that this score can be negative if the model is arbitrarily worse than random.

    (1)
  • Mean Absolute Error (MAE). An advantage of MAE is that it has a clear interpretation, the average absolute difference between the measured value and predicted value . This works well for cellular activity assay datasets, which use log normalized dose concentration value with similar concentration ranges across different assays. PK parameters are measured on different scales for some assays, which prevents comparison between assays with this metric.

    (2)
  • Mean Square Error (MSE). This is a risk metric corresponding to the expected value of the squared error (or loss).

    (3)

For classification models, we calculate:

  • Area Under the Receiver Operating Characteristics Curve (ROC AUC). The ROC curve plots the True Positive Rate versus the False Positive Rate as a binary classifier’s discrimination threshold is varied. The ROC AUC score is calculated by finding the area under the ROC Curve. This value can range from 0 – 1, where 1 is the best score.

  • Precision (Positive Predictive Value)

    (4)

    where TP = number of true positives and FP = number of false positives

  • Recall (True positive rate/ sensitivity)

    (5)

    where TP = number of true positives and TN = number of true negatives

  • Area under the precision-recall curve (PRC-AUC). The precision-recall curve plots precision versus recall as a binary classifier’s discrimination threshold is varied. It is a good measure of success of prediction when classes are very imbalanced. High scores show that the classifier is returning accurate results (high precision), as well as returning a majority of all positive results (high recall).

  • Negative Predictive Value (NPV)

    (6)

    where TN = number of true negatives and FN = number of false negatives

  • Cross entropy (log loss)

    (7)
  • Accuracy

    (8)

    where terms are defined as above.

2.6 Uncertainty quantification

Uncertainty quantification (UQ) attempts to measure confidence in a model’s prediction accuracy by characterizing variance in model predictions. Some common objectives for UQ are to use it to guide active learning or to weight model ensembles. AMPL generates UQ values for both random forest and neural network models.

2.6.1 Uncertainty quantification for random forest

Generating a value quantifying uncertainty is straightforward for random forest and is taken to be the standard deviation of predictions from individual trees. This quantifies how variable these predictions are, and thus how uncertain the model is in its prediction for a given sample.

2.6.2 Uncertainty quantification for neural networks

Our neural network-based UQ uses the Kendall and Gal methodkendall_what_2017 as implemented in DeepChem. This method combines aleatoric and epistemic uncertainty values. Aleatoric uncertainty cannot be reduced by adding more data but can be estimated. It is estimated by modifying the loss function of the model to predict both the response variable and the error of the model. Epistemic uncertainty arises because of limited data. It represents the uncertainty of the model. Normally this is calculated in a bootstrapped manner, as in the case of a random forest. However, since training neural networks is expensive, an alternate approach is to train one network to generate a set of predictions by applying a set of dropout masks during prediction. Prediction variability is then quantified to assess epistemic uncertainty.

2.7 Visualization and analysis

Plots generated by AMPL’s visualization and analysis module are shown in the Results section. Additional options include plots of predicted vs. actual values, learning curves, ROC curves, precision vs. recall curves, and 2-D projections of feature vectors using UMAP mcinnes_umap_2018. The module also includes functions for characterizing and visualizing chemical diversity. Chemical diversity analysis is crucial for analyzing domain of applicability, bias in dataset splitting, and novelty of de novo compounds. This module supports a wide range of input feature types, distance metrics, and clustering methods.

2.8 Hyperparameter optimization

A module is available to support distributed hyperparameter search for HPC clusters. This module currently supports linear grid, logistic grid, and random hyperparameter searches, as well as iteration over user-specified values. To run the hyperparameter search, the user specifies the desired range of configurations in a JSON file. The user can either specify a single dataset file or a CSV file with a list of datasets. The script generates all valid combinations of the specified hyperparameters, accounting for model and featurization type, and submits jobs for each combination to the HPC job scheduler. The module includes an option to generate a pre-featurized and pre-split dataset before launching the model training runs, so that all runs operate on the same dataset split. The user can specify a list of layer sizes and dropouts to combine, along with the maximum final layer size and a list of the numbers of possible layers for a given model, and the module combines these different options based on the input constraints to generate a variety of model architectures. The search module can check the model tracker database to avoid retraining models that are already available. It also provides users the option to exclude hyperparameter combinations that lead to overparameterized models, by checking the number of weight and bias parameters for a proposed neural network architecture against the size of the training dataset. Finally, the search module throttles job submissions to prevent the user from monopolizing the HPC cluster. Input parameters for hyperparameter search include:

  • Boolean flag indicating whether we are running the hyperparameter search script

  • UUID of hyperparam search run model was generated in

  • Comma-separated list of number of layers for permutation of NN layers

  • Comma-separated list of dropout rates for permutation of neural network layers

  • The maximum number of nodes in the last layer

  • Comma-separated list of number of nodes per layer for permutation of neural network layers

  • Maximum number of jobs to be in the queue at one time for an HPC cluster

  • Scaling factor for constraining network size based on number of parameters in the network

  • Boolean flag directing whether to check model tracker to see if a model with that particular param combination has already been built

  • Path where pipeline file you want to run hyperparam search from is located

  • Type of hyperparameter search to do. Options are grid, random, geometric, and user_specified

  • CSV file containing list of datasets of interest

2.9 Running AMPL

There are three ways to run AMPL:

  • Using a config file: Create a JSON file with desired model parameters and run full pipeline via command line

  • Using command line arguments: Specify model parameters via standard command line arguments

  • Interactively in a Jupyter notebook using an argparse.Namespace object or a dictionary

3 Results

Benchmark experiments were run to evaluate and validate components of the pipeline.

3.1 Data

Experimental datasets were made available by ATOM Consortium member GlaxoSmithKline from a variety of bioactivity and pharmacokinetics experiments. Selected datasets were used for training and evaluating models. These datasets are summarized in Table 1. Pharmacokinetic (PK) and safety datasets were curated separately, as they contain different types of experimental data and thus require different processing. The raw datasets were cleaned to remove rows with outlying, missing, and duplicate measurements, and processed to yield machine learning datasets with a single aggregate value per unique compound. These procedures informed the design of curation functions included in the pipeline. Curation of the PK datasets required the conversion of values to standard units, the removal of compounds with stability or recovery issues, and the exclusion of data that was generated using significantly different assay protocols. Subsequently, replicate experimental measurements were identified by matching duplicate canonical SMILES strings and averaged to produce a single value per compound. For the safety datasets, censored measurements were an additional concern. Since bioactivity assays are typically performed over a limited range of compound concentrations, IC50 or EC50 values may be reported as being above or below a maximum or minimum concentration, so that the measurements are censored. When all measurements for a compound are censored in the same direction, the user is given the option to either exclude the compound from the dataset, or include it with a relational operator indicating the direction together with the censoring threshold. In the case where some replicate measurements are censored and some are not, AMPL computes a maximum likelihood estimate for the mean activity, assuming a Gaussian distribution of measurements around the true mean. The distribution of response values is reported along with the mean and standard deviation.

Dataset Units Species Dataset Size Minimum Maximum Mean Median
Blood to Plasma Ratio Human 101 0.47 10.5 0.85 0.77
Blood to Plasma Ratio Dog 71 0.37 6.85 0.85 0.88
Plasma Protein Binding HSA fraction unbound Human 123734 0.0001 1 0.05 0.044
Plasma Protein Binding HSA fraction unbound Rat 2086 0.0001 1 0.036 0.033
Plasma Clearance (In Vivo) mL/min/kg Dog 1181 0.1 2946 12.6 15.2
Plasma Clearance (In Vivo) mL/min/kg Rat 10431 0.001 8763 30.2 38.2
Vd,ss L/kg Dog 1054 0.07 569 1.9 1.9
Vd,ss L/kg Rat 9681 0.01 2080 2.3 2.4
Hepatocyte Clearance mL/min/g liver tissue Human 1695 0.01 97 1.6 1.5
Hepatocyte Clearance mL/min/g liver tissue Dog 630 0.1 504 2 1.8
Hepatocyte Clearance mL/min/g liver tissue Rat 2098 0.02 878 2.9 2.9
Microsomal Clearance mL/min/g liver tissue Human 29162 0 156 2.8 2.4
Microsomal Clearance mL/min/g liver tissue Dog 2080 0.03 150 2.5 1.8
Microsomal Clearance mL/min/g liver tissue Rat 30563 0.01 198 3.9 3.7
LogD 27345 0.01 53703 258 407
Table 1: Pharmacokinetics datasets used to benchmark AMPL
Assay Target Primary Liability Experimental System Detection
BSEP pIC50 Bile Salt Export Pump Hepatic membrane vesicles
ADRA1B pIC50 Adrenergic 1B pIC50 CNS Intracell Ca
ADRA2C pIC50 2C Adrenoceptor CNS CHO K1
ADRB2 pEC50 2 Adrenoceptor CNS FRET
CHRM1 pEC50 Cholinergic Receptor Muscarinic 1 CNS CHO Intracellular Ca Fluorescence
CHRM1 pIC50 Cholinergic Receptor Muscarinic 1 CNS CHO Intracellular Ca Fluorescence
CHRM2 pEC50 Cholinergic Receptor Muscarinic 2 CNS CHO Intracellular Ca Fluorescence
DRD2 pEC50 Dopamine D2 CNS HEK293F Low Na GTPgS SPA
GRIN1 pIC50 GRIN1 GRIN2B NR2B NR1A 2B Subunit pIC50 CNS
HRH1 pIC50 Histamine Receptor H1 CNS Luminescence
HTR1B pIC50 5-hydroxytryptamine receptor 1B CNS 10ul LEADseeker GTPgS
HTR2A pEC50 5-hydroxytryptamine Receptor 2A CNS HEK Luminescence
HTR2A pIC50 5-hydroxytryptamine Receptor 2A CNS HEK Luminescence
HTR2C pEC50 5-hydroxytryptamine Receptor 2C CNS CHO Luminescence
HTR2C) pIC50 5-hydroxytryptamine Receptor 2C CNS CHO Luminescence
HTR3A pIC50 5-hydroxytryptamine Receptor 3A CNS FLIPR
KCNA5 (Kv1.5) pIC50 KCNA5 (Kv1.5) Cardiovascular CHO Electrophys
KCNE1 KCNQ1 (Kv7.1) pIC50 KCNE1 KCNQ1 (Kv7.1) Cardiovascular MinK Human Blocker CHO Electrophys
MAOA pIC50 Monoamine Oxidase A CNS FLINT
PDE3A pIC50 Phosphodiesterase 3A Cardiac SPA(cAMP Inhibition)
PDE4B pIC50 Phosphodiesterase 4B CNS SPA
Phospholipidosid pEC50 Phospholipidosis Induction Cellular Tox HEPG2 FLINT
PI3K pIC50 Phosphoinositide 3-kinase (pI3K) Cellular Tox TR FRET
COX 2 pIC50 Cyclooxygenase 2 Cardiovascular FLINT SAR
SCN5A (NaV1.5) pIC50 SCN5A (NaV1.5) Cardiovascular
SCL6A2 pIC50 Noradrenaline Transporter NET CNS BacMam Bind SPA
SLC6A4 pIC50 Seratonin Transporter (SERT) CNS BacMam binding SPA
OATP1B1 pIC50 Organic Anion Transport Polypeptide (SLCO1B1) Hepatic HEK Image
Table 2: Safety datasets used to benchmark AMPL

3.2 Experimental design for regression pharmacokinetic models

To evaluate AMPL’s performance, we built a total of 11,552 models on 15 pharmacokinetic datasets and 26 bioactivity datasets. These models include 9,422 regression models and 2,130 classification models. We evaluated a variety of deep learning model types and architectures and compared them to baseline random forest models. We explored the performance of four types of features: ECFP fingerprints, MOE descriptors, Mordred descriptors, and graph convolution-based latent vectors. For the neural network models, we searched over many combinations of learning rates, numbers of layers, and nodes per layer. For each combination of neural network hyperparameters, we trained for up to 500 epochs and used a validation set performance metric ( for regression, for classification) to choose an early stopping epoch for the final model. For random forest models, the only hyperparameter varied was the maximum tree depth, as previous experiments showed that other model hyperparameters had a minimal effect for our datasets. The complete set of hyperparameters varied was as follows:

  • Splitter Types: scaffold and random

  • Fraction for train set: 0.7

  • Fraction for validation set: 0.1

  • Fraction for holdout set: 0.2

  • Feature types: ECFP, MOE, mordred, and graph convolution

  • Model types: neural network and random forest

  • Neural network learning rates: 0.0001, 0.00032, 0.001, 0.0032, 0.01

  • Maximum number of epochs: 500

  • Number of layers: 1,2

  • Layer size options: 1024,256,128, 64, 32,16,8,4,1

  • Maximum final layer size: 16

  • Dropout rate: 0.1

3.3 Analysis of modeling performance

To identify which featurization type generated the most predictive models for each model type, models with the best validation set score were selected for each model/splitter/dataset combination. The number of “best” models for which each feature type yielded the highest test set score is plotted in Figure 2.

Figure 2: Number of times each featurization type produces the best model for the 15 PK datasets

Figure 2 shows that the chemical descriptors generated by the commercial MOE software outperformed those produced by the open source Mordred package in most cases. DeepChem’s graph convolution networks outperform all other feature types for neural network models. The model/featurization combination with the most accurate predictions on the holdout set is shown in Figure 3. MOE featurization with random forest models most frequently outperformed other featurization/model type combinations for both types of splitters.

Figure 3: Number of times each featurization type/ model type combination produces the best model for the 15 PK datasets

Figure 4 confirms that random forest models tend to outperform neural network models for the evaluated datasets.

Figure 4: Number of times each model type produces the best model for the 15 PK datasets

3.4 Investigation into neural network performance

Neural networks are known to perform more poorly on smaller datasets, so we wanted to examine the relationship between the size of a dataset and the test set values for the best random forest and neural network models for that dataset. Figure 5 shows the test set values for the best neural network and random forest models for each dataset, where best is defined as the model with the highest validation set value. The figure shows that as the dataset size increases, the score for the test set increases as well.

Figure 5: Plot of best test set values versus the dataset size for neural network and random forest models

The pattern is true for the overall best model, regardless of type, for both regression and classification, as shown in Figure 6. These results indicate that we will need to augment our datasets to further improve model performance. We plan to explore multiple avenues to address this requirement: conducting additional experiments, running simulations, sourcing public data, building multi-task models, and experimenting with transfer learning approaches.

Figure 6: Per-dataset model accuracy versus dataset size

We also examined the architectures that yielded the best model for each feature type for the neural network models. Our hypothesis was that larger datasets would perform better with larger networks. Figure 7 shows number of parameters in the hidden layers of the model versus the size of the dataset. The color indicates the dataset and the shape indicates the featurizer type. The number of parameters for the 2-layer networks was calculated by multiplying the first and second layers together. We can see a clear lower bound in the number of parameters for the best network for all featurizer types as the dataset size increases.

Figure 7: Number of hidden layer parameters versus number of samples for the best model for each dataset/featurizer combination

3.5 Summary of model performance

Figure 8 and Figure 9 show the full set of test set values for the best model for each molecular featurization representation and model type for random and scaffold splits respectively (picked as before by the best validation set value). Random sampling inflates the values of the holdout set, which is as expected since there is greater structural overlap between the set of compounds in the training and holdout set.

Figure 8: Performance accuracy for regression for random split

For scaffold split-generated holdout sets, there is a very clear pattern between dataset size and value, although the complexity of the predicted property and quality of the dataset also obviously has an effect.

Figure 9: Performance accuracy for regression for scaffold split

3.6 Model tuning results

To evaluate whether hyperparameter search improves model performance, the test set for a baseline model was compared with the test set from the best-performing model selected by looking at the validation set value. Small datasets and ECFP-based models, which showed poor neural network performance overall, showed little to no improvement, while better-performing datasets and featurizers showed greater improvement with hyperparameter search. This suggests that data augmentation will be necessary to improve prediction performance on the smaller problematic datasets, and that ECFP is a poor featurizer no matter the hyperparameters.

Figure 10: Histogram of improvement in values for the test set for the four featurizers for neural network models

3.7 Classification experiments

A set of classification model experiments were also conducted for a panel of 28 bioactivity datasets, without any hyperparameter tuning. In total 2,130 neural network and random forest models were generated. A dose concentration threshold was used to label active and inactive compounds on a per-dataset basis using thresholds provided by domain experts at GlaxoSmithKline. The classes were extremely unbalanced, which partially explains the high ROC-AUC scores.

Figure 11: Performance accuracy for classification

3.8 Uncertainty quantification

To explore the utility of the uncertainty quantification values produced by neural network and random forest models, a case study is presented for three representative PK parameter datasets: rat plasma clearance (in vivo), human microsomal clearance, and human plasma protein binding HSA. These datasets were selected to represent small, medium, and large sized datasets with low, medium, and high values.

3.8.1 Precision-recall plot analysis

Precision-recall curves measure the fraction of low error predictions made at varying UQ thresholds. Precision is defined as the fraction of predictions with UQ values less than the UQ threshold, with error less than some predefined threshold. For this analysis we use mean logged error and define “low-error” as samples with logged error below the mean (log served to normalize the distribution). Recall reports the fraction of low-error samples which pass the UQ filter threshold. Overall, we would like to use the UQ value as a threshold to identify low error samples at a higher rate than in the overall test set. Table 3 shows the percentage of low error samples in the test set as a whole for each dataset/model/featurizer combination.

Dataset Model and featurizer type Percent of total low error samples
Rat Plasma Clearance (In Vivo) Neural network + ECFP 41.4%
Rat Plasma Clearance (In Vivo) Neural network + GraphConv 41.8%
Rat Plasma Clearance (In Vivo) Neural network + MOE 42.9%
Rat Plasma Clearance (In Vivo) Neural network + Mordred 40.5%
Rat Plasma Clearance (In Vivo) Random forest + ECFP 42.5%
Rat Plasma Clearance (In Vivo) Random forest + MOE 41.7%
Rat Plasma Clearance (In Vivo) Random forest + Mordred 42.0%
Human Microsomal Clearance Neural network + ECFP 41.0%
Human Microsomal Clearance Neural network + GraphConv 41.0%
Human Microsomal Clearance Neural network + MOE 39.0%
Human Microsomal Clearance Neural network + Mordred 39.8%
Human Microsomal Clearance Random forest + ECFP 39.5%
Human Microsomal Clearance Random forest + MOE 38.5%
Human Microsomal Clearance Random forest + Mordred 39.6%
Human Plasma Protein Binding HSA Neural network + ECFP 43.4%
Human Plasma Protein Binding HSA Neural network + GraphConv 43.0%
Human Plasma Protein Binding HSA Neural network + MOE 43.1%
Human Plasma Protein Binding HSA Neural network + Mordred 43.5%
Human Plasma Protein Binding HSA Random forest + ECFP 42.0%
Human Plasma Protein Binding HSA Random forest + MOE 42.8%
Human Plasma Protein Binding HSA Random forest + Mordred 42.5%
Table 3: Percent of total low-error samples in the test set for the specified dataset, model/featurizer combinations
Figure 12: Precision-recall plot for rat plasma clearance (in vivo), varying UQ value
Figure 13: Precision-recall plot for human microsomal clearance, varying UQ value
Figure 14: Precision-recall plot for human plasma protein binding HSA, varying UQ value

In general, a low UQ threshold with accurate uncertainty would correspond to a precision of 1, which means confident predictions correspond to low-error predictions. To have the greatest utility, the curve should keep fairly high precision as the recall increases. UQ successfully filters out low confidence predictions in some cases but performance varies widely with the model/featurization type and the dataset. Figures 12, 13 and 14 show that precision drops quickly as recall increases and for some models precision is poor even when applying the lowest UQ threshold. Nevertheless, for each dataset there exists a UQ threshold for at least one model which could be used to increase the fraction of low error predictions over the baseline percentages shown in Table 3. For example, Figure 14 suggests that applying a UQ threshold could increase precision to 65% from around 42% with a recall of 10%. Later it is shown that for the human plasma protein binding HSA dataset, this could still yield a collection of compounds with a diverse range of response values.

3.8.2 Calibration curves

To further investigate how error changes as the uncertainty increases, we plotted calibration curves of mean error per uncertainty bucket, with the 95% confidence interval of error shown for each bucket as error bars. We would like uncertainty to serve as a proxy for error, so we would hope to see the mean error for the samples in a bucket to increase as the UQ threshold for that bucket increased. Results for neural network and random forest models built on MOE feature vectors and neural network graph convolution models are shown to demonstrate the variation in performance. For rat plasma clearance (in vivo), there is an overall upward trend for all three calibration curves, but it is not completely monotonically increasing for any of them. This is the smallest dataset of our case study, so increasing the bucket size may improve the choppiness of these curves, but overall UQ does not look like it would be a reliable proxy for error for this dataset.

Figure 15: Mean error per uncertainty bucket for rat plasma clearance (in vivo) neural network model with MOE features
Figure 16: Mean error per uncertainty bucket for rat plasma clearance (in vivo) random forest model with MOE features
Figure 17: Mean error per uncertainty bucket for rat plasma clearance (in vivo) neural network model with Graph Convolution features

Human microsomal clearance shows greater variation in the calibration curves. For MOE features with a neural network model, shows an inverse pattern where the error actually decreases as the uncertainty increases.

Figure 18: Mean error per uncertainty bucket for human microsomal clearance neural network model with MOE features

For MOE features with a random forest model, there seems to be no correlation, except for in the very highest bucket.

Figure 19: Mean error per uncertainty bucket for human microsomal clearance random forest model with MOE features

The graph convolution model, conversely, shows an upward trend, although it is not monotonically increasing.

Figure 20: Mean error per uncertainty bucket for human microsomal clearance neural network model with Graph Convolution features

These curves show that the featurizer and model type have a strong effect on the relationship between UQ and error. For human plasma protein binding HSA, which is the largest dataset with over 123,000 compounds, all calibration curves display the desired behavior: error increases as uncertainty increases and the 95 percent confidence intervals are small.

Figure 21: Mean error per uncertainty bucket for human plasma protein binding HSA neural network model with MOE features
Figure 22: Mean error per uncertainty bucket for human plasma protein binding HSA random forest model with MOE features
Figure 23: Mean error per uncertainty bucket for human plasma protein binding HSA neural network model with Graph Convolution features

3.8.3 Examining the relationship between UQ and predicted value

Since the UQ values quantify the variation in predictions, the relationship between UQ and the predicted values were checked for evidence of a correlation by examining plotted UQ versus predicted values. Rat plasma clearance (in vivo) shows a somewhat negative relationship, where the variation in predictions decreases as the magnitude of the predicted value increases. We found a similar though much less pronounced trend when examining error versus predicted value, so it looks like overall the model is predicting better for compounds with higher clearance values.

Figure 24: Uncertainty value versus Predicted for rat plasma clearance (in vivo) neural network model with MOE features
Figure 25: Uncertainty value versus Predicted for rat plasma clearance (in vivo) random forest model with MOE features
Figure 26: Uncertainty value versus Predicted for rat plasma clearance (in vivo) neural network model with Graph Convolution features

For human microsomal clearance, MOE feature vectors yield models where the UQ is strongly biased by the predicted value, especially for the neural network model, as seen in Figure 27. Error versus predicted value does not show this trend, so this is likely indicating that UQ contains no real information value for this model.

Figure 27: Uncertainty value versus Predicted for human microsomal clearance neural network model with MOE features

This trend exists for the MOE random forest model as well, although it levels off, suggesting slightly less biased UQ values.

Figure 28: Uncertainty value versus Predicted for human microsomal clearance random forest model with MOE features

The graph convolution model displays a more balanced relationship between UQ and predicted value, which mirrors what we saw in the previous two sub-sections, that this model’s UQ is more informative of error than the MOE models’ UQ.

Figure 29: Uncertainty value versus Predicted for human microsomal clearance neural network model with Graph Convolution features

Human plasma protein binding HSA, which showed the best calibration curves, also shows the least correlation between UQ and predicted value. UQ has a wide range of values for all predicted values.

Figure 30: Uncertainty value versus Predicted for human plasma protein binding HSA neural network model with MOE features
Figure 31: Uncertainty value versus Predicted for human plasma protein binding HSA random forest model with MOE features
Figure 32: Uncertainty value versus Predicted for human plasma protein binding HSA neural network model with Graph Convolution features

3.8.4 Correlation between UQ and error

While these plots provide useful methods for visualizing the behavior of uncertainty quantification, we wanted to identify a value that could summarize if we could trust a given model’s UQ results. Since we want the certainty of the model to be reflected in accurate predictions, we calculated the Spearman correlation coefficient between between binned prediction error and UQ. Results are shown in Figure 33. Correlations range from -0.088 to 0.33. While these correlations are fairly low, all p-values are , and all but one are , so there is significance to the weak correlations identified.

Figure 33: Spearman correlation coefficient between error and uncertainty values

4 Discussion

Key observations from the extensive series of model evaluations are summarized here:

  • Neural networks generally produced more accurate models only on the larger datasets.

  • The proprietary MOE descriptors outperformed the open-source Mordred descriptors for both random forest and neural networks. Among neural network representations, graph convolutions outperformed ECFP.

  • A range of neural network architectures performed best, depending on the dataset size. Small networks appear to be prominently featured in many datasets.

  • Model performance generally improved as dataset size increased, suggesting the need for public dataset integration or multi-task/ transfer learning approaches.

  • Hyperparameter tuning generally improved performance, in some cases dramatically.

  • Uncertainty quantification showed a weak correlation with error, and the efficacy of using UQ to filter predictions varied considerably between datasets and model types.

The differences in prediction accuracy show that the parameters needed for in silico drug discovery present a diverse set of data-driven modeling challenges. The extensive benchmarking suggests that there is no clear one best modeling approach for every predicted parameter. The differences in performance show the importance of having a rigorous model building pipeline that can be readily adapted and re-applied to build parameter specific models as new data becomes available.

5 Conclusions

In this paper, we present the ATOM Modeling PipeLine, or AMPL . This open-source software suite allows the user to build global and local models for a wide array of molecular properties that are needed for in silico drug discovery. Results of extensive benchmarking on a wide variety of pharmacokinetic and safety datasets were also presented, with an exploration of the effects of different featurization and model types on model accuracy. While the datasets used for developing and testing the pipeline are not publicly available, the software used to curate data and train, evaluate, and share new models is available as open source and benefits from having been tested on a wide array of pharmaceutically-relevant parameters. Additional public datasets are included with the pipeline release to support applying reproducible training and testing protocols that enable the broader modeling community to evaluate and improve modeling approaches over time.

6 Disclaimer

This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344.

7 Funding Sources

This work represents a multi-institutional effort. Funding sources include: Lawrence Livermore National Laboratory internal funds; the National Nuclear Security Administration; GlaxoSmithKline, LLC; and federal funds from the National Cancer Institute, National Institutes of Health, and the Department of Health and Human Services, under Contract No. 75N91019D00024.

8 Appendix

9 Benchmarking of AMPL on public datasets

AMPL is open source and available for download at
http://github.com/ATOMconsortium/AMPL. To support reproducibility of this pipeline, we provide model-building examples for three public datasets in AMPL’s open source repository. These datasets include:

  • Delaney et al. solubility dataset delaney

  • Wenzel et al. human liver microsome intrinsic clearance clearance

  • Drug Target Commons KCNH2 (hERG) inhibition assay herg

Since the data from our main benchmarking experiments are proprietary, we also benchmarked AMPL on these publicly-available datasets. Results are shown below.

{adjustbox}

max width= Dataset Model and featurizer type Train set Validation set Test set Delaney solubility Neural network + ECFP 0.66 0.21 0.29 Delaney solubility Neural network + GraphConv 0.76 0.55 0.54 Delaney solubility Neural network + Mordred 0.79 0.67 0.74 Delaney solubility Random forest + ECFP 0.91 0.27 0.37 Delaney solubility Random forest + Mordred 0.99 0.72 0.73 Wenzel microsomal clearance Neural network + ECFP 0.19 0.07 0.054 Wenzel microsomal clearance Neural network + GraphConv 0.11 0.064 0.067 Wenzel microsomal clearance Neural network + Mordred 0.40 0.21 0.13 Wenzel microsomal clearance Random forest + ECFP 0.90 0.17 0.21 Wenzel microsomal clearance Random forest + Mordred 0.92 0.15 0.12 KCNH2 (hERG) inhibition Neural network + ECFP 0.30 0.22 0.15 KCNH2 (hERG) inhibition Neural network + GraphConv 0.28 0.19 0.18 KCNH2 (hERG) inhibition Neural network + Mordred 0.24 0.20 0.19 KCNH2 (hERG) inhibition Random forest + ECFP 0.90 0.36 0.38 KCNH2 (hERG) inhibition Random forest + Mordred 0.94 0.39 0.36

Table 4: scores for public dataset models

References

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
398239
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description