Duality in Simplicity and Accuracy in QSPR: A Machine Learning Framework for Predicting the Solubility of Diverse Pharmaceutical Acids in Deep Eutectic Solvents

1. Introduction

Carboxylic acids and their derivatives play a hugely important role in pharmacy; they are not simple active medicinal substances but also key structural motifs used to design new biologically active molecules [1,2,3]. Due to the presence of a carboxyl group, they can participate in proton dissociation, form hydrogen bonds, undergo conjugation with plasma proteins, and influence absorption and distribution processes in the body. In pharmaceutical practice, aromatic and aliphatic acids are found as components of anti-inflammatory, analgesic, antioxidant, or diuretic drugs. Their bioavailability and effectiveness often depend on the form of salt, esters, or complexes in which they occur [4,5]. The experimental part of this study is focused on two pharmaceutically important compounds, namely mefenamic acid and niflumic acid, belonging to the broad class of non-steroidal anti-inflammatory drugs (NSAIDs) [6,7]. Mefenamic acid is widely used for treating mild to moderate pain, including menstrual pain (dysmenorrhea), and for inflammatory conditions requiring short-term pain relief [8,9,10]. Niflumic acid, in addition to its anti-inflammatory and analgesic properties, is also used as a research tool in pharmacology, for example, as a non-specific blocker of certain ion channels (e.g., chloride channels), which makes it an interesting subject for research beyond typical therapeutic applications [11,12,13]. Their physicochemical properties have direct practical consequences, particularly limited solubility in the aqueous phase, which complicates oral and parenteral formulations and necessitates the use of auxiliary agents such as solubilizers, surfactants, or salts [14,15,16]. Additionally, several other carboxylic acids were considered during the construction of a machine learning model. Flufenamic acid is used as an anti-inflammatory drug, although its clinical application is currently limited [17,18,19]. Ibuprofen and ketoprofen are some of the most widely used NSAIDs with analgesic and antipyretic properties [20,21,22,23,24]. They are widely available over-the-counter and are generally well-tolerated. Probenecid represents another category of acid use; it's a uricosuric drug used to treat gout and extend the half-life of antibiotics by inhibiting renal excretion [25,26]. In addition to synthetic drugs, there is also a broad group of natural acids, most often phenolic, found in plants and having significant biological potential [5,27,28]. Ferulic acid is a powerful antioxidant that protects against oxidative stress and is used in research on neurodegenerative diseases [29,30]. Caffeic acid shows anti-inflammatory, anti-cancer, and hepatoprotective activity [31,32,33]. p-Coumaric acid and syringic acid also belong to the phenolic acids, with antioxidant properties and potential significance in the prevention of lifestyle diseases [34,35,36,37,38].

The solubility of a particular substance is one of the most critical parameters in pharmacy, materials chemistry, and process engineering [39,40]. In the context of the pharmaceutical industry, the solubility of an active pharmaceutical ingredient (API) determines both pharmacokinetic properties and production and storage strategies. The poor solubility of many organic molecules limits bioavailability, necessitating the use of modified forms (e.g., amorphous formulations, nanocrystallization, and lipid carriers) and increasing product development costs [41,42]. Beyond the pharmacological aspect, solubility also influences the selectivity of extraction processes, the efficiency of crystallization, purification, and recycling of raw materials, as well as the environmental behavior of compounds [43,44]. It is also a multifactorial property affected by temperature, pH, ionic strength, pressure, and specifically the thermodynamic properties of the solid matter (lattice energy, polymorphism) and specific interactions with solvent molecules [45,46,47]. In light of this, precise solubility profiles are key not only for evaluating drug efficacy but also for predicting product stability and the risk of adverse physicochemical transformations during storage and administration.

Among various techniques used for solubility enhancement of APIs [41,48,49], deep eutectic solvents (DESs) are particularly interesting and promising. DESs are a flexible and increasingly studied class of liquid systems formed by mixing appropriate hydrogen bond donors and acceptors, which leads to a significant reduction in the melting point relative to the starting components [50,51,52,53]. DES are distinguished by a number of properties useful in a pharmaceutical context: low volatility (which promotes safety and reduces emissions), a wide spectrum of polarity and acidity, significant "tunability" through component selection, and the ability to solubilize compounds of various chemical natures [54,55]. In pharmacy, DES are being researched as extraction media for isolating natural products, as potential carriers for formulations that enhance the solubility and bioavailability of APIs, and as so-called THEDES (therapeutic deep eutectic solvents), i.e. systems in which the components themselves may have biological activity or facilitate drug stabilization and delivery [56,57,58,59]. At the same time, the vast combinatorial space of possible components and molar ratios necessitates the use of accelerated testing methodologies and predictive tools, as an empirical review of all variants is costly and time-consuming.

