An established literature-based data model was selected for the development and validation process. The data has been refined to include only verified records, focusing on features of OFDTs such as tablet hardness, thickness, friability, and punch size. To expand our database, we conducted a literature review using the Scopus, Web of Science, and Google Scholar databases. A keyword search strategy was employed, including terms like “oral disintegrating,” “fast disintegrating,” “rapidly disintegrating,” and “oral dispersible.” The formulation should specify the total quantity of all excipients. Additionally, tablet quality attributes such as hardness, thickness, friability, punch size, and disintegration time should be included (Hana et al. 2018).

A total of 248 articles were retrieved through a database search. Out of these, only 185 research articles were selected for data extraction. Upon further manual search, 93 articles did not meet the inclusion criteria and were excluded from the study. After thorough sorting, 92 articles yielded a total of 1076 formulations. The formulation data included the name of the active pharmaceutical ingredient (API), other excipients, and process details, all of which were documented in the dataset. The final dataset consisted of the following parameters for each formulation: API name, dose, excipient name, dose (each excipient displayed in a separate column), hardness, friability, thickness, punch size, and DT. This information was then used for modeling using ML techniques (Momeni et al. 2023).

According to the European Pharmacopoeia 10^th edition, orodispersible tablets should breakdown within 3 minutes. Therefore, any records from the database that exceed 180 seconds have been excluded from further analysis. A correlation study was conducted to investigate the relationship between the dependent variable (DT) and various independent factors such as API, process parameters, and composition (Szlek et al. 2022).

The process of developing the ML model is split into three phases: data pre-processing, modeling, and model interpretation, as shown in Fig. 1.

Figure 1. 

Schematic representation of an applied workflow.

After completing data collection, the data must be processed before building predictive models to ensure the robustness and effectiveness of ML models. Several commonly employed methods, such as data cleansing, dimension reduction, imbalanced data solutions, and data splitting strategies, are necessary for data analysis. Data cleaning is performed to identify missed observations and is done by replacing data points with median or mean values. However, there are limitations to replacing missing values, as a decrease in data size may impact the accuracy of the model. Dimensionality reduction is used to eliminate the least significant features in the dataset, reducing overfitting issues and simplifying the model’s complexity. Various approaches to dimensionality reduction, such as principal component analysis (PCA), high correlation filtering, and random forest feature selection, are commonly used in data processing. Imbalanced data solutions address the uneven distribution of different database classes, as using an unbalanced dataset in a prediction model can lead to poor performance. Data splitting is another crucial step in data processing, where the entire dataset is randomized and divided into three subcategories: training, validation, and testing. The training set is used to train the models; validation is for tuning hyperparameters and preventing overfitting; and the testing set is used to assess the prediction potential of unknown data. The recommended ratio for these categories is 70% for training, 20% for validation, and 10% for testing, though the ratios may vary depending on the data size. Therefore, data preprocessing and splitting strategies are essential steps before undertaking the task.

ML modeling tasks involve various techniques such as classification, regression trees, neural networks, and potentially many other algorithms. These models are trained using prepared databases, and their performance is evaluated using an error metric. Keeping track of different modeling methods and exploring various features can be challenging and computationally expensive. Therefore, AutoML (Automated Machine Learning) is utilized. AutoML approaches often use ensemble learning strategies, which combine several model types to produce predictions that are more reliable. In this case, TPOT AutoML employed the K-fold cross-validation technique to generate a definitive production model by selecting features based on a predefined threshold. Each fold consists of a distinct training-testing pair and a validation set, with 568 records randomly selected for training, 244 records for validation, and 348 records for testing.

