poor aqueous solubility 
 remains a major cause of attrition in the 
 drug development process . despite theoretical developments ,
the solubility 
 of druglike molecules still eludes truly quantitative computation . 
 in recent work
, we have shown that accurate 
 first - principles calculation is now becoming possible , provided that 
 both the crystalline and solution phases are described by accurate 
 theoretical models . before this
, energy terms from a computed thermodynamic 
 cycle ( see figure 1 ) had been used as descriptors 
 in a multilinear regression model for intrinsic solubility , delivering 
 accuracy much better than from direct computation and comparable with 
 the leading informatics approaches . since then , sophisticated machine 
 learning techniques have been 
 applied to many problems in the chemical sciences , while , as we have 
 shown , the accuracy of direct computation of hydration 
 energies and solubilities has improved significantly .
this led us 
 to revisit the idea of hybrid informatics - theoretical models for solubility .
cheminformatics methods have seen widespread use for property prediction , 
 particularly in the pharmaceutical industry where they have been applied 
 to ; aqueous solubility , melting point , boiling point , log p ( where p is the partition coefficient 
 between octanol and water ) , binding affinities , and toxicology predictions . such methods are usually much quicker than pure 
 chemical theory calculations , making high throughput virtual screening 
 ( htvs ) a possibility .
some methods have become accessible and easy - to - use 
 web - based tools .
however , informatics 
 methods suffer from the difficulty of decomposing the results into 
 intuitive , physically meaningful understanding and can not reflect 
 the physical details of the system . to understand the underlying physics 
 and chemistry , it is necessary to carry out an atomistic physics - based 
 calculation .
the exact nature of the theory varies between these 
 methods and the phase being studied .
crystal structures are often 
 modeled using one of the lattice energy minimizing simulation methods , plane - wave density functional theory ( dft ) methods , or periodic dft using atom - centered basis sets .
the simulation methods often contain empirical parameters , which 
 lowers the cost of these methods significantly , compared to dft .
popular solution - phase models include atomistic simulation methods 
 based on molecular mechanics and dynamics , quantum - mechanical implicit solvation methods ( such as the polarizable 
 continuum model ( pcm ) ) , and  hybrid  
 models ( such as the classical statistical mechanics - based reference 
 interaction site model ( rism ) or hybrid 
 quantum mechanics / molecular mechanics ( qm / mm ) methods ) .
these methods have the inherent problem for industrial 
 and drug discovery applications of being significantly more computationally 
 intensive than cheminformatics models , which makes high - throughput 
 computation infeasible .
the closest thing to an exception among contemporary 
 theoretical models may be 1d rism , which requires only a few minutes 
 of calculation time per compound and has been previously combined 
 with cheminformatics to build the 1d - rism / sdc method . by combining lower levels of theoretical chemistry 
 with cheminformatics , 
 we hope to produce results in good agreement with experiment , but 
 at a lower cost than higher - level theoretical methods , and with higher 
 accuracy than using cheminformatics descriptors alone .
a set of 100 broadly druglike 
 organic molecules was assembled with the prerequisites that each molecule 
 should have an available crystal structure in the cambridge structural 
 database ( csd ) and a well - documented 
 aqueous intrinsic solubility in the literature .
where possible , we 
 prefer experimental solubilities obtained with the cheqsol method , which has been shown to give reproducible results 
 with only small random errors .
the possibility of significant systematic 
 errors between different experimental methodologies remains an issue 
 and may possibly limit the accuracy with which modeling - based studies 
 can be validated .
a total of 122 potentially useful cheqsol 
 solubilities were obtained from the two solubility challenge papers and downloaded from the web
. while noting that several corrections had previously been 
 made , we also corrected or disambiguated the following names : amitriptyline , 
 5-bromogramine , 5,5-diphenylhydantoin , 4-hydroxybenzoic acid , nortriptyline , 
 and phenanthroline .
of the 122 compounds , 38 had corresponding crystal 
 structures and could be included in our dls-100 dataset . where a choice 
 existed , we selected the solubility and crystal structure of the least 
 soluble and , therefore , most stable polymorph . for druglike compounds 
 with known crystal structures ,