In this context, computational solubility prediction, especially using QSPR and advanced machine learning methods, is becoming a strategic tool [60,61,62,63]. QSPR models translate the relationship between molecular structure and macroscopic properties, while modern ML algorithms (including neural networks, graph models, and ensemble strategies) can capture complex, non-linear patterns and interactions between multiple variables. Neural networks offer particularly great flexibility: they can integrate different types of representations (physicochemical descriptors, fingerprints, molecular graphs, and features describing composition and properties) and learn relationships that are difficult to capture with traditional linear models [64,65,66]. Recently, such neural networks were successfully employed for solubility predictions of drugs [67,68,69]. At the same time, the question arises of balancing the simplicity and interpretability of models with their ability to accurately represent reality. In practice, this means the need for careful descriptor selection, cross-validation, estimation of prediction uncertainty, and techniques to prevent overfitting.

The purpose of this work was to create a predictive model for estimating the solubility of various pharmaceutically active carboxylic acids. The model was developed by combining COSMO-RS-derived molecular descriptors with machine learning methods, based on our newly established DOO-IT (Dual-Objective Optimization with ITerative feature pruning) framework. New experimental data for mefenamic acid and niflumic acid were obtained for this study, which were supplemented with solubility values found in literature for a number of acids used in the pharmaceutical realm. The constructed models were thoroughly validated, and their performance was discussed, highlighting their effectiveness and the potential for generalization.

2. Results and Discussion

2.1. Solubility Measurements of Mefenamic Acid and Niflumic Acid

Figure 1 summarizes the mole-fraction solubilities at 25 °C of mefenamic acid (x_MEF) and niflumic acid (x_NIF) in choline chloride- and menthol-containing DESs, with full numerical values in Table S1 (x_MEF) and Table S2 (x_NIF). Across all systems, x_MEF spans 1.38 × 10⁻⁴ to 1.40 × 10⁻² and x_NIF 2.38 × 10⁻⁴ to 2.11 × 10⁻², and menthol-based DESs generally afford higher solubility than their choline-chloride analogues for both compounds. The highest x_MEF is observed for Men/TRG 1:3 (1.40 × 10⁻²), with Men/TEG also giving elevated values, while the highest x_NIF appears for Men/GLY 1:1 (2.11 × 10⁻²) and remains high in Men/TRG 1:1–1:3. Within the choline chloride series, ChCl/DEG 1:1 provides the top x_MEF (6.49 × 10⁻³), and ChCl/TRG 1:1 gives the top x_NIF (1.47 × 10⁻²). Considering the hydrogen-bond donors (HBDs), TRG and TEG are associated with the largest solubilities in the menthol series, GLY is particularly favorable for x_NIF at 1:1, and DEG stands out for x_MEF in the choline chloride series; the effect of HBD fraction is system-specific; for example, x_MEF increases from 3:1 to 1:3 in Men/TRG, x_NIF maximizes at 1:1 in Men/GLY, and values decrease beyond 1:1 in ChCl/DEG (x_MEF) and ChCl/TRG (x_NIF).

2.2. Identifying an Optimal Predictive Model via the DOO-IT Framework

