Prediction of chemical compounds properties using a deep learning model

The discovery of new medications in a cost-effective manner has become the top priority for many pharmaceutical companies. Despite decades of innovation, many of their processes arguably remain relatively inefficient. One such process is the prediction of biological activity. This paper describes a new deep learning model, capable of conducting a preliminary screening of chemical compounds in-silico. The model has been constructed using a variation autoencoder to generate chemical compound fingerprints, which have been used to create a regression model to predict their LogD property and a classification model to predict binding in selected assays from the ChEMBL dataset. The conducted experiments demonstrate accurate prediction of the properties of chemical compounds only using structural definitions and also provide several opportunities to improve upon this model in the future.


Introduction
Deep learning [24] has been successfully applied in a number of problem domains from natural language processing [39], medical imaging analysis [32] to finance [26]. Deep learning architectures are also successfully used for many predictive tasks in chemistry and biology domains [55]. One application of deep learning in the chemistry domain is to predict important properties of chemical compounds. It allows for the assessment of chemical compounds before committing to an expensive synthesis process [7,8,20,43].
Deep learning theory is based on deep neural networks (DNN) which consist of many layers. Each layer is comprised of a number of neurons. Higher levels of the DNN represent more complex concepts. To improve network performance, layers are often implemented using different methodologies. The overall topology of a neural network is selected based on the problem to be solved and often is tuned during an experimental phase.
In this paper, the DNN learns to represent compounds by their chemical descriptors. This opens a wide range of opportunities to build sophisticated machine learning applications for predicting different properties of chemical compounds.
This research study investigates how the pretrained autoencoder can be used for building classification and regression models for predicting the LogD property and target binding of chemical compounds using data obtained from different sources.
The remainder of the paper is outlined as follows. Section 2 provides background information which includes a description of predicting properties of chemical compounds and other machine learning models used for their prediction. Section 3 introduces the mathematics behind the developed autoencoder together with visualizations of the learning mechanisms. Section 4 provides detail on the data sets used and experiments carried out, while Sect. 5 describes the obtained results, followed by discussion and conclusions in Sect. 6.

LogD
Lipophilicity is possibly one of the most important physicochemical properties of a potential drug. It plays a role in solubility, absorption, membrane penetration, plasma protein binding, distribution, CNS penetration and partitioning into other tissues or organs such as the liver and has an impact on the routes of clearance. It is important in ligand recognition, not only to the target protein but also CYP450 interactions, HERG binding, and PXR mediated enzyme induction. Most drugs entering a market are designed for oral administration. The absorption of drugs can either be via passive diffusion across membranes or via carrier mediated transport. Carrier mediated transport is energy dependent and requires a specific transporter protein. In contrast passive diffusion does not require the presence of a specific carrier transporter protein and is less structure specific than carrier-mediated transport; there is a general dependence on lipophilicity for structurally diverse compounds. However, the relationship with LogD is nonlinear with an optimum of LogD 1-2.
Measurement of LogP can be undertaken in a variety of ways, the most common is the shake-flask method, which consists of dissolving some of the solute in question in a volume of octanol and water, shaking for a period of time, then measuring the concentration of the solute in each solvent. This can be time-consuming particularly if there is no quick spectroscopic method to measure the concentration of the molecule in the phases. A faster method of logP determination makes use of high-performance liquid chromatography.
However, the majority of known drugs contain ionizable groups, as shown in Fig. 1, which shows the distribution of small molecule drugs with DrugBank [53] and are likely to be charged at physiological pH and LogP only correctly describes the partition coefficient of neutral (uncharged) molecules. LogD, the distribution constant is a better descriptor of the lipophilicity of a molecule. This can be determined in a similar manner to LogP but instead of using water, the aqueous phase is adjusted to a specific pH using a buffer. LogD is thus pH dependent, hence one must specify the pH at which the logD was measured. Of particular interest is the logD at pH = 7.4 (the physiological pH of blood serum).
Usually, it is not practical to determine the LogD of every compound made experimentally (and it may be of interest to calculate logD prior to synthesis) and so calculated results are used.

Binding
The majority of drug-like small molecules are specifically designed to bind to protein targets involved in disease related pathways. The activity of molecules in a biological assay may be captured by a variety of different measures IC50, EC50, Ki, % inhib, etc., but most are a measure of a binding event in some manner. The biological results varies considerably in quality from single point high-throughput screening (HTS) data [15] to full dose response curves. Much of these data are captured in ChEMBL, a database of bioactive drug-like small molecules and abstracted bioactivities.