one further cheqsol solubility was 
 available from palmer et al . and two from 
 narasimham et al .
crystal structures were selected on the basis 
 of stability , preferring the polymorph with the lowest literature 
 solubility or the lowest lattice energy according to our computations 
 where polymorph - specific experimental information was not available . 

we also applied the additional pragmatic selection criterion that 
 the asymmetric unit cell should contain only one molecule .
once structures 
 were identified , they were downloaded in either the shelx format ( .res ) 
 or csd legacy format ( .dat ) .
we chose to use chemistry development 
 kit ( cdk ) molecular descriptors in this 
 study , because these descriptors 
 do not require proprietary software and are applicable to solubility 
 prediction .
the cdk is an open source 
 cheminformatics java library . in order to use the cdk molecular descriptors
we thus decided to use one principal source 
 for smiles records , selecting the well - annotated database chemspider .
since we are modeling intrinsic solubility , 
 we wish to describe the neutral form of the druglike compound .
this 
 remains the case even if a protonated or deprotonated charged form 
 dominates at neutral ph or across the ph range of the cheqsol ( or 
 other ) experiment . to obtain a smiles string for each molecule in 
 the dls-100 dataset
, we wrote a taverna workflow , which uses web services provided by the chemspider database .
the workflow is freely available on the myexperiment repository at the following reference . in five cases
thus , we instead 
 took the smiles from the solubility challenge web site for cimetidine , pindolol , and phenobarbital , 
 and from wikipedia for griseofulvin and 
 glipizide . using the resulting 100 smiles , 
 we initially calculated all 268 available nonprotein cdk descriptors 
 for each compound .
we found that 145 of these descriptors were either 
 undefined for 2d structures , or had the same value for all 100 compounds ; 
 their deletion left 123 remaining descriptors .
we took 
 experimentally determined crystal structures of the compounds in our 
 dls-100 dataset as the initial input to our calculations .
dmacrys , a periodic lattice simulation program , was used 
 to perform the crystal structure minimizations and calculate vibrational 
 contributions arising from the crystal .
dmacrys works in conjunction 
 with the gdma2 and gausssian 09 ( g09 ) 
 programs .
the output of these calculations 
 gives us the enthalpy of sublimation and crystal portion of the entropy 
 of sublimation .
the selected crystal structures were input into 
 dmacrys , which was used to standardize the covalent bond lengths between 
 hydrogens and heavy atoms , as the experimentally determined bond lengths 
 are not accurate , because of the uncertainty in the hydrogen positions 
 obtained by x - ray diffraction , before any calculations were run .
electrostatic 
 interactions were calculated by multipole expansions ( obtained using gdma2 ) of molecular charge distributions 
 calculated at the mp2/6 - 31 g * * level using g09 .
dmacrys carries out a rigid - body minimization of the crystal structure , 
 hence arriving at minimized lattice energies .
this lattice energy 
 can be converted to an enthalpy of sublimation by the following formula : enthalpy of sublimation : 
 1where ulatt is 
 the lattice energy ( energy of the crystal assuming the crystal is 
 static and at 0 k relative to infinitely separated molecules ) and 
 the 2rt term arises from lattice 
 vibrational energy .
the entropy of sublimation 
 was calculated by : entropy of 
 sublimation : 
 2where srot is 
 the rotational entropy in the gas phase and strans is the entropy of translation in the gas phase .
the use of eq 3 makes these assumptions : 
 ( i ) the rotational and translational entropy of the crystal is minimal , 
 ( ii ) there is no change in electronic entropy between phases , and 
 ( iii ) the intramolecular entropy is constant between the two phases . 

the crystal entropy is calculated by locating the frequencies of the 
 phonon normal modes ( lattice vibrations ) at the gamma point .
this 
 is achieved using lattice dynamics , the results of which are used 
 to calculate the helmholtz free energy ( see eqs 
 s2 and s3 in the supporting information ) .