The DOO-IT (Dual-Objective Optimization with Iterative feature pruning) framework, aiming to find the most accurate and parsimonious machine learning model, was applied to an extended dataset of new solubility measurements and literature data for pharmaceutical acids in deep eutectic solvents. In total, N = 1020 points characterize the solubility of a variety of pharmaceutical acids in choline chloride and menthol-based deep eutectic solvents. The DOO-IT pipeline performs repeated dual-objective optimizations, minimizing at the same time MAE and model complexity. It combines an iterative backward feature pruning and candidate models selected based on the Pareto front are ranked by the corrected Akaike Information Criterion (AICc) to select solutions that balance fit and parsimony. This three-stage workflow allows for fully automated and unbiased selection of the optimal model, defined as both accurate and parsimonious. The repetition of the procedure with varying selection of the number and types of molecular descriptors allows for effective probing of the complex relationships between data and molecular descriptors. The result of so many trials of the protocol is provided in Figure 2. This is the central pillar of the model development, visualizing the outcome of the DOO-IT model selection workflow. The presentation of the AICc distribution as a function of the number of descriptors used for the model training enables direct quality model assessment. The red envelope curve, which traces the minimum AICc value achievable for a given number of descriptors, illustrates the fundamental trade-off between accuracy and simplicity inherent in QSPR modeling. A clear, non-linear trend is immediately apparent. The AICc values drop precipitously from a high baseline with 3-4 descriptors to a deep minimum at 6-7 descriptors. This initial decline signifies that each of these first few descriptors adds substantial explanatory power, capturing the primary physicochemical drivers of solubility and drastically improving the model's performance with minimal risk of overfitting. Beyond this point, the curve exhibits a more complex behavior, rising and then gradually descending to a second, even lower minimum. This pattern confirms that a simplistic feature selection approach is insufficient and a systematic, iterative pruning strategy is essential to navigate this rugged optimization landscape.

The most significant finding illustrated in Figure 2 is the existence of two distinct "basins of excellence," highlighted by the green boxes. This discovery directly refers to the "duality" in our article's title, emphasizing that not a single "best" model is presented, but a strategic choice between two high-performing, yet fundamentally different, modeling philosophies. The results of the AICc stability analysis reveal two distinct centers of excellent performance rather than a single minimum. One basin is centered around 6 descriptors (an ultra-parsimonious) and the other at 16 descriptors (the high-performance). The convex envelope of the best models shows that intermediate complexities are in many cases suboptimal or redundant for the present chemical space. This outcome indicates that distinct, scientifically meaningful descriptor combinations can achieve competitive performance via different mechanisms. Two separate models can therefore be regarded as the ones with optimal performance, depending on the intended task. A simple, general model (the 6-descriptor one) that captures the main drivers of solubility with high efficiency and a more complex, nuanced model (the 16-descriptor one) that adds subtle but important interaction terms to achieve the highest possible accuracy. This duality is a crucially important finding, once again proving the robustness and usefulness of the developed DOO-IT framework. It also highlights that the machine learning model space can be very complex and requires a thorough analysis to select the optimal model.