After completing the ML modeling process, it is necessary to evaluate the predictive performance of ML models to understand how well they generalize to new, unseen data. ML models are often prone to overfitting, which occurs when the model not only learns the underlying patterns in the training data but also the noise and random fluctuations that come with the data. To prevent overfitting and ensure model stability, the selected attributes indicate the significant impact they had on the model’s assumptions, highlighting the opaque nature of machine learning models. A small amount of data was removed during the model evaluation using the K-fold approach. This research involves training and validating models using a five-fold cross-validation method, followed by selecting features using a Python script. The training and validation procedures were repeated five times to thoroughly cover the input database and achieve the optimal model. After selecting the final input feature vector, the model was trained using a 10-fold cross-validation procedure. Root mean square error (RMSE), normalized root mean square error (NRMSE), coefficient of determination (R²), mean absolute error (MAE), and mean square error (MSE) are used to measure the robustness of the models. Seven algorithms from the TPOT AUTOML platform were utilized for feature selection and final model development: DTR, GBR, RFR, ETR, LASSO, SVM, and DL.

$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left({\text{ pred }}_i-ob{s}_i\right)}^2}{n}}$

$NRMSE=\frac{RMSE}{ob{s}_{max}-ob{s}_{min}}(100\%)$

$R^2=1-\frac{S{S}_{res}}{S{S}_{tot}}=1-\frac{{\sum }_{i=1}^n{\left(pre{d}_i-obs\right)}^2}{{\sum }_{i=1}^n{\left(ob{s}_i-obs\right)}^2}$

The variables in the equation are as follows: “obsi, predi” represent the practical and expected values, respectively; “i” is the data record number; “n” is the total number of records; “obs_max” is the highest experimental value; “obs_min” is the least observed value; “R²” is the coefficient of determination; “SS_res” is the sum of squares of the residual errors; “SS_tot” is the total sum of the errors; and “obs” is the arithmetic mean of observed values (Szlek et al. 2022).

The accuracy of ML results cannot be improved simply by fitting data into models. As the data becomes larger and more complex, better data handling techniques such as DTR, GBR, random RFR, ETR, LASSO, SVM, and DL become necessary to handle it.

DTR is a versatile algorithm used for classification and regression tasks simultaneously. It operates on the concept of breaking down complex problems into simpler, more manageable subproblems, making it an excellent choice for various applications. Decision trees (DT) have a hierarchical structure, with conditions applied from the tree’s root to its leaves. This structure allows for a step-by-step decision-making process. One of the key strengths of DT is its transparent and interpretable structure. The rules generated by DT are easy to understand. Once trained on a dataset, DT can produce logical rules that can be applied to new, unseen data by recursively dividing them into subgroups based on the conditions learned during the training phase.

GBR generates a series of decision trees, with each tree addressing the errors of the previous one. The model is generated iteratively, with each iteration adding a new decision tree to the ensemble and focusing on the errors or residuals of the combined model from prior iterations. The loss function is a crucial component in GBR as it determines the variation between the predicted and actual values of the desired variables. The algorithm minimizes this loss function during each iteration, ensuring that the model is continually improving. MSE is a commonly used loss function in GBR, where the average squared difference between the expected and actual values is calculated. Overall, GBR is a powerful algorithm for regression tasks and is widely used in practice due to its flexibility, high predictive accuracy, and ability to handle complex relationships in data (Ghazwani et al. 2023).

$f\left(x\right)=\sum _{m=1}^M{\beta }_m{h}_m\left(x\right)$

RFR is a machine-learning-based regression algorithm. It is built on bagging and random subspace algorithms. Due to its versatility, capability to handle uncertain data, and suitability for high-dimensional feature spaces (with many predictors), RFR is widely respected. In recent years, RFR has emerged as the most advantageous general-purpose algorithm. The “divide and conquer” approach, which involves bootstrapping data subsets, building decision trees on each subset, and then aggregating these results, best characterizes.

${\hat{f}}_{RF}^C\left(\mathbf{x}\right)=\frac{1}{C}\sum _{i=1}^C{T}_i\left(\mathbf{x}\right)$

The RFR employs a vector input variable x to generate an output. This is achieved by merging the predictions of the C decision trees. Ti(x) indicates a regression tree created from a subset of input variables and bootstrapped samples (Borup et al. 2023).

ETR is an enhanced method for addressing generalization (overfitting) concerns associated with random forest (RF). This approach is a recent advancement in the field of ML and can be viewed as an extension of the widely used RF. Its purpose is to minimize the risk of overfitting. Similar to RF, ETR trains each base estimator using a random subset of features. It does not select a feature and its corresponding value for use in node splitting (Hameed et al. 2021).

