Improving Chemical Autoencoder Latent Space and Molecular De Novo Generation Diversity with Heteroencoders

Chemical autoencoders are attractive models as they combine chemical space navigation with possibilities for de novo molecule generation in areas of interest. This enables them to produce focused chemical libraries around a single lead compound for employment early in a drug discovery project. Here, it is shown that the choice of chemical representation, such as strings from the simplified molecular-input line-entry system (SMILES), has a large influence on the properties of the latent space. It is further explored to what extent translating between different chemical representations influences the latent space similarity to the SMILES strings or circular fingerprints. By employing SMILES enumeration for either the encoder or decoder, it is found that the decoder has the largest influence on the properties of the latent space. Training a sequence to sequence heteroencoder based on recurrent neural networks (RNNs) with long short-term memory cells (LSTM) to predict different enumerated SMILES strings from the same canonical SMILES string gives the largest similarity between latent space distance and molecular similarity measured as circular fingerprints similarity. Using the output from the code layer in quantitative structure activity relationship (QSAR) of five molecular datasets shows that heteroencoder derived vectors markedly outperforms autoencoder derived vectors as well as models built using ECFP4 fingerprints, underlining the increased chemical relevance of the latent space. However, the use of enumeration during training of the decoder leads to a marked increase in the rate of decoding to different molecules than encoded, a tendency that can be counteracted with more complex network architectures.


Introduction
Autoencoders have emerged as deep learning solutions to turn molecules into latent vector representations as well as decode and sample areas of the latent vector space [1][2][3].An autoencoders consists of an encoder which compresses and changes the input information into a bottleneck layer and a decoder part which recreates the original input from the compressed vector representation (the latent space vector).After training, the autoencoder can be reassembled into the encoder which can be used to calculate vector representations of the molecules.These can be used as a sort of molecular fingerprints or GPS for the chemical space of the molecules.The decoder can be used to translate back from the latent representation to the molecular representation used during training, such as simplified molecular-input line-entry system (SMILES).This makes it possible to use the decoder as a steered solution for molecular de-novo generation, as the probability outputs of the decoder can be sampled, creating molecules which are novel but close to the point in latent space.Alternatively, the molecules of the nearby latent space can be explored by randomly permuting the latent vector.
Various encoder-decoder architectures have been proposed as well as the latent space has been reg- ularized and manipulated using variational autoencoders [1] and adversarial autoencoders [2].Both convolutional neural networks (CNNs), as well as recurrent neural networks (RNNs) have been used for the encoder part [1][2][3], whereas the decoder part has mostly been based on RNNs with either gated recurrent units (GRU) [4] or long short-term memory cells (LSTM) [5] to enable longer range sequence memory.
A famous painting of René Magritte, "The Treachery of Images", shows a pipe, and also has the text "Ceci n'est pas une pipe": This is not a pipe.The sentence is true as it is a painting of pipe, not the pipe itself, kindly reminding us that representation is not reality.Autoencoders based on SMILES strings [6] face the same fundamental issue.Is the latent space a representation of the molecules or is it a condensed representation of the SMILES strings representing the molecules?
Due to serialization of the molecular graph, a single molecule has multiple possible SMILES strings, which has been exploited as data augmentation with the SMILES enumeration technique [7].A simple challenge with different SMILES representations of the same molecule suggest that the latent space is more related to the SMILES string than to the molecule, which has also previously been noted [2]. Figure 1 shows the same three molecules after projecting into the latent space of an RNN to RNN autoencoder.The molecules end up being projected to very different areas of the latent space, although some clustering can be observed.The latent space is thus likely a mixture of SMILES representation information and chemical information representation.One way of solving this challenge could be to use special engineered networks and graph based approaches [8] for molecular generation.
As an alternative to engineering the outcome, it is here suggested that it is possible to use SMILES enumeration or chemception image embedding [9] to create chemical heteroencoders.The concept is illustrated in Figure 2. By translating from one format or representation of the molecule to the other, the encoder-decoder network is forced to identify the latent information behind both representations.This should in principle lead to a more chemically relevant latent space, independent of the formats or canonicalization used.
Here, the choice of representation and enumeration is explored for both training the encoder or decoder and the latent space similarity to SMILES and scaffold based metrics calculated.Moreover, it is tested if these changes influence the properties of the decoder when used for de-novo design of molecules.Further, an optimized and expanded heteroencoder architectures trained on ChEMBL23 datasets are used to extract latent vectors for subsequent use as input to QSAR models of five different molecular datasets.

