Many of the models in Table 1 work from the same data source, a compilation of α- and β-CD affinities published by Suzuki in 2001 (Suzuki, 2001) [additionally, the sources that cite a different paper by Katritzky ultimately use this same data, as the Katritzky paper cites Suzuki for observations (Katritzky et al., 2004)]. In addition to Suzuki, we also compiled complexes of α- and β-CD Rekharsky and Inoue and Suzuki (Rekharsky and Inoue, 1998). Complexes of γ-CD, missing from the Suzuki dataset and sparse in the Rekharsky and Inoue data, were collected from Connors (Connors, 1995). Once compiled, the data were cleaned for reliable information, one-to-one cyclodextrin complexes, a temperature of 298 ± 2 K, and a solvent of water with pH 7. To obtain structure-data files (SDFs) of the ligands, the names of the guest molecules were passed through the Chemical Identifier Resolver, a web interface provided by the National Cancer Institute’s Computer-Aided Drug Design Group (NCI/CADD). To handle the data, the R packages tidyverse, data. table, XML, RCurl, and Matrix were used (Bates et al., 2017; Dowle et al., 2017; Wickham, 2017; Duncan Temple Lang and the CRAN T and eam, 2018a; Duncan Temple Lang and the CRAN Team, 2018b).

Dataset splitting was performed using the R package caret (Kuhn and Quinlan, 2018). First, the cleaned data was split between α-, β-, and γ-CD. Structural and activity outliers in each category were removed. Structural outliers were detected using a statistical method relying on standard deviations of molecular descriptors (Roy et al., 2015). For activity outliers, molecules with reported ΔG values greater than 2.5 standard deviations from the mean were removed. Though traditional practice advises classifies outliers as values more than only two standard deviations away, in this case, retaining data points remained a priority and a larger margin was allowed. There were 9, 21, and 11 α-, β-, and γ-CD outliers, respectively. After removal, around 200, 250, and 100 α-, β-, and γ-CD observations remained.

The data was then split into training, testing data and external validation. For each separate cyclodextrin, an external validation set was created from a random 15% subset of the data. To create multiple training and test sets, the remaining modeling data was split with representative resampling of ΔG values into ten different 75:25 train to test data partitions. Though not as advanced as maximum dissimilarity algorithms, this method proved more practical due to the large number of descriptors (over 1,000) generated for each guest molecule. Furthermore, maximum dissimilarity algorithms, when implemented in this instance, had the unfortunate tendency to select highly similar training and test sets, defeating the purpose of creating multiple sets in the first place.

The process of docking is based on two processes: sampling and scoring (Jacob et al., 2012). Sampling refers to the capacity to search an active site on a protein, macromolecule or, in this case, affinity host. This can be performed with distance matrices, matching algorithms or incremental construction, multiple copy simultaneous searching, stochastic methods, or any combination of the aforementioned strategies. Scoring calculates the final binding affinity between the guest and host and can be dependent on force-field, empirical, or knowledge-based calculations. Docking generally involves the use of a host and a guest molecule which can be either rigid or flexible. Three types of conformation exist: rigid-rigid, rigid-flexible and flexible-flexible. In this paper we use AutoDock Vina, a version of AutoDock that uses Monte Carlo stochastic sampling coupled with a force field based scoring function from a resample of a drug like database to derive its weighted parameters. Vina in particular uses a flexible drug guest and a rigid cyclodextrin host, although it allows side chain mobility when docking ligands onto proteins.

The PyRx Virtual Screening Tool provides a variety of services, including molecular energy minimization, docking calculation, and visualization of molecules. PyRx version 0.8 was used here, as further editions require purchase (Dallakyan and J, 2015). Ostensibly, the source code of newer versions of PyRx is freely available, but actually implementing the code requires fairly advanced knowledge of Python, making public usage difficult. To begin, all guest molecules went through energy minimization to determine the most likely atomic configurations. AutoDock Vina, integrated within PyRx, calculated the change in Gibbs free energy (kcal/mol). We tested the effect on the docking process of changes in the search space, search exhaustiveness, and scoring force field type.

