the prediction of in vivo toxicity based on in vitro measurements 
 is a challenging task , which is in 
 the center of active development of modern computational toxicology .
the topcoder data science competition platform in collaboration 
 with environment protection agency ( epa ) organized the toxcast challenge 
 in 2014 .
lel is defined as  the lowest dose that 
 shows adverse effects in these animal toxicity tests .
the lel is then 
 conservatively adjusted in different ways by regulators to derive 
 a value that can be used by the agency to set exposure limits that 
 are expected to be tolerated by majority of the population . 
the total toxcast challenge included five 
 consecutive subchallenges , 
 which were executed over a seven month period and attracted 432 registrants 
 from 32 countries .
the first subchallenge was to identify software 
 libraries and/or methods to describe the chemical structure of various 
 compounds .
the second subchallenge was about identification of a specific 
 combination of in vitro assays , which could be used 
 to predict the in vivo systemic lel .
the third subchallenge 
 was executed privately and was entitled  predictive capability 
 tests  .
it had a goal  to build a prediction model ( algorithm ) 
 using data from high - throughput in vitro assays , 
 chemical properties , and chemical structural descriptors to quantitatively 
 predict a chemical s systemic lel .
 the final fifth subchallenge was about the documentation 
 of the results of the models . in this study , we present the 
 results of the rank - i model for the 
 prediction challenge . according to the challenge rules , the participants 
 were strictly forbidden to use any data other than the data that were 
 provided in this competition .
this was done in order to offer equal 
 conditions to all participants as well as to better evaluate the performance 
 of models developed using in vitro measurements . 
 indeed
, the use of any information outside of the data provided within 
 the challenge , e.g. , information about toxicological chemical pathways , 
 could potentially bias the comparison of algorithms .
the summary of 
 this prediction challenge was published by the epa . since september 2015
, this information has been available 
 from the web archives and included as supporting information ( see section  epa 
 toxcast lelpredictor marathon match results summary  ) to this 
 study .
this article analyzes the steps that were taken to develop 
 the 
 rank - i submission model by participant  novserj 
the supporting information also contains an extended 
 technical description , which was submitted to topcoder as part of 
 the contest .
data , rank - i model , and model predictions 
 are available at http://ochem.eu/article/68104 .
we believe 
 that this article will be interesting to both participants and organizers 
 of the challenge and will help them to better understand and interpret 
 the results of the challenge .
moreover , the model 
 reported in this study has the highest prediction ability for lel 
 end point as validated by the challenge organizers and thus can be 
 of potential value for people working on the risk assessment of chemical 
 compounds .
the in vitro measurements provided within the scope of the challenge included 
  a battery of more than 700 biochemical and cell - based in vitro assays to identify what proteins , pathways , and 
 cellular processes these chemicals interact with and at what concentration 
 they interact .  the total 
 data set used during the challenge incorporated 1,854 molecules .
the 
 experimental lel values were provided for 483 compounds that were 
 used as the training set .
the test set included lel values for 143 
 chemicals , which were kept secret and were split into provisional 
 ( 63 ) and final scoring ( 80 ) sets . during the challenge ,
topcoder participants 
 could submit predictions and receive the statistical results for the 
 provisional set .
the 80 compounds from the final scoring 
 set were used only once to rate the final model submissions .
the challenge organizers did not reveal information about the compounds , 
 which were used as the test sets .
after the challenge , the experimental 
 values for 143 compounds from the test set were kindly released by 
 the topcoder organizers to the authors of the article .
the detailed description of the rank - i submission is provided in supporting information ( see  technical 
 description  section ) .
ten different descriptor packages 
 implemented in the public platform ochem were used individually to create ten models for the resulting consensus . 

four descriptor packages , e - state , qnpr , isida fragmentor , and gsfrag , were based on the 2d representations 
 of the chemical structures .
the other six packages , inductive , chemaxon , adriana , mera / mersy , cdk , and dragon descriptors , 
 were calculated using 3d representations of the chemical structures . 

exactly the same procedure was used for all 
 models , and thus , the same number of descriptors and molecules were 
 utilized to develop each model with different machine learning methods . the selected descriptors count for each model after the 
 unsupervised filtering