GDB-8
The GDB-8 dataset [10,11] was downloaded and split randomly into a train and test set using a 0.9 to 0.1 ratio.

ChEMBL23
Structures were extracted from the ChEMBL23 database [12] and validated using in-house rules at Science Data Software LLC (salts were stripped, solvents removed, charges neutralized and stereo information removed).Maximum available length of the canonical SMILES string allowed for a molecule was 100 characters.10 thousand molecules were selected randomly for the held out test set.From the remainder of the 1.2 million molecules was randomly selected a training set of 400 thousand molecules and a validation set of 300 thousand molecules for used during training procedures.

QSAR datasets
Five experimental datasets were used, spanning physico-chemical properties as well as bioactivity.Four datasets (IGC50, BCF, MP, LD50) for QSAR modeling were downloaded from the EPA Toxicity Estimation Software Tool [13] webpage [14] and used as is without any additional standardization.The solubility was obtained from the supplementary information of [15].Information of the datasets are shown in Table 1.
Datasets from EPA's TEST suite were already split into train/test sets in a 75/25% ratio and were used accordingly.Molecules for the solubility dataset were obtained by resolving CAS numbers from the supporting info [15] and the dataset was randomly split using the same ratio as the other QSAR datasets.

1D and 2D Vectorization
SMILES were enumerated and vectorized with onehot encoding as previously described [7].
2D vectorization was done similar to the vectorization used in Chemception networks [9] with the following modifications: A PCA with three principal components was calculated on atomic properties from the mendeleev python package [19] (dipole_polarizability, electron_affinity, electronegativity, vdw_radius, atomic_volume, softness and hardness).The PCA scores were normalized with min-max scaling to be between zero and one to create the atom type encoding.PCA and scaling were performed with the Scikit-Learn python package [20].RDKit [21] was used to compute 2D coordinates and extract information about atom type and bond order.The normalized PCA scores of the atom types were used to encode the first three layers and the bond order was used to encode the forth layer.A fifth layer was used to encode the RDKit aromaticity perception.The 2D coordinates of the RDKit molecule were rotated randomly up to +/-180°around the center of coordinates before discretization into numpy [22] floating point arrays.
The first layer consisted of 64 LSTM cells [5] used in batch mode.The final internal memory (C) and hidden (H) states were concatenated and used as input to a dense layer (the bottleneck) of 64 neurons with the rectified linear unit activation function (ReLU) [25].Two separate dense layers with ReLU activation func-tions were used to decode the bottleneck outputs into the initial C and H states for the RNN based decoder.
The decoder consisted of a single layer of 64 LSTM cells trained with teacher forcing [26] in batch mode.The output from the LSTM cells was connected to a Dense layer with a softmax activation function matching the dimensions of the character set.A two-layer model was also constructed by increasing the number of LSTM cells to 128 and the number of LSTM layers to two in both the encoder and decoder.Accordingly, four separate dense networks were used to decode the bottleneck layer into the initial C and H states for the two LSTM layers in the decoder.
The networks were trained with mini-batches of 256 sequences for 300 epochs using the categorical cross entropy loss function and the Adam optimizer [27] with an initial learning rate of 0.05.The two layer model was trained with an initial learning rate of 0.01.The loss was followed on the test set and the learning rate lowered by a factor of two when no improvement in the test set loss had been observed for 5 epochs.
After training in batch mode, three models were created from the parts of the full model.A decoder model from the initial input to the output of the bottleneck layer.A model to calculate the initial states of the LSTM cells in the decoder, given the output of the bottleneck.Lastly, a stateful decoder model was constructed by creating a model with the exact same architecture as the decoder in the full model, except the LSTM cells were used in stateful mode and the input vector reduced to a size of one in the sequence dimension.After creation of the stateful model, the weights for the networks were copied from corresponding parts of the trained full model.
The image to sequence model CNN encoder was built from three different Inception-like modules [28].The first module consisted of a tower with a 1x1 2D convolutional layer (Conv2D) followed by a 3x3 Conv2D, a tower with a 1x1 Conv2D layer followed by a 5x5 Conv2D layer and a tower with just a single 1x1 Conv2D layer.The outputs from the towers were concatenated and sent to the next module.
The standard inception module was constructed with a tower of 1x1 Conv2D layer followed by a 3x3 Conv2D layer, a tower with a 1x1 Conv2D layer followed by a 5x5 Conv2D layer, an extra tower of a 1x1 Conv2D layer followed by a 7x7 Conv2D but with only half the number of kernels and a tower with a 3x3 Maxpooling layer followed by a 1x1 Conv2D layer.All strides were 1x1.The outputs from the four towers were concatenated and sent to the next module.
The inception reduction modules were similar to the standard module, except they had no 7x7 tower and used a stride of 2x2.
A standard inception module was stacked with a reduction inception module three times, giving 7 inception modules in total including the initial one.The number of kernels was set to 32.
The outputs from the last inception module were flattened and followed by a dropout layer with a dropout rate of 0.2.Lastly the output was connected to the bottleneck consisting of a dense layer with the ReLU activation function.The decoder part was constructed similar to the sequence to sequence models described above with one layer LSTM cells.The image to sequence model was trained similar to the sequence to sequence models for 200 epochs.
The models are named after the training data in a encoder2decoder naming scheme."Can" is training data with canonical SMILES, where "Enum" designates that the input or output was enumerated during the training."Img" shows that the data was the 2D image embedding.