The first basin, located at 6 descriptors, represents the pinnacle of model parsimony, which has important implications. The descriptors forming this model likely represent the most fundamental and universal properties of the organic acid molecules that dictate their behavior in DES. Indeed, this particular model utilized the following descriptors: log(x_API^COSMO), E_vdW,API, ΔHH1, ΔHH2, ΔHH4, and ΔHBA1. Hence, for the studied set of data, the dispersion contribution of API and relative hydrophobicity values are the most dominant contributions to the solubility. Also the relative value of acceptability at vicinal σ-potential regions is important. The second high-performance model utilizes 16 descriptors (log(x_API^COSMO), μ_API, E_tot,API, E_Misfit,API, E_HB,API, E_vdW,API, μ_DES, E_Misfit,DES, ΔHBD1(12), ΔHBD4, ΔHH1, ΔHH2, ΔHH4, ΔHBA1, ΔHBA2, and ΔHBA4. This set of descriptors reveals the complex nature of solute-DES interactions under saturated conditions and the necessity of extension of the core features of the former model with additional polarity hydrogen bonding capacity of solute and solvent. Hence, the extended 16-descriptor solution integrates a broader set of energetic and interaction terms, capturing solvation phenomena with higher resolution.

It is imperative to contextualize these findings within a crucial methodological framework. The dataset exclusively comprises organic acids of pharmaceutical relevance, including well-known compounds such as ketoprofen, ibuprofen, ferulic acid, probenecid, caffeic acid, p-coumaric acid, syringic acid, and flufenamic acid. A significant and unresolved challenge in modeling such systems is the unknown dissociation state (pKa) of acidic solutes within the complex, non-aqueous environment of deep eutectic solvents (DES), which are rich in both hydrogen bond donors (HBD) and acceptors (HBA). Lacking a reliable method to determine the precise ionization state of each acid in every unique DES, we adopted a necessary and consistent simplification: all molecular descriptors were calculated for the neutral, non-dissociated forms of the molecules. Despite this crude simplification, the remarkable accuracy of the resulting models is particularly noteworthy. It strongly suggests that the fundamental physicochemical properties of the parent molecule are the dominant drivers of solubility and that our modeling framework is robust enough to capture these governing relationships despite the simplification of the solute's ionization state.

2.3. Performance of the Optimal Solubility Models

As indicated above, the stability analysis based on the AICc criterion revealed the existence of two “basins of excellence”, which means that two specific models were selected for a detailed inspection. The results of the Dual-Objective Optimization for the 6-descriptor model are presented in Figure 3, while the 16-descriptor model is described in detail in Figure 4. The plots highlight the inherent trade-off between accuracy and model complexity. The collection of non-dominated solutions forms a clear Pareto front (dark purple points), which represents the best accuracy attainable at any given level of complexity. Models lying to the right of this front (grey points) are inferior since a simpler and more accurate alternative exists. Coloring of the points along the Pareto front according to the nu hyperparameter reveals a consistent pattern: lower nu values yield simpler models with reduced SV ratios, whereas higher nu values give rise to more complex models characterized by larger SV ratios.

The Pareto analysis and 1-SE rule for the 6-descriptor set favored the chosen model (trial 1793) that slightly trades the very best MAE for substantially lower complexity. The trial 1790 model offers a slightly lower MAE value but is more complex when taking into account the support-vector ratio. The parity plot shows tight clustering around y=x for both training and held-out test points, indicating low bias and excellent generalization. This final 6-descriptor model demonstrates strong predictive performance on the held-out test set, achieving a coefficient of determination (R^²) of 0.947, indicating that approximately 95% of the variance in the target variable is explained by the model. The mean absolute error (MAE) of 0.0746 and root mean squared error (RMSE) of 1.831 reflect reasonable average prediction errors. However, there is a small but noticeable gap between MAE and RMSE, suggesting that there are some outliers or systematic deviations. The mean absolute percentage error (MAPE) of 7.46% further corroborates the model’s high relative accuracy. Residual analysis reveals negligible bias (mean residual = 1.642) and a standard deviation (0.811) consistent with the RMSE, supporting the assumption of homoscedastic and unbiased errors. Collectively, these metrics indicate a robust, acceptable calibrated model with promising generalization capability.

In the case of the high-performance 16-descriptor solution, trial 955 was selected over trial 1597, being much less complex with only a minimal drop in the MAE value. The parity for this model is on an excellent level and the scatter is minimal, as evidenced in the right panel of Figure 3. The model significantly extends the descriptor set compared to the ultra-parsimonious 6-descriptor one and consequently provides a more complete description of the balance between solute–solvent and solvent–solvent interactions. The finalized model built on 16 descriptors shows strong performance on the holdout test set, with R² = 0.970, meaning it explains about 97% of the variance in the target. Error magnitudes are modest (MAE = 0.047; RMSE = 1.795), and the noticeable gap between RMSE and MAE suggests occasional outliers or mild systematic deviations. A MAPE of 4.67% further supports the model’s high relative accuracy. Residuals indicate limited bias (mean = 1.611) and variability (SD = 0.793) consistent with the RMSE, aligning with the assumption of homoscedastic, essentially unbiased errors. Collectively, these results point to a well-calibrated, robust model with promising generalization potential, although it is much more complex compared to an ultra-parsimonious one.

3. Materials and Methods

3.1. Materials

Mefenamic acid and niflumic acid (both ≥97%, Sigma-Aldrich, St. Louis, MO, USA) were used as received. The hydrogen-bond acceptors were choline chloride (ChCl, CAS 67-48-1, ≥99%) and menthol (Men, CAS 89-78-1, ≥98.5%), and the hydrogen-bond donors comprised ethylene glycol (ETG, CAS 107-21-1), diethylene glycol (DEG, CAS 111-46-6), triethylene glycol (TEG, CAS 112-27-6), tetraethylene glycol (TRG, CAS 112-60-7), glycerol (GLY, CAS 56-81-5), 1,2-propanediol (P2D, CAS 57-55-6), and 1,3-butanediol (B3D, CAS 107-88-0); all polyols/polyethers were obtained from Sigma-Aldrich with stated purity ≥99%. Methanol (analytical grade, CAS 67-56-1; Chempur, Piekary Śląskie, Poland) was used for sample handling where applicable. Unless otherwise specified, all chemicals were employed without additional purification.

3.2. Solubility Measurements Procedure

A similar methodology to [70] was employed, adapted here for mefenamic acid (MEF) and niflumic acid (NIF). Each DES was prepared at the target molar ratio by gentle heating with stirring until a clear single phase formed, then equilibrated to 25 °C. Pre-equilibrated aliquots were spiked with an excess of MEF or NIF, sealed, and agitated isothermally for 24 h at approximately 60 rpm. After equilibration, suspensions were maintained at 25 °C, supernatants were withdrawn, filtered through 0.22 μm PTFE syringe filters, and analyzed by UV–Vis. Spectra were collected from 200 to 500 nm in quartz cuvettes; analytical wavelengths were set at the absorption maxima (λ^MEF_max = 351 nm; λ^NIF_max = 344 nm). Concentrations obtained from UV–Vis were converted to mole-fraction solubilities (x_MEF or x_NIF) using molar masses and the density of each saturated solution; densities were determined gravimetrically at 25 °C. All measurements were performed in triplicate.

Calibration curves were established for each compound using methanolic stock solutions and serial dilutions. For MEF: calibration range 0.002 to 0.078 mg mL^-1, slope 28.265, intercept -0.010, linearity R² = 0.9993, LOD = 0.00261 mg mL^-1, LOQ = 0.00790 mg mL^-1. For NIF: calibration range 0.005 to 0.090 mg mL⁻¹, slope 18.808, intercept -0.006, linearity R² = 0.9994, LOD = 0.00272 mg mL^-1, LOQ = 0.00825 mg mL^-1.

3.3. COSMO-RS Computations

Application of the COSMO-RS framework [71,72,73,74,75] requires appropriate representation of molecular diversity. This is done by performing the conformational analysis prior to the determination of any thermodynamic properties. For this purpose, the default protocol was applied, taking advantage of the COSMOconf (version 2023, BIOVIA COSMOlogic) / TURBOMOLE (version 7.7, 2023, TURBOMOLE GmbH) tandem for the generation of the most representative structures for all solutes and solvent molecules. The applied protocol is consistent with previously published schemes [76,77,78]. For each molecule, up to ten low-energy conformations were determined for both gas and condensed phases, the latter accounting for solvent effects under the conductor-like screening model. The resulting "cosmo" and "energy" files were generated using the BP_TZVPD_FINE_24.ctd parameter set, essential for thermodynamic calculations in COSMOtherm, which requires application of the RI-BP/TZVP//TZVPD-FINE level of theory.

3.4. Molecular Descriptors

Two distinct sets of molecular descriptors were generated using the COSMO-RS theory. The first set of descriptors comprised interaction energies from solubility calculations [79,80,81]. While the standard iterative solubility protocol is typically sufficient, it frequently yields erroneous predictions of complete miscibility for highly soluble solutes in DES systems [87,88,89,90]. For these cases, complete solid-liquid equilibrium (SLE) calculations were mandated. Requisite thermodynamic fusion data for the solid solutes, including melting temperature, T_m, and enthalpy of fusion, ΔH_fus, were obtained by averaging available literature values [91]. The heat capacity of fusion was approximated as constant, ΔC_p,fus≈ΔS_fus≈ΔH_fus/T_m. The resulting Gibbs free energy of fusion values, ΔG_fus=ΔH_fus-TΔS_fus, utilized in the calculations are provided in the Supplementary Materials. The COSMO-RS output files yielded five primary descriptors for each solute: total intermolecular interaction energy, E_int,API, its constituent electrostatic misfit, E_misfit,API, hydrogen bonding, E_HB,API, and van der Waals, E_vdW,API, contributions, as well as the chemical potential and μ_API. Analogous descriptors for the DES solvent (E_int,DES, E_misfit,DES, E_HB,DES, E_vdW,DES, and μ_DES, were calculated as the sum of the individual DES components, weighted by their respective molar fractions in the solute-free mixture. Relative descriptors, defined as the difference between the solute and DES values, were also included. Additionally, the computed solubility values from COSMO-RS were similarly included, log(x_API^COSMO).

Apart from this, the final set of molecular descriptors was augmented with values derived from σ-potential distributions. The standard COSMO-RS output consists of 61 data points covering the charge density range of -0.03e/Å2÷+0.03e/Å2. Consistent with prior machine learning applications, this data was reduced by averaging values over 0.005 intervals. This process resulted in a 12-step function defining three characteristic regions of the σ-potential: hydrogen bond donor (HBD1÷4, -0.03e/Ų to -0.01 e/Ų), hydrophobicity, (HH1÷4, from -0.01e/Ų to +0.01 e/Ų), and acceptability (Hydrogen Bond Acceptor, HBA1÷4, from +0.01e/Ų to +0.03 e/Ų). Consequently, four descriptors were generated for each region, leading to 24 descriptors of this type for the solute, the solvent, and the relative difference between them.

3.5. Dataset

The values of molecular descriptors were computed for all components of studied systems, including pharmaceutical acids and DES constituents. The set of solutes included compounds for which new solubility measurements were included in this paper, namely mefenamic acid and niflumic acid. In addition, the values of already published solubility data of ketoprofen, ibuprofen, ferulic acid, probenecid, caffeic acid, p-coumaric acid, syringic acid, and flufenamic acid in DES were included. In total, the dataset comprised N=1020 mole fractions at saturated conditions in choline chloride, betaine, and menthol-based DES with variety of proportions of different HBA counterparts. Both dry and water-diluted systems were included if available. All data, including solubility values, fusion data, computed solubility, and all molecular descriptors, are available in the Supplementary Materials.

3.6. Machine Learning Protocol

3.6.1. Core Algorithm and Data Preprocessing

The machine learning workflow was centered on the nu-Support Vector Regression (nuSVR) algorithm [92], chosen for its demonstrated ability to effectively model complex, non-linear relationships often present in QSPR studies [93,94,95]. To handle these non-linearities, the Radial Basis Function (RBF) kernel was selected. The RBF kernel is a powerful and flexible choice, capable of mapping features into an infinite-dimensional space, which allows it to model intricate decision boundaries while requiring the tuning of only a single parameter, gamma. The optimization of the nuSVR hyperparameters was conducted as follows: the regularization parameter C and the nu parameter were directly optimized. The kernel coefficient gamma, which dictates the influence of each support vector, was optimized via a guided, data-driven strategy. For each optimization cycle, a baseline gamma_base value was heuristically determined from the median pairwise squared Euclidean distance of the training data subset [96,97]. The optimizer then refined this anchor by searching for an optimal logarithmic scaling factor. This approach focuses the search on a physically relevant scale, enhancing optimization efficiency. Prior to model training, two standard preprocessing steps were performed. First, the full dataset (N=1020) was partitioned into a training set (80%) and a held-out test set (20%) using a fixed random seed to ensure reproducibility. Second, all molecular descriptors in the training set were standardized by removing the mean and scaling to unit variance using the StandardScaler from scikit-learn [98,99]. As SVR algorithms are sensitive to feature scaling, this step ensures that no single descriptor disproportionately influences the model due to its magnitude. The same scaling transformation was subsequently applied to the test set.

3.6.2. Dual-Objective Optimization Protocol

To explicitly manage the inherent trade-off between model accuracy and simplicity, a dual-objective optimization (DOO) strategy was implemented using the Optuna framework (v. 3.2) [100,101,102]. The TPE sampler within Optuna was configured to simultaneously minimize two competing objectives, which were evaluated using a 5-fold cross-validation scheme on the training data. The first objective was predictive accuracy, quantified by the Mean Absolute Error (MAE). The second objective was model complexity, quantified by the mean support vector (SV) ratio. The SV ratio is calculated for each fold as the number of support vectors divided by the number of training samples in that fold, providing an intrinsic measure of complexity for nuSVR models. A model with a lower SV ratio is considered more parsimonious.

Hence, the outcome of a dual-objective optimization is a set of solutions forming a Pareto front. This front consists exclusively of non-dominated solutions. A solution is considered non-dominated if no other solution exists that is superior in one objective without being inferior in the other. In other words, to improve a non-dominated solution with respect to one objective, a trade-off in the form of a degradation in the other objective must be accepted. Conversely, a dominated solution is one for which at least one other solution exists that offers better performance in one objective while being no worse in the other, making it an objectively suboptimal choice.

3.6.3. Iterative Model Refinement and Candidate Selection

The framework employs an iterative backward pruning methodology to integrate feature selection directly into the optimization process. This automatic procedure relies, therefore, on both Dual-Objective Optimization and Iterative features pruning (DOO-IT). The procedure begins with the complete descriptor set. A full DOO is executed, producing a Pareto front of non-dominated models. From this front, a single candidate model for the current iteration is selected, governed by the 1-Standard Error (1-SE) rule [103,104]. This involves identifying the most accurate model on the front and defining a performance threshold based on its standard error; the simplest model (lowest SV ratio) within this threshold is then chosen. Once a candidate is selected, its features are ranked based on permutation importance with 10 repeats [105]. The least impactful descriptor is then eliminated, and the procedure repeats with a new, full DOO on the reduced feature set. This cycle continues until a specified minimum number of features is reached, generating a series of robust, parsimonious candidate models at each level of complexity.

3.6.4. Final Model Selection via Information Criterion

To objectively select the single best model from the family of candidates generated by the iterative procedure, an information-theoretic criterion was used: the corrected Akaike Information Criterion (AICc) [106,107,108,109]. AICc provides a principled method for model selection by balancing goodness-of-fit against model complexity, thereby penalizing overfitting. As nuSVR is a non-probabilistic model lacking an explicit likelihood function, a pseudo-AICc formulation was adopted, a common practice for algorithm-based models. This approach approximates the log-likelihood term by assuming Gaussian residuals and estimates it from the model's Residual Sum of Squares (RSS). The critical component in this calculation is the effective number of parameters, k. For nuSVR, k was estimated as the number of support vectors plus two additional terms representing the model's bias (intercept) and the estimated error variance. This pseudo-AICc formulation enabled a principled, quantitative comparison between models of varying complexity. The candidate model exhibiting the minimum pseudo-AICc value across all runs and iterations was chosen as the final, most justified solution for subsequent analysis.

The DOO-IT framework was implemented as a fully automated pipeline using Python 3.10 [110] with the scikit-learn [111], Optuna [102], and pandas [112] libraries. To rigorously assess solution stability, the entire procedure was repeated fifteen independent times. Each dual-objective optimization within this process was configured to run for 2000 trials, ensuring a comprehensive exploration of the solution space.

4. Conclusions

The study addresses the challenge of solubility prediction, a problem of central importance in pharmaceutical and green chemistry research. Accurate predictive models therefore provide a powerful tool to reduce experimental workload, accelerate drug development pipelines, and enable the rational design of novel solvent systems such as DESs that combine efficiency with environmental compatibility. The Dual-Objective Optimization with Iterative feature pruning (DOO-IT) framework was applied for this task.

The study demonstrated that stability analysis of the DOO-IT framework uncovered not a single global optimum but two distinct regions of predictive excellence for modeling solubility of pharmaceutical acids in deep eutectic solvents. On one side of the solution landscape lies an ultra-parsimonious six-descriptor model that combines high predictive performance with minimal computational cost. Importantly, this model does not rely on COSMO-RS solubility input, instead exploiting physically meaningful descriptors that capture key thermodynamic drivers such as solvent–solute chemical potential differences, hydrogen-bonding characteristics, and van der Waals interactions. As a result, it provides a reliable and interpretable tool for rapid prescreening of candidate systems, especially in contexts where extensive quantum-chemical calculations or fusion data are unavailable. At the other extreme, a sixteen-descriptor model emerges as the configuration with the lowest corrected Akaike Information Criterion and the highest overall predictive accuracy. This model integrates COSMO-RS logarithmic solubility as an anchor descriptor, which enables it to correct systematic deviations at solubility extremes and deliver near-perfect agreement with experimental values. Although more resource-intensive, this solution is particularly valuable for applications where the highest quantitative fidelity is required, such as in fine-tuning formulation parameters or guiding late-stage development decisions.

By revealing these two complementary “basins of excellence,” our analysis highlights the versatility of the DOO-IT framework in identifying multiple, scientifically meaningful optima that balance accuracy, parsimony, and interpretability. The findings also extend our previous works, where a single model was sufficient to describe a narrower chemical space. In the present, more diverse dataset, the appearance of dual optimal regimes underscores the importance of tailoring model complexity to the scope of the prediction task. Taken together, these results suggest a pragmatic two-tiered strategy for future studies of solubility in deep eutectic solvents and related systems. Initial high-throughput screening can be effectively performed with the COSMO-RS-free parsimonious model, while subsequent high-precision evaluation of promising candidates can benefit from the expanded descriptor set that incorporates COSMO-RS calculations. This workflow balances efficiency with accuracy, making it possible to explore broader chemical spaces without sacrificing predictive reliability. Looking forward, validating the transferability of both models to other classes of solvents, as well as developing ensemble or adaptive strategies that dynamically combine parsimonious and high-performance regimes, will further enhance the applicability of this approach in green chemistry and pharmaceutical design.