the coordinates of a single molecule were 
 extracted from the minimized 
 lattice and used as input for the gaseous optimization with g09 .
optimizations 
 were carried out at the m06 - 2x and hf levels of theory with a 6 - 31 g * 
 basis set .
all solution - phase calculations 
 were carried out with g09 using the self - consistent reaction field 
 ( scrf ) protocol .
we selected the smd ( solvation model based on density ) implicit solvent model based on previous work . although rism yielded more - accurate absolute 
 hydration energies than smd in our recent work , smd generated a higher correlation coefficient against 
 experimental results for hydration free energy prediction ( r = 0.97 vs r = 0.93 )
given the parametrized 
 nature of our present model , correlation is more important than absolute 
 agreement , and , hence , smd is a suitable solvation model .
solution - phase 
 calculations were carried out with the same methodologies as used 
 in the gas - phase calculations , m06 - 2x/6 - 31 g * and hf/6 - 31g*. geometry 
 optimization was again carried out , this time taking the gas - phase 
 optimized structure as the starting point .
smd solves for the free energy 
 of solution ( ghyd ) as a sum of 
 the electrostatic contributions and nonelectrostatic contributions . 

the electrostatic contributions are calculated by the solution of 
 the nonhomogeneous poisson equation ; this equation 
 is a second - order differential equation linking the electrostatic 
 potential , dielectric constant , and charge distribution .
the nonelectrostatic 
 contributions of cavitation , dispersion , and solvent structure are 
 calculated as a sum of atomic and molecular contributions using parameters 
 inherent to the smd method .
smd has been shown to provide significant 
 improvements over some other implicit solvent models for datasets 
 containing molecules similar to those used in this study .
the hydration free energy is given by eq 4,gibbs free 
 energy of hydration : 
 4where esolution is the total energy of 
 the system in the smd solvation model and egaseous is the total energy of the system in 
 a vacuum .
this scheme typically takes 
 a few hours of calculation time per molecule on two 2.8-ghz 6-core 
 intel xeon x5660 processors .
sublimation energies were calculated 
 in the 1 atm standard state , which is the conventional standard for 
 experimental sublimation energies to be quoted . however , solvation 
 free energies are usually quoted in the ben - naim standard state of 
 1 mol / l . in this work
, g corresponds 
 to the 1 atm standard state , while g * corresponds 
 to the ben - naim 1 mol / l standard state ( see figure 2 ) .
the difference between these 
 two standard states is a constant energy value of 1.89 kcal / mol ( 7.91 
 kj / mol ) . in this work ,
we calculate the sublimation free energy in 
 the 1 atm standard state and then apply the correction to 1
mol / l 
 in order to be consistent with the hydration free energy calculations ; 
 hence , gsolu is in the 1 mol / l 
 standard state for all predictions in this work .
the presence 
 of an asterisk ( * ) denotes that the values refer to the standard state 
 of 1
our final solution free - energy 
 prediction is then given as the sum of the predicted sublimation and 
 hydration free energies : gibbs free 
 energy of solution : 
 5 therefore , we 
 have two predictions 
 for each molecule : the first method couples dmacrys with g09 and the 
 smd solvation model at the hf/6 - 31 g * level of theory .
the second method is dmacrys coupled with 
 g09 and the smd solvation model at the m06 - 2x/6 - 31 g * level of theory . 

this will be referred to as smd(m06 - 2x ) . for convenience of 
 comparison with experimental values of solubility
, 
 we convert the free energy of solution to log s values , 
 and all experimental solubility values to log s values:6here , r is 
 the universal gas 
 constant and t is the absolute temperature ( in kelvin ) .
the conversion of experimental solubility to log s can be found in the supporting information ( eq 
 s7 ) .
values for the full dls-100 dataset , including smiles 
 and inchi , can be found in the supporting information 
 ( see zip file and dataset ) .
to model the data , we use linear 
 and machine learning regression models : partial least - squares regression , 
 random forest and support vector regression . for reporting the predictive 
 accuracy of these models
