Predicting Drug Solubility Using Different Machine Learning Methods --
  Linear Regression Model with Extracted Chemical Features vs Graph
  Convolutional Neural Network
Predicting the solubility of given molecules remains crucial in the
pharmaceutical industry. In this study, we revisited this extensively studied
topic, leveraging the capabilities of contemporary computing resources. We
employed two machine learning models: a linear regression model and a graph
convolutional neural network (GCNN) model, using various experimental datasets.
Both methods yielded reasonable predictions, with the GCNN model exhibiting the
highest level of performance. However, the present GCNN model has limited
interpretability while the linear regression model allows scientists for a
greater in-depth analysis of the underlying factors through feature importance
analysis, although more human inputs and evaluations on the overall dataset is
required. From the perspective of chemistry, using the linear regression model,
we elucidated the impact of individual atom species and functional groups on
overall solubility, highlighting the significance of comprehending how chemical
structure influences chemical properties in the drug development process. It is
learned that introducing oxygen atoms can increase the solubility of organic
molecules, while almost all other hetero atoms except oxygen and nitrogen tend
to decrease solubility.
John Ho, Zhao-Heng Yin, Colin Zhang, Nicole Guo, Yang Ha
Predicting Drug Solubility Using Different Machine Learning Methods - Linear Regression Model with Extracted Chemical Features vs Graph Convolutional Neural Network

###### Abstract

Predicting the solubility of given molecules is an important task in the pharmaceutical industry, and consequently this is a well-studied topic. In this research, we revisited this problem with the advantage of modern computing resources. We applied two machine learning models, a linear regression model and a graph convolutional neural network model, on multiple experimental datasets. Both methods can make reasonable predictions while the GCNN model had the best performance. However, the current GCNN model is a black box, while feature importance analysis from the linear regression model offers more insights into the underlying chemical influences. Using the linear regression model, we show how each functional group affects the overall solubility. Ultimately, knowing how chemical structure influences chemical properties is crucial when designing new drugs. Future work should aim to combine the high performance of GCNNs with the interpretability of linear regression, unlocking new advances in next generation high throughput screening.

## 1 Introduction

In the pharmaceutical industry, discovery of new drugs is both an expensive and time consuming endeavor. In order to accelerate the process and cut cost, scientists use early stage high throughput screening (HTS) to eliminate molecules without desired properties [1]. One such property is solubility, which governs the uptake, movement and metabolism of drugs in human bodies. Using theory or known experimental data to predict the solubility of molecules is a hot research field for decades. Back in 1968, Hansch et al. found that the octanol water partition coefficient (P) can be used to predict solubility [2]. Later, the Yalkowsky group came up with a general solubility equation (GSE) which included P and melting point (MP) [3]. Jorgensen and Duffy used Monte Carlo (MC) simulations to predict aqueous solubility with structural information such as molecular weight (MW), volume, solvent accessible surface area (SASA) and hydrogen bonds (HB) counts, as well as some other physical descriptors such as water Coulomb and van der Waals interactions (ESXL), hydrophobic and hydrophilic components. They achieved a reasonable prediction performance on a dataset of 150 organic molecules [4].

In recent years, with fast-growing computing power and new algorithms, researchers can afford to work with much larger data sets and more sophisticated machine learning (ML) models. There now exists multiple databases with
thousands of chemicals and their experimental solubility, such as the AQUASOL and PHYSPROP database used by Huuskonen et al. [5], the ESOL by Delaney [6], and some solubility handbooks [7]. AqSolDB is a newly developed database which merges multiple existing data sets [8]. From a methods standpoint, instead of traditional linear models and classic neural network (NN) models, the Barzilay group developed graph convolutional neural networks (GCNN) that convert the molecular structures into graphs which are input into a directed message passing neural network, and achieved state-of-the-art performance [9]. There has also been research going beyond drug solubility in aqueous solutions, like extending solute types to small proteins [10] or various organic solvents [11].

