Solid–Liquid Equilibrium in Co-Amorphous Systems: Experiment and Prediction

In this work, the solid–liquid equilibrium (SLE) of four binary systems combining two active pharmaceutical ingredients (APIs) capable of forming co-amorphous systems (CAMs) was investigated. The binary systems studied were naproxen-indomethacin, naproxen-ibuprofen, naproxen-probucol, and indomethacin-paracetamol. The SLE was experimentally determined by differential scanning calorimetry. The thermograms obtained revealed that all binary mixtures investigated form eutectic systems. Melting of the initial binary crystalline mixtures and subsequent quenching lead to the formation of CAM for all binary systems and most of the compositions studied. The experimentally obtained liquidus and eutectic temperatures were compared to theoretical predictions using the perturbed-chain statistical associating fluid theory (PC-SAFT) equation of state and conductor-like screening model for real solvents (COSMO-RS), as implemented in the Amsterdam Modeling Suite (COSMO-RS-AMS). On the basis of the obtained results, the ability of these models to predict the phase diagrams for the investigated API–API binary systems was evaluated. Furthermore, the glass transition temperature (Tg) of naproxen (NAP), a compound with a high tendency to recrystallize, whose literature values are considerably scattered, was newly determined by measuring and modeling the Tg values of binary mixtures in which amorphous NAP was stabilized. Based on this analysis, erroneous literature values were identified.


Introduction
The poor aqueous solubility of a significant number of newly developed active pharmaceutical ingredients (APIs) presents one of the most serious problems in the pharmaceutical industry. Therefore, strategies to enhance the aqueous solubility of APIs, and thus their bioavailability, have been developed intensively in recent years. A promising way to increase bioavailability is API amorphization. The amorphous forms of APIs exhibit a higher dissolution rate and apparent aqueous solubility in comparison with their crystalline counterparts [1][2][3], but are inherently thermodynamically unstable and tend to convert to their original crystalline forms. To stabilize API amorphous forms and improve their dissolution characteristics, formulation strategies based on glass solutions are widely explored. Glass solutions can be classified as polymeric and non-polymeric systems depending on the excipient used [4]. Co-amorphous systems (CAMs), whose phase behavior is studied in this work, belong to non-polymeric glass solutions and are defined as single-phase amorphous mixtures formed by low-molecular-weight components [5]. This work focuses on API-API CAMs, in which one of the APIs acts as an amorphous stabilizer for the second API through various mechanisms such as salt formation, hydrogen bonding, and π-π interactions [5,6]. Currently, only a small number of API-API CAMs have been reported, as it is challenging

Thermodynamic Fusion Properties and Glass Transition Temperatures of Pure APIs
First, the polymorphic forms of the APIs studied were identified by X-ray powder diffraction (see Figure S1 in the Supplementary Materials) and comparison to the Cambridge structural database (CSD) [26]. The melting temperatures T m and the enthalpies of fusion ∆ fus H of given polymorphs and the glass transition temperatures T g obtained in this work are listed in Table 1. The reported values for T m , ∆ fus H, and T g were measured using the same conditions as the data for binary phase diagrams to ensure the consistency between the two datasets. The fusion thermodynamic properties for IND, IBU, and NAP, as well as T g values for IND and IBU, are in close agreement with those reviewed and summarized by Štejfa et al. [27]. The T g value of NAP was determined in this work based by extrapolation from the measured T g values of binary mixtures with NAP (for details and comparison to the literature values, see Section 3.3.4). In the case of PAR and PRO, T m and ∆ fus H correspond to typical values reported for form I of these APIs, as collected by Acree and Chickos in their compendia [28,29]. The T g values obtained in this work for PAR and PRO are also close to those reported in the literature [13,30,31]. We note that the T g values depend on the thermal history of the sample, as well as on the experimental conditions under which they are measured, which may lead to differences in these values in the order of units of • C. Table 1. Thermodynamic fusion properties and glass transition temperatures for the APIs studied.  [27]. d This work. The value was obtained by extrapolation from the measured T g values of binary mixtures with NAP (see Section 3.3.4). e The value was determined based on the isobaric heat capacity data reported by Neau et al. [32]. f The value was determined in our laboratory by a combination of Tian-Calvet, power-compensated DSC, and relaxation calorimetry.