Similarity metrics
SMILES strings sequence similarities were calculated as the alignment score reported by the pairwise global alignment algorithm of the Biopython package [29].The match score was set to 1, the mismatch to -1, the gap opening to -0.5 and the gap extension to -0.05.The fingerprint similarity metric was calculated on basis of circular Morgan fingerprints with a radius of 2 as implemented in the RDKit library [21].The fingerprints were hashed to 2048 bits and the similarity calculated with the RDKit packages FingerprintSimilarity function.The latent space similarity between two molecules was calculated as the negative logarithm to the Euclidean distance of the vector coordinates.

Enumeration Challenge
The encoder was used to calculate the latent space of the test set, followed by a dimensionality reduction with standard principal components analysis (PCA) as implemented in the Scikit-Learn package [20].Three molecules were converted to different SMILES strings with SMILES enumeration [7].The latent space coordinates of the non-canonical SMILES were calculated with the encoder and transformed and projected onto the visualization of the principal components from the PCA analysis.

Error analysis of output
The percentage of invalid SMILES was quantified as the number of produced SMILES which could not be validated as molecules by RDKit.Subsequently the equality of the input and output RDKit molecules was checked.The similarity of the scaffold was checked by comparing the generalized murcko scaffolds [30] including side-chains.The atom equivalence was checked by comparing the molecular sum formulas.The number and nature of bonds was compared via a "bond sum formula" by counting the number of single, double, triple and aromatic bonds.

Multinomial sampling of decoder
Multinomial sampling was implemented as previously described [31].The sampling temperature was kept at 1.0.