Models for comparison
The developed DNN model provides an accurate prediction of LogD and binding properties. To gauge its performance, a total of ten machine learning (ML) techniques were used. Taking into consideration a large variety of different experimental setups, model implementations and evaluation metrics, authors tried to summarize the results from a number of sources and provide a 'single' performance metric to compare each of the techniques. 1 • Logistic regression (LR) [41] is a very straightforward and popular classification model, which has a long and successful history of been used in many ML applications and statistical modelling. It uses a logistic function for weighting a linear combination of input parameters. LR models have been used in a large number of publications involving chemical data types. • Kernel ridge regression (KRR) [57] uses a modified approach to find a regression function by adding a bias, which causes a drop in variance. In other words a better prediction can be achieved by considering a slightly less fit model during the training process. • Another very popular and well recognized model is random forests (RF) [56]. It can be applied for both classification and regression problems and uses an ensemble of decision trees, which are trained on different subsets of the original data.
• Extreme gradient boosting (XGB) [14] tree is similar to RF and is an ensemble method. During each training step it constructs a new tree model, which in combination with previous models, minimizes the overall prediction error. It is a very popular approach in the ML community, and has been successfully applied to many problems, consistently providing accurate results. • Multitask neural network (MNN) [49] is a special modified neural network architecture for solving simultaneously multiple problems. There are a number of different designs of these neural networks. However, they are all constructed based on the principle that some fully-connected layers are shared between different tasks. In this way training processes for solving one task can influence other tasks and vice versa. • Graph convolutional (GC) [18]   The previously reviewed graph-based models are using undirected graphs (where molecular bonds between atoms naturally do not have directions). However, the latest DAG design synthetically introduces an additional directional feature for the graph-based molecular representation. It identifies a central atom and creates a directional-structure of chemical compound from this central point. Generating the additional features may provide an improvement in prediction accuracy in comparison to other graph-based models. • WEAVE [31] supports graph-based model that utilizes both properties of chemical compounds: nodes (for atoms) and edges (for bonds). The constructed features matrix is processed by convolution-like filters. Similar to image convolution neural network, WEAVE provides more informative representations of chemical compound structures. • Influence relevance voting (IRV) [51] is not the most common approach, but it has certain important advantages against other ML models. Its prediction can be easily interpreted in a similar fashion to k-nearest neighbours (kNN). IRV tries to identify k-nearest neighbours using a neural network to compute a more complex similarity function.
The majority of these models rely on two key factors: availability of large volumes of training data and knowledge about the physical structure of chemical compounds (which can be used for converting them to graphs or fingerprints). Despite large collections of stock chemical compounds, an offering by many research and commercial entities, the number of compounds in individual assays investigating a specific problem, is relatively small. It may significantly impact the performance of many ML models.
Besides, if a model relies on the physical structure of chemical compounds the training process can become very computationally expensive. This study investigates transfer learning using the variational autoencoder. It reduces reliance on a large volume of training data from the specific assays. 2 This model also does not require any external knowledge about the physical structure of chemical compounds. All knowledge required for making accurate predictions derive from SMILES representation of chemical compounds.

Methodology
The proposed methodology uses a pretrained autoencoder to build classification and regression models for predicting properties of chemical compounds. A high-level overview of the proposed approach is presented in Fig. 2.
The majority of the proposed workflow stays the same for predicting LogD and binding properties. It starts with selecting a collection of chemical compounds from the ChEMBLv23 database [11]. These compounds are used for training variational autoencoder. Then, a specially designed process isolates encoder layers of the variational autoencoder. These layers are combined with an additional sub-network for performing classification or regression tasks (depending on the problem being solved). The constructed regression or classification neural networks are trained using screening data derived from HTS. It is important to mention that the encoder layers remain frozen throughout this training. The obtained model is used for evaluation, and if the assessment is successful, the model can be deployed in a production environment.
The rest of the methodology section is split into three parts. Section 3.1 describes a variational autoencoder. Section 3.2 focuses on classification and regression models for predicting desirable properties of chemical compounds. 2 The variational auto-encoder can be trained to generate chemical compounds fingerprints using a large database of compounds such as ChEMBL. Then, this model can use a small portion of training data (from assays), to predict specific chemical compounds properties of binding characteristics.