The open source software PaDEL-Descriptor calculated over 1,000 descriptors for the remaining molecules, including fingerprints, structural details, and physical properties (Yap, 2011). Additionally, PaDEL-Descriptor removed salts and minimized the energy of inputted files using an MM2 force field. To improve model interpretability, more abstract predictors, such as those related to eigenvalues for molecular matrices or autocorrelation, were excluded from calculation. The elimination of these descriptors did not produce any noticeable effect on final model accuracy and made feature selection less resource intensive.

Recursive feature elimination (RFE), implemented with caret, was used to subset the predictors used for model-building (Kuhn 2018). Using this method, a random forest model is created using all available descriptors. Once trained, the relative importances of the predictors are calculated and differently sized subsets (defined by the user) of variables are selected to create and evaluate new models. The best combination of predictors is then returned by the model. RFE was performed on each of the ten train-test splits. The predictors determined to be useful for all folds were saved and used for tuning and training the models. This resulted in 13 variables for ⍺-CD, 16 variables for β-CD, and 39 variables for γ-CD.

We investigated the accuracy of several models that appeared in previous attempts at cyclodextrin QSARS (Table 1), including Cubist models, generalized linear models (GLM or GLMNet), random forests, partial least squares models, and support vector machines. Additionally, two QSAR methods not previously published for cyclodextrin—multivariate adaptive regression splines (MARS) and gradient-boosted models—were created and evaluated. Model building was accomplished with R-packages Cubist, glmnet, randomForest, pls, e1071, earth, and gbm, respectively (Cutler and Wiener, 2015; Mevik and Liland, 2016; Friedman et al., 2017; Kuhn et al., 2017; Meyer et al., 2017).

Cross-validation was used to determine ideal tuning parameters for each QSAR. For faster QSARs—such as generalized linear models (GLM) and partial least squares (PLS)—tuning was performed using 10-fold cross validation. For more resource-intensive models or models with large parameter spaces—such as random forests, Cubist and support vector machines (SVM)– only five folds were used. Optimized models, QSARs built with the tuned parameters and trained on the entire training set, were used to predict the test for each combination of test and training set. Further fine tuning was also performed at this step. The model that produced the lowest root-mean square error (RMSE) and highest R ² (or an otherwise most ideal combination) on all the test sets became the final model, i.e., the model saved for future use. Furthermore, the models were evaluated according to Tropsha and Golbraikh standards for QSARs (Golbraikh and Tropsha, 2002). Although R ² and RMSE can be useful for generalizing predictive capacity, they may be misleading in certain cases, necessitating stricter additional standards of evaluation. As an additional test of reproducibility, the final models were used in ensemble to predict the values of the external validation set. Because this dataset was withheld from the entire model training process, the external validation set served to simulate model performance on new data.

Applicability domain describes the range of molecules where the model can be expected to generate reliable predictions. A new molecule outside of the applicability domain is structurally quite different from the set of data the model was trained on, and thus a prediction will rely on extrapolation and may not be accurate. The applicability domain of the models was determined with the same method used to detect outliers when cleaning the dataset (Roy et al., 2015).

Y-randomization was used to further verify the significance of the results. Many advanced QSAR methods are powerful enough to model data off of noise, so y-randomization ensures that the modelling process produces results significantly more accurate than what could be obtained by chance. Randomization can be achieved by permutation (randomly changing the positions of observed values) or random number generation (replacing observed values with completely new data). Different combinations of permutation and/or random generation yields five different modes of y-randomization to investigate: 1) original ΔG values vs randomly generated descriptors, 2) permuted ΔG vs original descriptors, 3) random ΔG vs original descriptors, 4) random ΔG vs random descriptors, and 5) permuted ΔG vs random descriptors. (Combinations including permutation of descriptors are not included because the large number of predictors in QSARs renders the effects of such a process virtually indistinguishable from random number generation.) However, because the y-randomization process is extremely resource-intensive (as each mode requires that several randomized iterations undergo the modeling process), only mode 1, the most common interpretation of y-randomization, was investigated (Rücker et al., 2007). The observed ΔG values were randomly assigned to guest molecules, and the entire model refitting process was re-done, from feature selection to external validation.

