Predicting complexation performance between cyclodextrins and guest molecules by integrated machine learning and molecular modeling techniques.

Most pharmaceutical formulation developments are complex and ideal formulations are generally obtained after extensive experimentation. Machine learning is increasingly advancing many aspects in modern society and has achieved significant success in multiple subjects. Current research demonstrated that machine learning can be adopted to build up high-accurate predictive models in drugs/cyclodextrins (CDs) systems. Molecular descriptors of compounds and experimental conditions were employed as inputs, while complexation free energy as outputs. Results showed that the light gradient boosting machine provided significantly improved predictive performance over random forest and deep learning. The mean absolute error was 1.38 kJ/mol and squared correlation coefficient was 0.86. The evaluation of relative importance of molecular descriptors further demonstrated the key factors affecting molecular interactions in drugs/CD systems. In the specific ketoprofen-CD systems, machine learning model showed better predictive performance than molecular modeling calculation, while molecular simulation could provide structural, dynamic and energetic information. The integration of machine learning and molecular simulation could produce synergistic effect for interpreting and predicting pharmaceutical formulations. In conclusion, the developed predictive models were able to quickly and accurately predict the solubilizing capacity of CD systems. Current research has taken an important step toward the application of machine learning in pharmaceutical formulation design.

Introduction

In pharmaceutical industry, the formulation developments intend to optimize the excipient composition and process parameters to acquire optimn class="Chemical">al formulations with satisfactory physicochemical and biological characteristics. However, the process of current formulation development highly depends on the trial-and-error experimentation by the experiences of individual formulation scientists. Thus, the multivariable and complex nature of formulation design seriously hinders the efficiency and success rate of product development2, 3. Moreover, the successful application of high-throughput experimentation in the field of drug discovery in the past decade has contributed to significant increase of active compounds, which put forward constant pressure for the productivity of formulation development4, 5. Although the statistical experimental designs and optimization techniques, such as orthogonal design or response surface methodology, to some extent simplify the formulation screening process, the extensive experimentation is still needed due to the unpredictable formulation performance. It is necessary to seek for a high efficiency and accuracy methodology for the prediction and virtual screening of pharmaceutical formulations.

Currently molecular simulation-based calculations in the formulation development provide an n class="Chemical">alternative technique to the traditional trial-and-error experiments. Numerous molecular modeling methods, like molecular dynamics (MD), docking, Monte Carlo and even quantum mechanics, have been extensively applied3, 7, 8. The molecular modeling approaches are highly capable of the explanation for the molecular mechanism for specific formulations and delivery systems, but have the limited ability to distinguish key index in some specific formulations. Another existed significant limitation is that the modeling systems with limited molecule number may not represent the real experiments due to the limited computing capability. These methods have difficulties in making balance between time consumption and accuracy. In any case, from the wide studies performed to date, it is highly necessary and urgent to build up predictive models which allow a high-accuracy prediction in the pharmaceutical formulations and contributing factors analysis.

Machine learning is able to learn from input data, automate analyticn class="Chemical">al model building and even update outputs as new data becomes available by computer algorithms9, 10. The latest machine learning algorithm—deep learning garnered significant attention, which is representation-learning methods with multiple levels of representation11, 12. These algorithms are able to make high-accuracy predictions for the classification and regression tasks, which lead to extensive application in the bioinformatics and computational biology. However, the application of machine learning in pharmaceutical research especially formulations development, is still very limited. Given commercial confidentiality and long research periods, the pharmaceutical datasets start out too small to guarantee implementation of these algorithms. Despite many challenges, two recent studies have advanced the application of machine learning in the pharmaceutical field14, 15, 16, 17. The remote liposome study has used machine learning approaches to distinguish the drug loading efficiency and to predict the drug loading on the 366 formulation dataset. Another attempt was the combination of machine learning and nondestructive ultrasonic to develop a predictive model for the tablet breaking force and disintegration time for advancing the quality-by-design paradigm. Because of the complex nature of formulation design process and multidimensionality of chemical structures, it is really necessary to explore whether machine learning could be further expanded into formulation prediction based on available data.

Current research selected the cyclodextrin (n class="Chemical">CD) complexation systems as the model system for the predictive model development because of its extensive application in pharmaceutical field, relatively controllable variables and specific evaluation index. CDs can accommodate varieties of non-polar molecules into the inner cavity to form non-covalent host-guest inclusion complexes, which has been successfully used in pharmaceutical formulations. The guest/CDs binding process is not fixed or permanent, but dynamic and reversible. Complexation free energy or equilibrium constant is a key index to evaluate the strength of host–guest complexation and disassociation/association stability of the complex, which can decide the usage of a specific drug–CDs inclusion complex in the formulation design20, 21. Many experimental techniques have been developed to measure the complexation free energy including phase solubility measurements, calorimetric titration, nuclear magnetic resonance chemical shift, freezing point depression, pH-metric methods and so on, among which phase solubility method described by Higuchi and Connors is the most widely used22, 23, 24, 25. However, these experimental measurements often bring serious challenges for the pharmaceutical scientists because of the limited synthetic amount of new compounds and the limited API's solubility at the early stage of drug discovery.