ad is the applicability 
 domain of 
 the model as defined by ochem ( see also 
 ref ( 20 ) ) .
the model 
 used in the challenge 
 was developed with associative neural networks ( asnn ) .
it can be considered in a simplified 
 way as a combination of k - nearest neighbor ( knn ) method applied in 
 the space of ensemble predictions .
the models developed with the asnn 
 were top - ranked in several benchmarking studies and that is why this method was selected for the epa challenge .
the default parameters for the asnn algorithm , as optimized during 
 previous studies and provided on the ochem web site were used .
they 
 included 64 neural networks in ensemble , 3 neurons in a hidden layer 
 trained by the supersab algorithm .
in addition to asnn , we also analyzed knn , support vector machines 
 ( as implemented in libsvm ) , and partial least squares methods . as with the asnn method , the default 
 parameters of these algorithms as provided on the ochem site were 
 used .
the unbiased estimation of the 
 models performance is critically important for selection and 
 decision making for development of models .
two protocols , cross - validation , 
 and bagging are frequently used to estimate validation accuracy for 
 the training set .
it uses n  1 subsets as the training set and predicts the remaining 
 chunk of the data .
bootstrap aggregation ( bagging ) is another 
 powerful approach to develop and validate models developed by leo 
 breiman .
it is based on the aggregation 
 of models , each one of which is developed with its own training set 
 (  bag  in the terminology of breiman ) .
each bag is formed 
 by random sampling with replacement from the initial training set 
 and has the same size as the initial set .
the molecules ( on average 
 37% ) , which do not participate in the respective training set , are 
 called  out - of - the - bag  .
the predictions for these molecules 
 are used to evaluate the predictive power of models .
, there were several theoretical 
 developments to identify 
 the most significant descriptors for neural networks .
some of these methods calculate the sensitivity ( importance ) of 
 input parameters according to derivatives of neural network weights 
 with respect to the error function , while 
 others provide such estimations based on the analysis of the magnitudes 
 of the neural network weights . for this 
 study
, we used a method from the second group , which provided the 
 best results in our previous studies .
the sensitivity si of a neuron i was calculated as1where wij were weights connecting neuron i and j , maxa was taken over all weights 
 ending at 
 the neuron j having sensitivity sj , and summation was taken over all weights 
 connecting the neuron i with the upper layer neurons . 

the sensitivity calculations were performed recursively starting from 
 the last layer neuron , which had sensitivity of 1 .
the pruning 
 procedure started once neural network training was completed . on each 
 step ,
the least significant descriptors with smallest si were eliminated , and the models were 
 retrained with the decreased set of descriptors .
the sets of descriptors , 
 which calculated the minimal errors , were considered as the optimal 
 ones . in order to avoid overfitting and overtraining ,
the neural networks were trained using the efficient partitioning 
 algorithm , which uses the early stopping 
 procedure . the root mean square error ( rmse ) 
 metric
the challenge organizers used the following scoring 
 function1ato rank the models . as a result , models with 
 lower rmse got a higher score and higher rank among the others . in 
 addition to rmse , the organizers also reported pearson correlation 
 coefficients and auc defined as  percentage of pairs where 
 predicted1 < predicted2 among those where ground_truth1 <
the workflow for the model 
 development used in the challenge ( see 
  technical description  in the supporting information ) was based on our previous expertise to develop 
 recently published models .
final statistical results for 
 the top - ranking models are summarized in table 2 .
below , we provide a detailed analysis of 
 the steps , which were used to develop the model .
prediction accuracy for the  out - of - the - bag  
 samples .
confidence intervals 
 were estimated 
 using the subsets , which were sampled from the training set , and each 
 had the same size as the respective test set ( see for more details 
 ref ( 23 ) ) .
best model based on the in vitro assay descriptors developed using the libsvm method 
 ( see also table s1 ) .
model based on molecular weight 
 ( mw ) and number of carbon atoms ( nc ) developed using the same approach 
 as the above in vitro model .
there were 37 molecules , including 
 11 molecules from the training set , for which descriptor generation 
 failed for different packages .
cdk descriptor package does not support 
 inorganic elements such as [ sn ] , [ hg ] , [ b ] , and [ as ] .
