Chemical Science View Article Online EDGE ARTICLE View Journal | View Issue DeltaDelta neural networks for lead optimization of This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. small molecule potency† Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. Cite this: Chem. Sci., 2019, 10, 10911 All publication charges for this article José Jiménez-Luna, a Laura Pérez-Benito, bc Gerard Martı́nez-Rosell,f have been paid for by the Royal Society Simone Sciabola,d Rubben Torella,e Gary Tresadern c and Gianni De Fabritiis *afg of Chemistry The capability to rank different potential drug molecules against a protein target for potency has always been a fundamental challenge in computational chemistry due to its importance in drug design. While several simulation-based methodologies exist, they are hard to use prospectively and thus predicting potency in lead optimization campaigns remains an open challenge. Here we present the first machine learning approach specifically tailored for ranking congeneric series based on deep 3D-convolutional neural networks. Furthermore we prove its effectiveness by blindly testing it on datasets provided by Janssen, Pfizer and Biogen totalling over 3246 ligands and 13 targets as well as several well-known openly available sets, representing one the largest evaluations ever performed. We also performed online Received 13th September 2019 Accepted 15th October 2019 learning simulations of lead optimization using the approach in a predictive manner obtaining significant advantage over experimental choice. We believe that the evaluation performed in this study is strong DOI: 10.1039/c9sc04606b evidence of the usefulness of a modern deep learning model in lead optimization pipelines against more rsc.li/chemical-science expensive simulation-based alternatives. considered a low error in other scenarios (e.g. below 1 Introduction 1 kcal mol1) is not a guarantee for successful ranking. In the lead optimization phase of drug discovery, the chemical It is therefore common to focus on relative binding free structure of a molecule is typically modied by a medicinal energy (RBFE) simulation methods,4–13 where the difference in chemistry team with the intent of improving its potency, affinity between two ligands is computed using a thermody- selectivity, and many other pharmacokinetic and toxicological namic cycle that alchemically perturbs only the small region parameters.1–3 These modications result in congeneric series, associated with the changing atoms. RBFE methods have shown a set of ligands with few atom changes between them, usually good results in several studies, with accuracy close to around a unique or small number of different scaffolds for 1 kcal mol1 and reasonable correlations. Despite this, these which there are experimental structures of the complex with the methods suffer from several issues, such as system preparation, target protein. Series range from few hundred to thousands of treatment of waters, force-eld selection, protein exibility and compounds and require considerable human, time and nan- computational cost, making their prospective application cial resources for synthesis and assays. It is therefore of great difficult in practice.14 On the other side, many empirical,15,16 value to have in silico predictive tools to accelerate this process. knowledge-based17,18 and machine learning19–24 scoring func- Series typically feature very small potency differences, which in tions have been designed for the task of predicting absolute turn is a challenge for predictors, as having what could be binding affinities. They mostly tackle the problem in a regres- sion setup, where the binding affinity is to be predicted using a Computational Science Laboratory, Parc de Recerca Biomèdica de Barcelona, a set of protein–ligand descriptors, modelling the interaction Universitat Pompeu Fabra, C Dr Aiguader 88, Barcelona, 08003, Spain. E-mail: among both. The fact that they model absolute affinities and are gianni.defabritiis@upf.edu trained on very chemically diverse bodies of data, such as iter- b Laboratori de Medicina Computacional, Unitat de Bioestadı́stica, Facultat de ations of the PDBbind25 database, limits their applicability Medicina, Universitat Autònoma de Barcelona, Spain c when predicting small structural differences between two Janssen Research and Development, Turnhoutseweg 30, 2340 Beerse, Belgium d ligands, such in the congeneric series case. While other Biogen Chemistry and Molecular Therapeutics, 115 Broadway Street, Cambridge, MA 02142, USA machine learning approaches have been presented for this e Pzer I&I, 610 Main Street, Cambridge, MA 02139, USA task,26–28 here we propose a modern 3D-convolutional-neural- f Acellera, Carrer del Dr Trueta, 183, 08005 Barcelona, Spain network-based continuous learning approach for relative g Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluis Companys binding affinity prediction in congeneric series and show strong 23, 08010 Barcelona, Spain predictive power using multiple blind benchmarks as well as † Electronic supplementary information (ESI) available. See DOI: public datasets at negligible computational costs. This study 10.1039/c9sc04606b This journal is © The Royal Society of Chemistry 2019 Chem. Sci., 2019, 10, 10911–10918 | 10911 View Article Online Chemical Science Edge Article serves as a very large evaluation of a modern machine-learning 2.2 Descriptor calculation pipeline for lead optimization in a real-life drug discovery We have recently reported a machine learning approach that scenario, thanks to the joint collaboration with several phar- can learn based on 3D features of the binding site interac- maceutical companies. tions.20 A similar encoding was used here that represents the protein–ligand binding by voxelizing both using a 24 Å pocket 2 Materials and methods centered box, with a density of 1 Å3 per voxel. The contribution 2.1 Datasets studied of each atom to each voxel is inversely proportional to their Euclidean distance r and the van der Waals radius rvdw of the The BindingDB protein–ligand validation sets29 were used to rst. We use several channels for both protein and ligand, in the This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. pretrain our models. For testing, we also extracted well-known sense that the atomic contribution to each voxel depends on Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. publicly-available literature test sets30 used for benchmarking their type. The contribution of each atom to each voxel is RBFE calculations. Furthermore we include a recent freely- assigned according to a pair correlation function dened by: available BRD4 bromodomain dataset.31 In regards to internal    pharmaceutical data, we tested on ve different congeneric rvdw 12 nðrÞ ¼ 1  exp  : (1) series from Janssen R&D. Three chemical series (sets 1, 2 and 3) r were phosphodiesterase 2 (PDE2) inhibitors with bioactivity versus PDE2, PDE3, and PDE10 (ref. 32 and 33) (publication We dene several channels for both protein and ligand, in the number WO2018083103A1), the fourth series were proto- sense that the atomic contribution to each voxel depends on oncogene tyrosine kinase (ROS1) inhibitors (publication their type. For the protein we dene eight pharmacophoric-like number WO2015144799A1) and the nal beta-secretase 1 descriptors, as detailed in ESI Table 3.† For the ligands we use (BACE1) inhibitors.34 We tested six congeneric series with Pzer, a simpler representation based on atom types contained in the three of which target a kinase, and the remaining an enzyme, set {C,N,O,F,P,S,Cl,Br,I,H}, for a total of 18 stacked channels. a phosphodiesterase (PDE) and an activator of transcription. We note that there is no particular reasoning behind this choice The sizes of these vary from 93 molecules up to 362, for a total of of descriptors other than they showed promising practical 955 tested compounds. Lastly, Biogen tested the proposed performance in previous studies. procedure on two different series, composed of 196 and 220 analogues targeting a tyrosine-protein kinase and a receptor- associated kinase, respectively. The size of the sets presented 2.3 Neural network architecture here allow, to the best of our knowledge the largest evaluation Regular feed-forward neural networks do not scale well when yet of a modern machine learning pipeline in lead optimization. the input is high dimensional (as in images, or in this case Out of the total 645 available congeneric series available in atomic interactions). CNNs are specically designed for BindingDB, 495 with IC50 affinity values were extracted and handling lattices, where local spatial information needs to be processed for further evaluation, as it was the unit with most preserved. While a feed-forward network would ignore such data available, containing a diverse set of targets. The majority interactions, a convolutional one arranges its neurons spatially, of these sets encompass a single protein–ligand crystal struc- and only connects locally to the output of the previous layer. The ture, the rest of the ligands modelled against the reference latter type of architectures have become the de-facto workhorse using the Surex docking soware.35 We then assign each in computer vision problems,39–41 providing state-of-the-art protein structure in the database to a family cluster using a 90% performance. Due to this success, many applications in bio- sequence similarity threshold, as per PDB conventions.36 For informatics and computational chemistry followed.42–51 each series in the same protein cluster we use a maximum The neural network we propose has a novel zero-symmetric common substructure (MCS) protocol as available in rdkit37 to architecture whose main building blocks are 3D-convolution remove identical ligands. This procedure ensures that the same operations. In this work we focus on predicting relative affini- ligand is not repeated against similar targets, avoiding potential ties for close analogues in lead optimization, therefore, our overtting problems and overoptimistic evaluations.38 Affinity approach is to build a network whose input is a pair of ligand values were log-converted to avoid target scaling issues (pIC50 ¼ binding voxelized representations belonging to the same series. log10 IC50). Ligands that could not be read by rdkit were A two-input convolutional neural network was designed, with removed. Histograms of the number of ligands and their affinity xed weights on both legs. The inputs are forwarded through range per series are provided in ESI Fig. 7,† with the average several convolution and pooling operations and then attened available number of ligands per series being 8.84. In the into a 192-dimensional latent representation. The symmetry Schrödinger and BRD4 sets, since only DG (per kcal mol1) property of relative binding affinity requires that inverting the information was available, we converted affinity values to the order of the input ligands should change the sign of the pre- pIC50 range assuming non-competitive binding. Descriptive dicted value, and we embed such symmetry in the network by information on these series is provided in ESI Table 1.† computing the difference between latent representations. A Compounds provided by Janssen were docked using a common nal linear layer with no bias is then applied to the result of this scaffold structure via the Glide soware. These congeneric difference, ensuring zero-symmetry by design and producing series range from 48 up to a 900 different compounds with the desired predicted difference in affinity. In contrast, calcu- varying affinity ranges (ESI Table 2†). lating relative affinities from an absolute predictor inevitably 10912 | Chem. Sci., 2019, 10, 10911–10918 This journal is © The Royal Society of Chemistry 2019 View Article Online Edge Article Chemical Science leads to the concatenation of errors from two separate the series are then sequentially added in a random order. It is predictions. well known that either a random53 or scaffold-based training A schema of our architectural choice is provided in Fig. 1 is test split produce overoptimistic results when testing machine- provided in the ESI.† It features two convolution operations learning algorithms on activity benchmarks. Since the indus- with a kernel size of 3 in each leg, followed by a max-pooling trial datasets in our study include a compound creation time- operation, and nally another convolution operation with the stamp, we also evaluate a more realistic temporal split,54 same kernel size for both before attening and performing the where at each training step we consider the rst n tested ligands latent difference between analogues. The ReLU activation and the differences of the posterior ones against the rst are function was used for all layers in the network except for the taken. The performance of the machine learned models is re- This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. last, which does not feature one. We include a dropout layer in ported as the root mean squared error (RMSE) and either Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. the end to control for overtting. Xavier initialization was used Pearson's correlation coefficient R or Spearman's r between for the weights. Training is performed using the Adam experimental and predicted affinity differences. We note that in stochastic gradient descent optimizer52 with standard hyper- all blind tests a single model was provided, and no explicit parameters (b1 ¼ 0.9, b2 ¼ 0.999, 3 ¼ 104) using a batch size of attempt to optimize hyperparameters in each set was made. 32 samples for 50 epochs. Data augmentation is performed during training by randomly rotating input coordinates to 3 Results mitigate the lack of rotational equivariance in CNNs. Further- more, given a set of relative binding predictions, its absolute We rst present results concerning our validation on the 495 counterparts can always be retrieved given a single experimen- protein-series datasets from the BindingDB, where the tally determined absolute reference, such as the one provided proposed model achieves an average correlation coefficient by a lead. If more than one is available, absolute affinities can be above 0.4 and an RMSE below 1.25 (pIC50 units) even when only computed towards each, in practice providing a predictive one binding-energy difference is taken per congeneric series absolute affinity distribution, whose average can then be (ESI Fig. 1†). This suggests that the method works reasonably interpreted as a maximum a posteriori (MAP) estimate of the well in the very low-data scenario, such as the beginning of absolute affinity and its standard deviation as a measure of its a lead optimization campaign. A noticeable performance boost uncertainty, given the current model state. is seen as more differences are included in training, with a correlation coefficient above 0.62 and an RMSE below 1.05 2.4 Continuous learning when another four different ligands from the same congeneric series are known in advance, with performance plateauing In the proposed continuous-learning approach, we explicitly beyond ve additional training ligands. A comparison against use the fact that congeneric series are sequentially generated in an absolute affinity model is also provided (i.e. one of the legs of a lead optimization campaign, and follow an incremental the architecture), where as expected it can be appreciated that it training and testing procedure. For each congeneric series at performs considerably worse than its relative counterpart. a given time the affinity of previously tested ligands is known Now we present results on the Wang et al.30 and BRD4 experimentally: differences for these are taken as training data, inhibitor datasets.31 In this and the rest of cases, we pretrained while for test data we predict differences between unknown and a model with all difference pairs available in the BindingDB known ones. While this approach is less ambitious than having database, which provides a prior for further ne-tuning. We a predictor for relative affinity with no experimentally tested then mixed new available data as training in each sequential data (such as a physical-based model), its applicability is iteration of each set with the rest of the BindingDB database for general, since it is the common scenario that medicinal only 3 epochs, signicantly reducing computational overhead. A chemists face in lead optimization campaigns. The training for FEP baseline provided by Wang et al.30 is used for comparison. the BindingDB sets starts with a reference structure in each The model efficiently interpolates differences for unseen series, for which we take the crystal structure ligand if available ligands, achieving considerably high correlation coefficients or the structure with the lowest average maximum common and low errors in all series with as few as 3–4 additional ligands substructure (MCS) distance to the rest. Ligands from the rest of and associated activity pairs, surpassing in many cases the much more expensive FEP baseline (ESI Fig. 2†). For instance, for the MCL1 target, aer testing 3 ligands, the correlation coefficient is above 0.8, surpassing the FEP baseline, and the RMSE is below 1.2 (pIC50 units). The same evaluation procedure was taken for the compounds available in the Janssen PDE sets (Fig. 2 and ESI Fig. 3†) for both a random and a temporal split, where a base- line against Glide score55 is also added. Excellent performance was seen on a random split given enough training data, and as Fig. 1 Architecture of the proposed model. A two-legged neural network with tied weights was constructed, and a pair of protein– expected, although the temporal split performance is lower, it is ligand voxelization is feed-forwarded through it to later perform still sufficiently high to be used in a real-life prospective lead a latent space difference. optimization scenario. For instance, for the rst PDE2 activity This journal is © The Royal Society of Chemistry 2019 Chem. Sci., 2019, 10, 10911–10918 | 10913 View Article Online Chemical Science Edge Article This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. Fig. 2Average Pearson's correlation coefficient R (1 standard deviation) based on 25 independent runs on different sets for the Janssen PDE2, PDE3 and PDE10 targets. set aer 20 ligands sorted by time, the Pearson's correlation results on ngerprint similarity for the rst PDE2 set (ESI coefficient R and RMSE were 0.77 and 1.35 (in pIC50 units) Fig. 5†), show that aer 20 ligands sorted by chemical similarity respectively. Results for the ROS1 and BACE sets, show a similar the R and RMSE were 0.83 and 1.12 (in pIC50 units). These trend and insights (Fig. 3 and ESI Fig. 4†). Interestingly, suggest better performance in this scenario than the proposed performance slightly decreases for the BACE target late in the temporal split, and closer to the random one. evaluation, suggesting that the network may have found an We then present the results provided by Pzer using unexplored activity cliff or that the last tested ligands are harder a temporal split in Table 1, where specic target names cannot to predict than others in the series. Furthermore, we also be disclosed. We compare such results with several baselines provide a type of split where only differences among the most such as molecular weight, clog P, a MM-GBSA pipeline56,57 and chemically close ligands are predicted, based on ECFP4 deep-learning absolute affinity predictor KDEEP,20 trained on the ngerprint similarity, as available in rdkit. That is, in each v.2016 iteration of the PDBbind database. The model proposed training step we predict from the remaining untested pool of here performs considerably better than the rest when given only ligands those that are closest to the ones in our training set, 10% of the training data, again highlighting the importance of with the intention of resembling a real-life lead optimization incrementally training these on the congeneric series of RBFE scenario, typically applied to close analogues. Split-based interest. An exception, however, is found in the Kinase #3 series, 10914 | Chem. Sci., 2019, 10, 10911–10918 This journal is © The Royal Society of Chemistry 2019 View Article Online Edge Article Chemical Science This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. Fig. 3 Average Pearson's correlation coefficient R (1 standard deviation) based on several independent runs on two sets for the Janssen ROS1 and BACE targets. for which no signicant improvement is observed when high potency compounds in a close to real-life lead optimization providing extra training data. This particular last case appears scenario, by retrospectively comparing them to the experimental to be particularly hard to predict, as all tested methods perform order of synthesis. With some of the large industrial datasets it poorly. We provide results using a temporal split for the last two was possible to test this and we used the most active compound congeneric series provided by Biogen, for which we also as a surrogate interesting lead molecule. The model is trained on compare against several baselines: (a) Glide score, (b) an MM- the rst experimentally tested compounds, and then is incre- GBSA pipeline, and (c) a standard QSAR approach using mentally trained by choosing from the remaining ones based on MACCS, ECFP4 and rdkit descriptors with a random forest a upper condence bound (UCB)-like criterion,58 dened as: model (Fig. 4 and ESI Fig. 6†). Our model reveals similar conclusions, signicantly outperforming all baselines. Curi- UCB ¼ m(x) + bs(x), (2) ously, it can also be seen that the proposed method does not perform signicantly worse than the aforementioned baselines where m and s are the average and standard deviation predicted in the second target when no training data is used. When some absolute affinities provided by the model for ligand x and b is is used, such as only 5 analogues, our proposed machine- a user-chosen factor controlling the balance between exploita- learning model signicantly outperforms all baselines. tion and exploration, that we x in our study to b ¼ 1.64. All the tests regarding internal pharmaceutical data were We stop the procedure once the model retrieves the analogue carried out blindly by providing fully-containerized soware to with the highest associated affinity, and compare this with its our collaborators, who executed the application and reported original synthesis experimental order in its corresponding corresponding results. Furthermore only one pretrained model series. We present results for this simulation-based benchmark was provided without any opportunity to overt to each specic in Table 2. In 4 out of 5 sets our proposed model is able to reach test set. the compound with the highest affinity faster than its experi- One aim of our study was to test whether machine-learning mental order or by random selection. Surprisingly, in all ten driven relative affinity predictions could efficiently identify key independent runs of the second set for the PDE2 target, the Table 1 Spearman's r performance results between experimental and predicted absolute affinities provided by Pfizer I&I, where other empirical, simulation, and machine-learning based affinity prediction methods are compared on several congeneric series. Performance is poor for most tested model except for the sequential approach proposed here, with Pearson correlations averaging over 0.5 with as few as 10% used analogues from the congeneric series at hand Mol. weight This work This work This work Target # ligands (r) clog Pa (r) MM-GBSA (r) KDEEP (r) (10% training, r) (20% training, r) (30% training, r) Kinase #1 362 0.19 0.06 0.56 0.42 0.49 0.64 0.73 Kinase #2 106 0.1 0.28 0.25 0.25 0.25 0.41 0.51 Kinase #3 95 0 0.04 0.25 0.27 0.3 0.3 0.31 Enzyme 93 0.43 0.24 0.01 0.49 0.43 0.26 0.59 Phosphodiesterase 100 0.37 0.36 0.67 0 0.49 0.64 0.73 Activator of transcriptions 199 0.13 0.08 0.66 0.29 0.72 0.84 0.94 Weighted avg. 0.19 0.14 0.47 0.25 0.49 0.59 0.69 Simple avg. 0.2 0.18 0.4 0.18 0.45 0.52 0.64 a Calculated log P as available in rdkit. This journal is © The Royal Society of Chemistry 2019 Chem. Sci., 2019, 10, 10911–10918 | 10915 View Article Online Chemical Science Edge Article This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. Fig. 4 Results over 5 runs on Biogen's tyrosine-protein kinase and receptor-associated kinase using a temporal split, and MM-GBSA and QSAR random forest pipelines as baselines. compound with the highest affinity was found aer only a single simulation-based benchmark. The accuracy of the method synthesis epoch. Furthermore, one would expect the average heavily depends on the amount of available data but can be affinity in the training set to increase at each synthesis epoch (as trained and applied in minutes on a single GPU, offering the model is tasked to pick compounds with increasingly higher a substantial improvement in performance compared with UCB). This is the case for 4 out of 5 sets again (Fig. 5), with the physics-based RBFE calculations which can take days for exception of the ROS1 target, which shows a non-monotonic a small number of analogues. While the results presented here trend, albeit its model reaches the compound with highest are encouraging, it is important to note that they remain affinity before its experimental order. In all tested cases, the retrospective: a proper prospective validation of the model, average training pool affinity for the ligands selected by the which would entail chemists synthesizing compounds accord- model is higher than experimental choice. Overall results are ing to the decisions taken by the trained model, remains a topic very promising and suggest that the proposed method could be of future study. In the long term, however, we expect that applied successfully in a prospective scenario. Particularly, in improving molecular simulations accuracy59,60 by the integra- the rst PDE2 set, we were able to reach potent compounds tion of physics and machine learning approaches would synthesizing up to six times less molecules than the baseline produce a more convenient approach for engineering drug method used by the medicinal chemistry team. discovery. In the meantime, methods such as the one proposed here provide accurate performance at a fraction of the compu- tational cost of other approaches. 4 Discussion In this work we have designed and tested a deep-learning based Code & data availability model for the task of predicting relative binding affinity predictions in congeneric series. This work provides evidence All the models here were developed using the PyTorch package that the method is able to efficiently rank compounds as shown for tensor computation and neural network training.61 Bind- by an evaluation on both publicly available and industrial data ingDB, Wang et al. and Mobley et al. set results are available and can be of use by computational and medicinal chemists in upon reasonable request. Python code for generating the early drug-discovery projects by providing informed choices of proposed featurization is available within the open-source future compounds to synthesize, as suggested by our HTMD soware.62 The code of the network architecture in Table 2 Simulation-based benchmark results over 10 independent runs for the different datasets. We show the amount of molecules the model is allowed to pick at each synthesis epoch, the experimental order of the compound with the highest affinity in the series, the average synthesis epoch our model found said molecule, the total necessary sampled ligands the proposed model has chosen before the target compound, and the sampling advantages over the experimental and random orders Chosen per synthesis Experimental Found at synthesis Total sampled Advantage over Advantage over Target Set # ligands epoch order epoch ligands experimental choice random choice PDE2 1 900 10 766 12.2 132 634 318 PDE2 2 303 10 61 1 20 41 131.5 PDE2 3 278 10 253 5.9 69 184 70 ROS1 — 165 10 73 3.1 41 32 41.5 BACE — 229 10 190 20.8 218 28 103.5 10916 | Chem. Sci., 2019, 10, 10911–10918 This journal is © The Royal Society of Chemistry 2019 View Article Online Edge Article Chemical Science This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. Fig. 5 Average model-picked training set affinity per number of compounds synthesized for the Janssen PDE2, ROS1 and BACE sets, as well as a baseline based on the actual experimental choice order of compounds. a PyTorch implementation is provided in the ESI.† An imple- 4 L. Wang, B. J. Berne and R. A. Friesner, Proc. Natl. Acad. Sci. mentation of this application is available through the http:// U. S. A., 2012, 109, 1937–1942. PlayMolecule.org repository of applications, where users can 5 E. B. Lenselink, J. Louvel, A. F. Forti, J. P. D. van Veldhoven, freely submit their protein in PDB format and two sets of the H. de Vries, T. Mulder-Krieger, F. M. McRobb, A. Negri, same congeneric series, for training and validation respectively J. Goose, R. Abel, H. W. T. van Vlijmen, L. Wang, in SDF format. Depending on the size of these last two, training E. Harder, W. Sherman, A. P. IJzerman and T. Beuming, and prediction time may vary, as the order of data for training ACS Omega, 2016, 1, 293–304. nðn  1Þ 6 S. Wan, A. P. Bhati, S. Skerratt, K. Omoto, increases by , and for testing nm factors, where n and m 2 V. Shanmugasundaram, S. K. Bagal and P. V. Coveney, J. are the number of training and testing instances respectively. At Chem. Inf. Model., 2017, 57, 897–909. the moment, predictions are limited to a default total of a 1000 7 D. A. Goldfeld, R. Murphy, B. Kim, L. Wang, T. Beuming, molecules per congeneric series, with runtimes averaging and R. Abel and R. A. Friesner, J. Phys. Chem. B, 2015, 119, 824–835. hour on a modern GeForce 1080Ti GPU. Larger experiments can 8 L. Pérez-Benito, H. Keränen, H. van Vlijmen and be arranged for users willing to run more computationally G. Tresadern, Sci. Rep., 2018, 8, 4883. demanding experiments. 9 M. Ciordia, L. Pérez-Benito, F. Delgado, A. A. Trabanco and G. Tresadern, J. Chem. Inf. Model., 2016, 56, 1856–1871. 10 C. Schindler, F. Rippmann and D. Kuhn, J. Comput.-Aided Conflicts of interest Mol. Des., 2017, 32, 1–8. G. D. F is a founder and current CEO of Acellera Ltd. J. J. receives 11 H. Keränen, L. Pérez-Benito, M. Ciordia, F. Delgado, funding from Acellera Ltd. and G. M. R. is an employee. T. B. Steinbrecher, D. Oehlrich, H. W. T. van Vlijmen, A. A. Trabanco and G. Tresadern, J. Chem. Theory Comput., 2017, 13, 1439–1453. Acknowledgements 12 G. Heinzelmann, N. M. Henriksen and M. K. Gilson, J. Chem. Theory Comput., 2017, 13, 3260–3275. The authors thank Acellera Ltd. for funding. G. D. F. acknowl- 13 M. Aldeghi, A. Heifetz, M. J. Bodkin, S. Knapp and edges support from MINECO (BIO2014-53095-P), MICINN (PTQ- P. C. Biggin, Chem. Sci., 2016, 7, 207–218. 17-09079) and FEDER. This project has received funding from 14 Z. Cournia, B. Allen and W. Sherman, J. Chem. Inf. Model., the European Union's Horizon 2020 research and innovation 2017, 57, 2911–2937. programme under grant agreement No 675451 (CompBioMed 15 Y. Cao and L. Li, Bioinformatics, 2014, 30, 1674–1680. project). 16 M. P. Brenner, L. J. Colwell, et al., Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 13564–13569. Notes and references 17 T. A. Halgren, R. B. Murphy, R. A. Friesner, H. S. Beard, L. L. Frye, W. T. Pollard and J. L. Banks, J. Med. Chem., 1 C. A. Nicolaou and N. Brown, Drug Discovery Today: Technol., 2004, 47, 1750–1759. 2013, 10, e427–e435. 18 O. Trott and A. Olson, J. Comput. Chem., 2010, 31, 455–461. 2 I. Kola and J. Landis, Nat. Rev. Drug Discovery, 2004, 3(8), 711. 19 P. J. Ballester and J. B. O. Mitchell, Bioinformatics, 2010, 26, 3 S. Ekins, J. D. Honeycutt and J. T. Metz, Drug discovery today, 1169–1175. 2010, 15, 451–460. This journal is © The Royal Society of Chemistry 2019 Chem. Sci., 2019, 10, 10911–10918 | 10917 View Article Online Chemical Science Edge Article 20 J. Jiménez, M. Škalič, G. Martı́nez-Rosell and G. De Fabritiis, 39 A. Krizhevsky, I. Sutskever and G. E. Hinton, Advances in J. Chem. Inf. Model., 2018, 58, 287–296. Neural Information Processing Systems, 2012, pp. 1–9. 21 E. N. Feinberg, D. Sur, Z. Wu, B. E. Husic, H. Mai, Y. Li, 40 O. Vinyals, A. Toshev, S. Bengio and D. Erhan, Proceedings of S. Sun, J. Yang, B. Ramsundar and V. S. Pande, ACS Cent. the IEEE Computer Society Conference on Computer Vision and Sci., 2018, 4, 1520–1530. Pattern Recognition, 2015, pp. 3156–3164. 22 M. Ragoza, J. Hochuli, E. Idrobo, J. Sunseri and D. R. Koes, J. 41 A. Radford, L. Metz and S. Chintala, ICLR, 2016, pp. 1–16. Chem. Inf. Model., 2017, 57, 942–957. 42 M. Skalic, J. Jiménez, D. Sabbadin and G. De Fabritiis, J. 23 D. D. Nguyen, Z. Cang, K. Wu, M. Wang, Y. Cao and Chem. Inf. Model., 2019, 59(3), 1205–1214. G.-W. Wei, J. Comput.-Aided Mol. Des., 2019, 33, 71–82. 43 B. Ramsundar, B. Liu, Z. Wu, A. Verras, M. Tudor, This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. 24 Z. Gaieb, C. D. Parks, M. Chiu, H. Yang, C. Shao, R. P. Sheridan and V. Pande, J. Chem. Inf. Model., 2017, 57, Open Access Article. Published on 16 October 2019. Downloaded on 3/12/2020 10:41:49 AM. W. P. Walters, M. H. Lambert, N. Nevins, S. D. Bembenek, 2068–2076. M. K. Ameriks, et al., J. Comput.-Aided Mol. Des., 2019, 33, 44 Z. Wu, B. Ramsundar, E. N. Feinberg, J. Gomes, C. Geniesse, 1–18. A. S. Pappu, K. Leswing and V. Pande, Chem. Sci., 2018, 9, 25 R. Wang, X. Fang, Y. Lu and S. Wang, J. Med. Chem., 2004, 47, 513–530. 2977–2980. 45 J. Jiménez, S. Doerr, G. Martı́nez-Rosell, A. S. Rose and G. De 26 W. Zhan, D. Li, J. Che, L. Zhang, B. Yang, Y. Hu, T. Liu and Fabritiis, Bioinformatics, 2017, 33, 3036–3042. X. Dong, Eur. J. Med. Chem., 2014, 75, 11–20. 46 M. Skalic, A. Varela-Rial, J. Jiménez, G. Martı́nez-Rosell and 27 A. Amini, P. J. Shrimpton, S. H. Muggleton and G. De Fabritiis, Bioinformatics, 2018, 35, 243–250. M. J. Sternberg, Proteins: Struct., Funct., Bioinf., 2007, 69, 47 M. Skalic, A. Varela-Rial, J. Jiménez, G. Martı́nez-Rosell and 823–831. G. De Fabritiis, Bioinformatics, 2019, 35(2), 243–250. 28 D. Zilian and C. A. Sotriffer, J. Chem. Inf. Model., 2013, 53, 48 C. Wehmeyer and F. Noé, J. Chem. Phys., 2018, 148, 241703. 1923–1933. 49 G. Derevyanko, S. Grudinin, Y. Bengio and G. Lamoureux, 29 T. Liu, Y. Lin, X. Wen, R. N. Jorissen and M. K. Gilson, Nucleic Bioinformatics, 2018, 34, 4046–4053. Acids Res., 2006, 35, D198–D201. 50 M. H. Segler, M. Preuss and M. P. Waller, arXiv preprint 30 L. Wang, Y. Wu, Y. Deng, B. Kim, L. Pierce, G. Krilov, arXiv:1708.04202, 2017. D. Lupyan, S. Robinson, M. K. Dahlgren, J. Greenwood, 51 R. Gómez-Bombarelli, J. N. Wei, D. Duvenaud, D. L. Romero, C. Masse, J. L. Knight, T. Steinbrecher, J. M. Hernández-Lobato, B. Sánchez-Lengeling, T. Beuming, W. Damm, E. Harder, W. Sherman, D. Sheberla, J. Aguilera-Iparraguirre, T. D. Hirzel, M. Brewer, R. Wester, M. Murcko, L. Frye, R. Farid, T. Lin, R. P. Adams and A. Aspuru-Guzik, ACS Cent. Sci., 2018, 4, D. L. Mobley, W. L. Jorgensen, B. J. Berne, R. A. Friesner 268–276. and R. Abel, J. Am. Chem. Soc., 2015, 137, 2695–2703. 52 D. P. Kingma and J. L. Ba, International Conference on 31 D. L. Mobley and M. K. Gilson, Annu. Rev. Biophys., 2017, 46, Learning Representations 2015, 2015, pp. 1–15. 531–558. 53 C. Kramer and P. Gedeck, J. Chem. Inf. Model., 2010, 50, 32 F. J. Rombouts, G. Tresadern, P. Buijnsters, X. Langlois, 1961–1969. F. Tovar, T. B. Steinbrecher, G. Vanhoof, M. Somers, 54 R. P. Sheridan, J. Chem. Inf. Model., 2013, 53, 783–790. J.-I. Andrés and A. A. Trabanco, ACS Med. Chem. Lett., 55 R. A. Friesner, J. L. Banks, R. B. Murphy, T. A. Halgren, 2015, 6, 282–286. J. J. Klicic, D. T. Mainz, M. P. Repasky, E. H. Knoll, 33 P. Buijnsters, M. De Angelis, X. Langlois, F. J. R. Rombouts, M. Shelley, J. K. Perry, D. E. Shaw, P. Francis and W. Sanderson, G. Tresadern, A. Ritchie, A. A. Trabanco, P. S. Shenkin, J. Med. Chem., 2004, 47, 1739–1749. G. VanHoof, Y. V. Roosbroeck and J.-I. Andrés, ACS Med. 56 S. Genheden and U. Ryde, Expert Opin. Drug Discovery, 2015, Chem. Lett., 2014, 5, 1049–1053. 10, 449–461. 34 F. J. R. Rombouts, G. Tresadern, O. Delgado, C. Martı́nez- 57 T. Hou, J. Wang, Y. Li and W. Wang, J. Chem. Inf. Model., Lamenca, M. Van Gool, A. Garcı́a-Molina, S. A. Alonso de 2010, 51, 69–82. Diego, D. Oehlrich, H. Prokopcova, J. M. Alonso, N. Austin, 58 N. Srinivas, A. Krause, S. M. Kakade and M. Seeger, arXiv H. Borghys, S. Van Brandt, M. Surkyn, M. De Cleyn, A. Vos, preprint arXiv:0912.3995, 2009. R. Alexander, G. Macdonald, D. Moechars, H. Gijsen and 59 K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko A. A. Trabanco, J. Med. Chem., 2015, 58, 8216–8235. and K.-R. Müller, J. Chem. Phys., 2018, 148, 241722. 35 R. Spitzer and A. N. Jain, J. Comput.-Aided Mol. Des., 2012, 26, 60 J. S. Smith, O. Isayev and A. E. Roitberg, Chem. Sci., 2017, 8, 687–699. 3192–3203. 36 C. Camacho, G. Coulouris, V. Avagyan, N. Ma, 61 A. Paszke, G. Chanan, Z. Lin, S. Gross, E. Yang, L. Antiga and J. Papadopoulos, K. Bealer and T. L. Madden, BMC Bioinf., Z. Devito, Advances in Neural Information Processing Systems, 2009, 10, 1. 2017, vol. 30, pp. 1–4. 37 G. Landrum, Online, http://www.rdkit.org, 2006. 62 S. Doerr, M. J. Harvey, F. Noé and G. De Fabritiis, J. Chem. 38 C. Kramer and P. Gedeck, J. Chem. Inf. Model., 2010, 50, Theory Comput., 2016, 12, 1845–1852. 1961–1969. 10918 | Chem. Sci., 2019, 10, 10911–10918 This journal is © The Royal Society of Chemistry 2019