This study was proposed to develop high-accurate predictive models of complexation free energy between CDs and guest molecules based on a dataset of 3000 date points by three machine learning techniques (e.g., light gradient boosting machine [LightGBM], random forest [RF] and deep learning [DL]). The relative contributions of molecular descriptors were ann class="Chemical">alyzed to infer underlying molecular interactions in the drug–CD complexation systems. Furthermore, ketoprofen (KTP) was used as a model drug and the KTP–CDs inclusion complexes were studied by the experimental methods and the MD simulations. The comparison of complexation free energy by experimental determination, simulation-based calculation and machine learning model was further conducted.

Materials and methods

Dataset preparation and distribution

The data for the complexation free energy between CDs and a diverse set of organic molecules were collected from published literatures from the year 1990 to 2018. The literature data were considered acceptable if they belonged to binary n class="Chemical">CD complexation system with the complexation free energy or the equilibrium constant by phase solubility study. When the experimental pH values were not clear or reported as aqueous solution, they were considered to be neural (pH = 7). The temperatures were assumed to be room temperature (T = 25 °C) if the specific temperature was not noted or reported. The available dataset was further processed to eliminate duplicates, to eliminate entries for the same compounds with different complexation free energies in the same experimental conditions, as well as to remove experiments under uncommonly used conditions (e.g., partially or completely organic solvent). The final dataset with 3000 formulations for the predictive model building contained 1320 guest molecules and 8 CDs, of which the proportions of each CD were shown in Fig. 1A. The complexation free energies between guest molecules and CDs were distributed over the range of 0 to −40 kJ/mol with the mean of −15.12 kJ/mol, as shown in Fig. 1B. The guest molecule set contained a large amount of structurally diverse organic compounds, such as phenols, ketones, alcohols, steroids, nitriles, barbitals, and hydrochloride.

Figure 1

Relative distribution of data points in the full dataset: (A) The percentage of experimental data for each cyclodextrin. (B) The distribution of complexation free energy between cyclodextrins and guest molecules. (C) The molecular weight distribution of guest molecules. (D) The XlogP3 distribution of guest molecules. The full dataset was comprised of 3000 data points with 1320 different structures. Hp-β-CD, 2-hydroxypropyl-β-cyclodextrin; RMCD, randomly methylated β-cyclodextrin; TMCD, (2,3,6-tri-O-methyl)-β-cyclodextrin; DMCD, (2,6-di-O-methyl)-β-cyclodextrin; SBE-β-CD, sulfobutylether-β-cyclodextrin.

Relative distribution of data points in the full dataset: (A) The percentage of experimental data for each n class="Chemical">cyclodextrin. (B) The distribution of complexation free energy between cyclodextrins and guest molecules. (C) The molecular weight distribution of guest molecules. (D) The XlogP3 distribution of guest molecules. The full dataset was comprised of 3000 data points with 1320 different structures. Hp-β-CD, 2-hydroxypropyl-β-cyclodextrin; RMCD, randomly methylated β-cyclodextrin; TMCD, (2,3,6-tri-O-methyl)-β-cyclodextrin; DMCD, (2,6-di-O-methyl)-β-cyclodextrin; SBE-β-CD, sulfobutylether-β-cyclodextrin.

Molecular descriptors and dataset slitting