the failed molecules 
 either included unsupported atoms for cdk or some groups , e.g. , [ n3 + ] . 

several other molecules , e.g. , rifampicin , -cyclodextrin , milbemectin , 
 emamectin benzoate , etc . were large chemical structures and failed 
 either due to time - out or structure conversion problems . according 
 to the challenge rules ,
as a simple 
 solution , we used an average value of all training set molecules , 
 loglel = 3.2602 log(m ) , as the predicted values for the failed 
 molecules the challenge organizers offered a provisional set of n = 63 compounds for the purpose of model analysis and selection . 

however , we decided to skip the testing on this set for the following 
 considerations .
thus , an attempt to rely on the models performance 
 for this set by , e.g. , submission of multiple predictions and selection 
 of a 
best  model using it , could contribute a higher 
 uncertainty and result in the selection of a nonoptimal model for 
 the final set .
indeed , the final model rmse was 0.88  
 0.04 for n = 472 training set molecules .
the provisional 
 set was not available , 
 and thus , we were unable to calculate the confidence intervals for 
 it .
we estimated the intervals by random sampling of n = 63 molecules from the training set , for each of which we calculated 
 the intervals .
thus , selecting the best model based on 
 the performance for the provisional test set is about twice as uncertain 
 compared to selecting based on the training set .
therefore , instead 
 of relying on the accuracy of models for the provisional test set , 
 a strategy to rely on the estimated validated results for the training 
 set is more reliable .
an even better strategy could be to select a 
 model based on the combined accuracies of the provisional and training 
 sets , but such an analysis was not implemented .
the confidence 
 intervals for n = 80 molecules 
 were about the same as that for n = 63 molecules . 

the wide confidence intervals for both provisional and final test 
 sets might have contributed to the fluctuations of ranks of challenge 
 models for both sets .
for example , the top final scoring model was 
 only ranked # 8 for the provisional submission , while the fourth model 
 was ranked # 27 .
vice versa , the models ranked top # 1 and # 4 for provisional 
 submissions were ranked as # 9 and # 34 ( out of 47 participants ) for 
 the final test set .
thus , indeed , the 
 provisional ranking score was not strongly predictive of the final 
 one : provisional and final models ranks were correlated only 
 with correlation coefficient r = 0.76 .
the 
 rmse of the eight top - ranked final models were in the range 
 1.12 to 1.16 and thus were within the confidence intervals of the 
 winning model .
thus , statistically speaking , these models had the 
 same performance , and their differences in performance were due to 
 chance .
the model 
 submitted to the topcoder challenge was a consensus of the bagging 
 models developed with the asnn method . in this section
, we briefly 
 describe the considerations that were used to develop and select this 
 model for the challenge .
below , we compare the performance 
 of different methods , which are described in the methods section . bagging vs cross - validation
table s1 in supporting information demonstrates 
 that models developed using bagging had consistently smaller validation 
 rmses as compared to the cross - validation results .
comparison of different machine 
 learning methods ( table s1 ) shows that 
 combinations of machine learning methods 
 and descriptors provided quite similar performances with rmse ranging 
 from 0.9 to 1.2 log units .
considering that 95% confidence intervals 
 of rmse were 0.4 log units , the majority of these models were 
 not significantly different .
the highest rmse = 1.2 ( i.e. , the lowest performance ) was 
 calculated using the pls method for mera / mersy descriptors . actually , 
 the failure of this method was due to several outlying molecules that 
 had predictions far beyond the range of the training set values
. they 
 may be due to the sparseness of descriptors used to develop models 
 and insufficient number of data points used to fit the coefficients 
 in pls .
if we limited the predicted values for all compounds to the 
 range of the training set values , the results of pls models became 
 similar to those of other methods .
it is interesting that the 
 model developed using in vitro assay measurements 
 consistently provided the lowest accuracy compared 
 to the other descriptors . considering 
 that models developed with different descriptor sets had approximately 
 similar performance , we decided to build our consensus models using 
 a simple average of individual models . for each machine learning approach , 

since the 
 individual models were calculated using the bagging approach , the 
 developed consensus models were also validated using the same protocol . 

