Solubility Characteristics of Acetaminophen and Phenacetin in Binary Mixtures of Aqueous Organic Solvents: Experimental and Deep Machine Learning Screening of Green Dissolution Media

The solubility of active pharmaceutical ingredients is a mandatory physicochemical characteristic in pharmaceutical practice. However, the number of potential solvents and their mixtures prevents direct measurements of all possible combinations for finding environmentally friendly, operational and cost-effective solubilizers. That is why support from theoretical screening seems to be valuable. Here, a collection of acetaminophen and phenacetin solubility data in neat and binary solvent mixtures was used for the development of a nonlinear deep machine learning model using new intuitive molecular descriptors derived from COSMO-RS computations. The literature dataset was augmented with results of new measurements in aqueous binary mixtures of 4-formylmorpholine, DMSO and DMF. The solubility values back-computed with the developed ensemble of neural networks are in perfect agreement with the experimental data, which enables the extensive screening of many combinations of solvents not studied experimentally within the applicability domain of the trained model. The final predictions were presented not only in the form of the set of optimal hyperparameters but also in a more intuitive way by the set of parameters of the Jouyban–Acree equation often used in the co-solvency domain. This new and effective approach is easily extendible to other systems, enabling the fast and reliable selection of candidates for new solvents and directing the experimental solubility screening of active pharmaceutical ingredients.