Both the guest molecules and CDs were characterized with an appropriate set of molecular descriptors. The molecular descriptors used in this study were computed using n class="Chemical">ALOGPS 2.1 program (VCCLAB, Virtual Computational Chemistry Laboratory, http://www.vcclab.org) and ChemAxon program (https://chemaxon.com/company). The total of 17 and 22 molecular descriptors were used to characterize the guest molecules and CDs, respectively. These calculated molecular descriptors were then combined with the corresponding experimental conditions (described as above, pH value and temperature) to assemble a hybrid set of molecular descriptors for each entry in the dataset.

Relative distributions of data points used in the dataset were shown in Fig. 1, which shows that the percentages of data points for each CD in the dataset vary greatly. Thus, the dataset was split by the method of stratified random sampling, which is a method of sampling that random samples extracted from each n class="Chemical">CD group was proportional to the size of the original group. The dataset was randomly split into training set (80%), test set (10%) and validation set (10%), among which the data points of each CD kept the same percentage as that in the original dataset. It guaranteed that the composition of subsets was representative and avoided the unbalance in these three datasets. The test set is completely independent from the training set and validation set. Running a training set and validation set through an algorithm teaches the model how to weigh different features, adjust their parameters according to their likelihood of minimizing errors in the results, while the test set serves to evaluate the generalization ability of the predictive model.

Machine learning methods

Currently there have been developed many predicting methods in machine learning. This research compared the predictive models from three widely used algorithms—LightGBM, RF and DL. LightGBM was the first time to be used in the pharmaceuticn class="Chemical">al system.

LightGBM uses an open source gradient boosting decision tree (GBDT) algorithm by Microsoft (http://lightgbm.apachecn.org/cn/latest/). LightGBM is a fast, distributed and high-performance gradient boosting framework with tree-based learning n class="Chemical">algorithm, which has been extensively used for both classification and regression tasks. It uses histogram-based algorithms, which buckets continuous feature values into discrete bins. Compared with other decision tree algorithms, it can reduce the calculation cost of the gain for each split and speed up training. Furthermore, LightGBM grows tree leaf-wise and vertically, while most decision tree learning algorithms grow trees horizontally by level-wise. In fact, LightGBM tends to obtain lower loss than level-wise algorithms by choosing the leaf with max delta loss to grow when holding leaf fixed. Leaf-wise may cause over-fitting when the dataset is small. The max_depth in LightGBM is an important parameter in handling with over-fitting. In current research, the max_depth for tree model was limited by set max_depth to 4. The complexity of the tree model was controlled by the parameter of num_leaves, which was set to 9. N_estimators was used to control the number of boosted trees to fit, it was set to 800. The number of boosting iterations was set to 1000. When the LightGBM was used on the sparse datasets, each parameter had a small adjustment.

Random forest uses sklearn.ensemble RandomForestRegressor (http://scikitlearn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html). RF is an ensemble method for classification or regression tasks by using bootstrap samples of the training data and random feature selection in tree construction. Each bootstrap samples grow a tree and the best split at each node is defined among a randomly selected subset of mtry descriptors. When the data is limited, random forest estimates the ensemble prediction performance by performing a type of cross-validation in parn class="Chemical">allel with the bootstrap procedure with the out-of-bag data. The RF prediction of new data is made by aggregating the prediction of the ensemble. In addition to built-in estimation of prediction accuracy, RF can also calculate the importance of each molecular descriptor based on the out-of-bag predictions with each descriptor permuted22, 30. During the data training, n_estimators controlled the number of trees in the forest was set to 300, while max_depth controlled the maximum depth of the tree was set to 14. When looking for the best split, the number of features (max_features) to consider was set to 32. Min_sample_split, represented the minimum number of samples required to split an internal node, was set to 2. Min_samples_leaf, controlled the minimum number of samples required to be at a leaf node, was set to 3.

Deep Learning using Keras deep learning library in python (http://keras.io/). DL is a form of machine learning method, which comprises of multiple processing layers to learn representations of data with multiple levels of abstraction. Unlike simple artificial neurn class="Chemical">al networks, DL has multiple hidden layers with different weights, and then passes the signals successively deeper in the network until the output layer. During training procedure, the deep network trains the first layer as an auto-encoder with the training sets as the input until the average of the objective function stops decreasing, and then trains the second layer as an auto-encoder taking the first layer's output as the input. The process is iterated for the desired number of layers and finally the output of the last layer as the input for the prediction layer. In this research, the neural network included three hidden layers—the first layer containing 512 neurons, the second layer containing 256 neurons, and the third layer containing 64 neurons. The first three hidden layers used Rectified Linear Units (ReLU) as the activate function, while the last layer used a linear function. The number of samples that going to be propagated through the network was controlled by the batch_size, which was set to 16. The parameter of epoch indicated the times of the entire dataset passed through the neural network, which was set to 3600. During the training process, it would shuffle data for each epoch to have different data for each batch.

Model selection and comparison

The predictive performance had been characterized by three statistics, which were the mean absolute error (MAE), the root means square error (RMSE) and the squared correlation coefficient (R2). They were defined as:where and were the experimentally observed and theoreticn class="Chemical">ally predicted complexation free energy between molecule i and corresponding CD, respectively. yobs,mean was the mean of experimental value and n was the number of data points.

Experimental study of KTP–CDs inclusion complexes

Materials and reagents

α-, β-, and γ-CD were purchased from J&K Scientific Co., Ltd. (Beijing, China). 2-Hydroxypropyl-β-n class="Chemical">cyclodextrin (Hp-β-CD; FW = 1541.54 g/mol, purity > 99%), dimethyl-β-cyclodextrin (DMe-β-CD, FW = 1331 g/mol, purity > 99%), and sulfobutylether-β-cyclodextrin (SBE-β-CD; FW = 2242.01 g/mol, purity > 99%) were purchased from Shanghai Aladdin Bio-Chem Technology Co., Ltd. (Shanghai, China). Ketoprofen (KTP, purity > 99%) was purchased from Dongkangyuan Technology Co., Ltd. (Wuhan, China). All other chemicals in the study were of analytical grade.

Phase solubility study was performed to determine the binding constant (K) according to the methods described by Higuchi and Connors. Excess amounts of KTP were added to 3 mL of various concentrations of aqueous n class="Chemical">CDs solutions. The concentration of three parent CDs were between the ranges of 0 to 15 mmol/L, while the concentration of the other three CD derivatives were 0 to 9 mmol/L. The mixtures were shaken at 37 ± 0.5 °C for two days until equilibrium. Afterward, the samples were filtered through a 0.45 μm Millipore filter and then assayed by UV spectrophotometer at 260 nm. The binding constant (K, L/mol) of the KTP–CD complex was calculated from the phase solubility diagrams using Eq. (4) and the complexation free energy (ΔG, kJ/mol) could be further calculated according to Eq. (5).where S0 was the intrinsic solubility of KTP in distilled water, and the slope was obtained from the phase solubility curve by linear regression between KTP and CD concentration. R is the gas constant, and T is Kelvin temperature.

Molecular dynamic simulation of KTP–CDs inclusion complexes

The molecular dynamic (MD) simulations of KTP–n class="Chemical">CDs inclusion complexes were carried out by using AMBER14 software package and Generalized Amber Force Field (GAFF). The initial structures were generated by docking KTP with various CDs by AutoDock Tools package and then these structures were loaded into LEAP model to construct the solvated systems by using the TIP3P water box with 20 Å radius. All systems were equilibrated in the NPT ensembles at T = 310 K and P = 1 bar, and then ran MD simulation for 100 ns with a time step of 2 fs. The molecular mechanics with a Possion–Boltzmann/surface area solvent (MM-PBSA) method was used to calculate the binding affinity “ΔGbinding”. The binding entropy (ΔS) was calculated by normal mode analysis using the ptraj program in AMBER Tools. The binding enthalpy (ΔE) comprising of electrostatic energy (ELE), van der Waals (VDW), non-polar and polar contributions to solvation (PBSOL), were calculated for the complexes according to Eq. (6). The binding free energy (ΔG) was further calculated according to Eq. (7).

Results and discussion

Predictive performance of three machine learning methods

With these 3000 data points available in hand, the predictive accuracies of three machine learning model were evaluated by using 80% of the data as training set and 10% of the data as vn class="Chemical">alidation set to predict the remaining 10% dataset (as test set). Fig. 2 showed the scatter plots of the predicted vs. observed complexation free energy on full dataset by three machine learning methods and the detailed statistical performance of the predictive models for CDs complexation free energy were presented in Table 1. For the LightGBM model, the MSE and RMSE of the training set was 1.41 and 1.94 kJ/mol, respectively, with a coefficient (R2) value of 0.86. Generally, the predictive accuracy of the built model highly relied on the similarity degrees of the molecules to be predicted and the molecules in the training set. In order to measure the built model's generalizability, an independent and new test set was performed to calculate. While applying the LightGBM model for calculations on the test set, MAE, RMSE and R2 were 1.38 kJ/mol, 1.83 kJ/mol and 0.86, respectively. Thus, the built LightGBM model showed reasonable statistical criteria for new dataset. Turning to random forest model, it did not show overfitting and also had achieved good predictive performance on both the training set and test set. However, the predictive performance of the random forest model did not exceed that of the LightGBM model. The statistical performance of the deep learning predictive models showed the MAE for the training set was 3.34 kJ/mol with a R2 value of 0.76, and the MAE for the testing set was 3.36 kJ/mol with a R2 value of 0.62. Compared with LightGBM model and random forest model, deep learning model showed a relatively poor predictive performance.

Figure 2

Scatter plot of predicted vs. observed complexation free energy on full dataset by three machine learning methods. The predicted results of (A) LightGBM model, (B) random forest model, and (C) deep learning model. For all the models, an 80/10/10 split of training, validation and test set was used to measure the predictive performance.

Table 1

Statistical performance of the predictive models for CDs complexation free energy.

Method	Statistical parameter	Training set	Validation set	Test set
LightGBM	R2	0.86	0.87	0.86
	MAE (kJ/mol)	1.41	1.28	1.38
	RMSE (kJ/mol)	1.94	1.63	1.83
RF	R2	0.83	0.83	0.81
	MAE (kJ/mol)	1.54	1.39	1.54
	RMSE (kJ/mol)	2.17	1.89	2.11
DL	R2	0.76	0.63	0.62
	MAE (kJ/mol)	2.53	2.41	2.56
	RMSE (kJ/mol)	3.34	3.12	3.36

RF, random forest; DL, deep learning.

Scatter plot of predicted vs. observed complexation free energy on full dataset by three machine learning methods. The predicted results of (A) LightGBM model, (B) random forest model, and (C) deep learning model. For all the models, an 80/10/10 split of training, vn class="Chemical">alidation and test set was used to measure the predictive performance.

Statistical performance of the predictive models for n class="Chemical">CDs complexation free energy.

Although three predictive models shown high prediction accuracy both in the training and test set, there still were large deviations for some data points between the observed vn class="Chemical">alues and predicted values in the test set, especially when the value of complexation free energy was less than −20 kJ/mol, as shown in Fig. 2. Significant proportion of the deviations may be caused by the value distribution of complexation free energy in the original dataset, as show in Fig. 1B. The value of complexation free energy in the original data set between −5 and −20 kJ/mol accounted for more than 90% data points in the full dataset. Therefore, the imbalance distribution of dataset had a serious influence on the model construction and the predictive accuracy. In addition, the size of the dataset also had a great impact on the model construction and its predictive accuracy, which had been verified by the following research on the sparse data.

By comparing with the statistical performance of the predictive models on the training and test set, the LightGBM model can deliver high predictive accuracy without succumbing to overfitting. Moreover, LightGBM has the unique advantages with a faster training speed, higher efficiency and better accuracy than any other boosting algorithm because its histogram algorithm can produce much more complex trees by leaf-wise split approach rather than a level-wise approach. But the biggest challenge for the LightGBM algorithm is that it is sensitive to overfitting, especially on the small dataset. Therefore, the parameters of the LightGBM model should be carefully tested, such as the maximum depth of tree, the minimum number of the records for a leaf, and the fraction of data to be used for each iteration, etc.

Random forest has become one of the most widely used machine learning algorithms in the pharmaceuticn class="Chemical">al researches because it can be used for both classification and regression tasks. In addition, random forest also has another two advantages, which attract the pharmaceutical scientists to bring this algorithm in their research. One advantage is that random forest can randomly sample the data and construct decision trees and aggregate many decision tress to limit overfitting, which is one of the biggest problem in machine learning. The other advantage is that random forest can rank the relative importance of each feature on the prediction. Recently, Ahneman et al. used various machine learning techniques to predict the reaction performance in C–N cross-coupling. The results showed that the random forest algorithm provided better predictive performance than the other five techniques. The random forest algorithm was also used to rank the important descriptors and the relative experiments verified its effectiveness.

Deep learning uses multi-layered artificial neurn class="Chemical">al network to characterize the input features at different levels and then optimizes the model performance through various training skills. In fact, deep learning algorithm can not only build the predictive models, but also have the potential to automatically learn representations from data without introducing hand-coded rules or human domain knowledge. Therefore, compared with simple artificial neural networks and other machine learning techniques, deep learning can learn directly from raw data and be more useful for more complex problems. However, many parameters in deep learning algorithms have to be tuned and large amount of data is needed to come up with somewhat generalizable models. Generally speaking, the predictive accuracy of the models increases with increased data. In the pharmaceutical formulation research, however, the available datasets are quite small, which greatly limit further application of deep learning. Our group compared the predictive performance of single-layer neural network model and deep learning model for predicting the disintegration time of oral disintegrating tablets. The raw data comprised of 145 formulations. Due to small dataset, the accuracy of the prediction was only 80% on the test set. Although deep learning has strong predictive potential, the challenges for pharmaceutical researchers is the limited data in the formulation prediction. Meanwhile, the black-box of deep learning models cause many criticisms because it is difficult to provide the explanation between input variables and the built model. Thus, there is still lots of work to successfully introduce deep learning into the pharmaceutical formulation research.

As the results shown in Table 1, all three models of LightGBM, random forest and deep learning showed high predictive accuracy for the complexation free energy between n class="Chemical">CDs and guest molecules. Based on the comparison of the MAE and R2, the predictive performance of LightGBM model was the best, followed by the random forest model, and then deep learning. Meanwhile, the LightGBM algorithm had the highest efficiency as well as the fastest training speed than another two algorithms. Furthermore, LightGBM algorithm can also directly calculate the feature importance, which helped to guide theoretical analysis. Thus, the LightGBM model was the optimal model for the following research.

Predictive performance of LightGBM model on sparse data

Machine learning is a sub-field of artificial intelligence and data science for new algorithms and the predictions. It is generally accepted that gathering as much and high-quality data as possible is the key to ensure high predictability of machine learning algorithms. However, it is really difficult to obtain a large quantity of formulation data in pharmaceutics because of long research cycle and high cost in formulation development. Therefore, it is significantly important to explore the exact relationship between the size of dataset and the predictive accuracy of a machine learning algorithm in the formulation prediction.

Current research has investigated the predictive performance of the LightGBM model with sparse data, varying from 3000, 2500, 2000, 1500, 1000 to 500. Each dataset was split by the method of stratified random sampling as described above. Predictive performance on sparse data by the method of LightGBM for the training and test set was shown in Fig. 3. In the training set, the value of MAE did not change significantly with the decrease of the dataset, fluctuating between 1.3 and 1.65 kJ/mol, and the R2 vn class="Chemical">alue had small fluctuation between 0.8 and 0.9. For the test set, however, it was clearly demonstrated that a gradual erosion in the predictive accuracy occurred from 3000 to 500. The value of MAE increased dramatically with the decrease of dataset, from 1.38 to 2.28 kJ/mol, while the value of R2 also decreased significantly with the decrease of dataset, from 0.86 to 0.58. Results showed that the LightGBM predictive model still had good performance on the training set even when the data points in dataset decreased to 500, which demonstrated that the molecular descriptors used in the current research could exactly captured the properties of input variables. But the predictive performance on the test set became significantly worse with the smaller dataset size, which showed that the LightGBM model had a poor generalization ability on small dataset.

Figure 3

Training and test set performance of the LightGBM model with sparse data. The data points were gradually reduced from 3000 to 500 and the smaller datasets were randomly extracted from the full datasets. A gradual increase of MAE and a significant drop of R2 occurred with the decrease of dataset size.

Training and test set performance of the LightGBM model with sparse data. The data points were gradually reduced from 3000 to 500 and the smn class="Chemical">aller datasets were randomly extracted from the full datasets. A gradual increase of MAE and a significant drop of R2 occurred with the decrease of dataset size.

In fact, machine learning tends to encounter predictive limitation when substantially different guest molecules are used in the testing set. Even though some predictive models for complexation free energy between n class="Chemical">CDs and guest molecules had been built in previously published papers, these models were limited to the dataset with no more than 300 data point20, 34. In current research, the full dataset with 3000 data points had about 1320 different chemical structures, which highly contributed to the high-accuracy predictive ability. Therefore, one method for maximizing the extrapolative ability of a model with limited formulation data maybe spread the training data across its interesting chemical space.

On the other hand, higher quality data contributed to higher predictive accuracy. In pharmaceuticn class="Chemical">al industry, formulation data can be gathered from the published literature and patents, pharmaceutical companies or self-measured data. However, very little data can be shared from pharmaceutical companies because of confidentiality. Moreover, very few applications of high-throughput screening methods in formulation development leads to very limited self-measured data. The CD formulation data in current research was collected from the published literature, but almost all datasets from the published literature are flawed, such as missed, repeated, experimental or analytical errors. Furthermore, the reproductivity crisis of scientific literature leads to worse situation of the reliability of data. Thus, the preparation of high-quality dataset is the key step in machine learning process.

Molecular descriptor contributions for LightGBM model

In order to reveal the molecular mechanism and facilitate the future formulation design of drug–n class="Chemical">CD complexes, the relative importance of molecular descriptors in the model construction were calculated. Fig. 4 listed the ranking of relative importance of the molecular descriptors in the LightGBM model for the prediction of complexation free energy, which based on the used times of the feature in the model.

Figure 4

The relative importance of the molecular descriptors in the LightGBM model. Molecular descriptor with the suffix of x represented the descriptor of guest molecules, while with the suffix of y meant the descriptor of cyclodextrins. pH and T (K) represented the experimental conditions in the phase solubility study. The relative importance of molecular descriptors was determined by the measuring the number of times of the features used in the LightGBM model.

The relative importance of the molecular descriptors in the LightGBM model. Molecular descriptor with the suffix of x represented the descriptor of guest molecules, while with the suffix of y meant the descriptor of cyclodextrins. pH and T (K) represented the experimentn class="Chemical">al conditions in the phase solubility study. The relative importance of molecular descriptors was determined by the measuring the number of times of the features used in the LightGBM model.

It can be clearly remarked that the major contribution to the complexation free energy stemmed from the guest molecules, followed by experimental conditions and n class="Chemical">CDs. The top five important molecular descriptors with the contribution value above 700 in the LightGBM model were minimum projection radius_x, solvent accessible surface area_x, complexity_x, XLogP3_x and maximum projection radius_x. The minimum projection radius_x and maximum projection radius_x both characterized the size of guest molecules, where the minimum projection radius was more important. The geometrical or steric parameters of guest molecules are decisive rather than chemical factors because the guest compounds with proper size are accommodated into the inner cavity of CDs. During the process of the molecular dynamic simulation of drug–CD complexes, the preferred binding orientation of guest molecules are that the smaller ends have priority to insert into the inner cavity of CDs. Therefore, the minimum projection radius played the most important role in the complexation process.

Solvent accessible surface area_x is the surface area of the guest molecules accessible to a solvent. In the CD complexation process, the solvent accessible surface area of guest molecules is closely related to the transfer free energy from the aqueous environment to the non-polar n class="Chemical">CD inner cavity. The process that guest molecules displace the polar water molecules from the apolar CD cavity is an energetically favorable interactions, which can help to shift the equilibrium to form the inclusion complexation. In general, one of the main driving force of CD complexation is the release of enthalpy-rich water molecules from the inner cavity, and to gain an apolar–apolar association between guest molecules and inner cavity.

XLogP3_x is the partition coefficient, which represents the ratio of concentrations of the guest molecules in a mixture of water and oil at equilibrium. The vn class="Chemical">alue of XLogP3 is proportional to the hydrophobicity of a molecule, which is related to the hydrophobic interaction during the drug–CD complexation process. Guest molecules with strong hydrophobicity can contribute to the stability of the CD complexation system because of the favored hydrophobic effect. The complexity_x of the guest molecules in current research is defined as the topological complexity, which can be characterized as a steric descriptor of guest molecules by the number of atoms or bonds, branching and cyclicity, etc. To some extent, it represents the flexibility or the rigidity of guest molecules contributing to the CD complexation. In a word, three of the five most important molecular descriptors were the steric descriptors of guest molecules, while another two were hydrophobic descriptors.

Apart from the five most important descriptors, molecular weight_x, topological polar surface area_x and van der Wan class="Chemical">als surface area_x of the guest molecules also correlated with the complexation free energy. pH value of the experimental conditions with the contribution value over 600 also showed the importance because most of the guest molecules were either weak acids or weak bases. The ionization of drug molecules has a great influence on the complexation process. For example, the acidic molecule has weak complexation capability in the alkaline solution than neutral environment. Therefore, the value of complexation free energy was strongly affected by pH or the temperature of the experimental conditions.

Experimental determination of KTP–CDs complexation free energy

Phase solubility study has been a widely applied technique to determine the complexation free energy for the drug–CD inclusion complexes. Phase solubility curves of n class="Chemical">KTP with different CDs in distilled water at 37 ± 0.5 °C clearly dedicated that the aqueous solubility of KTP linearly increased with the increasing concentrations of CDs except for α-CD with a low value of R2 (shown in Table 2), which could be identified as AL type. Concerning that the value of all the slopes in KTP–CD solubility curves was less than unity, the stoichiometry may be supposed to be 1:1 (except for α-CD). In order to avoid the deviations of intercepts, the stability constants were calculated based on the intrinsic solubility of KTP. The complexation free energy was calculated to −18.178 kJ/mol for Hp-β-CD, −18.021 kJ/mol for DMe-β-CD and −17.876 kJ/mol for SBE, which showed no big difference between KTP and these three β-CD derivatives. But the calculated binding free energy between KTP and three parent CDs showed an order of β-CD (−16.489 kJ/mol) < γ-CD (−11.986 kJ/mol) < α-CD. Therefore, the solubilization effect of CD derivatives on KTP was significantly stronger than that of native CDs.

Table 2

Phase solubility study of ketoprofen with different CDs in distilled water at 37 ± 0.5 °C.

System	Equation	Kca (L/mol)	R2b	ΔGc (kJ/mol)
KTP-α-CD	y = 0.007x + 1.157	6.408	0.578	N/A
KTP-β-CD	y = 0.397x + 1.038	598.523	0.993	−16.489
KTP-γ-CD	y = 0.103x + 0.999	104.388	0.983	−11.986
KTP-Hp-β-CD	y = 0.559x + 0.744	1152.340	0.985	−18.178
KTP-DMe-β-CD	y = 0.544x + 1.046	1084.530	0.999	−18.021
KTP-SBE-β-CD	y = 0.530x + 1.423	1025.145	0.993	−17.876

Stability constant.

Correlation coefficient.

Gibbs free energy at the temperature of 37 ± 0.5 °C.

Phase solubility study of ketoprofen with different n class="Chemical">CDs in distilled water at 37 ± 0.5 °C.

MD simulation calculation of KTP–CDs complexation free energy

The MM-PBSA is a widely used method to cn class="Chemical">alculate the complexation free energies for the CD complexation in the solvent by MD simulation. The calculated complexation free energy and energy components of KTP–CD inclusion complexes by MD simulations were shown in Table 3. The binding free energy by the MM-PBSA showed the order of Hp-β-CD < Me-β-CD < SBE-β-CD < β-CD < γ-CD < α-CD, which was well in agreement with the order determined by the phase solubility study. However, the absolute values by the MD simulation were relatively larger than those of the experimental determinations because MD simulations were performed in the ideal state without any consideration for the aggregation of CD molecules. In addition, KTP molecule could not be entrapped by the inner cavity of α-CD and the KTP-α-CD system did not obtain equilibrium. Thus, the binding free energy of the KTP-α-CD system could not be calculated.

Table 3

Calculated complexation free energy and energy components of KTP–CD inclusion complexes by MD simulation.

System	ELE	VDW	PBSOL	PBTOT (kcal/mol)	TΔS (kcal/mol)	ΔG (kcal/mol)	ΔG (kJ/mol)
α-CD	N/A	N/A	N/A	N/A	N/A	N/A	N/A
β-CD	−7.82	−25.91	16.95	−16.78	−13.21	−3.57	−14.94
γ-CD	−9.41	−26.38	18.71	−17.08	−13.81	−3.27	−13.68
Hp-β-CD	−10.33	−33.21	23.41	−20.14	−14.06	−6.08	−25.44
DMe-β-CD	−6.89	−29.06	16.7	−19.24	−14.28	−4.96	−20.75
SBE-β-CD	−10.25	−27.45	19.02	−18.68	−13.44	−4.94	−20.67

ELE, electrostatic energy as calculated by the MM force field; VDW, van der Waals contribution from MM; PBSOL, non-polar and polar contributions to solvation; PBTOT, binding enthalpy; TΔS, total Entropy. Temperature has been multiplied in as 310 K.

Calculated complexation free energy and energy components of n class="Chemical">KTP–CD inclusion complexes by MD simulation.

ELE, electrostatic energy as calculated by the MM force field; VDW, van der Wan class="Chemical">als contribution from MM; PBSOL, non-polar and polar contributions to solvation; PBTOT, binding enthalpy; TΔS, total Entropy. Temperature has been multiplied in as 310 K.

The final binding free energy was made up of four energy components: non-bonded ELE, VDW, n class="Chemical">PBSOL and total entropy. The VDW showed the highest contribution to the stability of the complexes, followed by TΔS and the ELE. In fact, the inclusion of a guest molecules in the CD cavity is the displacement of entrapped polar water molecules from the hydrophobic inner cavity by the polar guest molecules. Thus, the complexation equilibrium was driven by the favorable electrostatic energy and van der Waals, which was closely related to the geometric properties, ionized conditions, hydrogen bond donors or acceptors of the guest molecules. The dis-favorable solvation free energy stemmed from two factors: (1) the polar/apoplar interaction between entrapped water molecules and the CD inner cavity, and (2) the polar/apoplar interaction between free water molecules and guest molecules, which heavily depended on the hydrophobic properties of guest molecules. According to the analysis of relative importance of the molecular descriptors (shown in Fig. 4), six of top ten important descriptors in predicting the complexation free energy were the geometric properties, while three reflected the hydrophobic properties and one related to the ionized conditions of guest molecules. Thus, the analytical results of machine learning model further elucidated the factors and specific properties governing the binding free energy of the drug–CD complexation, which provided deep insight into the molecular interaction between the guest molecules and CDs.

Comparison of the experimental, theoretical and predicting results of KTP–CDs complexes

To demonstrate the predictive performance of built LightGBM model relative to typical MD simulation cn class="Chemical">alculation, the comparison of complexation free energy between KTP and CDs provided by the experimental determination, MD simulation calculation and LightGBM prediction were performed as shown in Fig. 5. Taking the experimental results as a reference, LightGBM model showed high accuracy predictability with the MAE value of 1.282 kJ/mol, while the MD simulation calculation presented a large deviation with the MAE value of 3.206 kJ/mol. Traditional development of CD formulations primarily relies on the experimental trial-and-error and limited understanding of the physicochemical properties of guest molecules, which lead to high cost and long research process. Simulation-based calculation and data-driven predictive model provide alternatives to experimental determination as implementation of pre-formulation screens. Due to the limitation of computing capacity and empirical force field, molecular simulations encounter the difficulty to mimic the real experimental process and thus generate larger deviation. However, molecular modeling is able to provide the structural, dynamic and energetic information of the CDs systems. Thus, among these three major strategies of formulation design, machine learning predictive model demonstrated significant advantages, not only with high predicting accuracy and efficiency, but also increasing understanding about the molecular interaction of the drug–CD complexation process. On the other hand, the main difficulty of data-driven machine learning model is large amount of reliable data. The integration of these three approaches (experiment, molecular modeling and data-driven machine learning) will be the main challenge in the next step.

Figure 5

Complexation free energy between ketoprofen and various cyclodextrins by the experimental determination (ΔG_Experiment), MD simulation calculation (ΔG_MD simulation) and LightGBM prediction (ΔG_LightGBM).

Complexation free energy between ketoprofen and various n class="Chemical">cyclodextrins by the experimental determination (ΔG_Experiment), MD simulation calculation (ΔG_MD simulation) and LightGBM prediction (ΔG_LightGBM).

Linear regression is the traditional tool to establish a predicting model in n class="Chemical">CD systems. Katritzky and co-workers built QSAR multiple regression models with a dataset of 218 compounds to predict the complexation free properties between diverse guest molecules and β-CDs. Many other published reports also applied multiple linear regression analysis to build predicting models between guests and α-CDs/γ-CDs36, 37, 38, 39. Even though some of these models showed good fit between input variables and complexation free energy, their application range were severely limited because their predictive ability was generally applicable to only guest molecules with similarly chemical structures10, 40, 41.

In addition, there were several prediction models by the machine learning methods, most of the predicting results for the test set were significantly poorer than those for the training set because of very limited dataset of no more than 300 data points, which made susceptible to overfitting and chance correlation. Meanwhile, previous efforts primarily focused on the native n class="Chemical">CDs, while there were very few applications of machine learning to predict the complexation free energy of CD derivatives. Merzlikine and co-workers applied cubist and random forest to build up the β-CDs and SBE-β-CDs complexation models, which was the first study to establish a predictive model for β-CDs derivatives. Furthermore, these predictive models were only applicable for one specific CDs, not general to all CDs. The dataset of current research comprised of 3000 data points with 1320 different structures of chemicals. Thus, the prediction models presented high accuracy both on the training and test set, which demonstrated the generalization ability of prediction models was significantly better than those of the past models.

Until now, great efforts have been made to elucidate the possible intermolecular interaction and molecular forces which may contribute to the stability of the inclusion complexes43, 44, 45. The binding forces for the formation and affinity of CD inclusion complexes include hydrophobic interaction, van der Wan class="Chemical">als interaction, hydrogen bonding, relief of high-energy water, and relief of conformational strain energy8, 46, 47. However, the relative contributions and even the nature of the different forces in the complexation process are still uncertain48, 49. Thus, although complexation free energy had been determined for thousands of drug–CD complexes, it is still unable to fast and accurately predict the actual complexation free energy of a complexation between a given CD and a given drug molecule in a given condition, and then to identify whether the drug–CD complexation could be feasible or successful in the formulations. In current research, the relative importance of each molecular descriptor in the LightGBM model for the prediction of complexation free energy not only help us to understand the molecular mechanism of CD complexation, but also to make a pre-formulation screening for the complexation free energy between CDs and guest molecules.

Conclusions

Vast resources and time are currently spending on the determination or calculation of complexation free energy between n class="Chemical">CDs and molecules. The research successfully established predictive models for the complexation free energy between CDs and guest molecules by three machine learning methods. The LightGBM model showed the highest prediction accuracy. The full dataset comprised of 3000 formulation with 1320 different structures. Therefore, the prediction models showed high predictive accuracy and good generalization ability without succumbing to overfitting. Furthermore, input variables included the molecular descriptors of compounds and experimental conditions, which indicated that the prediction models were not only applicable to specific CDs but to all CDs. Introducing experimental conditions into input variables to construct predictive models was also the first attempt in CD systems, which has significantly practical guidance on the CD formulation preparation. In addition, the predictive accuracy of LightGBM model on the sparse data decreased with the decreased size of dataset. The calculation of molecular descriptor contributions provided a clear framework about the factors that affect the complexation free energy between the CDs and guest molecules. Comparison with molecular modeling, the predictive models showed better agreement with experimental determination. However, currently research was only applicable to the drug–CD binary systems. Our future research will build the predictive models for the ternary or multiple-component CD systems and improved understanding of the molecular mechanism of these systems by integrated machine learning and molecular modeling techniques.

In a word, the inspiration from current research is that machine learning method can help to rule out useless avenues of formulation design and increase the efficiency of formulation researches, which will be especially important to early drug discovery with very scarce or expensive compounds and to the complicated formulation composition. Moreover, the prediction models can n class="Chemical">also help us to identify the key factors of formulation components and process parameters. The integrated methodology of machine learning and molecular modeling approaches will also help us to simplify formulation design and improve the productivity of other pharmaceutical formulation development.