the model based on the asnn method calculated the lowest rmse error 
 compared to consensus models developed using other machine learning 
 approaches .
the exclusion of the model based on in vitro descriptors did not change the accuracy of the asnn consensus model . 

the model based on a combination of both in silico and in vitro descriptors requires both sets of 
 descriptors .
this limits its application to compounds for which in vitro measurements are present , while the model based 
 exclusively on in silico descriptors can be applied 
 to any new compounds . therefore , we decided to submit the model developed 
 using only in silico descriptors to the challenge .
first , we 
 found that the models developed on the training set have large validation 
 rmse and that the provisional set statistics had a limited value for 
 model selection .
therefore , we followed the model development steps , 
 which were successful in our previous studies .
did the complexity of the final 
 model ( a consensus of several individual 
 models , each of which uses a different descriptor set and is a bootstrap 
 aggregation of multiple neural network submodels ) add any value , or 
 could we get some similar results using a simpler approach ? in order 
 to answer this question , we developed models using just two descriptors : 
 molecular weight and number of carbon atoms using linear regression . 

the rmse of this model on the training set was 1.0  0.04 log 
 unit .
the use of the same descriptors for the bagging approach decreased 
 rmse to 0.97  0.04 log unit for the libsvm method ( table 2 ) .
interestingly , 
 the best model calculated based on the in vitro assay 
 measurements had exactly the same accuracy ( rmse = 0.97  0.04 ) .
the topcoder organizers 
 kindly released information about the experimental values for the n = 143 test set molecules .
it allowed us to provide an 
 additional analysis of the results for this set and to better evaluate 
 the influence of the in vitro descriptors on the 
 prediction accuracy .
the model developed using in vitro descriptors ( see table 2 ) had higher rmse = 1.24 for the test set 
 as compared to that 
 of the consensus model , rmse = 1.08 , based on in silico descriptors .
the extension of the consensus model by inclusion of 
 the model based on in vitro descriptors increased 
 the rmse of the new consensus model to 1.10 log units .
we also explored 
 whether extension of in silico descriptors with in vitro descriptors can provide better prediction accuracy . 

for this study , we developed models using combinations of each descriptor 
 set with in vitro descriptors .
the rmse of the models 
 developed with in silico + in vitro sets were changed 
 in the range of 0.02 to 0.01 log units compared to rmses of 
 models calculated using only in silico descriptors . 

the rmse of the consensus model based on in silico + in vitro sets was 1.09 , i.e. , 0.01 log units higher compared to that of the 
 model based only on in silico descriptors .
the rmse calculated for the test 
 set was significantly higher than that for the training set compounds 
 ( table 2 ) .
an analysis 
 of extended functional groups ( efg ) was 
 done to identify whether both sets contained chemically different 
 compounds .
the efg consists of 583 manually curated functional groups , 
 which provide comprehensive coverage of heterocyclic compounds and 
 are relevant for medicinal chemistry .
the setcompare tool was used to determine statistically significantly 
 overrepresented chemical groups in training and in test sets using 
 hypergeometric distribution .
it was found that hydroxy compounds , 
 amines , saturated six - membered heterocycles containing one heteroatom , 
 etc .
, were overrepresented in the test set , while pnictogens , thiophosphoric 
 acid esters , halogen derivatives , etc . were overrepresented in the 
 training set .
thus , the chemical diversity 
 of molecules in both sets may have contributed to the observed differences 
 in the rmses for the training and test sets .
the 
 efg and setcompare tools were also used to analyze which chemical 
 features contributed to predictions with high errors .
a difference 
 of 1.5 log units between predicted and experimental values was used 
 to identify n = 62 compounds with high prediction 
 errors .
most of these compounds had extreme lel values , i.e. , either 
 low or high values , and only 7 compounds ( 10% ) were within [ 2 , 
 4 ] log(m ) interval as compared to 420 ( 75% ) of compounds that 
 were within this region for the remaining group .
thus , the model had 
 difficulties in predicting highly toxic and nontoxic compounds . according to the topcoder rules
, 
 the submitted results were scored using predictions for all test set 
 compounds . however , of course , some of the compounds from the test 
 set could be outside of the applicability domain ( ad ) and have lower 
 prediction accuracy .