, we averaged the rmse of log s over a 10-fold cross - validation of the dls-100 dataset . the cross - validation 
 fulfils two purposes in this study : parameter optimization and evaluation 
 of the accuracy of the models on unseen data . to ensure that each 
 test fold of data is truly unseen ,
the parameter optimization is carried 
 out in a separate layer of cross - validation within the training folds , 
 as we will discuss below . in order to avoid overfitting ,
the use of multivariate data presents 
 a danger of overfitting machine learning regression models ; moreover , 
 redundancy of attributes and correlation within the data add to the 
 risk of reaching misleading conclusions .
one 
 is the commonly used standardization method of variable scaling , equalizing 
 the distributions of the variables by normalizing the mean and standard 
 deviation of each column ( variable ) .
the 
 advantage of using this method is that it equalizes the prior importance 
 of all the attributes .
the second normalization method is principal 
 component analysis ( pca ) , transforming the data into a smaller subspace 
 where the new variables are uncorrelated with each other .
the pca data transformation method deals with 
 the redundancy of the data , and places emphasis on the variance of 
 the data .
the ability of each principal component to explain the data 
 is measured according to the variance accounted for .
third , we have 
 also fitted each model on the nonpreprocessed raw dataset , for comparison 
 with the results of the two different scaling methods . in this section , 

a summary of the regression models are presented ; detailed explanations 
 can be found in the supporting information .
the partial least 
 squares regression ( plsr ) model design is appropriate in a situation 
 where there is no limit to the x variables or predictors , 
 or where the sample size is small .
moreover , the plsr model is also 
 beneficial for analyzing strongly colinear and noisy data .
the goal 
 of a plsr model is to predict the output variable y from the input variables x and to describe the 
 structure of x. for this , plsr finds a set of components 
 from x that are relevant to y ; these 
 components are known as latent variables .
the intention of plsr is 
 to capture the information in the x - variables that 
 is most useful to predict y. a graphical representation is supplied in the supporting information figure s1(a ) .
random forest ( rf ) , a method 
 for classification and regression analysis , has very attractive properties 
 that have previously been found to improve the prediction of quantitative 
 structure  activity relationship ( qsar ) data .
an ensemble of many decision trees constitutes a random 
 forest , and each is tree constructed using the classification and 
 regression trees ( cart ) algorithm .
the 
 rf method is efficient in handling high - dimensional data sets and 
 is tolerant of redundant descriptors .
the main idea in support 
 vector regression ( svr ) is to minimize the risk factor based on the 
 structural risk minimization from structure 
 theory , to obtain a good generalization of the limited patterns available 
 in the given data .
first , the given data d are mapped 
 onto a higher dimensional feature space , using the kernel function k(xi , xj ) and then a predictive function 
 is computed on a subset of support vectors . here , we have used the 
 radial basis kernel function ( eq 7 ) to map the 
 data onto a higher dimensional space .
a graphical representation is 
 supplied in the supporting information ( figure 
 s1(b)).svr mapping 
 on radial basis kernel function : 
 7 to evaluate the performance of 
 various machine learning models , we report two statistics : the root 
 mean squared estimate ( rmse ) and squared pearson correlation coefficient r ( not to be confused with the coefficient of 
 determination ) . formulas for these are 
 given in the supporting information ( eq s5 ) .
we have also assessed statistical significance using menke and 
 martinez s method , which we have 
 used previously for similar analysis ( see supporting information ( eq s6 , tables s3s9 for r , and boxes s1s3 ) for statistical 
 significance ) .
we also analyzed the variable importance for the rf 
 method ( see table s17 in the supporting information ) .
variable importance was calculated in the cart program as implemented 
 in r. in order to compute and compare 
 the performance of the various regression models , we consider rmse 
 scores averaged over a 10-fold cross - validation . in the 10-fold cross - validation ,
the dataset is randomly 
 split into 10 partitions , where the training set consists of 90% of 
 the data and the test set consists of 10% of the data .
the predictivity on 
 the test fold is considered as an external measure to compute the 
 accuracy of the fitted model .
the entire process is repeated 10 times 
 in order to cover the entire dataset , with each fold forming the test 
 set on one occasion , and we record the average rmse .