Variational autoencoder
An autoencoder is a neural network which can be used to address representational learning problems [30]. It learns to reconstruct the original input using an informational bottleneck. Such neural networks have been successfully applied in various domains such as noise filtering in audiovideo content, translation, and compression. The proposed approach focuses on the autoencoder for reconstruction of chemical compounds. It takes simplified molecular-input line-entry system (SMILES) [52] and tries to reproduce these SMILES using latent space. SMILES are ASCII strings for describing the structure of chemical species. They can be imported by most molecule editors and converted into two-dimensional or three-dimensional structures of the molecules. For example Aspirin C 9 H 8 O 4 is represented by the following SMILES: CC(=O)OC1=CC=CC=C1C(=O)O, whose 2D structure is shown in Fig. 3.
A generic autoencoder is trained to minimize the reconstruction error L defined by Eq. 1: where x is the original input,x is the reconstructed output and R is a regularizer. The regularizer penalizes a model using large weights to prevent memorization and overfitting problems [17]. Ideally, a well-trained autoencoder should accurately reconstruct an input SMILES x such as already mentioned CC(=O)OC1=CC=CC=C1C(=O)O.
In this paper, we are focusing on variational autoencoder [22], a special type of autoencoder, which uses a probability distribution to reconstruct the original input x. Suppose z represents hidden variables (from the latent space) used to reconstruct the original input: The computation of the marginal distribution p(x) is very complex, since the following integral is intractable (in majority cases): There are two main approaches available to tackle this problem: Monte Carlo [38] and variational inference (used to build variational autoencoder) [27]. Let's approximate p(z|x) with another distribution q(z|x), where q can be chosen as a tractable distribution (such as Gaussian) [29]. Then, it is possible to find distribution parameters when q becomes close enough to p by minimizing the Kullback-Leibler divergence [47], which measures an amount of lost information for the chosen approximation: By replacing p(z|x) in Eq. 4 it is possible to derive to the following equation: and express p(x) as: where the second term is called a variational low bound: This allows us to rewrite Eq. 6 as: where p(x) in Eq. 8 can be considered as a constant, since x is given (as the original input). The rest of this equation represents a sum of two quantities, where KL-divergence needs to be minimized. The minimization of KL-divergence is effectively a maximization of the variational low bound L defined by Eq. 7. By substituting p(z, x) in Eq. 7 it is possible to derive to the following equation: where p(x|z) is the expectation with respect to q(z) and can be written as E qðzÞ log pðxjzÞ. The second term ÀKLðqðzjxÞjjpðzjxÞÞ represents the KL-divergence. This allows us to rewrite Eq. 9 as following: Let's build the variational autoencoder based on variational low bound Eq. 9. The observed distribution q is a function mapping x to z, which should match an another distribution p. The observed distribution p is a function mapping z tox, where p can be chosen. Both q and p are implemented as neural networks and for the further references are called encoder and decoder accordingly. A visualization of the variational autoencoder model is shown in Fig. 4. Let's select p to be the Gaussian distribution. This requires us to make the distribution q (in the latent layer) also similar to Gaussian. The cost function can then be expressed as: where N is a normal distribution defined by two parameters l-mean and r-variance. jx Àxj 2 in 11 has been derived from the definition of the reconstruction error for the Gaussian distribution pðxjxÞ ¼ e ÀjxÀxj 2 . As shown in Fig. 5, the encoder learns to represent the original input x as a set of attributes z in the latent space, where each attribute is defined as the probability distribution (with parameters l and r). The decoder learns to reconstructx close to the original input x using a set of attributes z from the latent space.
A set of attributes z in the latent space represents chemical compounds fingerprints and can be used for building classification and regression models. A typical approach for building classification and regression models is to join together trained encoder and problem-dependent prediction layers as shown in Fig. 6.
The training process in such architecture learns some function f(y|z) which predicts y (category value for a classification problem or real value for a regression problem) using chemical compounds fingerprints z generated by encoder. In the classical approach the encoder neural network q(z|x) only generates z and does not take part in the training process (so it weights remain frozen thorough out all training cycle). 3

Architecture
The variational autoencoder was implemented using a Convolution Neural Network [6] in combination with a few layers for supporting the variational training process. Its neural network topology is presented in Fig. 7.
It consists of two joined neural networks: encoder and decoder. The encoder neural network consists of nine layers. The Input layer takes SMILES transformed to the one-hot 150x78 matrix representation (each row represents a SMILES character and column its encoding). A visualization of SMILES transformation is shown in Fig. 8.
This figure demonstrates a process of one-hot encoding for the Aspirin SMILES, introduced earlier. If the input SMILES is less than 150 characters, the original sequence is padded with empty spaces on the right. When the input sequence is adjusted, the transformation process creates a zero-matrix with 150 rows (equals to the maximum number of characters in the input sequence) and 78 column (equals to SMILES vocabulary size). A transformation loop selects each character in the padded sequence and use a dictionary to identify the character code. A position of the selected    character and its code are used to set the according element to 1 in the one-hot matrix. 4 The input layer is followed by three 1D-convolution layers (Conv1D-1, Conv1D-2 and Conv1D-3) with 512, 256, 128 filters and 7, 5, 3 kernel sizes accordingly. These convolution layers perform a very similar role to 2D-convolution layers in image processing. They identify specific patterns of elements (atom and bonds) in chemical compounds and aggregate them in bigger substructures at each consequent layer. A useful insight of how these layers work is presented in the discussion Sect. 6.
The Flatten layer vectorizes convolution weights, so they can be processed by the following Dense-1 layer. This layer has 1024 neurons, which is exactly the same size as an output latent vector z. The size of this layer has been defined via a hyper-parameters tuning procedure. Dense-2 and Dense-3 layers implement variational learning, which computes l and r accordingly. The final Lambda layer combines l and r into a single latent vector (consisting of 1024 real values). Effectively this latent vector represents a chemical compound fingerprint.
The decoder consists of five layers. The input layer receives the latent vector from the encoder and passes it via two dense layers: Dense-4 and Dense-5. The Dense-5 scale-up original dimensionality from 1024 to 11700. This step is needed to transform a 1D into a 2D vector which is performed by the Reshape layer. This layer passes a 2D vector straight to the output. The final output of all these layers is 150x78 matrix, which can be reversed to the SMILES sequence.    4 The reverse process is applied to reconstruct the original sequence from the one-hot matrix representation.
A neural network topology of classification and regression models constructed based on pretrained encoder is shown in Fig. 9.
The encoder in classification and regression models has exactly the same topology structure as already described in variational autoencoder. However, the attached layers are designed to perform classification or regression tasks. The Dense-1 layer receives chemical compounds fingerprints and passes it to the next layer Dense-2. There is no difference in the current neural network topology whether it's applied for a classification or regression problem. The only difference is in an inactivation for the Dense-2 layer. 5 In case of the regression problem the activation function as relu and in case of the classification problem the activation function as sigmoid.
Due to the high level of complexity defined above methodology and architecture, it is helpful to review a simple example, which illustrates a prediction process.