Binary Solid-Liquid Phase Diagrams
The thermograms obtained for the binary mixtures studied were typical of eutectic systems (Figure 1a). For all binary systems and most of the compositions studied, CAMs characterized by a single glass transition temperature (T g ) were obtained by quenching the melt of initially crystalline binary physical mixtures during DSC analysis ( Figure 1b). As NAP is a poor glass former, i.e., it has a high tendency to recrystallize, the CAMs were successfully formed only up to x NAP = 0.7 in the mixture with IND and x NAP = 0.5 in the mixtures with IBU and PRO. For mixtures exceeding the given NAP mole fraction, recrystallization appeared during the cooling of the melt, while for mixtures with a lower NAP content, recrystallization was observed on the heating curve at temperatures above T g (e.g., NAP-IBU or NAP-PRO in Figure 1b) or no recrystallization from melt was detected during heating (e.g., NAP-IND in Figure 1b).
The eutectic temperatures (T E ) evaluated from thermograms as extrapolated onset temperatures showed a negligible variation with composition (in accordance with theory, see Figure 2) and the mean values are reported in Table 2. The standard deviation of the mean was significantly less than the uncertainty in the determination of phase transition temperatures. The experimentally obtained T L , evaluated as the top of the liquidus peak (as recommended by Höhne [9]), are listed in Table 3 and shown in Figure 2. The eutectic compositions (x E ) were estimated based on Tammann plots [33] (see Figure S2 in the Supplementary Materials). Close to the eutectic composition, the liquidus peak overlapped with the eutectic peak, which did not allow us to evaluate T L for these compositions and made the integration of eutectic peak, and thus the eutectic composition estimation using the Tammann plots, less reliable.