Using R’s “shiny” package, most of the process of running the QSAR could be implemented in a web app. The app was split into three main pages: Download, Upload, and Explore. “Download” accesses Chemical Identifier Resolver and obtain SDFs. The page also draws the obtained molecule using the package “ChemmineR,” allowing the user to check that the SDF is accurate. “Upload” implements the QSARs after the user provides the app with a CSV of the descriptors from PaDEL-descriptor. After calculating the affinity and analyzing the applicability domain of the molecules, the user is provided with both a graph and a table of the results. The third page “Explore,” stores the results of using the ensemble on FDA-approved drugs, as obtained from the annual publication “Orange Book: Approved Drug Products with Therapeutic Equivalence Evaluations.” (Food and Drug Administration, 2019)

Partial differential equations that model drug diffusion were solved using the R package deSolve (Soetaert et al., 2018) using the method of lines as previously done by Fu et al. (2011). In the model, the release media was assumed to be water and the delivery system was assumed to be flat, thin circular cyclodextrin disc. The boundary condition between the polymer and the media was approached as described by Wang and von Recum (2011). Diffusivities of drug molecules were calculated from molecular weight and viscosity using a modified Stokes-Einstein-Sutherland equation, as done by Vulic et al. (2015).

When predicting on the entire cleaned dataset (all modeling data, which includes the training, testing, and external validation set), AutoDock Vina yielded an R ² of 0.18 and a RMSE of 5.00 kJ/mol (Figure 2). Of the 547 cleaned complexes, docking provided calculations for 458, failing to provide data on 89 complexes. Adjusting settings in Vina, such as the minimization algorithm size of the steps in the calculation, did not yield significant differences in accuracy. In comparison, the affinity of only around 40 cleaned molecules could not be obtained by the ensemble QSAR. In these cases, the withheld molecules were determined to be outliers, and the actual QSAR model could still be used to predict a value.

The results of predicting on the test data for each QSAR and cyclodextrin type are markedly higher than docking (with the exception of γ-CD), reaching an R ² of around 0.5 to 0.7, as seen in Table 2 and Figure 3. The reported R ² for each QSAR type is calculated from an average of the performance of the model on all test splits. Additionally, Table 2 contains information on verification of all the QSAR types (3-VII Validation methods). Of the types investigated, only PLS and GLMNet failed to pass the salvo of verification criteria, both falling short of attaining an R ² of 0.6.

In terms of reliability, most models were able to handle the available data well, providing calculations for all provided molecules. Only the Random Forest and Cubist models failed to calculate the affinity of some molecules, possibly due to being based around decision-trees. The algorithm underlying both models attempts to draw predictions by categorizing entries based on their features. If they encounter a molecules entirely different from the data they trained on, the models may fail to create a prediction. Advantageously for our approach, the failure to calculate some values becomes less important where models are combined in an ensemble where the final prediction is averaged over many models.

The results of ensemble prediction (averaging the results of many different QSARs) can be seen in Figure 4. While ⍺- and β-CD models managed to reach moderately high predictive performance metrics, unfortunately, all γ-CD models lacked useable predictive power.

The models passing the verification in Table 2 were further verified using y-randomization. To ensure accuracy was not the result of the models building off of noise, 25 different permutations of ΔG values were created. All Q² values of the models created from the original data were calculated to lie well outside 3 standard deviations of the mean Q² of the randomized data. Additionally, the R ² values of the ensemble QSARs were significantly greater than the R ² values obtained from the ensemble models created from permuted data (means of 0.021 and 0.027 and standard deviations of 0.011 and 0.016 for ⍺- and β-CD, respectively).

Interpretability of a model provides a rough check if a model is calculating off of random noise or if the model is drawing logical calculations from physical properties to molecular behavior. Each model, due to differences in statistical algorithms and approaches, has differing levels of interpretability. GLM, being similar to linear models, have easily accessible coefficients associated with each predictor, so the relative impact of each factor can be compared with reasonable confidence. Cubist models, on the other hand, tend to be difficult to interpret as variables are processed through multiple levels of decision trees.