While these ML algorithms achieve outstanding performance, human scientists will find it difficult to gain mechanistic insights into the chemistry behind these solubility models. Such models are referred to as "black box": for example, it is almost impossible to understand a 20-layer deep learning NN or a GCNN with all molecules represented by huge matrices. From the perspective of chemists, we may have to step back from the focus on performance and seek more chemical insight. In this study, our goal is not to push the limits of predictive accuracy, but to take advantage of both classical models and modern sophisticated models to bring about a better understanding of the relationship between molecular structures and their chemical properties. With this knowledge, we may develop ML models in the future which are both highly accurate and human-interpretable.

## 2 Methods

Two ML models were applied in this study: A linear regression model and a GCNN model (Figure 1).

For the linear regression model, we used the molecular weight, total atom counts, and functional group counts as features to establish a multivariable regression with the experimental solubility values (logS). The features were obtained directly from the molecular structure using the RDKit module [12]. We also used L1 regularization with an alpha value of 0.01. For the GCNN, we used the Chemprop model [13]. Chemprop converts the atoms and bonds in the molecules into one-hot encodings, which are then concatenated into one tensor, representing each individual atom or bond. Next, Chemprop creates three tensors: one that maps each atom to its corresponding bonds (a2b); one that maps each bond to its corresponding atom (b2a), and one that maps each bond to its reverse bond (b2revb). It then concatenates each atom tensor into one large vector and each bond tensor into another large vector. With these five tensors, Chemprop finds the neighboring bonds of each bond and sums their vector representations. Finally, the model concatenates that sum to the vector representations of the bonds and atoms. The final vector representations of each bond are summed together to generate one feature vector for the entire molecule, which enters a standard feed-forward neural network with a single output (logS). We tested both models on three different datasets: the Delaney dataset, the Huuskonen dataset, and AqSolDB. We used 5-fold cross validation within each dataset, and used the root mean square error (RMSE) of the parity plots as the metric to evaluate the overall accuracy of the predictions.

Figure 1: **Setup of the linear regression model (top) and GCNN model (bottom), using the tyrosine molecule as an example.** The linear regression model uses human-engineered features such as molecular weight (MW) and number of functional groups to predict experimental solubility (logS), whereas the GCNN uses features learned via message passing along a graph.
## 3 Results

### Predicting Solubility

The parity plots for each model on different datasets are plotted in Figure 2, and the RMSEs are listed in Table 1.

For all three datasets, both the linear regression model and the GCNN model could make reasonable predictions, with the majority of predicted values less than 1 log unit away from the true values. This is comparable with other similar studies [14, 15]. The linear regression and GCNN model both showed the best performance (lowest RMSEs) on the Huuskonen dataset and the least-best performance on the AqSolDB dataset. Interestingly, the linear regression model shows relatively tight scatter on the Delaney and Huuskonen dataset but a very loose scatter on the AqSolDB dataset. This may be due to the large size of the AqSolDB dataset, as well as its relatively large solubility range.

The GCNN model overall did a better job than the linear regression model, especially when the dataset size is large. This is not surprising because the CNN is a more complicated model and has been shown to be a powerful method in many other fields, including computer vision. However, this does not suggest linear regression is inadequate. Considering the error due to experimental conditions such as pH and temperature, both models are sufficiently capable for drug design.

\begin{table}
\begin{tabular}{l l l l} \hline \hline Dataset & Size & RMSE, Linear Regression Model & RMSE, GCNN Model \\ \hline Delaney & 1127 & 1.13 & 0.59 \\ Huuskonen & 1282 & 1.08 & 0.49 \\ AqSolDB & 9982 & 1.83 & 0.76 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Performance of the Linear Regression Model and GCNN Model on Three Solubility Datasets

Figure 2: **Parity plots for the Delaney, Huuskonen, and AqSolDB datasets using linear regression model (top) and GCNN model (bottom).** Predictions are shown from the validation folds of 5-fold cross validation. Lines of best fit are shown in red.
### Understanding the relationship between molecular structure and solubility

Unlike the GCNN method that is a black box, the linear regression model can show us a straightforward relationship between the features we input and the solubility property we are interested in. We can run a feature importance analysis to see how each feature impacts the final results. The importance of each type of atom and functional group is shown in Figure 3 and Figure 4.