Prediction example
An example of neural network classifier built based on an encoder is shown in Fig. 10.
This is a high level of visualization intends to demonstrate some key concepts described above. Let us assume we have three types of compounds: red, green and blue. The red and green compounds have a very similar rectangular shape (despite the green compound has slightly rounded edges). The blue compound has the triangle shape. A constructed neural network should identify two classes: rectangle or triangle for the input compound.
According to the variational autoencoder methodology the trained encoder generates a latent vector representation for each input compound. This latent vector is a 'compact' representation 6 of the original input and in chemistry domain can be also refereed as a chemical compound fingerprint. It is likely that similar inputs will have a similar latent vector representation. It also means that similar compounds should be closely located in the latent space. The latent space is an abstract concept, which can be very useful to visualize a distributing of encoded compounds.
As it can be seen from Fig. 10 red and green crosses represent location of 'rectangular' instances in close approximation from each other. The blue cross represents a 'triangular' sample accordingly. The dashed line represents areas of distribution of compounds with similar shapes in the latent space.
MLP neural network can provide an efficient architecture to learn a distribution of compounds in the latent space. Since the latent vector has much smaller size in comparison to the original input, MLP does not require a complex topology. One or two hidden layers can be sufficient to handle prediction of properties of chemical compounds for the majority of cases.
The next section describes a series of experiments for evaluating the proposed solutions.

Data and experiment design
In order to evaluate the proposed methodology three studies are presented in this paper:

Experimental data and setup for SMILES reconstruction problem
ChEMBL is a chemical database of bio-active molecules with drug-like properties [16]. It is maintained by the European Bioinformatics Institute (EBI) of the European Molecular Biology Laboratory (EMBL), located at the Wellcome Trust Genome Campus in Hinxton, UK. ChEMBL data are widely used by pharmaceutical companies and research organizations around the World for creating screening libraries in drug discovery. ChEMBLv23 (version 23) has been selected for the current study. It includes approximately 1.7M chemical compounds. The initial study [21] already showed an accurate SMILES reconstruction using a variational autoencoder neural network. This work is focusing on an optimization of a training process. To train a neural network based on all chemical compounds containing in ChEMBL already requires a powerful architecture. However, to process collections such as ZINC [28] or Enamine [50] with hundreds of million chemical compounds requires a much more sophisticated approach. 5 In this publication we are referring to the binary-classification task. This topology requires a slight modification for multi-class classification problems. 6 A latent vector is a point representing the original input in the latent space. The latent space is a collection of vectors, generated by a complex compression function (encoder). It is extremely difficult (and more likely impossible) to demonstrate this on a practical example. This figure is an attempt to demonstrate the variational autoencoder concept in relation to chemical compounds. By introducing this figure we are trying to help our reader to understand how variational autoencoder works within the scope of this paper. If the reader wants to take a deeper dive into variational autoencoder theory, we recommend the following publication [33].
This experiment tries to define, what is an optimal size of a training data set for an accurate reconstruction of SMILES. This will help to scale down training data set and to preserve reconstruction accuracy at the same time. A design of proposed experiments is shown in Fig. 11.
All data taking part in the experiment were normalized and filtered using MolVS Open Source software [4]. SMILES exceeding 150 characters were removed. This had no significant impact on overall quality of experimental results, since only a very small percentage of these compounds were discarded. After filtering the data set comprises 1688073 samples. 7 This data set had been randomly split into training and evaluation partitions in proportions of 75% and 25% of samples respectively. The evaluation data (422,019 samples) remained unchanged throughout all experiments. It helped to score all produced models against the same benchmark. The size of training data varied from 10% to 100% of the original size (1,266,054 samples) depending on an experiment configuration. Generators were developed to feed data to a training model during a fitting processes. These generators streamed data directly from files using fixed size batches (equals to 1024), preventing any memory overflow.
Ten groups of experiments were carried out. Five tests were performed in each group, except for the last one. 8 Each test randomly selected a certain percentage of samples from the internal training partition. Then, these samples were divided into fitting and validation subsets in proportions of 75% to 25% respectively. The validation subset was used in each training epoch to assess model performance. This assessment was necessary to control learning rate, checkpoints and earlier stopping mechanisms. An approximate number of chemical compounds selected for each test is shown in Table 1.
The best setup (experiment 5) shown in Table 2 was selected to build a production autoencoder model (which was used in the further tests for building classification and regression models). A explanation for the selected data split (defined by the experiment 5) is provided in Sect. 5.1 and discussed in Sect. 6.   [21] is different, due to an outdated filtering algorithm. 8 The last group handled 100% of data, no random sampling was needed.