LASSO is a linear regression method that reduces the total sum of squares of residuals and the sum of the absolute values of the regression coefficients. The regression coefficient can be obtained using the given equation.

$\sum _{i=1}^m{(y^{\left(i\right)}-x^{\left(i\right)}b)}^2+\lambda \sum _{j=1}^n\left|b_j\right|$

Where m represents the number of samples, n represents the number of x variables - y(i) and x(i) ε R^{1x m} are the y and x values in the ith sample, respectively; bj is the jth regression coefficient; and λ is the hyperparameter. b is denoted by the following formula:

b = (b₁ b₂ ... b_n)^T

In LASSO, the regression coefficient, bj, can be reduced to zero, leading to the removal of the corresponding x. The study considered a range of values for λ from 2–15 to 2–14 …, 2–2, and 2–1 to find the value that maximizes the coefficient of determination, r², by 5-fold cross-validation. Scikit-learn was employed to estimate the LASSO (Kaneko 2021).

SVM is the most commonly used method in ML for classification, regression, and other tasks. SVM operates in high- or infinite-dimensional space and constructs a hyperplane or multiple hyperplanes. The hyperplane that maximizes the distance from the closest training data points in each class achieves significant separation. A larger margin typically results in a lower generalization error for the classifier. It is effective in high-dimensional spaces and can exhibit different behaviors based on mathematical functions like the kernel. SVM classifiers utilize various types of functions, such as linear, polynomial, radial basis function (RBF), and sigmoid, as kernel functions. However, if the dataset contains a higher level of noise, such as overlapping target classes, the performance of SVM is compromised (Gaye et al. 2021; Pérez and Bajorath et al. 2022).

DL is primarily used as a neural network, as shown in Fig. 2. DL can automatically extract features and transform basic representations into more complex abstraction layers without needing a separate feature extractor. DL is more sensitive to small and specific modifications in complex networks, leading to higher accuracy than typical ML techniques. DL algorithms have shown superior performance compared to other ML methods in effectively forecasting the in vitro performance of pharmaceutical formulations. Deep Neural Networks (DNN) have been utilized in pharmaceutical research, particularly in drug design, drug-induced liver toxicity, and virtual screening. Deep learning can create sophisticated and complex systems that represent various objects using chemical descriptors. This can greatly aid in the development of drugs and their prediction (Roggo et al. 2020).

Figure 2. 

The network architecture of DNN.

As ML models are inherently black boxes, efforts have been made to shed light on their prediction techniques. In our study, illustrated in Fig. 1, we used the Lundberg et al. SHAPLEY additive explanation (SHAP) approach to explain the relationship between input and output variables. The SHAP method is primarily based on cooperative game theory and has been applied in various domains, including the pharmaceutical industry. The concept of Shapley values is used to assess the contribution of each participant or individual to the overall team effort or outcome. In machine learning, Shapley values have been adapted to explain the impact of each feature in a predictive model. SHAP values offer a method to distribute the model’s prediction to each feature in a fair and consistent manner. The mathematical formula for SHAP values is provided below.

${\varphi }_j\left(val\right)=\sum _{S\subseteq \{1,\dots ,p\}\backslash \left\{j\right\}}\frac{\left|S\right|!(p-|S|-1)!}{p!}(val(S\cup \left\{j\right\})-val(S\left)\right)$

Where “S” is a subset of features in the model, “x” is the vector of feature values to be explained, and “p” is the number of features. The prediction “Valx(S)” is the result of estimating feature values in set “S” (Rozemberczki et al. 2022).

The Shapley value calculation method adheres to the axioms of efficiency, symmetry, dummy, and additivity, thus explaining the process for developing predictions. Random samples are used to replace the values of each attribute to assess their importance and impact. Computing the Shapley value can be computationally intensive due to the numerous potential coalitions of feature values that must be considered. Coalitions are carefully chosen to reduce repetitions, leading to decreased calculation time. However, this approach also leads to an increase in the variation of the Shapley value. The k-means method was utilized to reduce the number of repetitions needed to represent each feature’s impact. The k-means algorithm was set up with 12 centroids, each corresponding to a feature data domain in a cluster. A comprehensive SHAP matrix can be created by grouping the data domain of each characteristic. Displaying this matrix facilitates the understanding of the model’s predictions (Szlek et al. 2022).

