Quantification of Caffeine Interactions in Choline Chloride Natural Deep Eutectic Solvents: Solubility Measurements and COSMO-RS-DARE Interpretation

Solubility of active pharmaceutical ingredients is an important aspect of drug processing and formulation. Although caffeine was a subject of many studies aiming to quantify saturated solutions, many applied solvents suffer from not being environmentally friendly. This work fills this gap by presenting the results of solubility measurements in choline chloride natural deep eutectic solvents, ccNADES, comprising one of seven of the following polyalcohols: glycerol, sorbitol, xylitol, glucose, sucrose, maltose and fructose. The ratio of ccNADES components was optimized for maximizing caffeine solubility at room temperature. Additionally, temperature dependent solubility was measured for the first four systems exhibiting the highest solubility potential, both in their neat forms and in mixtures with water. Results were used for intermolecular interactions assessments using the COSMO-RS-DARE approach, which led to a perfect match between experimental and computed solubility values. An important methodological discussion was provided for an appropriate definition of the systems. Surprising linear trends were observed between the values of fitting parameters and water-ccNADES composition. In addition, comments on selection of the values of the fusion thermodynamic parameters were provided, which led to the conclusion that COSMO-RS-DARE solubility computations can effectively compensate for the inaccuracies of these important physicochemical properties.


Introduction
Caffeine (1,3,7-trimethylxanthine) is a stimulant of the central nervous system belonging to the class of methylxanthines. This group of compounds is widespread in the natural environment and can be found in tea and other plant leaves, coffee and cocoa beans, as well as cola seeds [1]. Methylxanthines can be classified as purine-based alkaloids, having a fused heterocyclic system in their moiety, which contains a pyrimidine and an imidazole ring [2]. The general properties of methylxanthines include rather poor solubility in water or ethanol [3] and low values of water-octanol partition coefficient [4]. Caffeine can be found in over sixty different plants, which usually contain other methylxanthines as well [5,6]. It has been used as a traditional beverage for thousands of years [7,8] and it is currently a very popular [9,10] physical [11,12] and cognitive [13,14] performance enhancer, although its immoderate consumption is related to some side effects [15][16][17]. Caffeine after oral intake is readily and completely absorbed from the small intestine [18]. It is then metabolized and in the form of mainly paraxanthine excreted with urine [19]. Similarly, as other methylxanthines, it binds with adenosine receptors [16] and induces psychomotor stimulant properties in the brain [18]. Caffeine is also used in different medical applications, namely in combination with painkillers, e.g., to relieve migraine [20,21], in bronchopulmonary dysplasia treatment [22], orthostatic hypotension treatment [23] and the treatment of apnea of prematurity [24].  The results of the DSC studies were consistent with the ones obtained through FTIR measurements. For clarity, only the results obtained for glycerol containing ccNADES and its water mixtures were presented in Figure 2. Both pure caffeine and the precipitates can be characterized by a melting point at circa 238 • C which is in satisfactory agreement with literature data for caffeine polymorphic form labeled as I. The mean value for the whole dataset is 238.24 • C with a standard deviation of 0.42 • C. When taking into account the position of the first endothermic signal there are however significant differences among the studied samples. Pure caffeine and the samples obtained after dissolution in neat ccNADES are characterized by a peak with a mean value of 157.98 • C and a standard deviation of 0.46 • C. This can be assigned as the transition from the polymorphic form II of caffeine to its form I, however compared to literature this peak is shifted toward higher temperatures. In the case of caffeine solubilized in mixtures containing both ccNADES and water the peak has a maximum located at a mean value of 151.58 • C with a standard deviation equal 0.59 • C, which once more is somewhat shifted compared to literature data. Nevertheless, it can be attributed to the loss of water experienced by the caffeine hydrate. It seems therefore that both FTIR and DSC measurements give a consistent picture, pointing out the fact that the presence of water in the mixtures is responsible for the formation of a hydrate, while the polymorphic form II of caffeine is the result of crystallization in the case of absence of water in the systems.