Neural Network Modelling for the ChEMBL dataset
The sequence to sequence autoencoder used for encoding the ChEMBL data and extraction of vectors for QSAR modelling was programmed in Python 3.6.1 [32] using Keras version 2.1.5[23] with the tensorflow backend [24].The encoder consisted of two bidirectional layers of 128 CuDNNLSTM cells in each oneway layer.The final C and H states were concatenated and passed as input to a dense layer with 256 neurons using the ReLU [25] activation function (the bottleneck).The output from the dense layer were decoded by four parallel dense layers with the ReLU activation function, whose outputs were used to set the initial C and H states of the decoder LSTM layers.The decoder itself consisted of two unidirectional layers of 256 CuDNNLSTM cells each .The decoder was trained under teacher forcing as described for the simpler networks above.Every non-linear activation was followed by Batch Normalization.No additional regularization was used.400k random structures from the CheMBL23 training set were pre-enumerated 50times for each SMILES string.The new 20M pairs were shuffled and used in both a canonical to enumerated and an enumerated to canonical setting and trained until model convergence.The same 400k canonical SMILES were also used to train an auto encoder from canonical to canonical SMILES.For the enumerated to enumerated training setting 50 pairs (when possible) of different SMILES strings were created for each molecule of the training set.The network was trained using mini-batches of 256 one-hot encoded SMILES strings, using the Adam optimizer with an initial learning rate of 0.005.The training was monitored and controlled by three callbacks.One callback monitored the loss of the validation set and lowered the learning rate by a factor two when no improvement had been observed for 2 epochs (Re-duceLROnPlateu).Another Callback stopped training when no improvement in the validation set loss had been observed for 5 epochs (EarlyStopping), and the last callback saved the model if the validation loss has improved (CheckPoint).Models typically converged after approximately 40 epochs, which usually took about 6 hours on a NVIDIA GTX 1080 Ti equipped server.

QSAR modelling
Subsequent QSAR modelling was performed using the machine learning capabilities of the Open Science Data Repository [33].An initial search for hyper parameters were performed after converting the molecules into ECPF4 fingerprints (radius 2, 1024 bits).The hyper parameter search for a neural network was performed using Tree of Parzen Estimators (TPE) algorithm [34] as implemented in Hyperopt [35] with the search space bounds listed in Table 2 .The performance on each dataset was optimized using 3-fold cross validation on the training set.The performance of the model with the final hyperparameters were subsequently tested on the held-out test set using an ensemble of 10 models build during 10-fold cross validation during training.The auto-/heteroencoders trained on the ChEMBL23 molecules were subsequently used to encode the QSAR datasets into vectors using the output from the bottleneck layer.The vectors were used as input to the QSAR models.The same hyper parameters were used as identified

GDB-8 dataset based models
All models were trained to full convergence and obtained close train and test losses, the later listed in Table 3 .The models trained on enumerated SMILES output have a markedly larger final loss, but all models show a low degree of malformed SMILES when sampling the latent space vectors calculated from the test set.

Molecular and sequence similarity
Using the same reference molecule, similarity metrics were calculated based on the latent space vectors of the test set, Morgan fingerprints and sequence alignment scores, followed by calculation of the correlation coefficients (R 2 ).Examples of SMILES alignments are shown in Figure 3 for two different alignments.Figure 4 shows an example scatter plot of the alignment scores and latent space similarity for the first molecule and the rest of the molecules in the test set.
The correlation between the same latent space similarity measurement and the Morgan fingerprint similarity was intended as a metric of the scaffold similarity independent of the SMILES strings, and an ex-   ample scatter plot is shown in Figure 5.Both the sequence alignment score and the fingerprint based similarity have correlation with the latent space similarity, which shows that the latent space is at least somehow related to our traditional understanding of similarities between molecules.The properties and the correlations of all the models trained on the GDB-8 dataset are listed in Table 3.The models with a decoder trained on canonical SMILES show a markedly larger correlation between the latent space and the SMILES sequence similarity metric than between the fingerprint based similarity and the latent space.In contrast, the fingerprint and sequence similarities correlations to the latent space similarity are more on the same level when the decoder is trained using enumerated SMILES.The heteroencoder based on the im-  Figure 6 and 7 show a result of similarity searching in the latent space of the test set using a query molecule.The molecules from the can2enum model seems qualitatively more similar than the ones that are most similar in the latent space produced by the can2can model.There is overlap between the two sets, so in some respects the two latent spaces are related.