Experimental data and setup for regression problem
LogD values for training and evaluation of a regression model were obtained from the ChEMBL database. To prevent influence outliers on the overall model only values in the interval between LogD 2 ½À20; þ20 were considered. Chemical compounds with normalization issues were also excluded during the pre-processing step. The total number of chemical compounds taking part in testing was 1669058.
A data set containing these values was obtained from ChEMBL using the following steps: 1. Data fields smi and val were retrieved from ChEMBL using the SQL statement, where 'smi' is a chemical compound SMILES and 'val' a LogD value. 2. MolVS Open Source software [4] was used to normalize each SMILES. 3. Records with lengthðsmiÞ [ 150 were removed; 4. Records with LogD values outside the specified interval val 2 ½À20; þ20 were removed; 5. All LogD values were normalized between 0 and 1. 10-fold cross-validation was conducted to assess the model performance. The design is presented in Fig. 12.
The data selected for cross-validation was split into 10 folds where nine folds (90% of data) were taken for training and one fold (10% of data) for evaluation accordingly. The data set allocated to training was also randomly split into fitting/validation partitions in the following proportion 75%/25%. The validation partition was used at each training epoch for assessing model quality. R2-score [10] metrics were recorded for each fold, and later generalized into the final result.
An additional 10-cross validation test was introduced to evaluate the regression model based on the experimental LogD data. Data were downloaded from latest ChEMBL version (on 30 Sept 2018) and pre-processed by industry experts in drug-discovery.
• The data was curated to remove results obtained with solvents other than octanol/aqueous buffer. 9 • Results derived from HPLC retention times were also removed. 10 • Results obtained from experiments conducted at pH other than pH7.4 were also removed (both high and low pH. 11 • Duplicates were removed by using of InChiKeys. 12  9 The vast majority of the historical data available was determined using octanol/aqueous buffer as the liquid phases, there is a body of data using cyclohexane/aqueous buffer but this was excluded to reduce the number of variables involved in the experiments. 10 Traditionally partition coefficients are determined experimentally using the Shake Flask Method [45]. Since this test can be timeconsuming, alternative methods have been explored; HPLC retention times [34], however, applying a regression equation derived from one chemical class to a second one may not be reliable. For internal consistency only data from experiments using the shake flask method were used. 11 The majority of the experiments were carried out at physiological pH (7.4) there is small amount of data from different pH but since this will affect the extent of ionization of the molecules alternative pH data was removed. 12 Some molecules had been tested many times, to give equal weight to every molecule a single record for each molecule was required. InChiKeys were generated for each molecule [3] and used to compare records, duplicate structures were removed.
After cleaning and pre-processing this data set consisted of 12413 samples which LogD property defined in the interval À12:0; þ12:0.
The second data set 'Lipophilicity' was obtained from the ML resource described in [55]. It consists of 4200 chemical compounds which LogD property is defined in the following interval ½À1:5; þ4:5. This data set was used to gauge the performance of developed system against other ML models.

Experimental data and setup for classification problem
ChEMBL 13 contains approximately 13.5M bio-activity measurements, where 1.1M assays are assigned to approximately 11K targets. The majority of available bioactivity data are highly unbalanced. More than 50% of assays have just a single measurement while others contain tens of thousands. On the other hand, a lot of targets belong only to a single assay, while others to hundreds. A large proportion of these data contain duplicate records. Such heterogeneity of data prevents clear identification, which measurements can be considered as active or inactive accordingly. A special protocol proposed in [40] helps to generate benchmark data sets for binary classification. It has the following six steps: A measurement was considered as inactive if com 2 I. Sets (A and I) of strings in comment-field are defined below. If a record is identified as inactive all further steps are discarded. A = ('active', 'note: corresponding ic50 reported as active') I =('inconclusive', 'not active', 'inactive', 'not active (inhibition \ 50% @ 10 um and thus dose-response curve)') 5. Removed all records passed the previous step where val ¼ ; j unt 6 ¼ 0 nM 0 j rel 6 2 f 0 [ 0 ; 0 ! 0 ; 0 \ 0 ; 0 ; 0 ¼ 0 ; 0 $ 0 g. 6. Assigned labels to each record according to the defined thresholds presented in Sect. 12: 7. All records with duplicates and contradictory measurements obtained during the previous steps were discarded.
56 assays were identified for this study. Each assay reflects in vitro measurements obtained during HTS. A break down between active or inactive chemical compounds together with associated target are presented in Tables 3 and 4. Three cross-validation experiments were carried out for each data set. A design of these experiment is shown in Fig. 12. This folds number was chosen due to the relatively small assay sizes. As it can be seen from Table 5, a large proportion of assays have approximately 800 compounds. Three folds provide a fair representation of active or inactive chemical compounds across training, validation and evaluations partitions.
The Maximum Unbiased Validation (MUV) [48] is another benchmark data set selected from PubChem BioAssay by applying a refined nearest neighbour analysis. The MUV data set contains 17 challenging tasks for around 90 thousand compounds and is specifically designed for validation of virtual screening techniques. The detail breakdown between active and inactive compounds is shown in Table 5.  For comparing prediction models the Receiver Operating Characteristics -Area Under The Curve (ROC-AUC) [12] was used. This metric is widely excepted by a ML community to assess performance of classification models.

Results
According to the experimental setup defined in Sect. 4, results are presented in three sub-sections. The first Sect. 5.1 shows results for variational autoencoder, the second Sect. 5.2 for prediction of LogD (regression problem) and the final Sect. 5.3 for prediction of compoundstargets binding (binary classification problem). The main rationale behind these experiments is to validate the versatility of latent vector based fingerprint, in other words, to prove that it works well regardless of the selected problem (classifier or regression).

SMILES reconstruction
Forty-six tests (9x5?1 for more details see 4.1) were conducted to assess the accuracy of SMILES reconstruction on different portion of training data. The results obtained are detailed in Table 6.
Changes in accuracy and editing distance for different size of training data sets are presented in Figs. 13 and 14. As it can be seen from Fig. 13, the reconstruction accuracy increases from 0.247 AE 0.027 for 10% of randomly selected samples to 0.877 AE 0.009 for 50% of samples accordingly. From 60% onwards accuracy stays around 0.8 on average with slight fluctuations. This is an expected result. However slight variations in accuracy starting from 60% of samples needs to be addressed. It is a difficult task to identify the exact reason of what is influencing these According to the carried out experiment the best result (accuracy 0.877 AE 0.009) was obtained for 50% of randomly selected samples. This configuration was selected for building a production variational autoencoder, which can be later used for training classification and regression models with SMILES input. Training was conducted using 40 epochs, and produced a model with accuracy 0.872. The editing distances for the production model were 0.682 and 0.677 for Hamming and Levenshtein respectively.

LogD prediction
Two cross validation experiments described in Sect. 4.2 were carried out on ChEMBL and Lipophilicity data sets. Results of these experiments are presented in Table 7.
As it can be seen from the obtained results the best performance is achieved for the ChEMBL data set, with an average coefficient of determination of 0.907 AE 0.008. A scatter plot with an alignment of true and predicted values is shown in Fig. 15. The training process was carried out with LogD values normalized in the interval ½À20; þ20 14 and all training cycles took 20 epochs. The high prediction accuracy is a very much expected result, since the majority of LogD data points in ChEMBL are computed using the ACD/Labs software [1]. In this scenario the regression model is simply learning to predict an outcome of another computational algorithm (such as ACD/Labs). It makes a learning task much more straightforward. This assumption is vindicated by results obtained in other studies. For example experiments with ChEMBL data set described in [1] also show a high accuracy using a SVM [19] model.
Much more interesting results are obtained for the Lipophilicity data set. The average R2 score equals to 0.542 AE 0.021, which is noticeably less in comparison to ChEMBL data. The training process was carried out with the LogD normalized interval ½À1:5; þ4:5 and all training cycles took 30 epochs. A longer training cycle reflects the complexity of building a predictive model on real experimental data. A scatter plot with an alignment of true and predicted values for the Lipophilicity data set is shown in Fig. 16.
Eight ML models were investigated on this data. A comparison chart for all these models is presented in Fig. 17. The highest score 0.697 is obtained for MPNN, which uses a generalized model [23]. It is very suitable for processing graph structured data, which makes it efficient in predicting properties of chemical compounds (since chemical compounds can be easily represented as an undirected graph). MPNN is closely followed by GC model [31], with R2 score equals to 0.662. GC utilities principles of circular fingerprints described in [35] representing molecular structures by atom neighbourhoods. Similar to GC, Weave [31] (0.636) is another graph-based model which processes chemical compounds as a undirected graph using a convolution approach. In contrast to the previous three models, XGB [14] provides a different approach for making predictions. It is an ensemble approach which combines predictions of individual decision trees. XGB coefficient of determination for Lipophilicity equals to 0.577 and is closely followed by the model developed with in this study (0.542 AE 0.021). Directed Acyclic Graph (DAG), Kernel Ridge Regression (KRR) and Random Forests (RF) are the lowest performing in this evaluation with R2 scores 0.507, 0.496 and 0.483 respectively.
The key challenge of developing ML models for predicting the properties of chemical compounds is to encode molecules into fixed-length strings or vectors representation [54]. Despite SMILES providing unique representations of molecules, the majority of ML models are also relaid on additional information such as electronic or topological, profiles of chemical compounds. To derive these features, the models, we are gauging against, applied different factorizations: Extended Connectivity Fingerprints (ECFP), Coulomb matrix, Grid features, etc. These approaches are computationally expensive, and there will always be a trade-off between speed, accuracy and expense. Because of this, it is not surprising that some models provide better performance compared to our approach. As has been already mentioned, our approach is purely data-driven and inspired by -et al. [25]. It should provide alternatives to replace crafted featurization methods with the learning ability of DNN. In the future development, we are planning to improve the encoding method and directly compare it to existing featurization approaches.
Open-source research in AI always provides a solid benchmark for assessing in-house models. However, it would be interesting to compare the developed model against a commercial application. The following results show comparison of predictions made by ChemAxon software [5] against the developed regression model. This  14 The deviation of predicting LogD values base on the interval ½À20; þ20 may be seen to be high. However, predicting LogD is a particularly challenging problem since we have to predict lipophilicity (LogP) and the ionization (pKa) of a molecule. The distribution can range considerably but the hope is that by identifying outliers we can identify those structural classes that provide unique challenges to the algorithm. In some cases this will be because the structural class is not well exemplified within the training set, in other cases, it may be that very minor structural changes have very profound effects and the more coarse-grained models don't give sufficient accuracy.
work was carried out in collaboration with Cambridge MedChem Consulting [2]. Data were selected and preprocessed by industry experts in drug discovery. All undertaking steps for preparing this experiment are described in Sect. 4.2. 10-folds cross-validations were performed on selected data set to compare models. 15 The obtained results are shown in Table 8. Results for the developed regression are shown in the 'ARM' column and for the commercial software in the 'ChemAxon' column accordingly.
An average R2 score for our developed regression model equals to 0.695 AE 0.013, which is nearly twice that obtained by ChemAxon with R2 of 0.338 AE 0.034. It is hard to comment on the underlining ChemAxon algorithm for making prediction without available source code. However, from the description presented on the company website [5] it possible to assume, that it is a deterministic algorithm which approximates a chemical compound structure into a specific property value. Two scatter plots represented in Figs. 18 and 19 show an alignment of true and predicted values based on our developed regression model and ChemAxon respectively 16 The experiments clearly demonstrate the validity of our proposed model to predict the LogD property of chemical compounds. However a large proportion of tasks required simple classification, for example whether a compound binds to the specified target. An evaluation of the classification model constructed based on variational autoencoder is presented in the next Sect. 5.3.

Binding prediction
Two cross validation experiments described in Sect. 4.3 were carried out on 56 ChEMBL data sets. Results of these experiments are presented in Table 9.
Considering the large volume of obtained results, they were split into five groups based on ROC-AUC metric (see Fig. 20). The first group combines 10.7% of assays with Despite competent results obtained on vast majority of tested assays, the authors carried out additional investigation to rank the developed classifier against other ML algorithms. Similar to the regression problem, the main objective here is not a direct comparison of different ML algorithms, since it requires different experimental setup. This study projects prediction accuracy observed for the developed model on results already described by -et al. [55]. The ROC-AUC characteristics for 17 binding-assays are presented in Table 10.
An average ROC-AUC of a 10-folds cross-validation experiment was recorded for each assay. The summary field represents an average ROC-AUC across all 17 experiments. It was used for ranking the developed classifier against 6 ML algorithms. A visual representation of this ranking is shown in Fig. 21.
Similar to the regression problem, GC -a graph based model, scored 0.775 the best result for binding classification. 17 It closely followed by BYPASS, representing a multitask neural network and LOGREG representing logistic regression model, scoring 0.764 and 0.749

Discussion and conclusion
During the last few years deep neural networks de facto have become an industry standard for creating sophisticated AI models. A significant amount of effort have been devoted to improve classical ML algorithms. Often XBG and RF produce even better results than DNNs. Such a variety of approaches makes it a very difficult task to choose the right technique. This has also made a huge impact on scientific publications. Researchers are forced to show rigorous testing, with comparison of every proposed technique against the recognized leaders. It expected that the developed technique must outperform others, which in practice leads to compromising with experiments design and cherry-picking phenomenon.
It can be seen from our comparison of published models, the majority show very similar performance. However, it is also important to score each model in term of practical application. For example a model can deliver impressive accuracy, but when it needs to be deployed in a production environment, scalability and efficiency diminish all advantages gained in perfecting the quality of predictions.
This work has been inspired by a research effort of using variational autoencoder to generate chemical compounds fingerprints, using these fingerprints for predicting specific properties of chemical compounds. Considering the high complexity of chemical compounds to train quality variational autoencoder requires a large data set of SMILES. A typical size of HTS assays consists of several hundred, maybe thousands of chemical compounds. Such volume of data samples does not deliver a sufficient variety of chemical compounds structures, which makes it impractical to train a variational autoencoder based on assay data.
A decision was made to use ChEMBL data to obtain a large representation of chemical compounds structures. It worked very well, with the developed variational autoencoder model reconstructing nearly 90% of SMILES using 1024 latent vector. Taking into the account Hamming and Levenshtein editing distances the final model in average has one misplaced or incorrect atom or bond. Obviously such error is unforgivable in chemistry, but the primary objective of variational autoencoder is to produce latent space where similar structures are crumbled together. Let's demonstrate this on simple example.
The search space for similar chemical compounds was reduced to 10,000 samples (out of 674,040 initially allocated for evaluation, for more detail see Sect. 4.1) to make computation more efficient. K-nearest neighbour [36] identified five closest chemical compounds defined in Table 11, and illustrated in Fig. 23.
The last column in Table 11 represents Tanimoto similarity score [9]. It is widely used in the chemistry domain to assess similarity between two chemical compounds. Two compounds can be considered similar if Tanimoto score is greater than 0.85 (for Daylight fingerprints). As it can be seen from the obtained results in Table 11 three compounds retrieved using latent space are also similar This example demonstrates that latent vector based fingerprints can be used to define similarity between two chemical compounds. It also clearly shows that selected chemical compounds are closely located in latent space.
The closest approximation of similar compounds in the latent space provides potential capability for the developed model to be applied to generating new chemical compounds and forecasting the desired properties. Such ML models become a hot topic in pharmaceutical domain, which can be witnessed by increasing the number of highquality research publications in this space [25,46]. In this paper, the main focus was on the transfer learning, to use the trained encoder as a base for classification or regression networks which can predict properties of chemical compounds. However, the decoder can be potentially used for generating novel chemical compounds. By introducing a small modification into the latent representation of a target chemical compound, it is possible to generate a novel structure. Despite this simple idea the implementation of such a model is very complex and outside the scope of this publication.
The trained variational autoencoder forms a solid base for creating different classification and regression models. An interesting pattern was observed during a training process. The majority of trained models converge to a stale state (where no longer improvement observed) during 2-3 epochs. However when the same topology of neural network was trained without pretrained variational autoencoder, the training process continue up to 100 epochs. Longer training is not a problem for relatively small data set, where full learning cycle can be complied in the matter of hours. However in case of such collection volume as ChEMBL, training may go on for days. Also, the experiments did not reveal any degradation in accuracy with shortening the training cycle.
The question is why are classification and regression models built based on variational autoencoder so efficient in the training process? To answer this question, let us come back to the methodology described in Sect. 3. Constructed classification and regression models consists from two parts. The first part is an encoder isolated from a trained variational autoencoder. The second part is MLP, which performs actual predictions. Effectively, the MLP is trained to make predictions based on chemical compound fingerprints. Since the 'hard work' has been already done by variational autoencoder, only a few cycles are required to learn differentiation rules (to solve classification or regression problem).
A further observation which was called 'latent space drift' was also noted. A number of publications using a similar approach do not clearly reveal their mechanisms of using an encoder. It can be used with frozen and unfrozen layers. Using an encoder with frozen layers makes a lot of   Table 11 sense. If it is already trained to encode SMILES, then the attached layers can only be trained to utilize the obtained chemical compounds fingerprints. However, a preliminary study showed that if an encoder is left unfrozen, the training result is generally better. An initial investigation found that fingerprint points in latent space become adjusted according to the target (predicting) property. This process is explained in Fig. 24. All this internal analysis of neural network behaviour becomes possible due the specialized AUROMIND software. It provides a set of tools for 'debugging' a training process. The authors plan to describe some of the core principles behind developed tools in upcoming publications.