Figure 4: **The linear regression weights of each type of functional group feature for the Delaney dataset.** Positive weights indicate features contributing to a relative increase in solubility, whereas negative weights indicate features which contribute to a relative decrease in solubility.

Figure 3: **The linear regression weights of each type of atom feature for the Delaney dataset.** Positive weights indicate features contributing to a relative increase in solubility, whereas negative weights indicate features which contribute to a relative decrease in solubility.
From chemistry background knowledge, we know that it is intermolecular forces between solute and solvent molecules which determines solubility. Overall, more polar molecules with more hydrogen bonds (either donors or acceptors) tend to have higher solubility in aqueous solutions. The feature analysis results here give a more quantitative view of these effects. For example, O atoms have a very positive effect on solubility, because they both increase the overall polarity of the molecules and can form hydrogen bonds with solvent water molecules. Halogens, on the other hand, have known negative effects on solubility [16]. Although they can also increase polarity, they are not hydrogen bond donors and generally form weaker hydrogen bonds with water. This suggests hydrogen bonds play a major role in aqueous solubility and are even more significant than polarity. Furthermore, halogenated molecules tend to be soluble in hydrophobic solvents [17]. This characteristic holds significant implications for drug delivery across cell membranes, and this extended solubility topic is worth further investigation within the pharmaceutical domain.

Interestingly, adding functional group counts can improve the RMSE significantly compared with only considering the MW and atom counts (Table 2). This suggests that the same type of atom, when incorporated into distinct functional groups, can exhibit different effects on solubility (this can also be seen in Figure 4). Some atoms such as N, S and P can form different functional groups which can have either positive or negative effects on solubility. Furthermore, a protonated amine tends to have much higher solubility than a non-protonated amine, and thus environmental conditions such as pH also need to be taken into consideration. Additionally, the effects of temperature should be taken into account, as temperature not only impacts solubility but also influences reactivity and stability. Finally, it is important to note that the aforementioned ML models can only predict a single solubility value for a given molecular structure. However, in reality, we must deal with high-dimensional data comprising various solubility values under different conditions for each compound. Addressing this challenge would require significant effort in terms of data collection, data cleaning, and algorithm development. Eventually, we believe that a complex NN-based model, coupled with interpretable feature analysis, will emerge as the tool of choice, rather than a simple linear regression approach.

### From solubility prediction to drug design

As discussed above, the predictive performance of simple solubility models is adequate for high throughput screening, even using the GSE that was developed many years ago. The other side of the coin for such solubility studies lies in their potential to aid future drug design. This entails a reverse approach: If I want to design a drug molecule with a certain solubility value (or other physical properties), what functional groups should I include?

Using the feature importance analysis results in this research, scientists can have a general idea of what functional groups to include. For example, in order to increase aqueous solubility, one can introduce an OH group on a side chain, while to increase the ability to cross cell membrane, a halogen atom might be the best choice. However, in reality, when multiple factors need to be simultaneously considered, the capabilities of human decision-making can be quickly exceeded. This is where the GCNN model can become invaluable. Through a well-trained neural network that establishes connections between defined molecular substructures and their corresponding properties, coupling the GCNN with a molecular generative model [18] could allow viable drug candidates with desired properties to be generated at scale. We believe that this approach will power the next generation of HTS in the pharmaceutical industry.

## 4 Conclusion

This study explores the prediction of aqueous solubility using two models: a GCNN with learned features and a linear regression model with human-engineered features. Both models show reasonable prediction accuracy on a variety of datasets, with the GCNN performing better overall. However, the linear regression model can provide valuable insights into the relationship between certain features and solubility, highlighting the importance of some atoms, functional groups, and hydrogen bonds. We believe that a GCNN model, coupled with feature analysis, may be a valuable direction for future research.

\begin{table}
\begin{tabular}{l c c c} \hline \hline Dataset & Size & RMSE, Atom Features Only & RMSE, Atom and Functional Group Features \\ \hline Delaney & 1127 & 1.13 & 0.96 \\ Huuskonen & 1282 & 1.08 & 0.89 \\ AqSolDB & 9982 & 1.83 & 1.73 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Performance of the Linear Regression Model with Only Atom Feature and with Atom and Functional Group Features on Three Solubility Datasets