Error analysis
The models in general produce large percentages of valid SMILES (Table 3).However, using enumeration in the input and output significantly increases the percentage of the outputs where the decoded molecule is not the same as the encoded molecule.The input and output molecules were further compared with regard to scaffold, molecular sum formula and equality of bonds.Figure 8 shows the error types and overlap for the can2enum model.494 molecules out of 1000 tested was valid SMILES but not the same as the input molecule.Further 220 had the wrong scaffold, 42 the wrong sum formula and 14 the wrong bondtypes or number of bonds.The majority (251) had the right scaffold, the right atoms and the correct bonds, but had seemingly assembled the molecule in a wrong order.The bond types and atoms are in principle simple accounting operations independant of the SMILES enumeration, whereas the models struggle more with the scaffold reconstruction and the atom order, which are influenced by the SMILES enumeration.The results for the other heteroencoder models are qualitatively similar (not shown).

Enumeration Challenge
The encoders capabilities to handle different SMILES from the same molecule were tested by projection to a PCA reduction of the latent space (c.f. Figure 1).Figure 9 shows the improvement that can be obtained by training the heteroencoders with enumerated SMILES for either the encoder and decoder.Training the encoder with enumerated SMILES strings gives the tightest clustering (enum2can), showing that the encoder has learned to recognize the same molecule independent on actual serialization of the SMILES string.By showing multiple different SMILES strings to the encoder during training, the encoder can produce the latent space coordinate most suitable for recreating the SMILES form of the decoder, irrespective of the SMILES form shown to the encoder.The enum2enum model has a similar tight clustering as the enum2can model (not shown).The can2enum model also show more tight clustering than the can2can model from Figure 1, indicating that the heteroencoding itself changes the latent space although the encoder itself was not trained on different SMILES forms.Alternatively, the model is doing a more complicated task which could work as regularization leading to better generalization.

Sampling using probability distribution
Figure 10 illustrates the difference between probability sampling of the can2can and can2enum model.The decoder outputs a probability distribution at each step, which can be sampled randomly according to the probabilities (Multinomial sampling).For the can2can model, there is little difference between this sampling strategy and the simple selection of the most probable next character.The can2enum model instead show a lot more uncertainty in the next characters in the beginning of the sampling.The first character is most likely "C", but also "N" and "F" are possibilities.As the model samples "C", the next character is either a "C", a branching "(", or start of a ring "1".Because it then samples "C", it has to choose a ring start next.Towards the end of sampling, the decoder gets completely certain with the last 6 characters, probably because there is only one way to finish the molecule with the already sampled characters.Table 4 shows some statistics on the sampled molecules using the latent coordinates from a single molecule.The model trained on canonical SMILES in both encoder and decoder are very sure about the SMILES it want to recreate, as only one SMILES form and one molecule is sampled.In contrast, the decoders trained with the enumerated SMILES create different SMILES forms of the correct molecule, but also creates other molecules as well.The more complex model with two LSTM layers handles the task a bit better than the single layer version, which only produce the molecule presented for the encoder 20% of the times.Examples of the sam-

enumerated SMILES and projected to the latent space reduced to two dimensions with principal components analysis (PCA). Using enumerated SMILES for training of the encoder leads to the tightest clustering, but also training
with the enumerated SMILES in the decoder improves the clustering (c.f. Figure 1).Small blue dots are the test set used for the PCA reduction.

A B Sampling of Can2Can model Sampling of Can2Enum model
pled molecules from the two layer model are shown in Figure 11.

QSAR modelling using ChEMBL trained heteroencoders
Training of the models used on the ChEMBL datasets, resulted in final losses of approximately 0.001, 0.01, 0.10, 0.11 for the can2can, enum2can, can2enum and enum2enum configurations of the training sets, respectively Reconstruction performance of the different encoder/decoder configurations is presented in the Table 5.The performance of the QSAR modelling done based on ECFP4 fingerprints and the different vectors obtained from the bottleneck layers of the three different neural network models are shown in Table 6.There is some improvement from the ECFP based baseline models to the can2can vector based models, with further improvement for the models based on the latent vectors extracted with the heteroencoders.The three different heteroencoders seem to produce latent vectors which perform very similar to each other in the QSAR modelling, with a tendency for the average performance to rise from enum2can to enum2enum over can2enum.
An interesting observation was that approximately 40% of the neurons of the bottleneck layer for each configuration are never activated.This is likely related to the use of the ReLU activation function and diminishing or enlarging the bottleneck layer resulted in almost the same percentage of inactive neurons (results not shown).The latent vector is thus even denser than the chosen number of neurons.