Experimental Solubility of Caffeine
As it was said before, caffeine was the subject of numerous solubility studies in neat and mixed solvents and authors also contributed to this field [77] with already published results which confirm the validity of the approach applied here. This was established by comparing the obtained solubility values and the ones provided in literature [62][63][64] with Table S1 in Supplementary Materials showing detailed results of this comparison. Since there are no reported data regarding caffeine solubility in the ccNADES systems studied here, measurements in water and some organic solvents were used. Our protocol was also successfully applied for the characteristics of other compounds dissolved in ccNADES including theophylline [90], theobromine [91], curcumin [92] and sulfonamides [93], which further confirms the adequacy of the used methodology.
Although it is expected that in general the solubility of many organic solutes in ccNADES is higher compared to traditional organic solvents and their mixtures [90][91][92][93], the selection of the optimal composition is not a straightforward task. Hence, in the first stage, the experimental pre-screening was performed using binary solvents comprising choline chloride and one of the seven studied polyols in three different molar ratios. The results of this preliminary testing performed at 25 • C were collected in Figure 3. Three different molar ratios were tested, namely a unimolar proportion of both constituents, as well as 2:1 excess amounts of either of them. the selection of the optimal composition is not a straightforward task. Hence, in the first stage, the experimental pre-screening was performed using binary solvents comprising choline chloride and one of the seven studied polyols in three different molar ratios. The results of this preliminary testing performed at 25 °C were collected in Figure 3. Three different molar ratios were tested, namely a unimolar proportion of both constituents, as well as 2:1 excess amounts of either of them. From the obtained data it can be inferred that in most cases the 1:1 molar ratio of choline chloride and the polyolic constituent yields the highest solubility of caffeine. However, for the ccNADES system including glycerol the composition with an excess amount of the polyol is responsible for the most efficient caffeine solubilization. In this case, the mole fraction of caffeine in the solvent was found to be xC = (297.97 ± 0.70)·10 −4 , which stands for almost a 15-times increase in solubility compared to water. All the other solvents were most efficient in their unimolar composition. The following concentrations were found for saturated systems of caffeine in ccNADES comprising sorbitol (xC = (220.77 ± 0.23)·10 −4 ), xylitol (xC = (207.36 ± 0.82)·10 −4 ) and glucose (xC = (171.22 ± 0.58)·10 −4 ). This yields an increase of solubility compared to water from 8.5-to 11-times. Detailed values of the mole fractions of caffeine in the above solvents together with their standard deviations can be found in Table S2 in Supplementary Materials. From the obtained data it can be inferred that in most cases the 1:1 molar ratio of choline chloride and the polyolic constituent yields the highest solubility of caffeine. However, for the ccNADES system including glycerol the composition with an excess amount of the polyol is responsible for the most efficient caffeine solubilization. In this case, the mole fraction of caffeine in the solvent was found to be x C = (297.97 ± 0.70)·10 −4 , which stands for almost a 15-times increase in solubility compared to water. All the other solvents were most efficient in their unimolar composition. The following concentrations were found for saturated systems of caffeine in ccNADES comprising sorbitol (x C = (220.77 ± 0.23)·10 −4 ), xylitol (x C = (207.36 ± 0.82)·10 −4 ) and glucose (x C = (171.22 ± 0.58)·10 −4 ). This yields an increase of solubility compared to water from 8.5-to 11-times. Detailed values of the mole fractions of caffeine in the above solvents together with their standard deviations can be found in Table S2 in Supplementary Materials.
The four optimal compositions were chosen for the next stage of experiments in which the ccNADES were treated as cosolvents and mixed with water in varying proportions. A total of twelve different compositions were studied including both neat ccNADES and water. Four different temperatures from 25 • C to 40 • C with 5 • C intervals were studied in this series of measurements. The obtained results were collected in Figures 4-7 and the corresponding detailed values of caffeine mole fractions can be found  in Supplementary Materials, Tables S3-S6. When taking into account the obtained solubility profiles, a similar picture emerges for all four studied ccNADES systems. It turned out that a relatively small amount of water added to the NADES promotes the increase of caffeine solubility. For all of the studied systems, the highest solubility was obtained when the amount of ccNADES in the solution was equal x* ccNADES = 0.8. The order of increased solubility of caffeine found for pure ccNADES remains the same in their mixtures with water and the increase in temperature was responsible for elevated solubility of caffeine. As mentioned beforehand, the NADES based on choline chloride and glycerol resulted in the highest solubility of caffeine, i.e., x C = (559.49 ± 2.06)·10 −4 at 40 • C, which was 19% higher than the pure ccNADES. For sorbitol and glycerol, the most effective compositions responsible for caffeine solubility were equal to x C = (413.65 ± 2.09)·10 −4 and x C = (394.23 ± 2.43)·10 −4 at 40 • C, respectively. The least effective composition involving glucose as the second constituent yielded caffeine solubility equal to x C = (356.94 ± 1.17)·10 −4 also at 40 • C, which amounted to an almost 30% increase in solubility compared to pure ccNADES. While the solubility increase in caffeine in these mixtures is not very pronounced, it is worth noting that an addition of small amounts of water to the system is a factor improving the handling of the systems due to decreased viscosity and density. Moreover, problems with high hygroscopicity of ccNADES systems become omitted by such optimal formulations.
The four optimal compositions were chosen for the next stage of experiments in which the ccNADES were treated as cosolvents and mixed with water in varying proportions. A total of twelve different compositions were studied including both neat ccNADES and water. Four different temperatures from 25 °C to 40 °C with 5 °C intervals were studied in this series of measurements. The obtained results were collected in Fig    The four optimal compositions were chosen for the next stage of experiments in which the ccNADES were treated as cosolvents and mixed with water in varying proportions. A total of twelve different compositions were studied including both neat ccNADES and water. Four different temperatures from 25 °C to 40 °C with 5 °C intervals were studied in this series of measurements. The obtained results were collected in Fig    When comparing the results obtained in this study with solubility of caffeine in organic solvents [77], it is evident that ccNADES perform generally better, both in their neat form and when considering mixtures with water. For example, at 25 • C the solubility of caffeine in NADES containing choline chloride with glycerol (the most effective among studied neat ccNADES) is 165% of caffeine solubility in pure DMSO, which was found to be the best performing organic solvent. Among the organic solvents only DMSO and DMF can compete in solubilization efficiency with the studied ccNADES. Since both of these solvents are related to serious toxicity issues [94,95], the choice of NADES as solubilizing media is even more evident.  When taking into account the obtained solubility profiles, a similar picture emerges for all four studied ccNADES systems. It turned out that a relatively small amount of water added to the NADES promotes the increase of caffeine solubility. For all of the studied systems, the highest solubility was obtained when the amount of ccNADES in the solution was equal x*ccNADES = 0.8. The order of increased solubility of caffeine found for pure ccNADES remains the same in their mixtures with water and the increase in temperature was responsible for elevated solubility of caffeine. As mentioned beforehand, the NADES based on choline chloride and glycerol resulted in the highest solubility of caffeine, i.e., xC = (559.49 ± 2.06)·10 −4 at 40 °C, which was 19% higher than the pure ccNADES. For sorbitol and glycerol, the most effective compositions responsible for caffeine solubility were equal to xC = (413. 65   When taking into account the obtained solubility profiles, a similar picture emerges for all four studied ccNADES systems. It turned out that a relatively small amount of water added to the NADES promotes the increase of caffeine solubility. For all of the studied systems, the highest solubility was obtained when the amount of ccNADES in the solution was equal x*ccNADES = 0.8. The order of increased solubility of caffeine found for pure ccNADES remains the same in their mixtures with water and the increase in temperature was responsible for elevated solubility of caffeine. As mentioned beforehand, the NADES based on choline chloride and glycerol resulted in the highest solubility of caffeine, i.e., xC = (559.49 ± 2.06)·10 −4 at 40 °C, which was 19% higher than the pure ccNADES. For sorbitol and glycerol, the most effective compositions responsible for caffeine solubility were equal to xC = (413. 65

Intermolecular Interactions of Caffeine with Solvents Components
The observed high solubilizing nature of studied ccNADES indicates that caffeine can strongly interact with solvent constituents. These interactions can be quantified by inspection of the structure and energetics of pairs identified via conformational analysis. Among many potential intermolecular complexes, the most energetically favorable ones are shown in Figure 8. The electron density distributions are presented in the form of socalled "metafiles" representing caffeine conformers involved in the intermolecular complex formation in ccNADES solutions. These caffeine conformers are not considered during typical conformational analysis and are added as a result of clusters examination. Hence, the optimization of pairs leads to eventually new geometries with alternative electron density distributions compared to monomeric caffeine. The way of presentation offered in Figure 8 emphasizes that interactions of caffeine with solvent molecules extend the pool of potential diversity of solute structures in the solutions. It is interesting to note that the most favorable caffeine dimer adopts a stacking conformation. In addition, the alternative hydrogen bonded pair has been found but it happens to be less stable by 1.6 kcal/mol. However, all the most stable pairs formed by caffeine with ccNADES constituents are stabilized by a strong hydrogen bond formed between the N-acceptor center of caffeine located on imidazole ring and the hydrogen atom of hydroxyl group of solvent molecule. The values of the most important energetic contributions to the stabilization of caffeine-ccNADES pairs are collected in Table 1. The last column comprises the final affinity values including all corrections, as enumerated in the methodology section. The most obvious conclusion drawn from the presented data is that polyalcohols have the highest affinity toward caffeine while the interactions of caffeine with choline chloride or water are much less attractive. Dimerization of caffeine is characterized by the lowest affinity among studied intermolecular complexes. From the fact of affinity discrepancies arises a methodological question of the number of structures to be included for adequate representation of caffeine diversity in the pool of conformers. It is important to emphasize that the origin of COSMO-RS-DARE [96] is the observation that interactions occurring in the solutions and leading to new species can alter electron density distributions due to new compound formation.    Hence, the σ-profiles were generated for all pairs and compared to the one obtained for monomeric caffeine. The results are presented in Figures 9 and 10. The former exemplifies detailed σ-profiles plots for glycerol containing system. It is quite visible that dimerization is the only case for which no significant changes in σ-profiles are observed. This is particularly evident on the top panel of Figure 9 presenting differential σ-profiles obtained after subtracting one of the caffeine monomers from the one obtained for pairs. For the rest of complexes, a non-trivial alteration of charge distributions takes place. The dotted lines split the region of σ-profiles into three sections typically associated with hydrogen donating and accepting capabilities separated by a hydrophobicity region. It is worth noting that all three regions are strongly affected upon pair formation which suggests that such conformers should not be removed from solubility computations. The changes of HB accepting capability of caffeine after binding to ccNADES is quite expected due to modification of densities on the imidazole ring after binding with the hydroxyl group of ccNADES constituents. Less obvious are however, the other alterations. Negative values on differential plots of σ-plots indicate decreasing electron density of caffeine involved in pairs formation compared to monomeric caffeine. First, the hydrogen atom connected to the carbon of the imidazole ring is apparently quite acidic and can be directly involved in hydrogen bonding [77]. This is supposed to be typical for all methylxanthines, as it was already documented for theophylline [90] and theobromine [91]. Hence, the pair formation with proton donors also reduces the donating character of this acidic center. A bit surprising is the increase in hydrophobicity of caffeine bound to ccNADES constituents or water, which can be attributed to the overall withdrawing of the electron density from the rest of the caffeine molecule not directly involved in the intermolecular complex formation. On the other hand, caffeine dimerization introduced no significant change in the charge distribution, so adding it to the overall pool of conformers is supposed to be superfluous. On the contrary, the pairs of caffeine with water or ccNADES constituents should be included in COSMO-RS DARE computations as introducing new density distribution characteristics of caffeine. However, these σ-profiles, although different from monomeric caffeine, are quite similar to each other, as documented by RMSD values presented in Figure 10. This suggests the possibility of serious simplifications of COSMO-RS DARE computations. In principle, each new conformer resulting from new complex formation can be characterized by introducing a new set of interaction parameters, as defined in Equation (1). However, due to similarities of σ-profiles associated with pairs formation, in all further computations only one set of adjustable parameters was used for the representation of intermolecular interactions with new caffeine conformations. It happened that this simplification does not reduce the accuracy of the computed values of solubility and provides averaged interaction parameter for a given caffeine-ccNADES-water system. The above section provided important methodological clues for efficient COSMO-RS-DARE computations. To complete this part, the exemplary input files were included and described in the Supplementary Materials for ensuring the reproducibility of the results presented here (see Scheme S1 and Table S7).  Comparison of σ-profiles similarities of studied complexes of caffeine with ccNADES constituents and water. As a measure of similarity, the RMSD (root-mean-square deviations) values were computed for σ-profiles presented in Figure 9 of pairs ("mcos" file) with respect to monomeric caffeine ("cosmo" file).

Results of COSMO-DARE Computations
Successful solubility computations within the COSMO-RS framework require additional characteristics of solid species. It is known [97,98] that although the COSMO-RS approach is formulated for liquid bulk systems, solids can also be treated if missing thermodynamic data are provided from external sources [99]. According to details provided in the methodology section, the solubility computations are performed as an itera-  Comparison of σ-profiles similarities of studied complexes of caffeine with ccNADES constituents and water. As a measure of similarity, the RMSD (root-mean-square deviations) values were computed for σ-profiles presented in Figure 9 of pairs ("mcos" file) with respect to monomeric caffeine ("cosmo" file).

Results of COSMO-DARE Computations
Successful solubility computations within the COSMO-RS framework require additional characteristics of solid species. It is known [97,98] that although the COSMO-RS approach is formulated for liquid bulk systems, solids can also be treated if missing thermodynamic data are provided from external sources [99]. According to details provided in the methodology section, the solubility computations are performed as an itera- Figure 10. Comparison of σ-profiles similarities of studied complexes of caffeine with ccNADES constituents and water. As a measure of similarity, the RMSD (root-mean-square deviations) values were computed for σ-profiles presented in Figure 9 of pairs ("mcos" file) with respect to monomeric caffeine ("cosmo" file).

Results of COSMO-DARE Computations
Successful solubility computations within the COSMO-RS framework require additional characteristics of solid species. It is known [97,98] that although the COSMO-RS approach is formulated for liquid bulk systems, solids can also be treated if missing thermodynamic data are provided from external sources [99]. According to details provided in the methodology section, the solubility computations are performed as an iterative procedure seeking convergence of the computed values of chemical potentials until convergence is achieved, as defined in Equation (4). Unfortunately, such default computations are often far from accurate, although for some systems one can observe a fairly quantitative correspondence. For example, solubility of benzamide and salicylamide [100] can be quite accurately predicted without additional parametrization in neat and binary solvents. However, estimations made for a very similar compound, ethenzamide, are much less accurate [101]. It was already documented that the DARE extension of COSMO-RS computations can be quite effective, as it was in the case of nicotinamide, ethenzamide, sulfonamides and methylxanthines [91,[101][102][103]. Unfortunately, for the majority of cases, default COSMO-RS is able to predict solubility at most qualitatively. This is also the case for caffeine. In Figure 11, the results of solubility computations of this solute for series of neat and binary solvents are presented and compared to available experimental data. The list of considered here systems is provided in the Supplementary Materials (see Tables S8 and S9). Comparison of computed and experimental data taken from the literature definitively exposes the inadequacy of the COSMO-RS model to describe the solubility of caffeine in neither neat nor binary solvents. In addition, the solubility values computed for caffeine in ccNADES both via traditional COSMO-RS and with COSMO-RS DARE are also provided in Figure 11. Table 2. The values of parameters defining linear relationship G int as a function of ccNADES composition derived from data presented in Figure 12. tive procedure seeking convergence of the computed values of chemical potentials until convergence is achieved, as defined in Equation (4). Unfortunately, such default computations are often far from accurate, although for some systems one can observe a fairly quantitative correspondence. For example, solubility of benzamide and salicylamide [100] can be quite accurately predicted without additional parametrization in neat and binary solvents. However, estimations made for a very similar compound, ethenzamide, are much less accurate [101]. It was already documented that the DARE extension of COSMO-RS computations can be quite effective, as it was in the case of nicotinamide, ethenzamide, sulfonamides and methylxanthines [91,[101][102][103]. Unfortunately, for the majority of cases, default COSMO-RS is able to predict solubility at most qualitatively. This is also the case for caffeine. In Figure 11, the results of solubility computations of this solute for series of neat and binary solvents are presented and compared to available experimental data. The list of considered here systems is provided in the Supplementary Materials (see Tables S8 and S9). Comparison of computed and experimental data taken from the literature definitively exposes the inadequacy of the COSMO-RS model to describe the solubility of caffeine in neither neat nor binary solvents. In addition, the solubility values computed for caffeine in ccNADES both via traditional COSMO-RS and with COSMO-RS DARE are also provided in Figure 11.  Table 2.
A perfect match between experimental values and COSMO-RS-DARE computed solubility of caffeine in ccNADES is clearly visible. On the contrary, the results obtained for neat and binary solvents using the default COSMO-RS approach are rather non-acceptable. Errors are very high and even qualitative solubility analysis seems to be dubious. From this perspective, the advantage of DARE extension is unquestionable. The only problem with these kinds of computations is the necessity of using two fitting pa-  Table 2.
parison of intermolecular stabilization of caffeine complexes with water and polyalcohols. As it was documented in Table 1, the caffeine affinity to any sugar predominates in all systems. That is why in the case of sufficient concentration of polyalcohol in the system these kinds of contacts determine the G int values. On the other hand, increased dilution with water reducing the concentration of sugars systematically diminishes their stabilizing role in ccNADES. This observation is very useful since it enables a significant reduction of the number of adjustable parameters with quite small loss of accuracy.  A perfect match between experimental values and COSMO-RS-DARE computed solubility of caffeine in ccNADES is clearly visible. On the contrary, the results obtained for neat and binary solvents using the default COSMO-RS approach are rather non-acceptable. Errors are very high and even qualitative solubility analysis seems to be dubious. From this perspective, the advantage of DARE extension is unquestionable. The only problem with these kinds of computations is the necessity of using two fitting parameters for every solution. Indeed, the observed significant improvement in the case of COSMO-RS DARE is granted from optimizing values of intermolecular interaction parameters, G int = H int + T·S int . In principle, variations of ccNADES + water composition alters the fundamental physicochemical properties of the resulting ternary solvents and each one is to be treated as a separate solvent system. Hence, for every solute-solvent system there are introduced two fitting parameters including the correction for the energetics of intermolecular interactions of those conformers which were added to the model as a consequence of intermolecular complex formation in the solutions. The values of these parameters were optimized by fitting to experimental data. Resulting optimized values are collected in Figure 12. It is interesting to see that in pure water the values of adjustable parameters adopt the lowest values and systematically and nearly linearly increase with amounts of added ccNADES up to the molar fraction equal to 0.2. Then an almost constant trend is noticed irrespectively of the amount of water. Observed trends can be quantified by linear function of the NADES composition. This was illustrated by gray lines in the case of a system comprising glycerol or glucose. The other two systems behave similarly. Additionally, parameters characterizing observed linear functions are collected in Table 2.
The accuracy of such a simplified model was tested by additional computations using G int values computed from data provided in Table 2. The result of the simplified COSMO-RS-DARE were provided in Figure 11. It is interesting to notice that the loss of accuracy is not very high and correspondence with experimental data is acceptable.
The splitting of the linear trend for x * ccNADES around 0.2 can be explained by comparison of intermolecular stabilization of caffeine complexes with water and polyalcohols. As it was documented in Table 1, the caffeine affinity to any sugar predominates in all systems. That is why in the case of sufficient concentration of polyalcohol in the system these kinds of contacts determine the G int values. On the other hand, increased dilution with water reducing the concentration of sugars systematically diminishes their stabilizing role in ccNADES. This observation is very useful since it enables a significant reduction of the number of adjustable parameters with quite small loss of accuracy.

Role of Fusion Data in Solubility Model
In final comments the influence of accuracy of fusion data on solubility computations is discussed. It was already mentioned that fusion properties, including melting temperature, fusion enthalpy or change of heat capacities upon melting, are added to the COSMO-RS model as external parameters for determining the values of Gibbs free energy, ∆G fus . This value is used for definition of the reference state for computed chemical potentials of the solute in the bulk phase of a given composition. From general thermodynamics consideration [104], the chemical equilibrium is obtained in the case of equal values of chemical potentials of the pure solid solute and dissolved state. There is a long discussion in the literature [105,106] about the proper way of computation of this important thermodynamic quantity. The problem is additionally complicated by the fact that in many cases fusion data are not directly measurable due to solid sublimations [107], thermal instabilities, polymorphism and formation of a variety of solvates, as for example curcumin [108], fullerene C60 [109], nitrofuran antibiotics [110], L-norvaline [111], β-Alanine [112] among many other cases.
Lack of fusion data is a direct obstacle in using COSMO-RS and available predictive QSPR [99] or additive [113,114] models are quite inaccurate and do not help resolving the problem. Even in the case of known experimental data characterizing fusion properties, their variety obtained by different authors often results in differing values. The most spectacular example is benzoic acid, for which one can find tens of different values of heat of fusion and melting temperatures [115]. Since solubility computations in the COSMO-RS approach are very sensitive to the values of ∆G fus , the knowledge of these values is mandatory. This problem is also important for COSMO-RS-DARE computations. In this work the values of fusion data were assumed as a mean value of published data, which are collected in Table 3. Table 3. The values of fusion data characterizing solid caffeine.

Tm [K] H fus [kJ/mol]
509.05 1 510.2 [116], 510.3 [88], 508.8 [117], 508.7 [118], 505.4 [119], 510.0 [120], 507.7 [121], 509.6 [122], 509.1 [62], 510.8 [123] 20.66 1 20.08 [116], 19.52 [88], 20.37 [117] It is interesting to quantify the sensitivity of G int to the values of ∆G fus . For this purpose, COSMO-RS-DARE was used with different values of Gibbs free energies for the most representative binary systems. The obtained results are provided in Figure 13. From these results very interesting conclusions can be drawn. Although in general the values of intermolecular parameters are sensitive to fusion data, there are quite large ranges of shifting values, δG fus , for which almost identical results are obtained. This is a very important conclusion suggesting a direct applicability of COSMO-RS-DARE to systems for which fusion data are not available or inaccurate. It is then possible to compute solubility for such solids for which either precise fusion data are missing or there are available data with large discrepancies by rational guesses of melting temperature and fusion enthalpy. Even estimates coming from predictive models can be used as such guess values. For the generalization of intermolecular interaction parameters resulting from the COSMO-RS-DARE approach, it is necessary to add a relationship of the kind provided in Figure 13. The rational guess might be understood as both reasonable values and those belonging to the region of relatively low impact on intermolecular interaction parameters. As it was documented in Figure 13, for caffeine the range of such a rational guess is quite large. It happened that reducing the values of ∆G • fus by at least of 0.2 kcal/mol is a reasonable presumption for caffeine fusion in the studied systems. It seems to be very interesting to see if this is a general conclusion, which immediately provokes forthcoming studies. which fusion data are not available or inaccurate. It is then possible to compute solubility for such solids for which either precise fusion data are missing or there are available data with large discrepancies by rational guesses of melting temperature and fusion enthalpy. Even estimates coming from predictive models can be used as such guess values. For the generalization of intermolecular interaction parameters resulting from the COS-MO-RS-DARE approach, it is necessary to add a relationship of the kind provided in Figure 13. The rational guess might be understood as both reasonable values and those belonging to the region of relatively low impact on intermolecular interaction parameters. As it was documented in Figure 13, for caffeine the range of such a rational guess is quite large. It happened that reducing the values of ΔG°fus by at least of 0.2 kcal/mol is a reasonable presumption for caffeine fusion in the studied systems. It seems to be very interesting to see if this is a general conclusion, which immediately provokes forthcoming studies. Figure 13. The correlation between the COSMO-RS DARE interaction parameters and fusion data ΔGfus = ΔG°fus + δGfus, used for solubility computations. The value of ΔG°fus corresponds to experimental fusion data provided in Table 3. Figure 13. The correlation between the COSMO-RS DARE interaction parameters and fusion data ∆G fus = ∆G • fus + δG fus , used for solubility computations. The value of ∆G • fus corresponds to experimental fusion data provided in Table 3.  119 56-81-5). Methanol (CAS: 67-56-1) was purchased from Avantor Performance Materials (Gliwice, Poland). All of the above compounds had a ≥99% purity and were used without initial procedures, apart from choline chloride which was dried before use.

Preparation of the Calibration Curve
The calibration curve used for the determination of caffeine solubility was the same as in the previous study [77]. The preparation of the curve involved obtaining a stock solution in methanol, successful dilutions and measuring the absorbance values at 272 nm wavelength. Three different curves were prepared in this manner and averaged. Detailed values of absorbance and concentrations are provided in the Supplementary Materials in Table S10, while such statistical parameters as the determination coefficient R 2 , limit of detection (LOD) and limit of quantification (LOQ) of the calibration curve are shown in Table S11.

Preparation of Samples in NADES and Their Mixtures with Water
Choline chloride was the main compound involved in NADES preparation. The second NADES constituent was one of the seven polyols described in the materials section. Choline chloride was mixed in different molar ratios with the second constituent in sealed test tubes and placed in a 90 • C water bath until a uniform solution was formed. The resulting ccNADES were used in their neat form or in mixtures with water to which they were added in different ratios.
An excess amount of caffeine was then added to test tubes containing either the neat ccNADES or the ccNADES-water mixture so that a saturated solution could be formed. The samples were then incubated for 24 h at different temperatures in an Orbital Shaker Incubator ES-20/60 from Biosan (Riga, Latvia). The temperature was adjusted with a 0.1 • C accuracy and varied by no more than ±0.5 • C during the daily cycle. The samples were simultaneously mixed at 60 rev/min during heating. Because of the increased viscosity and density of ccNADES, the samples were centrifuged at 1000 rev/min for 5 min with the use of an EBA 20 centrifuge from Hettich (Tuttlingen, Germany) so that the undissolved precipitate would remain on the bottom of the test tube. Next, the samples were filtered using a PTFE syringe filer with 0.22 µm pore size. The test tubes, syringes, filters, etc. were initially heated at a temperature corresponding to the one of the samples. Finally, a defined amount of the obtained filtrate was transferred to test tubes containing methanol and samples diluted in such way were measured spectrophotometrically. Furthermore, the density of the samples was measured for the determination of the caffeine mole fractions.

Solubility Measurements
The determination of caffeine concentration in the samples was conducted using a spectrophotometric method. The A360 spectrophotometer from AOE Instruments (Shanghai, China) was employed for this task and the spectra were recorded in the wavelength range from 190 nm to 700 nm with a resolution of 1 nm. The spectrophotometer was calibrated using methanol, which was also used for the dilution of the measured samples, so that their absorbance values would fit within the linearity limit of the method. Based on the calibration curve prepared earlier and absorbance values measured at 272 nm it was possible to determine the concentration of caffeine in the samples. The concentration values were averaged from three measurements and expressed as mole fractions.

FTIR-ATR Measurements
The FTIR-ATR analysis of solid caffeine precipitates was conducted based on the sediments collected from the test tubes obtained from the described above solubility experiments. After drying, the samples were measured with a FTIR Spectrum Two spec-trophotometer from Perkin Elmer (Waltham, MA, USA) utilizing an ATR (attenuated total reflection) device. The wavenumber range was set to be 450-4000 cm −1 .

Differential Scanning Calorimetry Measurements
The differential scanning calorimetry (DSC) measurements were also performed in order to characterize caffeine solid residues. The procedure involved co-grinding 0.4 g of caffeine samples with appropriate amounts of solvents. A Retsch (Haan, Germany) MM 200 mill was utilized for this task with a frequency of 25 Hz and 30 min duration. The specimens were ground using 5 mL stainless steel jars with two stainless steel balls and obtained powders were analyzed with a DSC 6000 from PerkinElmer (Waltham, MA, USA). The heating rate was set to 10 K/min and the inert atmosphere was provided by 20 mL/min nitrogen flow. The samples were placed in standard aluminum pans and before the measurements the apparatus was calibrated with the use of indium and zinc standards.

Affinities and Intermolecular Interactions
The solute-solvent affinities are characterized as a global measure by providing values of the Gibbs free energy of corresponding reactions, ∆G r , and also as decomposed sets of meaningful energetic contributions. The former value is computed as concentration independent quantity by assessments of activities in addition to the values of mole fraction equilibrium constant. For this purpose, COSMOtherm software [125] was invoked with the set of two default commands, namely "reaction" (for concentration-dependent values, ∆G r (x) ) and "K_activity" (for concentration-independent values, ∆G r (a) [126]. As a result, for all considered solutions, the ∆G r (a) values are the same irrespectively of concentrations provided that the same qualitative composition is kept, whereas ∆G r (x) are allowed to be different for distinct binary compositions of the same ccNADES. Illustration of both types of affinities is provided in the Supplementary Materials ( Figure S1). In this paper, the ∆G r (a) quantity is used as a fundamental affinity measure. Apart from this, the computed values of interaction energies in the gas phase were augmented with the corrections accounting for zero-point vibrations, ∆E ZPE , electron correlation, ∆E corr , and basis superposition error, ∆E BSSE . It is important to emphasize that these contributions characterize weighted sums over all conformers included in the analysis. The weighting factors correspond to Boltzmann probability of a given conformer. Hence, the affinities were represented by the following reaction scheme: {X i } + {Y j } = {XY k } where X, Y stand for monomers and XY represents resulting pairs. In the case of X = Y a dimer is obtained, and otherwise hetero-molecular dimers are formed. In principle, each reagent is repressed by different conformers and in the above scheme i-th, j-th and k-th conformers are denoted in corresponding subscripts. The interaction energies account for weighted contribution of every conformer simply by the following sum: where E connotes the energy of the compound after averaging over individual energetic contributions, ε, of a given conformer.
The type of contributions is denoted by ω = {int, ZPE, corr}, and N i , N j , N k are the numbers of conformers of a particular species. Hence, this formula defines the way of computations of interactions energy (int), electron correlation (corr) and zero-point vibrational (ZPE) corrections. The necessity of extending the system characteristics beyond DFT was emphasized by Hellweg and Eckert [127]. The values of BSSE, ω = {BSSE}, were determined solely on the geometrical basis of the pair using reasonable simplification [128]. The final corrected intermolecular interactions energies in the gas phase can be computed by assuming the independence of all contributions in an additive way: It is important to stress that each component of the above equation is computed on a different level of theory. The intermolecular interaction energies in the gas phase correspond to the results of full geometry optimization of both monomers and pairs using RI-DFT with BP97 GGA functional using def2-SVPD basis set and Grimme D3-BJ dispersion corrections [129]. The resulting geometries were used further for zero-point vibrational energy computations at the same level. Much more sophisticated computations were indispensable for accounting for the electron correlation. For this purpose, the single point energy computations were conducted at RI-MP2/def2-QZVPP level. Finally, BSSE was accounted on the same level as geometry optimization in the gas phase using gCP [128]. Since COSMOtherm is parametrized for BP86 (B88-VWN-P86) functional, it is necessary to generate "cosmo" and "energy" files which are coherent with BP_TZVPD_FINE_22.ctd parametrization. Hence, RI-DFT BP86/def2-TZVPD//BP86/def-TZVP computations were performed. This means that geometry was re-optimized using a smaller basis set, and single point energy was followed with a more extended one. These types of computations include the fine grid tetrahedron cavity and van der Waals dispersion term based on the "D3" method of Grimme et al. [129]. Completing this stage allows for computations of the values of Gibbs free energies of reactions involving solute and solvents molecules. The final values of affinity were computed as follows: All these contributions are collected in Table 1.

Solubility Computations
Theoretical solubility predictions were performed using the Conductor-like Screening Model for Real Solvents (COSMO-RS) approach. It relies on quantum chemistry computations applying the continuum solvation models based on which statistical thermodynamics analysis of electrostatic surface interactions is conducted for assessment of chemical potentials of components in the bulk phase. Although it is defined for liquid bulk systems, solids can also be considered if fusion thermodynamic characteristics is provided. There is no need for providing the details here, as they are available in original publications [130][131][132].
It is of crucial importance to mention, that the microscopic properties resulting from quantum chemical optimizations are used for generation of histograms of electron charge density, which are termed σ-profiles. Based on this information the thermodynamic macroscopic properties of components in the liquid systems are quantified by means of statistical thermodynamics through analyzing the densities probability distributions leading to socalled σ-potentials. This is sufficient for assessment of the combinatorial contribution to the chemical potential, leading eventually to such fundamental thermodynamic properties, as activity coefficients, excess properties, partition coefficients, phase diagrams, as well as solubility. Irrespectively of the computed property, a correct representation of the molecular structure is mandatory for adequate representation of statistics of segments surface electron densities. Hence, exploration of the conformational space is necessary. Here, the COMSOconf [133] software was used for generation of the most energetically favorable structures of all considered molecules. Taking advantage of BIOVIA TmoleX [134], computations were conducted using TURBOMOLE [135] software with RI-DFT BP86 (B88-VWN-P86) functional and def-TZVP basis set for geometry optimization and def2-TZVPD basis set for single point calculations with inclusion of the fine grid marching tetrahedron cavity. Additionally, the parameter sets with hydrogen bond interaction were included and van der Waals dispersion terms were quantified using the "D3" method of Grimme et al. [129]. All of the calculations of solubility were performed using COSMOtherm (version 22.0.0, revision 5273 M) [125] with BP_TZVPD_FINE_20.ctd parametrization.
Solubility computations of solid solutes were performed by iteratively solving the following equation: where µ i o stands for the chemical potential of a pure compound i, µ s i x s,(k) i is the solute chemical potential generated after k-th iteration with estimated solubility value equal to x s,(k) i and ∆G fus represents the solute Gibbs free energy of fusion. The starting point, for which k = 0, represents the zero-order approximation and characterizes solubility at infinite dilution conditions of the solute estimated using values of limiting activity coefficient, (4) is solved iteratively until self-consistency is obtained by computing solubility using chemical potential estimated for mole fraction of the previous iteration. The form of solubility definition in Equation (4) requires the knowledge of the exact or at least approximated value of fusion Gibbs free energy, which must be provided from external measurements or estimations.

COSMO-RS DARE and Metafiles
Application of COSMO-RS for predicting physicochemical properties of any bulk system requires a proper representation of the chemical diversity of all components by defining a representative pool of conformations. This is not only related to the flexibility of the molecules adopting potentially many geometries but also to direct intermolecular interactions eventually leading to new compounds. This can originate from reactions occurring in the systems, for example, polymerization or hydrogen bonded complexes between solute solvent molecules. These species are outside of routine conformational analysis but might be important due to charge densities alterations affected by the very fact of complex formation. This problem is addressed by COSMO-RS with Dimerization, Aggregation and Reaction Extension (COSMO-RS-DARE) [96].
In this approach, the list of conformers of a given compound is extended by inclusion of products of any reaction, in which σ-surface includes only one molecule and weighting out all the rest of product. This ensures that the number of atoms and total formal charge suites monomeric forms not involved in the product formation. This is conducted via so-called "mcos" files, which have atomic weights set to one for atoms of considered species and zero to atoms of the second interacting component. Hence, such pseudo-conformers can be used for enriching the overall pool of conformers of a given compound. However, these interaction σ-surface segments of these pseudo-conformers are not captured with default settings of COSMO-RS and a new descriptor must be added ensuring that all new surface segments are included in computations. This is performed by introducing two adjustable parameters in the formal form of the interaction Gibbs free energy between i-th and j-th compounds: Hence, the COSMO-RS-DARE methodology computes interaction energy between surface segments including misfit, hydrogen bonding and van der Waals forces using the above formula of two interacting pseudo-conformers. The parameters of Equation (5) are estimated by fitting the computed solubility values to experimental ones for every solvent.

Contacts Conformational Analysis
The application of COSMO-RS DARE requires the definition of pseudo-conformers necessary for "mcos" files. Hence, conformational screening of bi-molecular systems is necessary, in which clusters of potentially probable structures are generated and optimized. The COSMOtherm program offers such kind of computation by performing the contact statistics based on the probability of interactions between molecule surface segments. Practically it is conducted by using "CONTACT = {1 2} ssc_probability ssc_weak ssc_ang = 15.0" command. This leads to automatic generation of series of geometries of high probability of contacts by alteration of the value of the dihedral angle between the two contacting molecules in the cluster with a step size of 15 • . Weak interactions are also included in the probability statistics. The generated geometries were further optimized adopting the same scheme as for conformations analysis conducted for monomers. The number of pairs used for further analysis was restricted by the energy window of 2.5 kcal/mol, and pairs with stabilization energy above this threshold were discarded.

Conclusions
The solubility of caffeine was studied in several choline chloride-based natural deep eutectic solvents (ccNADES) and their mixtures with water. The results of performed optimization of the NADES component ratio suggested that unimolar proportion of cationic and anion counterparts had the highest ability of caffeine dissolution. The only exception was found for the system comprising choline chloride and glycerol in 1:2 molar ratio. Such ccNADES formulation was found to be the most effective among all seven polyolic components studied in 2:1, 1:1 and 1:2 molar ratios. Other ccNADES, particularly the ones utilizing sorbitol, xylitol and glucose, were also very effective in terms of caffeine dissolution. It was found that water and ccNADES exhibit co-solvency since the addition of a small amount of water to the ccNADES resulted in an increased caffeine solubility compared to the water-free ones. The composition involving a 0.8 molar ratio of ccNADES in the mixture turned out to be the most effective one. In general, it can be concluded that ccNADES perform better than standard organic solvents studied earlier, both in the water-free form and in aqueous mixtures. The usage of NADES is justified not only by their effectiveness but also by their beneficial health and environmental properties, especially when compared to classic organic solvents usually experiencing serious toxicity issues.
Such high solubilizing character ccNADES is granted from strong intermolecular interaction of caffeine with solvent constituents resulting in stable bi-molecular clusters formation. The highest contributions come from affinities of caffeine to sugars followed by interactions with choline chloride and water. Caffeine dimers are also allowed but these pairs are least stable among all considered intermolecular complexes. This is most probably related to stacking interactions rather than hydrogen bonding, although the latter is typical for all the rest of clusters. It was found that in all cases, the nitrogen center of the imidazole ring of caffeine plays the role of an accepter, while the hydroxyl group of the second counterpart of the hydrogen bonded complex acts as the hydrogen donor. It is interesting to note that not including these interactions in the default solubility computation performed with an aid of COSMO-RS approach results in caffeine solubility values which are at most in qualitative agreement with experimental data. This holds not only for saturated solutions of caffeine in neat solvents, but also in binary solvent mixtures and ccNADES. Fortunately, extending the system definition by adding interaction parameters dramatically increases accuracy leading to almost an exact match between measured and computed data. Unfortunately, the enthusiasm resulting from such effectiveness is tempered by the fact that fitting parameters cannot be obtained without prior experiments. However, the interaction parameters provide interesting characteristics of the studied systems and can be used for precise computation of solubility values for other formulations and temperatures. Moreover, fitting parameters used in COMSO-RS-DARE computations exhibit quite simple linear trends with ccNADES concentrations in aqueous mixtures. However, there is a clearly identifiable splitting of such linear trends around x ccNADES ≈ 0.2. This observation was rationalized by diversities of intermolecular interactions between caffeine and water compared to the affinities to the rest of ccNADES constituents.
Finally, there has been an established correlation between caffeine fusion data and interaction parameters. This is a very interesting and valuable observation suggesting that COSMO-RS-DARE is able to compensate for errors of solubility data originating not only from inappropriate conformers definition in default COSMO-RS solubility computations but also inaccuracies or lack of thermodynamic data of fusion. This suggests that solubility computations can be extended also for systems for which the fusion data were not measured or are immeasurable. The latter case is quite common and originates from the fact that many organic solids are unstable in the melting temperatures, which results from decomposition below melting, high sublimation and solid forms transformations. Bearing in mind this feature of COSMO-RS-DARE, not reported so far, some important notes on a proper definition of "meta" files were provided. It was suggested that detailed analysis of σ-profiles similarities can effectively reduce the size of the optimization problem.