the complete 
 design of the workflow is represented in a flowchart ( scheme 2 ) ; similar workflows have been used for classification 
 in other studies .
the complete workflow of this 
 analysis was written in r using the caret 
 package ; all scripts are available in 
 the supporting information . in out - of - bag validation ,
one evaluates the performance 
 of the 
 model by separating training and test data through bootstrap sampling ; 
 this is convenient only for the rf method .
it is not appropriate to 
 compare rf out - of - bag predictions with other models such as pls and 
 svr , which are not based on bootstrap sampling .
so , we used 10-fold 
 cross - validation to evaluate the performance of our various models . for each 
 model
, we use 90% of the total data designated as the training set 
 in order to find the optimum values for these parameters .
we selected 
 a range incorporating 20 different possible values for each model 
 parameter , in order to select its best value .
for each parameter , 
 a further level of 10-fold cross - validation is carried out in order 
 to retrieve the rmse of the models using each possible parameter value . 
 here ,
the training portion of 90% of the original data is further 
 split into 10 new folds of 9% , with nine ( 81% of the original data ) 
 being used to build each model and one ( 9% ) as an internal validation ; 
 this process of model building and internal validation is repeated 
 to predict each of the 10 possible internal validation folds .
this 
 internal cross - validation step is repeated 20 times , once for each 
 possible value of the parameter being assessed .
then , based on the 
 value giving the lowest average rmse score in the internal validation 
 folds , the optimum parameter value is selected . finally , the model 
 is fitted on the complete training set of 90% of the original data 
 using the selected parameter values .
the given 90%:10% split of the 
 data into training and test sets was 
 used to fit the final model for each fold of the main 10-fold cross - validation , 
 once the optimum parameter values have been selected .
the average 
 rmse and r values over the 10 folds were 
 considered in order to compare the usefulness of different descriptor 
 sets and to evaluate the performance of the fitted models .
the full dls-100 dataset , with the experimental 
 log s values , can be found as supporting information or downloaded from the mitchell group 
 web server ( http://chemistry.st-andrews.ac.uk/staff/jbom/group/informatics_solubility.html ; see the supporting information ( csv_smiles_si.csv 
 and table s1 ) ) , which is consistent with the excellent suggestions 
 from walters .
the dataset includes csd 
 refcodes , chemspider numbers , smiles , experimental log s values and inchi for all molecules .
the log s values 
 in this work come from refs ( 2 , 16 , 18 , and 19 ) .
where 
 possible , we have selected data obtained from the cheqsol method ; 
 where this was not available , we have selected reliable sources using 
 different determination techniques .
a good solubility prediction can 
 be considered as a prediction of approximately the same error as that 
 of the experiment .
the experimental values have been shown in a number 
 of previous papers to vary considerably . here
, we consider the experimental accuracy limit to be between 
 0.6 and 1 log s unit ( where 1 log s unit represents 5.7 kj / mol at 298 k ) .
previous work has reported 
 the experimental error in solubility prediction to be as great as 
 1.5 log s units and , on average , the error to be 
 at least 0.6 log s units . in 2006 ,
dearden noted , as was later 
 reiterated in the solubility challenge , that models with rmse predictions 
 of < 0.5 log s units are likely to be overfitted . for a prediction to be useful
, it must have an rmse within the standard 
 deviation of the experimental data ; otherwise , a trivial prediction 
 using the mean of the experimental data is a more accurate prediction 
 of the log s value . for 
 the dls-100 dataset ,
first , a purely theoretical prediction , in which no machine learning 
 is used and where predictions are made using only physics - based calculations . 

second , theoretical energies are used as the sole descriptors in machine 
 learning models .
third , cheminformatics descriptors , calculated using 
 the cdk , are used as the sole input to machine learning methods .
finally , 
 cheminformatics descriptors and theoretically computed energies are 
 combined as input to machine learning methods . for each of these methods , 