Discussion
Changing the representations used for training autoencoders (here called heteroencoders) have a marked influence on the properties and organization of the latent space.Although a perfect correlation to the standard fingerprint similarity is not wanted or expected, it is more reassuring that the dependence to the SMILES sequences is at a similar level to the fingerprint based similarity, than the situation where the correlation to the SMILES sequence is much larger than the correlation to the fingerprint metric.The greater balance between the two correlations strongly indicates that the latent space is just as relevant for the molecular scaffold as it is to the SMILES sequence in itself.
The dataset used in the first part of this study was of a limited size and molecular complexity (only 8 atoms).Additionally, as the dataset is fully enumerated the same graph structures are very likely present in both training and test set, which could be the basis of the excellent reconstructions of the test sets.The model could in principle memorize all graph structures instead of learning the rules behind the graph scaffolds, and then simply assign a specific sequence of atoms to the memorized graph.Even though the dataset was somewhat simple, the models trained on enumerated data may have struggled because of low neural network fitting capacity.This indeed seems to be the case as the 2-layer enum2enum model has much lower final loss (Table 3) and also much better reconstruction statistics when reconstructing (Table 3) and sampling the molecules (Table 4) than the single layer enum2enum model.Even though there is a rough correlation between the SMILES validity rate and the molecule reconstruction error rate, with heteroencoders the molecule reconstruction rate becomes a more relevant term to measure than the SMILES validity rate, as the former can diverge a lot, while the SMILES validity error rate is still low.
Models employed in other studies are more complex with larger and multiple layers [1-3, 8, 9, 31].The heteroencoder concept was thus further expanded to also handle ChEMBL datasets.The expansion of the networks to two layers in both the encoder and decoder, use of bidirectional layers in the encoder and a larger number of LSTM cells allowed to fit the larger molecules, although the uncertainty in the reconstruction of the molecule is still present (c.f.Table 5).It is likely that even more complex architectures with three LSTM layers or a further enlargement of the number of LSTM cells would be needed to lower the molecule reconstruction error further.
The image to sequence model seems to be an outlier in comparison with the SMILES based models, in the respect that the latent space don't show much correlation with neither the SMILES sequence to be decoded or the molecular graph.The model produces a very low percentage of invalid SMILES and also has a low error rate with respect to molecule reconstruction, but is also decoding to canonical SMILES, which is an easier task then decoding to enumerated SMILES.The various other tests showed no big difference or benefit when compared to the much simpler use of different SMILES representations.On the other hand, the heteroencoder architecture may be useful for ar- The failure of the image to SMILES heteroencoder to produce significantly better latent representation fits with the observation that the latent space is mostly influenced by the decoding procedure, not the encoding procedure.The various encoders, whether based on images, canonical SMILES or trained on enumerated SMILES, seem to learn to recognize the molecules anyway and create a latent space that is best suitable for recreation of the decoders form.It thus seems that using enumeration techniques or other formats for the decoder will influence the latent space the most.
Training autoencoders on enumerated or different data further seems to improve the latent space with respect to its relevance for QSAR modelling.This is encouraging as it suggests that the extracted vectors are not only relevant for reconstruction of the molecular scaffold in itself, but additionally capture the variations underlying biological as well as physico-chemical properties of the molecules.It seem that already the encoder independence of the SMILES form for the enum2can leads to a more smooth latent space (c.f. Figure 9 panel A), which increases the relevance for QSAR modelling.This is in contrast to the results in Table 3, where a less skewed correlation to the decoded SMILES serialization in the encoder part is forced by training on enumerated data in the output, which however only leads to marginal gains in QSAR model performance.
The improvement seems quite marked and larger than what other studies have found.Winter et al. also used the heteroencoder approach in parallel to our work and found improvements for QSAR modelling [36].However, the improvements over baseline models were not as marked as in our results.The differences in network architectures (our use of bidirectional layers, LSTM vs. GRU's and batch normalization as example) and maybe also the choice of training data (Drug like molecules of ChEMBL) could be possible explanations.Future benchmarking on common datasets will likely show the way to the best network architecture and what unlabelled datasets to use for specific tasks.
On the other hand, the solubility dataset we used have previously been carefully modelled with chosen features and topological descriptors, resulting in a R 2 of 0.92 and a standard deviation of prediction of the test set of 0.6 [15].Likewise, a carefully crafted QSAR model of BCF obtained a R 2 of 0.73 and an RMSE of 0.69 [37], which is on par with our model using the can2enum derived latent vectors.However a later benchmark showed better performance for the CORAL software for prediction of BCF (R 2 : 0.76, RMSE: 0.64) [38], suggesting that further improvements are possible.
Thus the QSAR models based on heteroencoder derived latent vectors seem to almost match the performance of highly optimized QSAR models from selected features (c.f.Table 6), and it may rather be the ECFP4 and can2can model derived latent vectors that are mediocre for the tested type of QSAR tasks.Further, the ECFP fingerprints and auto-/heteroencoder derived latent vectors are of different dimension and nature.The fingerprints are 1024 dimensional, but binary, where the latent vectors are 256 dimensional and real valued.To make sure that the improvements were not due to different optima of the model hyper parameters for the different data, the neural network architectures for the QSAR models were optimized based on the ECFP4 fingerprint input.Some improvement of the fingerprint based models were observed, but reusing the ECFP4 hyper parameters for the latent vector based modelling still resulted in a large improvement in model performance for these input types.Further tuning of the hyper parameters of the models based on the latent vectors could likely further increase the performance to some degree (not tested).On the other hand, the denser dimensionality (256 << 1024) could help protect against over fitting and make the choice of hyper parameters less critical for these models.Either way, the use of heteroencoder derived latent vectors seem to be the better choice.
Feature generation for a dataset of chemical structures using the ChEMBL trained auto-/heteroencoders described in the publication is publicly available and hosted on the OSDR platform [39], where it is possible to encode molecular datasets into the latent vector space for subsequent uses, such as in QSAR modelling.
The increased relevance of the latent space with respect to bioactivity and physico-chemical properties are likely to increase the relevance and quality of the de-novo generated libraries where the neighborhood of as example lead compounds are sampled on purpose.However, the use of enumeration for training the decoder comes at the cost of greater uncertainty in the decoding, at a marginal improvement to the relevance of the latent space for QSAR modelling when compared to the enum2can model.On the other hand, the greater uncertainty and "creativity" in decoding could be beneficial and further help in creating more diversity in the generated libraries, but if this is the case has yet to be investigated.The choice of enumeration for decoder and/or encoder will thus likely depend on the intended use-cases.

Conclusion
The pilot study using a fully enumerated train and test set with 8 atoms showed that the latent space representation is sensitive to the chosen formats of the input and output in the training.Using canonical SMILES for the decoder gives a latent space representation which seem closer correlated to the SMILES strings than the molecular graphs.In contrast, training the encoders on input and output from different representations in chemical heteroencoders (image or enumerated non-canonical SMILES), gives a latent space with a better balance between SMILES similarity and a traditional molecular similarity metric.Forcing the encoder-decoder pair to trans-code between different formats or SMILES representation indeed seem to give a latent space closer to an abstract idea or encoding of the underlying molecule.The latent space properties were mostly found to be influenced by the choice of training data and representations used as the decoder targets.The multinomial sampling of the decoder shows higher variance for the decoders trained to predict enumerated data, both with regard to the SMILES form of the correct molecule but also by sampling other molecules.The changed properties of the decoders have broken the dependence on producing canonical SMILES, and may make them more relevant in de-novo design approaches in drug discovery, where a balance between similarity and variance is the goal.
The improved performance when using the latent space vectors from heteroencoders for QSAR modelling, further emphasizes their increased relevance, not just being a more SMILES independent representation of the molecule, but also for a better description of the chemical space relevant for biological as well as physico-chemical properties.This should hopefully lead to more drug-discovery relevant de-novo generated libraries.The increased relevance however comes at the price of greater uncertainty in the decoding, although more complex decoders seem to perform better at that metric.

Contributions
Esben Jannik Bjerrum developed the concept of heteroencoders, performed the training and tests using the GDB-8 datasets and prepared the initial [15] and final manuscripts.Boris Sattarov optimized the network architecture for ChEMBL data, prepared the QSAR datasets, trained and tested the performance of the latent vectors for QSAR application and helped prepare the final manuscript.

Figure 1 :
Figure 1: Enumeration challenge of a sequence to sequence model trained on canonical SMILES.The non-canonical SMILES of the same molecule is projected to different parts of the latent space reduced to two dimensions with principal components analysis (PCA).The small blue dots are the test set used for fitting the PCA.Some clustering of the enumerated SMILES can be observed.

Figure 3 :
Figure 3: Examples of optimal SMILES alignments of a molecule with two other molecules.The score is +1 for character match, -1 for mismatch, gap openings -0.5 and gap extension -0.05.Gaps are show with dashes, "-", and are not SMILES single bonds.

Figure 4 :
Figure 4: Scatter plot of the latent space similarities and the alignment scores of the SMILES strings.

Figure 5 :
Figure 5: Scatter plot of the latent space similarities and the circular fingerprint similarities.

Figure 6 :
Figure 6: Molecules similar in latent space using the can2can model.The reference molecule is in the upper left corner and similarity drops row-wise in normal reading direction.

Figure 7 :
Figure 7: Molecules similar in latent space using the can2enum model.The reference molecule is in the upper left corner and similarity drops row-wise in normal reading direction.

Figure 8 :
Figure 8: Venn diagram of the errors encounted during molecule reconstruction of 1000 molecules for the GDB-8 can2enum model.

Figure 9 :
Figure 9: SMILES enumeration challenge of the GDB-8 dataset based Enum2can and Can2enum encoders.The same three molecules were encoded from 10enumerated SMILES and projected to the latent space reduced to two dimensions with principal components analysis (PCA).Using enumerated SMILES for training of the encoder leads to the tightest clustering, but also trainingwith the enumerated SMILES in the decoder improves the clustering (c.f.Figure1).Small blue dots are the test set used for the PCA reduction.

Figure 2 :
Chemical heteroencoders are similar to autoencoders but translates from one format of the information or representation of the information to the other.The molecule toluene can be represented as a canonical SMILES strings, in different enumerated SMILES or via a 2D embedding.The autoencoder converts the canonical SMILES string to the latent space and back again (blue arrow), whereas many more possibilities exists for heteroencoders (green arrows).

Table 1 :
The datasets used for QSAR modelling.

Table 2 :
Hyper parameter search space using the Tree of Parzen estimator method inHyperopt.

Table 3 :
Properties of the models trained on different different input and output representations of the GDB-8 dataset.All values were calculated using the test dataset.

Table 4 :
Statistics on molecule generation with multinomial sampling at t=1.0, n=1000, GDB-8 dataset based models.Examples of different sampled molecules using multinomial sampling with the decoder from the two layer LSTM model (enum2enum 2-layer).The one in the upper-left corner is the reference molecule used to encode the latent space coordinates.

Table 5 :
Reconstruction performance on the ChEMBL datasets of the different encoder/decoder configurations

Table 6 :
Performance of the QSAR models on the held out test-set for different input data.The best performance for each metric and dataset is highlighted in bold.R 2 is the squared correlation coefficient (closer to one is better), RMSE is the root mean square error of prediction on the test set (lower is better).