The pre-processed database contained 92 direct compression OFDTs (new data entries), including 28 unique APIs and 50 variable coding compositions (excipients were topologically coded). There were 5 variables encoding formulation dimensions such as thickness [mm], hardness [N], friability [%], punch size [mm], and DT [s]. Descriptive statistics (Table 1) show that the variables did not follow a normal distribution and the formulations were significantly positively skewed (right-skewed distribution), as shown in Fig. 3. The database was divided using a 10-fold cross-validation method that was calibrated to ensure input variables were classified fairly across the splits.

Figure 3. 

Box and violin plots of specific features from the database. The graph illustrates the interquartile range (IQR) using boxes, which include the first quartile (Q1), median (horizontal line), and third quartile (Q3). The lower whisker is calculated as Q1–1.5*IQR and the upper whisker as Q3 + 1.5*IQR. These graphs demonstrate the dispersion of numerical data through the kernel density function.

The raw data and curated data are available at (Raw database) https://doi.org/10.6084/m9.figshare.25880377. (Curated database) https://doi.org/10.6084/m9.figshare.25880560. (accessed on 22 May 2024)

The selection of features and development of the final model were conducted using an automated approach using the TPOT AUTOML method. The dimensions of the TPOT AUTOML are provided in Table 2. This table displays the accuracy of the models that were developed. The values of RMSE, NRMSE, R², MAE, and MSE are also analyzed to assess the accuracy and precision of the model’s output.

As expected, ML techniques, especially ETR, have proven to be the best pipeline in the TPOT AutoML analysis for the curated dataset. However, it is clear that the deep learning model, when hyper-tuned, performed on par with its ML counterparts, achieving great R² values and NRMSE percentages. After the initial evaluation of the deep learning models trained using a five-fold cross-validation scheme, it was found that when the DL was trained with three hidden layers of 100 neurons each, 56 input neurons, and one output neuron, using a tanh activation function with 2200 epoch values combined with a 10-fold cross-validation scheme, the model accuracy significantly improved. This improvement was reflected in the NRMSE and R² values, as shown in Fig. 4.

Figure 4. 

Scatter plot between actual and predictive values for the disintegration time of the DL model.

The input variables were categorized into two main groups: composition and manufacturing parameters. Features below the variable importance level were eliminated, except for those in the composition section. The final input vector consisted of 18 inputs. (Table 3) displays the chosen features and their scaled significance.

Table 3 shows that the quantities of disintegrants (Crospovidone and croscarmellose sodium) and binder (Microcrystalline cellulose) have varying levels of significance. The data suggests that the number of disintegrants will have the most significant impact on the predicted outcomes. The lower significance of components like lubricants and fillers may be due to the positive skewness of the variable distribution.

The SHAP summary graph in Fig. 5 illustrates the influence of features on the model’s predicted output. The colored line shows actual feature values and their impact on a prediction along the x-axis. A SHAP summary plot helps identify overall impacts and underlying assumptions.

Figure 5. 

SHAP dependence diagram for the ML models for the top 20 attributes.

A higher disintegration time (DT) is predicted with a greater concentration of disintegrants such as crospovidone and croscarmellose sodium. Conversely, when it comes to fillers like Mannitol and Avicel, a higher quantity of Mannitol and Avicel leads to lower DT. For binders like Microcrystalline Cellulose (MCC), Hydroxypropylmethylcellulose (HPMC), and Polyvinyl Pyrrolidine (PVPK30), two distinct effects are observed. A higher amount of HPMC and PVPK30 results in higher DT, while MCC decreases DT at higher concentrations. Lubricants like magnesium stearate (MgSt) and sodium stearyl fumarate (SSF) also play a role. SSF tends to increase DT, likely due to its hydrophilic properties. On the other hand, a higher amount of MgSt, which is more lipophilic, lowers DT due to the occlusion effect (Szlek et al. 2022).