we present the results and discussion , with comparison between the 
 methods made on the basis of rmse and r ( correlation coefficients for cheminformatics and combined models 
 can be found in the supporting information ( tables 
 s3s9 ) ; rmse values can be found in the supporting information ( tables s10s16 ) ) .
in addition to these results , we have replicated the solubility 
 challenge using 2d molecular descriptors alone .
the theoretical methodologies 
 described earlier utilize a thermodynamic cycle to access the free 
 energy of solution .
table 1 shows the r correlation coefficient and the rmse for the 
 predictions made by these methods .
chart 1 shows 
 the linear fit to the data from the smd(hf ) method , which has the 
 lower rmse and the higher r correlation 
 coefficient of the two purely theoretical methods .
the rmse for the smd(hf ) method is nearly three 
 times the suggested criterion of 1 log s unit of 
 error .
the situation is even worse for the smd(m06 - 2x ) method for 
 which the rmse is just over four times this criterion ( see charts s4s6 in the supporting information ) .
both methods produce results outside the useful prediction criterion 
 of 1.71 log s units . from these results
first , it is clear that the given methodologies 
 do not adequately quantify the physics occurring in the solution process 
 ( i.e. , solid to solution ) .
second , we can conclude that , if it is 
 possible to explain the underlying structure of these data using a 
 general model , based on the predicted log s values , 
 such a model will be inherently nonlinear . compared with our 
 previous work , in 
 which theoretical models provided a good prediction of log s ,
our theoretical methodology here differs only marginally , 
 in the use of mp2 multipoles , and still produces good results ( see supporting information ( chart s1 and table s2 ) ) for the same 25 molecules in this work ( dataset dls-25 ) .
the predictions 
 for the additional 75 molecules alone show worse predictions than 
 for the full 100-molecule set presented above ( see charts s2 and s3 in the supporting information ) .
the additional 
 75 molecules therefore appear to form a more difficult dataset to 
 predict .
it is likely that improved results can be obtained from purely 
 theoretical calculations , if some of the approximations made here 
 are improved ; for example , improved modeling of the solvated phase 
 to more accurately describe the solvent and its effects on the solute 
 could increase accuracy . also , we note that the intramolecular degrees 
 of freedom are neglected in the dmacrys calculations , and further 
 assumptions are made by using eqs 2 and 3 in the methods section .
we subsequently applied machine learning methods to the theoretical 
 energies in order to carry out nonlinear regression analysis .
the 
 average rmse scores over 10-fold cross - validation ( see the methods section for details ) is represented as two - dimensional 
 ( 2-d ) column charts ( see charts 2 and 6 ) .
different grayscale column bars represent the 
 different machine learning methods used in this study .
the standard 
 deviation is shown as an error bar ( black line ) . the use 
 of the calculated energies as descriptors in the machine 
 learning models yields considerably improved results , compared to 
 those from the predictions made without machine learning .
the results 
 now , while still missing the 1 log s unit error criterion , 
 do make useful predictions in which the rmse is within the standard 
 deviation of the experimental data ( 1.71 log s units ) . 

the rf and svr models produce notably better results than pls . charts 2 and 3 show that the method 
 minimizing the rmse ( 1.21 log s units ) is rf with 
 hf when scaled with pca .
an additional point of interest is that the chemical descriptors 
 alone using rf or svr can provide a marginally better prediction of 
 log s than the machine learning methods with only 
 the energies as descriptors .
in particular , we noticed that fitting 
 the rf model on data that are scaled to a given mean and standard 
 deviation produces a statistically significant improvement in its 
 prediction with cheminformatics descriptors alone rather than theoretical 
 energies ( see the supporting information ( boxes 
 s1s3 ) ) . in all other cases ,
this suggests that slightly more useful information about the molecules 
 log s values is conveyed by the cheminformatics descriptors 
 than by the theoretical energies alone ( see chart 4 ) .
when the descriptors and energies are 
 combined as input for the machine learning methods , we obtain results 
 that are generally only very slightly better than those obtained from 
 cheminformatics descriptors alone .
this implies that the theoretical 
 energies contain very little extra useful information not already 
 present in the descriptors .