ochem uses standard deviation of the predictions 
 of models , which contribute to the consensus model , as distance to 
 the model ( abbreviated as std - cons ) .
std - cons 
 was found as the best definition of the distance to the model .
ochem defines ad of a model as the value of std - cons that covers 
 95% predictions from the training set . in the previous benchmarking 
 study of 12 definitions of distances to the model applied to 11 models , 
 we found that std - cons provided the best separation of molecules with 
 large and low prediction errors notwithstanding the used model .
therefore , we concluded that ad of models is 
 determined by the composition of the training set of molecules rather 
 than by the used descriptors or machine learning methods . in the current 
 study ,
seven compounds from the test set were outside of the ad of 
 the consensus model .
the consensus model as well as individual submodels 
 calculated significantly higher rmses , which were in the range of 
 2.1 to 2.7 log units , for these seven compounds ( table 1 ) .
this result supports the 
 previous conclusions about the universal nature of the ad of models 
 and good discriminating power of std - cons distance to the model .
it 
 also indicates that taking into consideration the ad of the models 
 is important to avoid predictions with high errors .
it is interesting 
 that four out of seven compounds had lel < 3 , while other three 
 compounds had lel > 5.5 .
thus , the used ad has identified compounds 
 that had experimental toxicity values in the ranges that are difficult 
 to predict .
the final consensus model was based on 10 submodels , which were 
 developed with n = 4036 descriptors ( table 3 ) .
these descriptors were selected 
 from the initial set following the unsupervised filtering procedure . 

we explored whether the performance of this model can be further improved 
 by using a supervised descriptor selection procedure based on neural 
 network pruning of the least sensitive descriptors .
the application 
 of this procedure decreased the numbers of descriptors in 5 to 100-fold 
 ( table 3 ) .
models developed 
 using these descriptors had on average lower training set rmses as 
 compared to those based on descriptors selected by the unsupervised 
 filtering , while the opposite result was calculated for the test set 
 rmses ( table 3 ) .
thus , 
 selection of descriptors optimal for the training set introduced variable 
 selection bias .
indeed , during supervised 
 selection of descriptors we evaluated the performance of models for 
 the training set molecules multiple times .
this resulted in selection 
 of descriptors with improved fit for this set but at the same time 
 decreased the prediction accuracy for test set compounds , which have 
 different chemical diversities .
the neural networks are very efficient 
 methods to work with high - dimensional data and can be also efficiently used without a need of supervised variable 
 selection .
rmse is the root mean squared error 
 calculated for the training ( n = 483 ) and full test 
 set ( n = 143 ) .
the final consensus model was based on 10 submodels , which were 
 developed with n = 4036 descriptors ( table 3 ) .
these descriptors were selected 
 from the initial set following the unsupervised filtering procedure . 

we explored whether the performance of this model can be further improved 
 by using a supervised descriptor selection procedure based on neural 
 network pruning of the least sensitive descriptors .
the application 
 of this procedure decreased the numbers of descriptors in 5 to 100-fold 
 ( table 3 ) .
models developed 
 using these descriptors had on average lower training set rmses as 
 compared to those based on descriptors selected by the unsupervised 
 filtering , while the opposite result was calculated for the test set 
 rmses ( table 3 ) .
thus , 
 selection of descriptors optimal for the training set introduced variable 
 selection bias .
indeed , during supervised 
 selection of descriptors we evaluated the performance of models for 
 the training set molecules multiple times .
this resulted in selection 
 of descriptors with improved fit for this set but at the same time 
 decreased the prediction accuracy for test set compounds , which have 
 different chemical diversities .
the neural networks are very efficient 
 methods to work with high - dimensional data and can be also efficiently used without a need of supervised variable 
 selection .
rmse is the root mean squared error 
 calculated for the training ( n = 483 ) and full test 
 set ( n = 143 ) .
in this study , we highlighted the steps used 
 to develop the rank - i 
 submission for the epa toxcast challenge , which was organized by the 
 topcoder community .
we have shown how to consider limitations of the 
 training and test data sets and that following  occam s 
 razor  principle helps to provide a top - entry to the challenge . 