Evaluation for the relative importance of variables are shown in Figure 5. Random forest was the only QSAR type with a pre-packaged importance function for variable analysis. PLS variables were analyzed using a function obtainable from Mevik et al. (2007). Max Kuhn’s caret package was used to evaluate GLMNet, the two SVM kernels, and Cubist. Unfortunately, caret was unable to process the final models for GLMNet and SVM, and the reported variable importance values were actually derived from models created within caret’s “train” function, and are thus slightly different from the models saved in the ensemble. To determine importance, the “train” function removes a variable, rebuilds the model, and analyzes the effect on accuracy. The more important a variable, the larger the drop in accuracy. After each variable has been tested, the function can then rank the importance of the descriptors.

For β-CD, XLogP, a measure of lipophilicity, appears to be important for all models, consistent with how the structure of cyclodextrin allows for easier complexation with small hydrophobic drugs. The same reasoning can be extended to LipoAffinityIndex and MLogP, additional approaches to quantifying lipophilicity. The number of carbons, nC, is also consistently important, possibly due to a relationship with molecule size. WTPT-2 is the PaDEL weighted path descriptor divided by the number of atoms, and also may be important due to encoding information on molecular size.

However, not all the variables can be linked to set chemical properties. SpMax and SpMin relates to eigenvalues of a modified connectivity matrix, a numerical representation of atomic and molecular bonds, and may not be associated with any interpretable physical property (the same analysis can also be used for GATS predictors). To aid interpretability, building a model with predictors easily attributed to physical or chemical properties may be advised. The extent to which interpretability should trade off with accuracy remains in question. Our findings go in accordance to those in the general literature were lipophilicity tends to highly influence model output.

After collecting a list of FDA-approved drugs and drug combinations from the Orange Book, an annual publication listing all approved pharmaceuticals, the names were cleaned for individual active compounds. In total, 1,401 unique molecules could be extracted. Of these, 1,116 could be downloaded from Cactus and 1,031 could be processed by PaDEL. Many of the molecules that could not be analyzed by PaDEL would have proven impractical for cyclodextrin delivery, such as simple ionic salts (e.g., potassium chloride), or large molecules made of more than 100 atoms. Running the remaining guests through applicability domain analysis yielded 638 molecules, 45.5% of the original set. While less than half of FDA-approved drugs could pass through the model, the 600 available guests spans a wide range of properties and uses, allowing the page to be useful for candidate selection (Figure 6).

Though the app could be uploaded online through shinyapps. io, server time limitations on the account hosting the app make it impractical for usage by a large number of individuals simultaneously. In order to run the app for more than a few hours, such as with screening a large dataset of molecules, the code would have to be downloaded through GitHub. In addition, the user would need to download the R libraries and the IDE RStudio, potentially negating the goal of creating an accessible, intuitive interface. The “Explore” page partially alleviates this obstacle, as it allows the user to perform a quick search of a pre-predicted affinity rather than spend time downloading the structure file, launching PaDEL, and running the QSAR.

To demonstrate the extensibility of the shiny app and its value in the design of drug delivery strategies the results of the QSAR predictions can be fed into a mechanistic model of drug delivery (Figure 7). The results demonstrate the ranges of values expected from the strongest affinity binding predictions and from the weakest. As expected, strong predictions produce much slower release profiles and weak predictions produce faster release profiles. Conservation of mass was verified as the sum of all mass within the system from the polymer and media compartment across all times as a test of the implementation. Notably, the implementation in R required a modification of the method of lines for appropriate modeling of the polymer to liquid media interface. (Linge and Langtangen, 2016). Future efforts in creating new modules could explore substructure searching to identify alternative strategies for weak binders or drugs that demonstrate unsuitable release profiles. It is expected that not all drugs will follow this simplistic model of drug release. For those cases our lab has built a whole suite of approaches to alter elution rates such as a wide range of formulations, supramolecular interactions, Schiff-base formation and multi-arm PEG substitutions.