NAP (1)-PRO (2) IND (1)-PAR (2)
The combined expanded uncertainty U c (0.95 level of confidence) in the determination of T L values is estimated to be 0.3 • C. b Deviation between the calculated and the experimentally determined liquidus temperature, The results may be affected by the approximative set of PC-SAFT parameters for PRO (for details, see Section 3.2.2).
The experimental SLE data were used to evaluate the performance of two computational models, PC-SAFT EOS and COSMO-RS-AMS, to predict the phase diagrams for these systems. To calculate the solubility, i.e., the liquidus curves, two types of thermodynamic data are needed (see Equation (1)): (i) thermodynamic fusion properties of pure APIs and (ii) the activity coefficients (γ L API ) of a given API in the liquid solution. The required melting temperatures, fusion enthalpies, and differences in the liquid and crystalline heat capacities are listed in Table 1. To predict γ L API , PC-SAFT EOS and COSMO-RS-AMS were employed. As the fusion thermodynamic data used in Equation (1) were determined experimentally (as is common because it is known that their prediction is burdened with large uncertainties and their determination is rather straightforward), the assessment of the performance of the two models reduces to a comparison of the quality of the γ L API .prediction. To make this comparison as fair as possible, PC-SAFT EOS was applied with the binary interaction parameters k ij set to 0 (for details, see Section 3.3.2), i.e., the γ L API prediction was carried out based on solely pure-component parameters for given APIs without any experimental input from the binary systems studied. COSMO-RS-AMS is a quantum-chemistry-based model that requires only the molecular structure as input for the γ L API prediction. The liquidus curves calculated using the two computational models are plotted together with the experimental data in Figure 2. Ideal solubility calculations are also shown in the phase diagrams to assess the ability of the two models to predict the direction of deviations from ideality. As shown in Figure 2, the experimental liquidus curves in the binary systems NAP-IND ( Figure 2a) and NAP-IBU (Figure 2b) are very close to those predicted based on the ideal solubility assumption, i.e., using γ L API equal to 1. For these two systems, COSMO-RS-AMS closely captures the observed trend in experimental T L data, while PC-SAFT EOS performs well only for the NAP-IND system, for which the predicted behavior by PC-SAFT EOS is also close to the ideal system. For the NAP-IBU system, PC-SAFT EOS predicts significant positive deviations, i.e., γ L API > 1, which is in disagreement with the experimental observation. For the IND-PAR system (Figure 2d), the trend in the experimental T L data is well described by both models. Based on the experimental T L data for NAP-PRO (Figure 2c), the system exhibits significant positive deviations from ideality, which is remarkably well captured by COSMO-RS-AMS, while the PC-SAFT EOS predicts negative deviations from ideality. As discussed in Section 3.3.2, the PC-SAFT parameter set for PRO was obtained using an approximate procedure owing to the unavailability of experimental thermodynamic data for pure PRO, which may be the reason PC-SAFT EOS does not provide satisfactory results for systems containing PRO. This situation points to the substantial limits of PC-SAFT EOS, a relatively highly parametrized model, for the initial screening of suitable excipients for a given API during which a high number of pairs of API-excipients are considered, including newly proposed or developed APIs or excipients. In such situations, it is highly probable that suitable thermodynamic data for PC-SAFT parametrization will not be available for all considered or preselected materials, and one would have to opt for approximative parametrization procedures, as in the case of PRO, which may lead to unreliable phase diagram predictions and, subsequently, API-excipient compatibility ranking.
The comparison of T E and x E obtained as an intersection of liquidus curves calculated using the two models and using the assumption of ideal solubility with the experimentally determined values are provided in Table 2. As mentioned above, we consider x E determined based on the Tammann plots to possess higher uncertainty owing to overlapping eutectic and liquidus peaks for compositions close to the eutectic composition, which makes the integration of eutectic peak less reliable. Table 3 presents the comparison between the calculated and experimental T L values. On the basis of these quantitative results, it can be concluded that COSMO-RS-AMS provides overall more reliable predictions of SLE for the systems studied than PC-SAFT EOS. Given the fact that PC-SAFT EOS is a more parametrized model compared with COSMO-RS-AMS and its parametrization requires the thermodynamic data (e.g., solubility in organic solvents, liquid density, or vapor pressure) for each pure component forming the mixture, which may not be easily accessible, especially for newly developed APIs, COSMO-RS-AMS seems to be a more suitable model, allowing for rapid initial screening of suitable excipient candidates for an API under consideration, formulation composition, or processing conditions. To conclude, in addition to a superior performance in predicting the binary phase diagrams for the systems studied, the lower input data requirement of COSMO-RS-AMS (only optimized molecular structure must be provided) presents a significant advantage of this model compared with PC-SAFT EOS, especially when screening of the compatibility of a high number of API-excipient pairs is required.
The phase behavior of the NAP-IND system was previously studied by Rades and co-workers [2,10,11]. Phase diagrams are presented only in graphical form and modeled as ideal systems in [2,10]. In addition, the differences between the crystalline and liquid heat capacities were neglected in the solubility calculations. In [11], only the eutectic composition evaluated based on the dependence of the enthalpy of the eutectic peak on composition is reported (x NAP = 0.60), which is in close agreement with the values obtained in this work (see Table 2). Experimentally determined T E = 127.9 • C [10] agrees well with the value obtained in this work (T E = 128.4 ± 0.3 • C), while a slightly higher T E of about 130 • C can be read from the graphical representation of the phase diagram in [2]. The eutectic composition reported in [2] (x NAP = 0.55) was obtained as an intersection of liquidus curves calculated based on the assumption of ideal solubility, while a slightly higher value (x NAP = 0.60) was determined experimentally by the authors in [10]. However, it is not clear by which method this value was derived from the experimental data (the experimental data points seem to be only connected by connecting lines and not described by correlation or computational model). Despite certain small differences and unclear data treatment in [10], it can be concluded that the phase diagrams presented in [2,10] and the eutectic composition reported in [11] agree well with the results of this work.
For the IND-PAR system, Fael and Demirel [18] reported that the physical mixture exhibited a melting peak at 142 • C, which can be associated with the eutectic temperature. This value is slightly higher than that determined in this work (T E = 138.7 • C).