a similar consensus 
 approach was used to achieve the overall best balanced accuracy for 
 12 end points for another toxcast challenge organized by nih .
the consensus modeling 
 was also successfully used in the cerapp project to identify potential 
 endocrine disruptors .
it is rather 
 surprising that the rank - i model did not involve the in vitro descriptors .
the lel 
 is defined as lowest effect level dose across multiple animal studies . 

this can contribute to considerable differences in the determined 
 quantitative toxicity thresholds due to interspecies variations as 
 well as differences in the experimental protocols .
these factors could 
 contribute to the biological noise of the measured values and make 
 their prediction a difficult task .
the relatively 
 weak performance of this model and all others in the challenge can 
 also point out the limitation of the brute - force machine learning 
 approach to this problem .
the in vitro assay data 
 may need to be treated as more than just a table of numbers , and one 
 will need to incorporate biological knowledge into the structure of 
 the model .
indeed , pharmacokinetic and pharmacodynamics properties 
 of the analyzed molecules could be essential for their toxicity .
thus , 
 we can expect that the use of systems biology methods can contribute 
 to more accurate predictions of the lel .
it should be mentioned that 
 the use of external data was explicitly forbidden for the purpose 
 of the toxcast challenge .
we can not exclude the possibility that some of the 
 currently used in vitro assays could be insufficient 
 for the analyzed end 
 point .
for example , if toxicity is caused by metabolites of the analyzed 
 compound the in vitro assays ignoring metabolic activation 
 may not correctly report toxicity .
currently , it is not clear whether 
 such problems frequently occur , but recent studies suggest that taking 
 into consideration the metabolic activation was an important factor 
 for prioritization of potentially emerging contaminants .
of course , the same problem can also contribute 
 to difficulties with prediction of toxicity based on in silico descriptors . which of these factors contributed to the low 
 accuracy of the model ?
such analysis is beyond the scope of the article 
 and will hopefully be answered in the future with new computational 
 studies by the scientific community .
importantly , the public availability 
 of the test set compounds released in this article will help other 
 users to develop and benchmark new approaches to predict lel and benchmark 
 their results against the rank i model of the toxcast challenge .
moreover , 
 since the model is publicly available and does not use in 
 vitro descriptors , it can be used to predict the lel of new 
 compounds in prospective studies and can be benchmarked using new 
 measurements , which may be available in the future .
we believe publishing 
 models online in a usable and reproducible manner will become an integral 
 part of future computational chemistry . in summary , we have described the protocol for developing 
 the rank - i 
 model of the epa toxcast challenge . the model is based only on in silico descriptors , and we were not able to increase 
 its prediction ability using in vitro measurements 
 in a postmarathon study presented in this article .
the relatively 
 low accuracy of this model indicates high complexity of the lel and 
 suggests that pure brute - force machine - learning approaches may not 
 be sufficient to accurately predict such a complex biological end 
 point .
possibly , systems biology approaches can help to develop better 
 models for the prediction of lel using the available in vitro measurements . at the same time
, we can not exclude the possibility 
 that the currently used in vitro assays may not be 
 sufficient to correctly characterize this end point .
the developed 
 model and used data are publicly available at http://ochem.eu/article/68104 and can be used by interested users to answer these questions as 
 well as to benchmark new ideas , methods , or approaches .
 the 
 toxcast epa challenge was managed by topcoder in spring 2014 . 
 
 the goal of the challenge was to develop a model to predict the lowest 
 effect level ( lel ) concentration based on in vitro measurements and calculated in silico descriptors . 
 
 this article summarizes the computational steps used to develop the 
 rank - i model , which calculated the lowest prediction error for the 
 secret test data set of the challenge . 
 the model was developed using 
 the publicly available online chemical database and modeling environment 
 ( ochem ) , and it is freely available at http://ochem.eu/article/68104 . 
 surprisingly , this model does not use any in vitro measurements . 
 the logic of the decision steps used to develop the 
 model and the reason to skip inclusion of in vitro measurements is described . 
 we also show that inclusion of in vitro assays would not improve the accuracy of the model . 
Introduction
Data
Methods
Results
Development of Models Using Descriptors Optimized with Pruning
Discussion