the joint results do present a statistically 
 significant improvement for pls and rf , once scaled by the mean / standard 
 deviation , compared to those for the theoretical energies alone . in 
 light of this , and given that the descriptors alone produce a marginally 
 improved result compared to chemical theory , it is fair to say the 
 cheminformatics descriptors are seen to contain a modest amount of 
 additional information not incorporated in the theoretical energy 
 terms .
this suggests that the 123 descriptors of the cheminformatics 
 descriptors and the 10 theoretical energy descriptors convey similar 
 information , with only a small amount of additional information being 
 conveyed by adding the descriptors to the energies and almost no information 
 gained by adding the energies to the set of descriptors .
interestingly , the best result in terms of rmse is from the descriptors 
 with the m06 - 2x energies , which , on their own , produced the worse 
 of the two pure theory results in this work ( see charts 5 and 6 ) .
the rf model performs 
 particularly well over all descriptor sets , even without any type 
 of scaling , the best rmse result being only 0.13 units outside the 
 1 log s unit target . the best single prediction , 
 in terms of the rmse , was made by the pls model , using descriptors 
 and the m06 - 2x energies scaled by the standard deviation and the mean , 
 with an rmse of 1.11  0.04 log s units .
all 
 of these methods make predictions inside the standard deviation of 
 the experimental data ; therefore , all of the predictions are useful . 

we also note that the rf model shows small but statistically significant 
 improvements with all scaling methods ( using the theoretical energies 
 and cheminformatics descriptors combined ) when compared to some models 
 trained on the theoretical energies only ( see supporting information ( boxes s1s3 ) ) .
this is the 
 only model to show such improvements with all scaling methods in the 
 present work .
we analyzed the relative variable importance ( see table s17 in the supporting information ) and 
 found that x log p ( from ref ( 57 ) ) was consistently rated 
 as the most important feature .
x log p is a computed estimate of the base-10 logarithm of the octanol : water 
 partition coefficient ( the ratio of concentrations of solute solvated 
 in the two different solvents ) .
this has been seen in many previous 
 studies and is not so surprising given that it provides information 
 specifically about the solvated phase.x log p uses an atom additive model for the prediction 
 of log p. in the supporting information , we include tables ( table s17 in the supporting 
 information ) displaying the 10 most important descriptors ; 
 here , we will briefly comment upon these .
we find kier and hall s 
  path and chain indices to be 
 of importance ; these quantify the degree of bonding to heavy atoms 
 within a given path or chain length .
broto 
 autocorrelation , which describes the 
 charge and mass distribution along a given path length , is found in 
 the top 10 .
finally , we also note randic s weighted path descriptors , which are used to account for molecular branching . 

once the theoretical energies are added to the descriptor set , the 
 free energies of hydration and solution are ranked in the top 10 , 
 along with the theoretical log s prediction .
explanations 
 of the molecular descriptors used in this work can be found in ref ( 61 ) .
one may expect that molecular branching would play an important role , 
 because it gives information on the extent and flexibility of the 
 molecule , hence contributing some entropic information .
coupling this 
 descriptor with the kier hall descriptors , information can be acquired 
 on the composition of such chains , in terms of heavy atoms .
again , 
 here , information is imparted concerning the distribution of heavy 
 atoms and electronic factors .
for example , the degree of charge separation 
 across a molecule and the localization of charges are important factors 
 in determining particularly enthalpic but also entropic contributions . 
 the theoretical energies in the top 10 are all closely related quantities ; 
 it is not surprising that the ( purely theoretical ) prediction of log s is found in the top 10 : since this is the quantity we 
 are trying to predict , it is expected to provide sufficient information 
 to the model to be found in the top 10 .
the free energies of solution 
 and hydration provide direct information from electronic structure 
 theory and statistical thermodynamics on the interactions of a given 
 molecule , in a given conformation , within its environment , and on 
 the energetics of phase transitions . as a benchmark , we also 
 present our method s predictions 
 of the solubility challenge set based solely on cheminformatics descriptors 
 ( see table 2 ) . as