Kinetic Stabilization of CAMs and Glass Transition of Binary Mixtures
As mentioned above, all initially crystalline physical mixtures were transformed to the amorphous state for compositions up to x NAP = 0.5 for the NAP-IBU and NAP-PRO and x NAP = 0.7 for the NAP-IND system. The mixture IND-PAR formed CAMs in the whole concentration interval. Although IBU, PRO, and IND have all been classified as good glass formers (Class 3 according to the classification proposed by Baird et al. [13]), IND offers the highest stabilization for NAP, an API with a high tendency to recrystallize (Class 1). It is also important to note that the mixtures close to the eutectic composition were successfully transformed into CAMs for all systems studied, which is a prerequisite step for the subsequent monitoring of their physical stability.
The glass transition temperature of binary mixtures plays an important role in the kinetic stabilization of CAMs. The physical stability of CAMs is typically proportional to the difference between its T g and the storage temperature. The T g values obtained for the CAMs studied are listed in Table S1 in the Supplementary Materials. Based on the physical stability study of three CAMs, including NAP-IND, Kissi et al. [10], and Beyer et al. [11], CAMs corresponding to the eutectic composition form the most stable CAMs. T g values at x E are close to typical storage temperature of 25 • C for NAP-IND (T g (x E ) ≈ 23.8 • C) and NAP-PRO (T g (x E ) ≈ 23.3 • C), significantly below for NAP-IBU (T g (x E ) ≈ −40.7 • C, x E was assumed to be x NAP = 0.1), and slightly higher for IND-PAR (T g (x E ) ≈ 30.7 • C). Although the kinetic stabilization of NAP-IND and IND-PAR CAMs possessing the eutectic composition derived on the basis of their T g values is rather limited, Kissi et al. [10] found that CAMs did not show any sign of recrystallization in about 35 days (the range was 31 to 37 days) when stored at room temperature under dry conditions. The prolonged physical stability for NAP-IND CAMs was also reported by Löbmann et al. [2] (x NAP = 0.5, storage temperatures of 4 and 25 • C, dry conditions, stability of at least 21 days) and Beyer et al. [11] (x NAP = 0.6, storage temperature of 21 • C, dry conditions, monitoring period of 56 and 112 days). Fael and Demirel [18] reported that the physical stability of IND-PAR CAMs for 2:1, 1:1, and 1:2 molar ratios was up to 7 months (various storage conditions were examined: 4, 25, and 40 • C under dry conditions and 29 • C under mild humid conditions and relative humidity of 55%). The optimum composition in terms of physical stability was found to be a molar ratio of 2:1, followed by a molar ratio of 1:1. The literature findings summarized above, along with the T g values of NAP-IND and IND-PAR CAM, which are close to typical storage temperatures, suggest that intermolecular interactions play a significant role in stabilizing CAMs.
It is important to mention that, in both studies [10,11] indicating the blends corresponding to eutectic composition as the most stable CAM, the eutectic composition was close to that of the equimolar mixture. Therefore, as stated by Kissi [10], the observed relation between CAM physical stability and the eutectic composition should be investigated for mixtures whose eutectic composition is located further away from the equimolar mixture. Such a system can be, e.g., NAP-PRO, whose phase behavior was studied in this work (it forms CAM at eutectic composition, which significantly differs from the equimolar composition).
Significant discrepancies in the T g values for NAP were identified in the literature [2,[14][15][16]. As reliable T g values for pure components (APIs and excipients) that form amorphous formulations represent key information in the evaluation of their kinetic stability, we attempted to clarify the situation regarding the T g of pure NAP in this work based on measuring and modeling T g of binary amorphous mixtures in which NAP was stabilized in the amorphous state. The experimental values on T g of the three mixtures containing NAP, i.e., NAP-IND, NAP-IBU, and NAP-PRO, listed in Table S1 in the Supplementary Materials, were correlated by the Gordon-Taylor equation, Equation (7), and the Kwei equation, Equation (8), with T g of pure NAP (T g, NAP ) as a fitted parameter. The modeled T g curves are shown in Figure 3 and the obtained T g, NAP values are summarized in Table 4. The mean value of T g, NAP = 6.4 ± 1.4 • C (the uncertainty quoted is the standard deviation of the mean) is in excellent agreement with the values determined by Paudel et al. [14] (T g, NAP = 6.2 • C, amorphous NAP prepared by spray drying) and by Löbmann et al. [2] (T g, NAP = 5.0 • C, amorphous NAP prepared melt quenching), but differs significantly from the value reported by Blaabjerg et al. [15] (T g, NAP = 56.1 • C, sample prepared by melt quenching method) and adopted, for example, by Kawakami [16]. Some variation in the measured T g values is expected because of their dependence on the thermal history of the amorphous material and, in general, on the experimental conditions using which they are determined. However, such a large deviation of approximately 50 • C cannot be explained by these phenomena and the value reported by Blaabjerg et al. [15] can be considered erroneous. In this work, we provide clear evidence that T g, NAP is about 6 • C, in accordance with the two previous studies [2,14].

Samples Description
The APIs studied in this research are listed with their basic characteristics in Table 5. APIs were used as received from the manufacturer without further purification. The thermodynamic fusion properties and T g values determined in this work are provided in Table 1.

Differential Scanning Calorimetry
A Q1000 differential scanning calorimeter (TA Instruments, Inc., New Castle, DE, USA) was used to determine melting temperatures and enthalpies of pure APIs, phase diagrams, and T g values. The physical mixtures containing two different APIs in different molar ratios were prepared by grinding with a mortar and pestle for 10 min. Approximately 5-10 mg of sample was then hermetically sealed in aluminum pans and analyzed by DSC. The experiments consisted of two heating cycles. The heating rate of 2 • C min −1 was applied for SLE measurements during the first heating run. Subsequently, the melt was quenched using a cooling rate of 10 • C min −1 and T g values were measured during the second heating cycle with a heating rate of 10 • C min −1 .

X-ray Powder Diffraction
X-ray powder diffraction analysis was performed to identify the polymorphic forms of the APIs studied using a θ-λ powder diffractometer X Pert3Powder in Bragg-Brentano parafocussing geometry using wavelength CuK radiation (λ = 1.5418 Å, U = 40 kV, I = 30 mA). Data were gathered using an ultrafast detector 1D PIXcel angular range 5-50 • (2θ) with a step size of 0.039 • (2θ) and 0.7 s for each step. HighScorePlus 4.0 software was used to analyze the obtained diffractograms. The polymorphic forms were identified based on the comparison to the CSD [26].

Modeling of Solid-Liquid Equilibria
The API solubility (mole fraction x L API ) was calculated according to the following equation: where ∆ fus H is the fusion enthalpy of pure API, T m is its corresponding melting temperature (in K), T is absolute temperature (in K), ∆ fus C p is the difference between the isobaric heat capacity of the liquid and the crystalline phase, and R is the universal gas constant. γ L API is the activity coefficient of one of the APIs in the liquid API-API mixture. γ L API was calculated using the PC-SAFT equation of state, COSMO-RS-AMS, or set as 1 in the case of ideal solubility calculations.

PC-SAFT Equation of State
According to the PC-SAFT equation of state, the residual Helmholtz energy (a res ) is commonly calculated as a sum of three different contributions resulting from repulsion (hard chain), van der Waals attraction (dispersion), and hydrogen bonding (association) [19,20]: From Equation (2), other thermodynamic properties of a system can be calculated, including the activity coefficient γ L API [19]. PC-SAFT considers molecules to be chains constituted by spherical segments. Materials are then characterized using the following set of pure component parameters: the number of segments within a chain (m i ), the diameter of the segment (σ i ), the dispersion energy parameter (ε i /k, k is the Boltzmann constant), the association energy parameter (ε i assoc /k), the association volume (κ i assoc ), and the number and type of association sites per molecule (N i assoc ). Given by the semi-empirical nature of the model, values of these parameters are routinely fitted to experimental data. Values of the PC-SAFT parameters for the studied APIs are given in Table 6. The parameters for IND, IBU, NAP, and PAR were taken from the literature [34]. For each API, they were obtained by fitting them to the properties of the pure liquid API along with the API solubility data in a row of pure solvents. However, the PC-SAFT parameters for PRO were not found in the available literature. At the same time, the literature lacked experimental solubility data of PRO in pure solvents. Therefore, an alternative parametrization approach was applied to PRO. First, the parameter κ i assoc was set to a constant value of 0.01, as usual for APIs in the literature [34]. For the association energy parameter, ε i assoc /k, we used the value typical for the phenolic OH group (1650 K) [35,36], while the structural parameters m i and σ i were estimated using a group contribution approach [37,38]. Finally, the dispersion energy parameter, ε i /k, remained the only adjustable parameter and was fitted to the only experimental solubility data point of PRO in a pure solvent (ethanol) available in the literature reported by Yagi et al. [39] (for completeness, the PC-SAFT parameters for ethanol were taken from [20]). Therefore, owing to this crudely approximative nature of the PC-SAFT parameter set for PRO, the PC-SAFT results for systems with this API should be taken with caution and considered to be only an illustration of what to expect from the model when parametrized in an alternative way because of the inaccessibility of experimental data. To emphasize this, all results from PC-SAFT for PRO are denoted with "!".  [34]. b Values determined in this work using an alternative approach that combined fitting, group contribution method, and structural similarity (see the text for details).
For the calculation of the thermodynamic properties of mixtures, the combination rules for the cross parameters σ ij and ε ij /k between components i and j are applied as follows [19]: where k ij is the binary interaction parameter, which can be calculated by fitting the experimental solubility data or set to 0. In the case of k ij = 0, SLE data are calculated based only on the PC-SAFT parameters of the pure components. In this work, the PC-SAFT was used solely with k ij = 0, i.e., no experimental data for the binary systems studied were employed. Combination rules for the association parameters can be found elsewhere [20]. γ L API is determined from the following equation: where ϕ L i is the fugacity coefficient of API in the liquid API-API mixture and ϕ L 0,i is the fugacity coefficient of the pure liquid. The fugacity coefficients are obtained from the general equation: where a res is the reduced residual Helmholtz energy obtained from PC-SAFT, Z is the compressibility factor, and ρ is the system molar density.

COSMO-RS-AMS
COSMO-RS represents an efficient and successfully used methodology for predicting the thermodynamic properties of fluid systems. It combines quantum chemical calculations of molecular properties and a statistical mechanical procedure to obtain the macroscopic properties of a solution [21]. The quantity that bridges these two steps is the sigma-profile of a molecule, which is a surface histogram with respect to the charge density calculated quantum chemically using density functional theory (DFT) and the COSMO model to imitate the solvent environment. As such, the COSMO-RS methodology allows for a priori predicting of phase equilibria without any experimental data, also including SLE in pharmaceutical systems (e.g., [40,41]). Multiple implementations of the COSMO-RS methodology are available in the literature (e.g., [22,42,43]). In this work, we used COSMO-RS as implemented in the Amsterdam Modeling Suite (AMS), version 2022.101 [22,44] (COSMO-RS-AMS). The sigma-profiles for IND, IBU, NAP, and PAR used in this work were taken from Klajmon [41], while that for PRO was determined in this work using AMS. Molecular geometry of PRO, which is the only input in producing the sigma-profile, was taken from the HAXHET01 crystal structure of form I [45] available in the CSD [26], which was considered to be appropriately representing predominant molecular geometries of PRO in the condensed phase [41].

Modelling of the Glass Transition Temperature Curve
The Gordon-Taylor equation [46] and the Kwei equation [47] were used to model T g values of the binary mixtures studied and to estimate the T g value for pure NAP. The Gordon-Taylor equation is defined as follows [46]: where T g is the glass transition temperature of a binary mixture, T g1 and T g2 are the glass transition temperatures of pure APIs, and x 1 and x 2 are their molar fractions. k is a parameter that is determined by fitting experimentally measured T g values. By adding the second fitting parameter q to the Gordon-Taylor equation, Equation (7), the Kwei equation is obtained, defined as follows [47]: where k has the same meaning as in the Gordon-Taylor equation, Equation (7), and q is the second fitted parameter.

Conclusions
In this work, the performance of the two computational models, PC-SAFT EOS and COSMO-RS-AMS, to predict the phase diagrams for four binary systems combining two APIs was evaluated based on the comparison to experimental SLE data. Overall, COSMO-RS-AMS outperformed PC-SAFT EOS and, given the fact that it is a significantly less parametrized model compared with PC-SAFT EOS, it can be considered as a more suitable computational tool for initial screening of phase diagrams of API-excipient pairs, and thus their compatibility.
Melting of the initial binary crystal mixtures and their subsequent quenching lead to the formation of CAMs for all binary systems and most of the studied compositions. NAP, an API with a high tendency to recrystallize, was successfully stabilized in its amorphous form in mixtures with IND, IBU, and PRO (APIs with good glass-forming ability), which allowed us to determine its T g and reconcile the values published in the literature.

Supplementary Materials:
The following supporting information can be downloaded at https:// www.mdpi.com/article/10.3390/molecules28062492/s1. Figure S1: XRPD patterns of APIs studied; Figure S2: Tammann plots for the systems studied; Table S1: T g values obtained for binary systems studied.