Introduction
Acetaminophen (N-(4-hydroxyphenyl)acetamide, CAS: 103-90-2), also known under the name paracetamol, was synthesized as early as 1878 [1] and was later used clinically [2]. Since then, acetaminophen has become one of the most often used and prescribed drugs in the world. Paracetamol mainly acts as an analgesic and antipyretic, has a very weak antiinflammatory effect and has a low risk of side effects compared to aspirin and nonsteroidal anti-inflammatory drugs [3]. The importance of this drug and its widespread occurrence are the reasons that it is still the subject of many research papers dealing with all aspects of its usage [4,5].
Phenacetin (N-(4-ethoxyphenyl) acetamide, CAS: 62-44-2) is structurally similar to acetaminophen but contains an alkoxy group instead of the hydroxyl substituent. It was widely used together with acetaminophen as a pain-relieving and fever-reducing drug until its withdrawal from medicinal usage. Even though phenacetin was, for a long time, used more often than acetaminophen, it was banned by different organizations due to its adverse effects on human health, including carcinogenic and kidney-damaging properties [6,7]. Interestingly, nowadays, it has been found to be one of the main cocaine adulterants worldwide [8,9].
One of the main concerns of the pharmaceutical industry is the utilization of appropriate solvents for all stages of creating a pharmaceutical formulation. This is not a trivial task, since the limited solubility of various active pharmaceutical ingredients (APIs) is widely recognized as a major obstacle [10][11][12]. One of the strategies used to overcome this problem is the selection of a proper solvent system for the studied pharmaceutical formulation, which can be a tedious task. Therefore, a screening stage can be very helpful in order to obtain preliminary data and limit resource-consuming experiments. In this context, the utilization of machine learning can be a very useful approach for solubility predictions and solvent design. The growing awareness in the scientific community about the nonprecedential value of this methodology is associated with expanding the range of free programming ecosystems for direct use, including Scikit-learn [13], TensorFlow [14], PyTorch [15], Keras [16], OpenNN [17] or DeterminedAI [18], to mention only a few representative examples. Bearing this in mind, it is not surprising that many different approaches were used for inferring solubility from different representations of molecular structures. Indeed, recently, Panapitiya et al. [19] assessed the current deep learning architectures for solubility predictions and molecular representations used for depicting patterns between structural molecular properties and measured molecular solubility. Despite great efforts made in the development of deep machine learning models (DMLM) devoted to solubility predictions, there are still serious limitations prohibiting the general use of already developed models [20]. First of all, they concentrate mainly on aqueous solubility [19,[21][22][23][24][25][26][27], with few exceptions [21,[28][29][30] devoted to studying organic solvents. An additional limitation comes from the fact that the temperature dependence is often not explicitly defined. Moreover, the applicability domain is not always defined, which makes difficult generalizations and predictions for new systems. On the other hand, DMLM requires a huge number of data, which are unavailable for many potentially interesting but sparingly explored solvents or their mixtures. Hence, screening for new and green solvents is limited.
The aim of this paper is threefold. First of all, the theoretical tool for screening effective and environmentally friendly solvents for acetaminophen and phenacetin dissolution was developed. For this purpose, a nonlinear DMLM was developed using a newly developed ensemble of neural networks (ENN). For learning purposes, the curated experimental solubility data taken from the literature and new measurements were used. Molecular descriptors computed within the COSMO-RS [31] framework were used as a molecular representation of the physicochemical properties of the studied systems in saturated solutions. Finally, the applicability of the developed model for screening for new effective, neat and binary solvents mixtures was documented and discussed.

Solubility Data Collection and Curation
The results of solubility measurements of acetaminophen (A) and phenacetin (P) were searched in the literature. Both active pharmaceutical ingredients were the subjects of quite intensive investigations in an extended range of neat solvents and multicomponent solvent mixtures. Unfortunately, due to differences in the applied methodology, in some cases, quite large deviations could be noticed between the reported values. Since nonlinear models are very sensitive to the quality of the dataset used for training, the observed discrepancies need to be resolved prior to the application of DMLM. Hence, in the very first step, the data have been meticulously analyzed for removing incongruences. This step was done by a careful inspection of the back-computed data obtained after fitting to the Buchowski-Ksiazczak equation [32] and the Jouyban-Acree model [33] in the case of neat and binary solvents, respectively. The former model, also known in the literature as the λh-model, is a two-parameter equation that is often used for fitting of the solidliquid equilibrium data. Formally, the mole fraction of the solute is defined and computed as follows: where λ and h are adjustable parameters, and T m stands for the melting temperature. This equation is very popular and is often used for solubility data interpretation due to its high flexibility and the ability to adequately represent temperature-related solubility. Very often, the accuracy of the back-computed mole fraction is within the experimental uncertainty. The popularity of the λh-model is also related to the fact that the parameters have physical meaning, since h represents the energetics of the solubilization, and λ is an indicator of the solute association. The curation procedure was performed in a two-step manner. Initial fitting was done using the whole set of data reported for the given solution.
Then, outliers were identified by finding cases with a percentage error deviating higher than 10% from the computed values. After removing these points, the final fitting was done, and the obtained parameters of the applied theoretical models were collected in the Supplementary Materials S2 in "neat solvents". In all cases where data needed to be cured, the ones back-computed with the aid of Equation (1) were used as experimental points. If only single measurements were available, the original experimental data were used. As an illustration of this step, the solubility of A and P in neat water is provided in Figure 1. It is visible that a cloud of points was collected for saturated water solutions of A. The points annotated with open circles are those that were excluded from the pool of data due to not fulfilling the inclusion criterion. The details of all the considered systems are provided in the Supplementary Materials S1 in Tables S1 and S2.
the case of neat and binary solvents, respectively. The former model, also known in the literature as the λh-model, is a two-parameter equation that is often used for fitting of the solid-liquid equilibrium data. Formally, the mole fraction of the solute is defined and computed as follows: where λ and h are adjustable parameters, and Tm stands for the melting temperature. This equation is very popular and is often used for solubility data interpretation due to its high flexibility and the ability to adequately represent temperature-related solubility. Very often, the accuracy of the back-computed mole fraction is within the experimental uncertainty. The popularity of the λh-model is also related to the fact that the parameters have physical meaning, since h represents the energetics of the solubilization, and λ is an indicator of the solute association. The curation procedure was performed in a two-step manner. Initial fitting was done using the whole set of data reported for the given solution. Then, outliers were identified by finding cases with a percentage error deviating higher than 10% from the computed values. After removing these points, the final fitting was done, and the obtained parameters of the applied theoretical models were collected in the Supplementary Materials S2 in "neat solvents". In all cases where data needed to be cured, the ones back-computed with the aid of Equation (1) were used as experimental points. If only single measurements were available, the original experimental data were used. As an illustration of this step, the solubility of A and P in neat water is provided in Figure 1. It is visible that a cloud of points was collected for saturated water solutions of A. The points annotated with open circles are those that were excluded from the pool of data due to not fulfilling the inclusion criterion. The details of all the considered systems are provided in the Supplementary Materials S1 in Tables S1 and S2.  [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48] and phenacetin (P) [49][50][51] in neat water. Open circles represent data points with percentage deviations from back-computed values higher than 10% and were excluded from the analysis.
In principle, the Buchowski-Ksiazczak model is also applicable to multicomponent systems, provided that every solvent composition is considered as a new solvent. However, the solubility of A and P in binary solvent mixtures was determined for different compositions, which makes it difficult to fit the data coming from different authors. Thus, another very popular and accurate model was used instead. In the co-solvency literature, the Jouyban-Acree model [33] was proven to be quite accurate. This approach models the excess solubility by fitting with an aid of a polynomial in the following form: * • * • * • * • • * * (2) Figure 1. The example of the solubility data curation of acetaminophen (A) [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48] and phenacetin (P) [49][50][51] in neat water. Open circles represent data points with percentage deviations from backcomputed values higher than 10% and were excluded from the analysis.
In principle, the Buchowski-Ksiazczak model is also applicable to multicomponent systems, provided that every solvent composition is considered as a new solvent. However, the solubility of A and P in binary solvent mixtures was determined for different compositions, which makes it difficult to fit the data coming from different authors. Thus, another very popular and accurate model was used instead. In the co-solvency literature, the Jouyban-Acree model [33] was proven to be quite accurate. This approach models the excess solubility by fitting with an aid of a polynomial in the following form: The three parameters, J 0 , J 1 and J 2 , are regressed against experimental data. The curation procedure was similarly applied as for neat solvents. Details of the studied systems are provided in the Supplementary Materials S1 in Tables S3 and S4. Additionally, the values of the obtained parameters with statistical measures of accuracy were collected in the Supplementary Materials S2 in "binary solvents". The following analytical-grade compounds were used without any initial procedures. Acetaminophen (A, CAS: 103-90-2) and phenacetin (P, CAS: 62-44-2) were purchased from Sigma-Aldrich (Poznań, Poland), as well as dimethyl sulfoxide (DMSO, CAS: 67-68-5), N,N-dimethylformamide (DMF, CAS: 68-12-2) and 4-formylmorpholine (4FM, CAS: 4394-85-8). The 5.0-grade nitrogen was obtained from Linde (Warsaw, Poland) and used during the DSC experiments.

Solubility Measurement Protocol
First, an excess amount of either of the APIs was placed in a test tube, which was then filled with 10 mL of solvent, which enabled to obtain saturated solutions. The binary mixtures were prepared by mixing the organic solvent with water in appropriate molar fractions. For each system, three samples were prepared. Samples prepared in such a manner were placed in an Orbital Shaker ES-20/60 incubator from Biosan (Riga, Latvia) and incubated for 24 h at four different temperatures ranging from 20 • C to 40 • C with intervals equal to 5 • C. The temperature was adjusted with ±0.1 • C accuracy, and its daily variance was ±0.5 • C. Shaking of the samples at 60 rev/min was applied together with heating. The next step involved the filtration of the samples using syringes with 0.22 µm pore size PTFE filters. In order to avoid precipitation, the test tubes, syringes and filters were heated before the filtration at appropriate temperatures, matching the ones of the solutions. Next, a small volume of the filtrate was transferred using an automatic pipette to a tube filled with methanol in order to conduct spectroscopic measurements. For the calculation of mole fractions of the solute, the density of the samples was also determined by weighing a fixed amount of the solution in a 10 mL volumetric flask. The calibration curve was prepared with the use of stock solutions of acetaminophen and phenacetin in methanol with concentrations equal to 2.00 · 10 −5 mg/mL and 1.89 · 10 −5 mg/mL, respectively. Small amounts of these stock solutions were successfully diluted with methanol in 10 mL volumetric flasks in order to obtain solutions with decreasing concentrations. The absorption of these solutions was measured using an A360 spectrophotometer from AOE Instruments (Shanghai, China). The analytical wavelengths were selected as 249 nm for both acetaminophen and phenacetin. Three separate curves were prepared, and the final curve used for solubility determination was the result of their averaging. The spectra of the prepared samples were measured using the same device with a wavelength range from 215 to 500 nm with a 1 nm resolution. In order to ensure that the absorbance values do not exceed 2.5, the samples were diluted with methanol accordingly. The solubility of acetaminophen and phenacetin in the samples was determined based on the calibration curves, and their mole fractions in the solutions were calculated.

FTIR and DSC Characteristics of Solid Residues
The Fourier-transform infrared spectroscopy (FTIR) and differential scanning calorimetry (DSC) techniques were utilized to analyze the solid residues left after solubility measurements. Before the actual analysis, the samples were left to dry by air. A FTIR Spectrum Two spectrophotometer from Perkin Elmer (Waltham, MA, USA) was used with a diamond attenuated total reflection (ATR) device. For the DSC measurements, the DSC 6000 calorimeter, also from Perkin Elmer (Waltham, MA, USA), was used. The parameters were set as the heating rate of 5 K/min and the nitrogen flow of 20 mL/min. Zinc and indium reference standards were used for calibration, and the samples were placed in standard aluminum pans.

Deep Machine Learning Approach (DML)
Among many available algorithms developed for DML, artificial neural networks (ANNs) are one of the most universal, flexible and effective ways of model development.
For the purpose of this project, the simplest form of neural network architecture was implemented, which comprised just one hidden layer in-between the input and output ones. Despite its simplicity, such ANNs can be very effective [28] if an ensemble of neural networks model (ENNM) is constructed. The effectiveness of the network is boosted by tuning the ANN hyperparameters for the most accurate reproduction of the experimental data. Such properties as the number of neurons in the hidden layer, the form of activation functions and the error function were optimized. The training procedure was done after splitting the whole dataset into three subsets, namely training (70%), test (15%) and validation (15%), with the proportions indicated in parentheses. An ensemble of neural networks was constructed by collecting the best ANNs fulfilling the inclusion criteria: (1) accuracy (RMSD < 0.06), (2) precision (number of outliers ≤ 2%) and (3) reliability (predicted solubility within the formal range of log(x) between 0 and 1 for at least 99% of computed values). The applicability domain (AD) [52][53][54] was analyzed for each ANN and the whole ensemble.

Affinity Descriptors
Despite the very large number of available molecular descriptors, their application to the temperature-related solubility models is not straightforward, and many developed parameters are not intuitive [55]. It seems to be preferable to use such physicochemical properties, which are directly related to the chemical structures in multicomponent systems at defined external conditions. Fortunately, the COMOS-RS framework offers a very intuitive set of chemical system characteristics that might be used for DMLM development. For example [28], the σ-chemical potential trends, µ s (σ), are rich in information characterizing the chemical diversity of a given compound in varying environments. It represents the summarized σ-profiles with the inclusion of the mixture composition. It is derived [56] by iteratively solving the exact equation: where σ represents the charge density; P s (σ) is the so-called σ-profile representing the histogram of charge densities; a eff is an average molecular contact area; e(σ, σ ) is the sum of the three (misfit, hydrogen bonding and dispersion) contributions to the intermolecular interaction and RT is the multiplication of the temperature and the gas constant. In Figure 2, there are presented exemplary plots of descriptors derived using µ σ of A and P in pure water at room temperature. For data reduction purposes, the descriptors are defined as averaged values in the ranges defined in Figure 2 for three meaningful subregions quantifying fundamental affinity properties. It is assumed [56] that the electronegative charge distribution region, σ ∈ <−0.03, −0.01>, characterizes the affinity for HB donors, HBA (hydrogen-bonding acceptability). On the opposite scale, there is the electropositive polarity interval σ ∈ <0.01, 0.03> to which the affinity for HB acceptors is associated, HBD (hydrogen bond donicity). The intermediate range, σ ∈ <−0.01, +0.01>, characterizes contributions to nonpolar interactions and is regarded as a measure of hydrophobicity, HYD. In Figure 2, there are the lines plotted for the pure solute and the solute at an infinite dilution in a given solvent or solvent mixture. Additionally, there are provided plots of relative values of an inverse plot of the solute and solution [28]. These values represent the area between the former plots and were used as molecular descriptors. Hence, the following measures of mutual affinity were used for model training: 1. HBD = HBD solute (inversed) − HBA solvent -the relative solute donicity with respect to solute acceptability; 2. HBA = HBA solute (inversed) − HBD solvent -the relative value of solute acceptability with respect to solvent donicity D; 3. HYD = HYD solute (inversed) − HYD solvent -the relative hydrophobicity measure.
BIOVIA, San Diego, CA, USA). The RI-DFT BP86 (B88-VWN-P86) level of theory was used. The geometry optimization utilized the def-TZVP basis set, while the single point calculations were conducted using the def2-TZVPD basis set. Additionally included were the fine grid tetrahedron cavity and the parameter sets with a hydrogen bond interaction and van der Waals dispersion term based on the "D3" method of Grimme et al. [57]. The final solubility calculations were done with COSMOtherm (version 22.0.0, BIOVIA, San Diego, CA, USA) [58], and the BP_TZVPD_FINE_21.ctd parametrization was used.

Estimated Solubility as a Molecular Descriptor
The COSMO-RS (Conductor-like Screening Model for Real Solvents) [59][60][61][62] is a commonly applied approach for the theoretical characteristics of liquid multicomponent systems. It takes advantage of both quantum chemistry and post-quantum statistical thermodynamics for computing many fundamental properties, including chemical activity. Although this approach cannot deal with solids, the solubility can be still computed provided that fusion properties are known. Fortunately, for both studied aromatic amides, these quantities were measured and reported several times in the literature [63].
where ∆ ̅ stands for the partial molar Gibbs free energy of melting at conditions corresponding to the solubility measurements. By definition, this term for the pure solute The above definition can be used for deriving the temperature-dependent values of the molecular descriptors for any solute in any solvent or solvent mixture. However, prior to the application of Equation (3), it is necessary to properly represent the molecular structure. This is done after performing a full conformational analysis and identifying the most energetically favorable conformations. For these calculations, COMSOconf software (version 20.0.0, BIOVIA, San Diego, CA, USA) utilized TURBOMOLE (revision V7.5.1, TURBOMOLE GmbH, Karlsruhe, Germany) interfaced with TmoleX 2021 (version 21.0.1, BIOVIA, San Diego, CA, USA). The RI-DFT BP86 (B88-VWN-P86) level of theory was used. The geometry optimization utilized the def-TZVP basis set, while the single point calculations were conducted using the def2-TZVPD basis set. Additionally included were the fine grid tetrahedron cavity and the parameter sets with a hydrogen bond interaction and van der Waals dispersion term based on the "D3" method of Grimme et al. [57]. The final solubility calculations were done with COSMOtherm (version 22.0.0, BIOVIA, San Diego, CA, USA) [58], and the BP_TZVPD_FINE_21.ctd parametrization was used.

Estimated Solubility as a Molecular Descriptor
The COSMO-RS (Conductor-like Screening Model for Real Solvents) [59][60][61][62] is a commonly applied approach for the theoretical characteristics of liquid multicomponent systems. It takes advantage of both quantum chemistry and post-quantum statistical thermodynamics for computing many fundamental properties, including chemical activity. Although this approach cannot deal with solids, the solubility can be still computed provided that fusion properties are known. Fortunately, for both studied aromatic amides, these quantities were measured and reported several times in the literature [63].
where ∆ f us G m i stands for the partial molar Gibbs free energy of melting at conditions corresponding to the solubility measurements. By definition, this term for the pure solute at its melting point equals zero. The actual solubility computations conducted using COSMO-RS rely on iteratively solving the following equation: is achieved. The computation completes itself when the computed solubility difference is below a certain threshold value.
It is worth admitting that, for many systems already studied, the solubility computations within the COSMO-RS framework are often only qualitatively accurate [64][65][66]. Nevertheless, despite the discrepancies between measured and computed solubility values, they still capture an important part of system information, since the computed values of solubility are quite well correlated with the experimental ones (R 2 = 0.907). However, in the studied population, the value of the mean average percentage error (PE) exceeds 53% of the mole fraction solubility, and there are cases with PE > 1000%. This means that the accuracy of the computed solubility is not sufficient for screening purposes but can be very effective if used as molecular descriptors for nonlinear modeling. There is also another important issue related to solubility computations using the COSMO-RS approach. For some systems, the complete miscibility of the solute with solvent is obtained. This is illustrated in Figure 3 for acetaminophen. Dotted lines with crosses represent results obtained using COSMOtherm. Unexpectedly, at higher concentrations of DMSO (x 2 * > 0.6), the program fails in solubility computations, indicating complete miscibility. These artificial results were observed also for DMF and 4-formylmorpholine. For the consistent application of computed solubility, extrapolation was used in such cases, as indicated in Figure 3 by plots with open symbols. This prevents from using nonphysical values of solubility as molecular descriptors. Despite these shortcomings, it happened that the selection of the computed solubility as a molecular descriptor is appropriate, since neural networks associate a high contribution to this parameter. Superscripts, i and i + 1 in the above equation correspond to the values obtained i the two subsequent iterations. The repetition of the iterative cycle lasts until convergenc is achieved. The computation completes itself when the computed solubility difference i below a certain threshold value.
It is worth admitting that, for many systems already studied, the solubility compu tations within the COSMO-RS framework are often only qualitatively accurate [64][65][66] Nevertheless, despite the discrepancies between measured and computed solubilit values, they still capture an important part of system information, since the compute values of solubility are quite well correlated with the experimental ones (R 2 = 0.907 However, in the studied population, the value of the mean average percentage error (PE exceeds 53% of the mole fraction solubility, and there are cases with PE > 1000%. Thi means that the accuracy of the computed solubility is not sufficient for screening pur poses but can be very effective if used as molecular descriptors for nonlinear modeling There is also another important issue related to solubility computations using th COSMO-RS approach. For some systems, the complete miscibility of the solute wit solvent is obtained. This is illustrated in Figure 3 for acetaminophen. Dotted lines wit crosses represent results obtained using COSMOtherm. Unexpectedly, at higher concen trations of DMSO (x2* > 0.6), the program fails in solubility computations, indicatin complete miscibility. These artificial results were observed also for DMF an 4-formylmorpholine. For the consistent application of computed solubility, extrapolatio was used in such cases, as indicated in Figure 3 by plots with open symbols. This pre vents from using nonphysical values of solubility as molecular descriptors. Despite thes shortcomings, it happened that the selection of the computed solubility as a molecula descriptor is appropriate, since neural networks associate a high contribution to this pa rameter.

Intermolecular Interactions as Molecular Descriptors
The final set of molecular descriptors comes from the inspection of the solubility files and extraction of the three major contributions to the intermolecular interactions. Within the COSMO-RS framework, bulk systems are modeled as an association of closely packed molecules that are ideally screened and enclosed within a virtual conductor. The Coulomb interactions are the results of the screening of two contacting segments by their surroundings and, in turn, the back-polarization of the solute molecule. This "misfit" results in a specific interaction energy per unit area that represents the electrostatic part of the total energy: where a eff stands for the effective contact area between two surface segments, with α being the adjustable parameter. Similarly, hydrogen bonding (H-Bond) can be described by the charge density of the screening of neighboring strongly negative and positive centers for donors and acceptors, respectively. Such interactions are assumed to take place when contact occurs between two sufficiently polar surface pieces with opposite polarity, which is formally defined as: E HB σ, σ = a e f f c HB min(0; min(0; σ don + σ HB ); min(0; σ acc − σ HB )) where σ HB and σ HB are adjustable parameters. Additionally, the Van der Waals (vdW) interactions occurring between surface segments have to be included in the COSMO-RS model, which can be defined by the following relation: where τ vdW , τ vdW and c vdW are element-specific adjustable parameters. The Van der Waals energy is related only to the type of element present in the atoms experiencing surface contact. Any prediction made within the COSMO-RS framework results in obtaining these contributions, which makes them directly available from output files of a given computed property. Therefore, such energetic information can be utilized as a set of universal descriptors characterizing the properties of different systems. Here, the values of the total energy, E tot , of A or P in the given system were used as a molecular descriptor. In addition, the relative values of the contributions to the interaction energies were used. These energetic molecular descriptors were computed for the ith components as follows: where j corresponds to misfit, HB or vdW contributions. The descriptor values for multicomponent systems were computed as sums weighted by molar fractions. This set of molecular descriptors, taking advantage of intermolecular interactions, was already demonstrated [67] as adequate for solubility modeling.

Predictive Models Used for Data Curation and the Presentation of the Final Predictions
Although DML models offer very valuable help in predicting physicochemical properties, the direct application of the obtained model is prohibited by the necessity of demonstrating at least some elementary programming skills in Python or another high-level programming language. For overcoming this limitation, the final predictions were provided in a more intuitive way of application in terms of the parameters of the Jouyban-Acree model. For this purpose, the final ENNM was used for solubility data predictions, and the obtained values were fitted with an aid of Equation 2 (see Supplementary Materials S2 in "screening A" and "screening P").

Statistical Measures
The three following commonly used statistical measures were used in this paper. Firstly, the percentage error, PE, with the following general formula: where x exp is the experimental value, and x est is the estimated value. Secondly, the mean averaged percentage error, MAPE, described as: where x 1 est is the estimated mole fraction value, x 1 exp is the experimental mole fraction value and N is the number of experimental points. Finally, the root means squared deviation, RMSD, defined as:
ics 2022, 14, x FOR PEER REVIEW The distributions of the molecular descriptors are types of molecular descriptors characterize the solut pressed by computed solubility, the energetics of inte density distributions interpreted in terms of solute-so The distributions of the molecular descriptors are documented in Figure 5. The three types of molecular descriptors characterize the solute activity in saturated systems expressed by computed solubility, the energetics of intermolecular interactions and charge density distributions interpreted in terms of solute-solvent affinities.

Developed Model of Solubility
The performed training of the ANNs resulted in an ensemble model, the performance of which is documented in Figure 6. An acceptable accuracy was obtained with R 2 = 0.996. The number of outliers was less than 2%, defined by three normalized standard deviations.

Solubility Screening
For screening purposes, hypothetical aqueous binary mixtures of green organic solvents were tested. As the measure of environmental friendliness, the values of the indices defined by Harten et al. [83] were used. Solvent substitution is possible using PARIS III (Program for Assisting the Replacement of Industrial Solvents III, Version 1.4.0) software supported by the U.S. Environmental Protection Agency (EPA). PARIS III utilizes a series of computed indices assessing human toxicity potential upon ingestion or inhalation, as well as potentials arising from terrestrial, aquatic, ozone depletion, global warming, photochemical oxidation and acid rain toxicity. In particular, the environmental index (EI) was used as an indicator of the overall relative measure of hazardous properties. For screening purposes, only such solvents were selected that are characterized by EI < 0.5 [83]. This tight criterion unexpectedly excludes DMSO, despite that it is generally considered a green solvent. The reason for this is a very high EI = 11.7 value of DMSO coming from the photochemical oxidation potential contribution. Hence, the selection was done by arbitrarily excluding this contribution, which resulted in the reduction of the EI value of DMSO down to 0.26. It is worth noting that the final list of potential green solvents does not include some commonly used alcohols, starting from methanol and ending on 1-octanol. Additionally, many often-used esters are not included. 4-formylmorpholine is characterized by EI = 0.51, so it can be considered almost green using the criterion adopted for the purpose of this study. The final list of solvents used for screening was shortened from all 5200 available in PARIS III down to 50. The results of the screening are provided in the Supplementary Materials S2 in "screening A" and "screening P". The distributions of the molecular descriptors are documented in types of molecular descriptors characterize the solute activity in sat pressed by computed solubility, the energetics of intermolecular inter density distributions interpreted in terms of solute-solvent affinities. Figure 5. The collection of histograms of normalized values of molecular DMLM development. The following notation is used: HBA and HBD deno affinity and donicity, respectively. HYD denotes hydrophobicity. Ej describ energies of the solute with respect to the solvent, where j = tot stands for the to Figure 5. The collection of histograms of normalized values of molecular descriptors used for DMLM development. The following notation is used: HBA and HBD denote hydrogen-bonding affinity and donicity, respectively. HYD denotes hydrophobicity. ∆E j describes relative interaction energies of the solute with respect to the solvent, where j = tot stands for the total energy, j = HB for the hydrogen bonding, j = vdW for the dispersion contribution and j = misfit for the electrostatic contribution.

Developed Model of Solubility
The performed training of the ANNs resulted in an ensemble model, the perfor mance of which is documented in Figure 6. An acceptable accuracy was obtained with R = 0.996. The number of outliers was less than 2%, defined by three normalized standar deviations.

Solubility Screening
For screening purposes, hypothetical aqueous binary mixtures of green organi solvents were tested. As the measure of environmental friendliness, the values of the in dices defined by Harten et al. [83] were used. Solvent substitution is possible usin . Among these top-best solvents, only DMSO and 4FM can be regarded as green ones, according to the tight criteria set for the purpose of this study. Theoretical screening confirmed the above sequence and additionally enabled the identification of greener alternatives. It happened that the highest solubility of acetaminophen was found for DMSO, which was followed by neat ethyltriglycol and aqueous solutions of methyltriglycol and 2-pyridin-2-ylethanol. The exemplary solubility profiles are collected in Figure 7. Electron density scales are the same as in Figure 2. The red line represents solubilities predicted based on parameters fitted with Equation (2), while the dashed one stands for average ENNM solubility data.

Conclusions
The study was aimed at finding green and effective solvents adequate for acetaminophen and phenacetin processing in pharmaceutical practice. This goal was achieved via both the experimental and theoretical screening of binary aqueous solvent mixtures. The latter was achieved by formulating the ensemble neural network model, ENNM, which was derived after representing the molecular features of the studied compounds by using COSMO-RS-derived molecular descriptors. However, the experimental pool of data was cured prior to model development and extended with newly measured data. The considered aqueous binary solutions in which A and P were measured followed the previously observed high solubilizing potential of 4-formylmorpholine, which can also be considered a green solvent (EI = 0.51). This intuition was correct, since acetaminophen In general, the solubility of phenacetin is slightly lower compared to acetaminophen, and the experimentally observed sequence includes the following neat solvents: N,N-dimethylacetamide (log(x P exp,25 • C ) = −0.67, EI = 1.6), DMF (log(x P exp,25 • C ) = −0.93, EI = 2.2) and 4FM (log(x P exp,25 • C ) = −1.21, EI = 0.5). Additionally, 1,4-dioxane + water (x 2 * (opt) = 0.65, log(x A exp,25 • C ) = −1.33, EI = 0.8) has a relatively high potential of P dissolution. Unfortunately, the solvents that have been used so far are hardly green, with the exception of 4FM. Hence, it is interesting to inspect the results of the theoretical screening.
In the case of phenacetin, the screening procedure pointed to aqueous solutions of oleic acid, linoleic acid and octanoic acid as the potentially most-suited solvents. Unfortunately, these solvents are immiscible with water and cannot be used in practice. The aspect of miscibility was not a part of the model development. Probably the replacement of water with DMSO might help in obtaining the homogeneous mixture, which would also result in a green solvent (EI ≈ 0.2). Fortunately, the next solvents found on the prediction list are the same as the ones identified for acetaminophen. Hence, as documented in Figure 7, ethyltriglycol, methyltriglycol and pyridine-2-ethanol are suggested as first-choice solvents for future experiments with acetaminophen and phenacetin.

Conclusions
The study was aimed at finding green and effective solvents adequate for acetaminophen and phenacetin processing in pharmaceutical practice. This goal was achieved via both the experimental and theoretical screening of binary aqueous solvent mixtures. The latter was achieved by formulating the ensemble neural network model, ENNM, which was derived after representing the molecular features of the studied compounds by using COSMO-RS-derived molecular descriptors. However, the experimental pool of data was cured prior to model development and extended with newly measured data. The considered aqueous binary solutions in which A and P were measured followed the previously observed high solubilizing potential of 4-formylmorpholine, which can also be considered a green solvent (EI = 0.51). This intuition was correct, since acetaminophen is highly soluble in neat 4FM, although the solubility in such solvents as DMF and DMSO is slightly higher. 4-formylmorpholine was also identified as a fairly effective solvent for phenacetin dissolution, although the absolute solubility values are much lower compared to acetaminophen. A closer inspection into the type of solvents used so far for solubility studies of both considered aromatic amides revealed that environmentally friendly solvents were sparingly used. This was the direct motivation for the extensive screening of greener solvents, which are purchasable at a reasonable price. After shortening the list of candidates included in the PARIS III database, comprising 5200 solvents, by applying these criteria, the final screening was done for 50 solvents and their mixtures with water. It was found that both aromatic amides studied here can be effectively dissolved in such green solvents as ethyltriglycol, methyltriglycol and pyridine-2-ethanol. It is very tempting to perform appropriate measurements in the future for confirming these findings. Finally, it is worth adding that the newly developed model for the solubility prediction of acetaminophen and phenacetine relies on three types of criteria for the selection of the neural networks into an ensemble. The first one addresses the model accuracy, which is evaluated using the RMSD measure of every ANN. The second criterion is the model precision, which is meant as a minimization of the number of outliers lying beyond three times the normalized standard deviation. The third criterion addresses reliability, including restrictions on the formal range of the predicted values. In the case of solubility expressed as logarithmic values of mole fractions, this range restricts the predicted values to the span between 0 and 1. Hence, this triple combination of accuracy-precision-reliability criteria is an essential part of the developed ENNM.