suitable crystal structures 
 are not available for all molecules in the solubility challenge , we 
 could not calculate the theoretical energies .
tables 2 and 3 , and 
 chart 7 , demonstrate that our method can make 
 predictions for the solubility challenge dataset within the coveted 
 1 log s unit rmse error and , in fact , makes predictions 
 that are consistent with some commercially available methods and deep - learning 
 methods .
a recent publication reported 
 rmse scores of 0.95 log s units for the commercially available package mlr - sc and 0.90 log s units for a 
 deep - learning method .
however , these 
 results are not directly comparable with ours , for two reasons .
first , 
 our results have been calculated for a 10-fold cross - validation and 
 for the canonical training : test split ( see tables 2 and 3 ) .
second , the deep - learning 
 result ( rmse = 0.90 ) given by lusci et al . is contingent on correcting eight putative errors in the cheqsol 
 solubility data , the most substantial of which is for indomethacin , 
 a compound that has been shown to hydrolyze under alkaline conditions .
while we have corrected names and smiles for 
 the solubility challenge set , we have not adjusted any solubility 
 values therein .
it is also reasonable to suggest that , using the solubility 
 challenge set as a benchmark , our 100-molecule set could be considered 
 as a 
, given the improved prediction 
 offered by our method when the solubility challenge set is used instead .
our current work shows that accurate solution 
 free energies are 
 not calculable via the simple theoretical procedure that we present 
 here . a significant portion of the important physics in the solution 
 process is not captured using the approximate methodologies that we 
 utilize in this work .
this reaffirms that , currently , qspr methodologies 
 are the most - accurate and time - efficient methods for accurate solution 
 free energy predictions . in addition , we show that state - of - the - art 
 machine learning methods , with a modest number of cheminformatics 
 descriptors , are capable of making solution free - energy predictions 
 that are consistent with those of commercially available programs 
 and newer deep - learning approaches . here
, theoretical energies and 
 cheminformatics descriptors are generally shown to not be complementary 
 for such predictions .
since both sets of descriptors ( theoretical 
 energies and cheminformatics descriptors ) produce a similar level 
 of accuracy when used alone in the machine learning methods , and little 
 improvement is seen when they are combined , we can conclude that the 
 information conveyed is of a similar nature and that the theoretical 
 energies are , for this reason , a more efficient form of information 
 storage , as 10 descriptors contain equivalent information to 123 molecular 
 descriptors .
however , in terms of time , the molecular descriptors 
 are much less expensive to calculate and their use is therefore more 
 time - efficient . additionally , we note that the rf method has produced 
 promising predictions in this work , with relatively low rmse
. this 
 method has consistently produced good results and would be our recommended 
 method to make solubility predictions .
 we 
 present four models of solution free - energy prediction for druglike 
 molecules utilizing cheminformatics descriptors and theoretically 
 calculated thermodynamic values . 
 we make predictions of solution free 
 energy using physics - based theory alone and using machine learning / quantitative 
 structure  property relationship ( qspr ) models . 
 we also develop 
 machine learning models where the theoretical energies and cheminformatics 
 descriptors are used as combined input . 
 these models are used to predict 
 solvation free energy . 
 while direct theoretical calculation does not 
 give accurate results in this approach , machine learning is able to 
 give predictions with a root mean squared error ( rmse ) of 1.1 
 log s units in a 10-fold cross - validation for our 
 drug - like - solubility-100 ( dls-100 ) dataset of 100 druglike molecules . 
 
 we find that a model built using energy terms from our theoretical 
 methodology as descriptors is marginally less predictive than one 
 built on chemistry development kit ( cdk ) descriptors . 
 combining both 
 sets of descriptors allows a further but very modest improvement in 
 the predictions . 
 however , in some cases , this is a statistically significant 
 enhancement . 
 these results suggest that there is little complementarity 
 between the chemical information provided by these two sets of descriptors , 
 despite their different sources and methods of calculation . 
 our machine 
 learning models are also able to predict the well - known solubility 
 challenge dataset with an rmse value of 0.91.0 log s units . 